4 PROB ScoreToProb(SCORE Score)
\r
6 if (MINUS_INFINITY >= Score)
\r
8 return (PROB) pow(2.0, (double) Score/INTSCALE);
\r
11 static const double log2e = log2(exp(1.0));
\r
13 double lnTolog2(double ln)
\r
18 double log2(double x)
\r
21 return MINUS_INFINITY;
\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
29 SCORE ProbToScore(PROB Prob)
\r
32 return MINUS_INFINITY;
\r
33 // return (SCORE) floor(INTSCALE*log2(Prob));
\r
34 return (SCORE) log2(Prob);
\r
37 WEIGHT DoubleToWeight(double d)
\r
40 return (WEIGHT) (INTSCALE*d);
\r
43 double WeightToDouble(WEIGHT w)
\r
45 return (double) w / (double) INTSCALE;
\r
48 SCORE DoubleToScore(double d)
\r
50 return (SCORE)(d*(double) INTSCALE);
\r
53 bool ScoreEq(SCORE s1, SCORE s2)
\r
55 return BTEq(s1, s2);
\r
58 static bool BTEq2(BASETYPE b1, BASETYPE b2)
\r
60 double diff = fabs(b1 - b2);
\r
63 double sum = fabs(b1) + fabs(b2);
\r
64 return diff/sum < 0.005;
\r
67 bool BTEq(double b1, double b2)
\r
69 return BTEq2((BASETYPE) b1, (BASETYPE) b2);
\r
72 const double dLn2 = log(2.0);
\r
75 double pow2(double x)
\r
77 if (MINUS_INFINITY == x)
\r
82 // lp2(x) = log2(1 + 2^-x), x >= 0
\r
83 double lp2(double x)
\r
85 return log2(1 + pow2(-x));
\r
88 // SumLog(x, y) = log2(2^x + 2^y)
\r
89 SCORE SumLog(SCORE x, SCORE y)
\r
91 return (SCORE) log2(pow2(x) + pow2(y));
\r
94 // SumLog(x, y, z) = log2(2^x + 2^y + 2^z)
\r
95 SCORE SumLog(SCORE x, SCORE y, SCORE z)
\r
97 return (SCORE) log2(pow2(x) + pow2(y) + pow2(z));
\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
103 return (SCORE) log2(pow2(w) + pow2(x) + pow2(y) + pow2(z));
\r
106 SCORE lp2Fast(SCORE x)
\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
116 for (int i = 0; i < iTableSize; ++i)
\r
117 dValue[i] = (SCORE) lp2(i*dScale);
\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
129 // SumLog(x, y) = log2(2^x + 2^y)
\r
130 SCORE SumLogFast(SCORE x, SCORE y)
\r
132 if (MINUS_INFINITY == x)
\r
134 if (MINUS_INFINITY == y)
\r
135 return MINUS_INFINITY;
\r
138 else if (MINUS_INFINITY == y)
\r
143 dResult = x + lp2Fast(x-y);
\r
145 dResult = y + lp2Fast(y-x);
\r
146 assert(SumLog(x, y) == dResult);
\r
150 SCORE SumLogFast(SCORE x, SCORE y, SCORE z)
\r
152 SCORE dResult = SumLogFast(x, SumLogFast(y, z));
\r
153 assert(SumLog(x, y, z) == dResult);
\r
157 SCORE SumLogFast(SCORE w, SCORE x, SCORE y, SCORE z)
\r
159 SCORE dResult = SumLogFast(SumLogFast(w, x), SumLogFast(y, z));
\r
160 assert(SumLog(w, x, y, z) == dResult);
\r
164 double VecSum(const double v[], unsigned n)
\r
167 for (unsigned i = 0; i < n; ++i)
\r
172 void Normalize(PROB p[], unsigned n)
\r
176 for (i = 0; i < n; ++i)
\r
179 Quit("Normalize, sum=0");
\r
180 for (i = 0; i < n; ++i)
\r
184 void NormalizeUnlessZero(PROB p[], unsigned n)
\r
188 for (i = 0; i < n; ++i)
\r
192 for (i = 0; i < n; ++i)
\r
196 void Normalize(PROB p[], unsigned n, double dRequiredTotal)
\r
200 for (i = 0; i < n; ++i)
\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
209 bool VectorIsZero(const double dValues[], unsigned n)
\r
211 for (unsigned i = 0; i < n; ++i)
\r
212 if (dValues[i] != 0.0)
\r
217 void VectorSet(double dValues[], unsigned n, double d)
\r
219 for (unsigned i = 0; i < n; ++i)
\r
223 bool VectorIsZero(const float dValues[], unsigned n)
\r
225 for (unsigned i = 0; i < n; ++i)
\r
226 if (dValues[i] != 0.0)
\r
231 void VectorSet(float dValues[], unsigned n, float d)
\r
233 for (unsigned i = 0; i < n; ++i)
\r
237 double Correl(const double P[], const double Q[], unsigned uCount)
\r
239 double dSumP = 0.0;
\r
240 double dSumQ = 0.0;
\r
241 for (unsigned n = 0; n < uCount; ++n)
\r
246 const double dMeanP = dSumP/uCount;
\r
247 const double dMeanQ = dSumQ/uCount;
\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
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
262 const double dCorrel = dSum1 / sqrt(dSum2*dSum3);
\r
266 float Correl(const float P[], const float Q[], unsigned uCount)
\r
270 for (unsigned n = 0; n < uCount; ++n)
\r
275 const float dMeanP = dSumP/uCount;
\r
276 const float dMeanQ = dSumQ/uCount;
\r
281 for (unsigned n = 0; n < uCount; ++n)
\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
291 const float dCorrel = dSum1 / (float) sqrt(dSum2*dSum3);
\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
300 for (unsigned n = 0; n < uCount; ++n)
\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
311 else if (v < dValue)
\r
316 assert(uNumberEqual >= 1);
\r
317 assert(uNumberEqual + uNumberLess + uNumberGreater == uCount);
\r
318 Ranks[n] = (float) (1 + uNumberLess + (uNumberEqual - 1)/2.0);
\r
322 void Rank(const double P[], double Ranks[], unsigned uCount)
\r
324 for (unsigned n = 0; n < uCount; ++n)
\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
335 else if (v < dValue)
\r
340 assert(uNumberEqual >= 1);
\r
341 assert(uNumberEqual + uNumberLess + uNumberGreater == uCount);
\r
342 Ranks[n] = (double) (1 + uNumberLess + (uNumberEqual - 1)/2.0);
\r
346 FCOUNT SumCounts(const FCOUNT Counts[])
\r
349 for (int i = 0; i < 20; ++i)
\r