infer.cc

// $Id$

#include <iostream>
#include <cstdlib>
#include <cmath>

#include "infer.h"

extern "C" { 
#include "libprob.h"
};

#define POSIX_SOURCE 1

#if !defined(M_PI)
#define M_PI (4 * atan(1.0))
#endif

template <class T> T sqr(const T& value) { return value * value; }

Infer::Infer(float udLh, float ut, unsigned int uV) {
  // the following bounds are checked to ensure that the exponent
  //  does not underflow, which is assumed to occur for results
  //  of less than 1e-37  => abs(t)<13.0

  // assign to copies
  dLh = udLh;
  t = ut;
  V = uV;
  if (V<=0.0) V=1.0;

  if (fabs(t)<13.0) {
    Em_ = V * pow(2*M_PI,-2) * dLh * (sqr(t) - 1) * exp(-sqr(t)/2.0);
  } else {
    Em_ = 0.0;  // underflowed exp()
  }

  if (fabs(t)<8.0) {
    B_ = pow((gamma(2.5)*Em_)/(V*(0.5 + 0.5*erf(-t/sqrt(2)))),(2.0/3.0));
  } else {
    // the large t approximation
    float a1 = V * dLh / (4 * M_PI * M_PI);
    float a3 = pow((gamma(2.5)  / V ),(2.0/3.0));
    float tsq = t*t;
    float Em_q = a1 * (2.0 * M_PI) * (tsq - 1.0)*t 
                     / ( 1.0 - 1.0/tsq + 3.0/(tsq*tsq)) ;
    B_ = a3 * pow(Em_q,(2.0/3.0));
  }
  

//    cout << "E{m} " << Em_ << endl;
//    cout << "Beta = " << B_ << endl;
}
  
float Infer::operator() (unsigned int k) {
  // ideally returns the following:
  //    return 1 - exp(-Em_ * exp(-B_ * pow( k , 2.0 / 3.0)));
  // but in practice must be careful about ranges
  // Assumes that exp(+/-87) => 1e+/38 is OK for floats
  float arg1 = -B_ * pow(k , 2.0 / 3.0);
  if (fabs(arg1)<87.0) {
    float exp1 = exp(arg1);
    float arg2 = -Em_ * exp1;
    if (fabs(arg2)<87.0) {
      float exp2 = exp(arg2);
      return 1.0 - exp2;
    }
  }
  // now must invoke the approximation instead
    float a1 = V * dLh / (4 * M_PI * M_PI);
    float a3 = pow((gamma(2.5)  / V ),(2.0/3.0));
    float tsq = t*t;
    float c = 1.0 / ( 1.0 - 1.0/tsq + 3.0/(tsq*tsq));
    float arg = -tsq*(0.5 + a3*pow(a1,2.0/3.0)*pow(2.0*M_PI,1.0/3.0)*c*pow(k,2.0/3.0));
    float exp1 = 0.0;
    if (arg>-87.0) {
      exp1 = exp(arg);
    }
    float p = a1*(tsq-1.0)*exp1;
    return p;
}