binaries/src/MSAProbs-0.9.7/MSAProbs/ScoreType.h

   1 /////////////////////////////////////////////////////////////////
   2 // ScoreType.h
   3 //
   4 // Routines for doing math operations in MSAPROBS
   5 /////////////////////////////////////////////////////////////////
   6
   7 #ifndef SCORETYPE_H
   8 #define SCORETYPE_H
   9
  10 #include <cmath>
  11 #include <algorithm>
  12 #include <cfloat>
  13 #include <assert.h>
  14
  15 typedef float ScoreType;
  16
  17 const float LOG_ZERO = -2e20;
  18 const float LOG_ONE = 0.0;
  19
  20 /////////////////////////////////////////////////////////////////
  21 // LOG()
  22 //
  23 // Compute the logarithm of x.
  24 /////////////////////////////////////////////////////////////////
  25
  26 inline ScoreType LOG(ScoreType x) {
  27         return log(x);
  28 }
  29
  30 /////////////////////////////////////////////////////////////////
  31 // EXP()
  32 //
  33 // Computes exp(x).
  34 /////////////////////////////////////////////////////////////////
  35
  36 inline ScoreType EXP(ScoreType x) {
  37         //return exp(x);
  38         if (x > -2) {
  39                 if (x > -0.5) {
  40                         if (x > 0)
  41                                 return exp(x);
  42                         return (((0.03254409303190190000 * x + 0.16280432765779600000) * x
  43                                         + 0.49929760485974900000) * x + 0.99995149601363700000) * x
  44                                         + 0.99999925508501600000;
  45                 }
  46                 if (x > -1)
  47                         return (((0.01973899026052090000 * x + 0.13822379685007000000) * x
  48                                         + 0.48056651562365000000) * x + 0.99326940370383500000) * x
  49                                         + 0.99906756856399500000;
  50                 return (((0.00940528203591384000 * x + 0.09414963667859410000) * x
  51                                 + 0.40825793595877300000) * x + 0.93933625499130400000) * x
  52                                 + 0.98369508190545300000;
  53         }
  54         if (x > -8) {
  55                 if (x > -4)
  56                         return (((0.00217245711583303000 * x + 0.03484829428350620000) * x
  57                                         + 0.22118199801337800000) * x + 0.67049462206469500000) * x
  58                                         + 0.83556950223398500000;
  59                 return (((0.00012398771025456900 * x + 0.00349155785951272000) * x
  60                                 + 0.03727721426017900000) * x + 0.17974997741536900000) * x
  61                                 + 0.33249299994217400000;
  62         }
  63         if (x > -16)
  64                 return (((0.00000051741713416603 * x + 0.00002721456879608080) * x
  65                                 + 0.00053418601865636800) * x + 0.00464101989351936000) * x
  66                                 + 0.01507447981459420000;
  67         return 0;
  68 }
  69
  70 /*
  71  /////////////////////////////////////////////////////////////////
  72  // LOOKUP()
  73  //
  74  // Computes log (exp (x) + 1), for 0 <= x <= 7.5.
  75  /////////////////////////////////////////////////////////////////
  76
  77  inline ScoreType LOOKUP (ScoreType x){
  78  //return log (exp(x) + 1);
  79  if (x < 2){
  80  if (x < 0.5){
  81  if (x < 0)
  82  return log (exp(x) + 1);
  83  return (((-0.00486373205785640000*x - 0.00020245408813934800)*x + 0.12504222666029800000)*x + 0.49999685320563000000)*x + 0.69314723138948900000;
  84  }
  85  if (x < 1)
  86  return (((-0.00278634205460548000*x - 0.00458097251248546000)*x + 0.12865849880472500000)*x + 0.49862228499205200000)*x + 0.69334810088688000000;
  87  return (((0.00059633755154209200*x - 0.01918996666063320000)*x + 0.15288232492093800000)*x + 0.48039958825756900000)*x + 0.69857578503189200000;
  88  }
  89  if (x < 8){
  90  if (x < 4)
  91  return (((0.00135958539181047000*x - 0.02329807659316430000)*x + 0.15885799609532100000)*x + 0.48167498563270800000)*x + 0.69276185058669200000;
  92  return (((0.00011992394456683500*x - 0.00338464503306568000)*x + 0.03622746366545470000)*x + 0.82481250248383700000)*x + 0.32507892994863100000;
  93  }
  94  if (x < 16)
  95  return (((0.00000051726300753785*x - 0.00002720671238876090)*x + 0.00053403733818413500)*x + 0.99536021775747900000)*x + 0.01507065715532010000;
  96  return x;
  97  }
  98
  99  /////////////////////////////////////////////////////////////////
 100  // LOOKUP_SLOW()
 101  //
 102  // Computes log (exp (x) + 1).
 103  /////////////////////////////////////////////////////////////////
 104
 105  inline ScoreType LOOKUP_SLOW (ScoreType x){
 106  return log (exp (x) + 1);
 107  }
 108
 109  /////////////////////////////////////////////////////////////////
 110  // MAX()
 111  //
 112  // Compute max of three numbers
 113  /////////////////////////////////////////////////////////////////
 114
 115  inline ScoreType MAX (ScoreType x, ScoreType y, ScoreType z){
 116  if (x >= y){
 117  if (x >= z)
 118  return x;
 119  return z;
 120  }
 121  if (y >= z)
 122  return y;
 123  return z;
 124  }
 125
 126  /////////////////////////////////////////////////////////////////
 127  // LOG_PLUS_EQUALS()
 128  //
 129  // Add two log probabilities and store in the first argument
 130  /////////////////////////////////////////////////////////////////
 131
 132  inline void LOG_PLUS_EQUALS (ScoreType &x, ScoreType y){
 133  if (x < y)
 134  x = (x <= LOG_ZERO) ? y : LOOKUP(y-x) + x;
 135  else
 136  x = (y <= LOG_ZERO) ? x : LOOKUP(x-y) + y;
 137  }
 138
 139  /////////////////////////////////////////////////////////////////
 140  // LOG_PLUS_EQUALS_SLOW()
 141  //
 142  // Add two log probabilities and store in the first argument
 143  /////////////////////////////////////////////////////////////////
 144
 145  inline void LOG_PLUS_EQUALS_SLOW (ScoreType &x, ScoreType y){
 146  if (x < y)
 147  x = (x <= LOG_ZERO) ? y : LOOKUP_SLOW(y-x) + x;
 148  else
 149  x = (y <= LOG_ZERO) ? x : LOOKUP_SLOW(x-y) + y;
 150  }
 151
 152  /////////////////////////////////////////////////////////////////
 153  // LOG_ADD()
 154  //
 155  // Add two log probabilities
 156  /////////////////////////////////////////////////////////////////
 157
 158  inline ScoreType LOG_ADD (ScoreType x, ScoreType y){
 159  if (x < y) return (x <= LOG_ZERO) ? y : LOOKUP(y-x) + x;
 160  return (y <= LOG_ZERO) ? x : LOOKUP(x-y) + y;
 161  }
 162  */
 163
 164 /*
 165  /////////////////////////////////////////////////////////////////
 166  // LOG()
 167  //
 168  // Compute the logarithm of x.
 169  /////////////////////////////////////////////////////////////////
 170
 171  inline float LOG (float x){
 172  return log (x);
 173  }
 174
 175  /////////////////////////////////////////////////////////////////
 176  // EXP()
 177  //
 178  // Computes exp(x), fr -4.6 <= x <= 0.
 179  /////////////////////////////////////////////////////////////////
 180
 181  inline float EXP (float x){
 182  assert (x <= 0.00f);
 183  if (x < EXP_UNDERFLOW_THRESHOLD) return 0.0f;
 184  return (((0.006349841068584 * x + 0.080775412572352) * x + 0.397982026296272) * x + 0.95279335963787f) * x + 0.995176455837312f;
 185  //return (((0.00681169825657f * x + 0.08386267698832f) * x + 0.40413983195844f) * x + 0.95656674979767f) * x + 0.99556744049130f;
 186  }
 187  */
 188
 189 const float EXP_UNDERFLOW_THRESHOLD = -4.6;
 190 const float LOG_UNDERFLOW_THRESHOLD = 7.5;
 191
 192 /////////////////////////////////////////////////////////////////
 193 // LOOKUP()
 194 //
 195 // Computes log (exp (x) + 1), for 0 <= x <= 7.5.
 196 /////////////////////////////////////////////////////////////////
 197
 198 inline float LOOKUP(float x) {
 199         assert(x >= 0.00f);
 200         assert(x <= LOG_UNDERFLOW_THRESHOLD);
 201         //return ((-0.00653779113685f * x + 0.09537236626558f) * x + 0.55317574459331f) * x + 0.68672959851568f;
 202         if (x <= 1.00f)
 203                 return ((-0.009350833524763f * x + 0.130659527668286f) * x
 204                                 + 0.498799810682272f) * x + 0.693203116424741f;
 205         if (x <= 2.50f)
 206                 return ((-0.014532321752540f * x + 0.139942324101744f) * x
 207                                 + 0.495635523139337f) * x + 0.692140569840976f;
 208         if (x <= 4.50f)
 209                 return ((-0.004605031767994f * x + 0.063427417320019f) * x
 210                                 + 0.695956496475118f) * x + 0.514272634594009f;
 211         assert(x <= LOG_UNDERFLOW_THRESHOLD);
 212         return ((-0.000458661602210f * x + 0.009695946122598f) * x
 213                         + 0.930734667215156f) * x + 0.168037164329057f;
 214
 215         //return (((0.00089738532761f * x - 0.01859488697982f) * x + 0.14415772028626f) * x + 0.49515490689159f) * x + 0.69311928966454f;
 216 }
 217
 218 /////////////////////////////////////////////////////////////////
 219 // LOOKUP_SLOW()
 220 //
 221 // Computes log (exp (x) + 1).
 222 /////////////////////////////////////////////////////////////////
 223
 224 inline float LOOKUP_SLOW(float x) {
 225         return log(exp(x) + 1);
 226 }
 227
 228 /////////////////////////////////////////////////////////////////
 229 // MAX()
 230 //
 231 // Compute max of three numbers
 232 /////////////////////////////////////////////////////////////////
 233
 234 inline float MAX(float x, float y, float z) {
 235         if (x >= y) {
 236                 if (x >= z)
 237                         return x;
 238                 return z;
 239         }
 240         if (y >= z)
 241                 return y;
 242         return z;
 243 }
 244
 245 /////////////////////////////////////////////////////////////////
 246 // LOG_PLUS_EQUALS()
 247 //
 248 // Add two log probabilities and store in the first argument
 249 /////////////////////////////////////////////////////////////////
 250
 251 inline void LOG_PLUS_EQUALS(float &x, float y) {
 252         if (x < y)
 253                 x = (x == LOG_ZERO || y - x >= LOG_UNDERFLOW_THRESHOLD) ?
 254                                 y : LOOKUP(y - x) + x;
 255         else
 256                 x = (y == LOG_ZERO || x - y >= LOG_UNDERFLOW_THRESHOLD) ?
 257                                 x : LOOKUP(x - y) + y;
 258 }
 259
 260 /////////////////////////////////////////////////////////////////
 261 // LOG_PLUS_EQUALS_SLOW()
 262 //
 263 // Add two log probabilities and store in the first argument
 264 /////////////////////////////////////////////////////////////////
 265
 266 inline void LOG_PLUS_EQUALS_SLOW(float &x, float y) {
 267         if (x < y)
 268                 x = (x == LOG_ZERO) ? y : LOOKUP_SLOW(y - x) + x;
 269         else
 270                 x = (y == LOG_ZERO) ? x : LOOKUP_SLOW(x - y) + y;
 271 }
 272
 273 /////////////////////////////////////////////////////////////////
 274 // LOG_ADD()
 275 //
 276 // Add two log probabilities
 277 /////////////////////////////////////////////////////////////////
 278
 279 inline float LOG_ADD(float x, float y) {
 280         if (x < y)
 281                 return (x == LOG_ZERO || y - x >= LOG_UNDERFLOW_THRESHOLD) ?
 282                                 y : LOOKUP(y - x) + x;
 283         return (y == LOG_ZERO || x - y >= LOG_UNDERFLOW_THRESHOLD) ?
 284                         x : LOOKUP(x - y) + y;
 285 }
 286
 287 /////////////////////////////////////////////////////////////////
 288 // LOG_ADD()
 289 //
 290 // Add three log probabilities
 291 /////////////////////////////////////////////////////////////////
 292
 293 inline float LOG_ADD(float x1, float x2, float x3) {
 294         return LOG_ADD(x1, LOG_ADD(x2, x3));
 295 }
 296
 297 /////////////////////////////////////////////////////////////////
 298 // LOG_ADD()
 299 //
 300 // Add four log probabilities
 301 /////////////////////////////////////////////////////////////////
 302
 303 inline float LOG_ADD(float x1, float x2, float x3, float x4) {
 304         return LOG_ADD(x1, LOG_ADD(x2, LOG_ADD(x3, x4)));
 305 }
 306
 307 /////////////////////////////////////////////////////////////////
 308 // LOG_ADD()
 309 //
 310 // Add five log probabilities
 311 /////////////////////////////////////////////////////////////////
 312
 313 inline float LOG_ADD(float x1, float x2, float x3, float x4, float x5) {
 314         return LOG_ADD(x1, LOG_ADD(x2, LOG_ADD(x3, LOG_ADD(x4, x5))));
 315 }
 316
 317 /////////////////////////////////////////////////////////////////
 318 // LOG_ADD()
 319 //
 320 // Add siz log probabilities
 321 /////////////////////////////////////////////////////////////////
 322
 323 inline float LOG_ADD(float x1, float x2, float x3, float x4, float x5,
 324                 float x6) {
 325         return LOG_ADD(x1, LOG_ADD(x2, LOG_ADD(x3, LOG_ADD(x4, LOG_ADD(x5, x6)))));
 326 }
 327
 328 /////////////////////////////////////////////////////////////////
 329 // LOG_ADD()
 330 //
 331 // Add seven log probabilities
 332 /////////////////////////////////////////////////////////////////
 333
 334 inline float LOG_ADD(float x1, float x2, float x3, float x4, float x5, float x6,
 335                 float x7) {
 336         return LOG_ADD(x1,
 337                         LOG_ADD(x2, LOG_ADD(x3, LOG_ADD(x4, LOG_ADD(x5, LOG_ADD(x6, x7))))));
 338 }
 339
 340 /////////////////////////////////////////////////////////////////
 341 // ChooseBestOfThree()
 342 //
 343 // Store the largest of three values x1, x2, and x3 in *x.  Also
 344 // if xi is the largest value, then store bi in *b.
 345 /////////////////////////////////////////////////////////////////
 346
 347 inline void ChooseBestOfThree(float x1, float x2, float x3, char b1, char b2,
 348                 char b3, float *x, char *b) {
 349         if (x1 >= x2) {
 350                 if (x1 >= x3) {
 351                         *x = x1;
 352                         *b = b1;
 353                         return;
 354                 }
 355                 *x = x3;
 356                 *b = b3;
 357                 return;
 358         }
 359         if (x2 >= x3) {
 360                 *x = x2;
 361                 *b = b2;
 362                 return;
 363         }
 364         *x = x3;
 365         *b = b3;
 366 }
 367
 368 #endif