IT++ Logo Newcom Logo

llr.h

Go to the documentation of this file.
00001 
00035 #ifndef LLR_H
00036 #define LLR_H
00037 
00038 #include <itpp/base/vec.h>
00039 #include <itpp/base/mat.h>
00040 #include <itpp/base/specmat.h>
00041 #include <itpp/base/matfunc.h>
00042 #include <limits>
00043 
00044 
00045 namespace itpp {
00046   
00050   typedef signed int QLLR;
00051 
00055   typedef Vec<QLLR> QLLRvec;
00056 
00060   typedef Mat<QLLR> QLLRmat;
00061   
00065   const QLLR QLLR_MAX = (std::numeric_limits<QLLR>::max() >> 2);
00066   // added some margin to make sure the sum of two LLR is still permissible
00067 
00117   class LLR_calc_unit {
00118   public:
00120     LLR_calc_unit();
00121 
00126     LLR_calc_unit(short int Dint1, short int Dint2, short int Dint3);
00127 
00144     void init_llr_tables(short int Dint1 = 12, short int Dint2 = 300, 
00145                          short int Dint3 = 7);
00146     // void init_llr_tables(const short int d1=10, const short int d2=300, const short int d3=5);
00147 
00149     inline QLLR to_qllr(const double &l);
00150 
00152     QLLRvec to_qllr(const vec &l);
00153 
00155     QLLRmat to_qllr(const mat &l);
00156   
00158     inline double to_double(const QLLR &l) const { return ( ((double) l) / ((double) (1<<Dint1))); };
00159 
00161     vec to_double(const QLLRvec &l);
00162 
00164     mat to_double(const QLLRmat &l);
00165 
00170     inline QLLR jaclog(QLLR a, QLLR b);
00171     // Note: a version of this function taking "double" values as input
00172     // is deliberately omitted, because this is rather slow. 
00173 
00180     inline QLLR Boxplus(QLLR a, QLLR b);
00181 
00185     QLLR logexp(QLLR x);
00186 
00188     ivec get_Dint();
00189   
00191     void operator=(const LLR_calc_unit&);
00192 
00194     friend std::ostream &operator<<(std::ostream &os, const LLR_calc_unit &lcu);
00195     
00196   private:
00198     ivec construct_logexp_table();
00199     
00201     ivec logexp_table;
00202     
00204     short int Dint1, Dint2, Dint3;
00205 
00206   };
00207 
00211   std::ostream &operator<<(std::ostream &os, const LLR_calc_unit &lcu);
00212 
00213   
00214 
00215   // ================ IMPLEMENTATIONS OF SOME LIKELIHOOD ALGEBRA FUNCTIONS =========== 
00216 
00217   inline QLLR LLR_calc_unit::to_qllr(const double &l)
00218     { 
00219       it_assert0(l<=to_double(QLLR_MAX),"LLR_calc_unit::to_qllr(): overflow");
00220       it_assert0(l>=-to_double(QLLR_MAX),"LLR_calc_unit::to_qllr(): overflow");
00221       return ( (int) std::floor(0.5+(1<<Dint1)*l) ); 
00222     };
00223   
00224   inline QLLR LLR_calc_unit::Boxplus(QLLR a, QLLR b)
00225     {
00226       /* These boundary checks are meant to completely eliminate the
00227          possibility of any numerical problem.  Perhaps they are not
00228          strictly necessary.  */
00229       if (a>=QLLR_MAX) { 
00230         return b;
00231       }
00232       if (b>=QLLR_MAX) {
00233         return a;
00234       }
00235       if (a<=-QLLR_MAX) {
00236         return (-b);
00237       }
00238       if (b<=-QLLR_MAX) {
00239         return (-a);
00240       }
00241 
00242       QLLR a_abs = (a>0 ? a : -a);
00243       QLLR b_abs = (b>0 ? b : -b);
00244       QLLR minabs = (a_abs>b_abs ? b_abs : a_abs);
00245       QLLR term1 = (a>0 ? (b>0 ? minabs : -minabs) : (b>0 ? -minabs : minabs));
00246       QLLR apb = a+b;
00247       QLLR term2 = logexp((apb>0 ? apb : -apb));
00248       QLLR amb = a-b;
00249       QLLR term3 = logexp((amb>0 ? amb : -amb));
00250 
00251       QLLR result = term1+term2-term3;
00252       if (result>=QLLR_MAX) {         
00253         return QLLR_MAX; 
00254       } else if (result<=-QLLR_MAX) {
00255         return (-QLLR_MAX);
00256       } else {
00257         return result;
00258       }
00259     }
00260 
00261   inline QLLR LLR_calc_unit::logexp(QLLR x)
00262     {
00263       it_assert0(x>=0,"LLR_calc_unit::logexp() is not defined for negative LLR values");
00264       int ind = x>>Dint3;
00265       if (ind>=Dint2) {
00266         return 0;
00267       }
00268 
00269       it_assert0(ind>=0,"LLR_calc_unit::logexp() internal error");
00270       it_assert0(ind<Dint2,"LLR_calc_unit::logexp() internal error");
00271 
00272       // With interpolation 
00273       // int delta=x-(ind<<Dint3);
00274       // return ((delta*logexp_table(ind+1) + ((1<<Dint3)-delta)*logexp_table(ind)) >> Dint3);
00275 
00276       // Without interpolation
00277       return logexp_table(ind);
00278     }
00279 
00280   inline QLLR LLR_calc_unit::jaclog(QLLR a, QLLR b )
00281     {
00282       QLLR x,maxab;
00283   
00284       if (a>b) {
00285         maxab = a;
00286         x=a-b;
00287       } else {
00288         maxab = b;
00289         x=b-a;
00290       }
00291  
00292       if (maxab>=QLLR_MAX) {
00293         return QLLR_MAX; 
00294       } else {
00295         return (maxab + logexp(x));
00296       }
00297     };
00298 
00299 }
00300 
00301 #endif
SourceForge Logo

Generated on Sat Aug 25 23:40:55 2007 for IT++ by Doxygen 1.5.2