binaries/src/mafft/extensions/mxscarna_src/ScoreType.hpp

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