ExpEvaluator.h

Go to the documentation of this file.

#ifndef EXPEVALUATOR_H_
#define EXPEVALUATOR_H_

#ifdef __cplusplus
#define _USE_MATH_DEFINES
#include "log.h"
#include "Quadrature.h"
#include <malloc.h>
#include <math.h>
#include "exponentials.h"
#endif


class ExpEvaluator {

private:

  bool _interpolate;

  bool _linear_source;

  FP_PRECISION _exp_table_spacing;

  FP_PRECISION _inverse_exp_table_spacing;

  FP_PRECISION _sin_theta_no_offset;

  FP_PRECISION _inverse_sin_theta_no_offset;

  int _table_size;

  FP_PRECISION* _exp_table;

  Quadrature* _quadrature;

  FP_PRECISION _max_optical_length;

  FP_PRECISION _exp_precision;

  int _azim_index;

  int _polar_index;

  int _num_exp_terms;

  int _num_polar_terms;


public:

  ExpEvaluator();
  virtual ~ExpEvaluator();

  void setQuadrature(Quadrature* quadrature);
  void setMaxOpticalLength(FP_PRECISION max_optical_length);
  void setExpPrecision(FP_PRECISION exp_precision);
  void useInterpolation();
  void useIntrinsic();
  void useLinearSource();

  FP_PRECISION getMaxOpticalLength();
  FP_PRECISION getExpPrecision();
  bool isUsingInterpolation();
  FP_PRECISION getTableSpacing();
  int getTableSize();
  FP_PRECISION* getExpTable();
  int getExponentialIndex(FP_PRECISION tau);
  FP_PRECISION getDifference(int index, FP_PRECISION tau);
  FP_PRECISION convertDistance3Dto2D(FP_PRECISION length);
  FP_PRECISION getInverseSinTheta();

  void initialize(int azim_index, int polar_index, bool solve_3D);
  FP_PRECISION computeExponential(FP_PRECISION tau, int polar_offset);
  FP_PRECISION computeExponentialF1(int index, int polar_offset,
                                    FP_PRECISION dt, FP_PRECISION dt2);
  FP_PRECISION computeExponentialF2(int index, int polar_offset,
                                    FP_PRECISION dt, FP_PRECISION dt2);
  FP_PRECISION computeExponentialH(int index, int polar_offset,
                                   FP_PRECISION dt, FP_PRECISION dt2);
  void retrieveExponentialComponents(FP_PRECISION tau, int polar_offset,
                                     FP_PRECISION* exp_F1,
                                     FP_PRECISION* exp_F2,
                                     FP_PRECISION* exp_H);

  FP_PRECISION computeExponentialG2(FP_PRECISION tau);
  ExpEvaluator* deepCopy();
};


inline int ExpEvaluator::getExponentialIndex(FP_PRECISION tau) {
  return int(tau * _inverse_exp_table_spacing);
}


inline FP_PRECISION ExpEvaluator::getDifference(int index, FP_PRECISION tau) {
  return tau - index * _exp_table_spacing;
}


inline FP_PRECISION ExpEvaluator::convertDistance3Dto2D(FP_PRECISION length) {
  return length * _sin_theta_no_offset;
}


inline FP_PRECISION ExpEvaluator::getInverseSinTheta() {
  return _inverse_sin_theta_no_offset;
}


inline FP_PRECISION ExpEvaluator::computeExponential(FP_PRECISION tau,
                                                     int polar_offset) {

#ifndef THREED
  FP_PRECISION inv_sin_theta = 1.f / _quadrature->getSinThetaInline(_azim_index,
                                                   _polar_index + polar_offset);
#else
  FP_PRECISION inv_sin_theta = _inverse_sin_theta_no_offset;
#endif
  FP_PRECISION exp_F1;
  expF1_fractional(tau * inv_sin_theta, &exp_F1);

  return inv_sin_theta * exp_F1;
}


inline FP_PRECISION ExpEvaluator::computeExponentialF1(int index,
                                                       int polar_offset,
                                                       FP_PRECISION dt,
                                                       FP_PRECISION dt2) {

  /* Calculate full index */
  int full_index = (index * _num_polar_terms + polar_offset) * _num_exp_terms;

  if (_interpolate) {
    return _exp_table[full_index] + _exp_table[full_index + 1] * dt +
        _exp_table[full_index + 2] * dt2;
  }
  else {
    int polar_index = _polar_index + polar_offset;
    FP_PRECISION tau = index * _exp_table_spacing + dt;
    FP_PRECISION inv_sin_theta = 1.0 / _quadrature->getSinTheta(_azim_index,
                                                                polar_index);
    return (1.0 - exp(- tau * inv_sin_theta)) / tau;
  }
}


inline FP_PRECISION ExpEvaluator::computeExponentialF2(int index,
                                                       int polar_offset,
                                                       FP_PRECISION dt,
                                                       FP_PRECISION dt2) {
  /* Calculate full index */
  int full_index = (index * _num_polar_terms + polar_offset) * _num_exp_terms;

  if (_interpolate)
    return _exp_table[full_index + 3] + _exp_table[full_index + 4] * dt +
        _exp_table[full_index + 5] * dt2;
  else {

    int polar_index = _polar_index + polar_offset;
    FP_PRECISION tau = index * _exp_table_spacing + dt;
    FP_PRECISION inv_sin_theta = 1.0 / _quadrature->getSinTheta(_azim_index,
                                                                polar_index);
    FP_PRECISION tau_m = tau * inv_sin_theta;
    FP_PRECISION F1 = (1.0 - exp(- tau_m)) / tau;
    return 2.0 / tau * (inv_sin_theta - F1) - inv_sin_theta * F1;
  }
}


inline FP_PRECISION ExpEvaluator::computeExponentialH(int index,
                                                      int polar_offset,
                                                      FP_PRECISION dt,
                                                      FP_PRECISION dt2) {
  /* Calculate full index */
  int full_index = (index * _num_polar_terms + polar_offset) * _num_exp_terms;

  if (_interpolate)
    return _exp_table[full_index + 6] + _exp_table[full_index + 7] * dt +
        _exp_table[full_index + 8] * dt2;
  else {
    int polar_index = _polar_index + polar_offset;
    FP_PRECISION tau = index * _exp_table_spacing + dt;
    FP_PRECISION inv_sin_theta = 1.0 / _quadrature->getSinTheta(_azim_index,
                                                                polar_index);
    FP_PRECISION tau_m = tau * inv_sin_theta;
    FP_PRECISION F1 = (1.0 - exp(- tau_m)) / tau;
    FP_PRECISION G1 = 1.0 / tau + 0.5 * inv_sin_theta - (1.0 + 1.0 / tau_m) * F1;
    return 0.5 * inv_sin_theta - G1;
  }
}


inline FP_PRECISION ExpEvaluator::computeExponentialG2(FP_PRECISION tau) {

  FP_PRECISION exp_G2;
  expG2_fractional(tau, &exp_G2);

  return exp_G2;
}


inline void ExpEvaluator::retrieveExponentialComponents(FP_PRECISION tau,
                                                        int polar_offset,
#ifdef SWIG  //FIXME Find out how to use restrict with SWIG
                                                        FP_PRECISION* exp_F1,
                                                        FP_PRECISION* exp_F2,
                                                        FP_PRECISION* exp_H) {
#else
                                              FP_PRECISION* __restrict__ exp_F1,
                                              FP_PRECISION* __restrict__ exp_F2,
                                              FP_PRECISION* __restrict__ exp_H) {
#endif

#ifndef THREED
  FP_PRECISION inv_sin_theta = 1.f / _quadrature->getSinThetaInline(_azim_index,
                                                   _polar_index + polar_offset);
#else
  FP_PRECISION inv_sin_theta = _inverse_sin_theta_no_offset;
#endif

  /* Limit range of tau to avoid numerical errors */
  tau = std::max(FP_PRECISION(1e-8), tau * inv_sin_theta);

  /* Compute exponentials from a common exponential */
   FP_PRECISION exp_G;
   expG_fractional(tau, &exp_G);
   *exp_F1 = 1.f - tau*exp_G;
   *exp_F1 *= inv_sin_theta;

   exp_G *= inv_sin_theta;
   *exp_F2 = 2.f*exp_G - *exp_F1;
   *exp_F2 *= inv_sin_theta;

   *exp_H = *exp_F1 - exp_G;

    /* Quadratic exponential interpolation tables */
//     __builtin_assume_aligned(_exp_table, VEC_ALIGNMENT);
// 
//     tau /= inv_sin_theta;
//     int exp_index = getExponentialIndex(tau);
//     FP_PRECISION dt = getDifference(exp_index, tau);
//     FP_PRECISION dt2 = dt * dt;
//     int full_index = (exp_index * _num_polar_terms + polar_offset)
//       * _num_exp_terms;
//     *exp_F1 = _exp_table[full_index] + _exp_table[full_index + 1] * dt +
//         _exp_table[full_index + 2] * dt2;
//     *exp_F2 = _exp_table[full_index + 3] + _exp_table[full_index + 4] * dt +
//         _exp_table[full_index + 5] * dt2;
//     *exp_H = _exp_table[full_index + 6] + _exp_table[full_index + 7] * dt +
//         _exp_table[full_index + 8] * dt2;
}


#endif /* EXPEVALUATOR_H_ */