binaries/src/mafft/extensions/mxscarna_src/probconsRNA/ScoreType.h

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