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 
00043 
00044 namespace itpp {
00045   
00049   typedef signed int QLLR;
00050 
00054   typedef Vec<QLLR> QLLRvec;
00055 
00059   typedef Mat<QLLR> QLLRmat;
00060   
00064   const QLLR QLLR_MAX=(INT_MAX>>2) ;
00065  // added some margin to make sure the sum of two LLR is still permissible
00066 
00116   class LLR_calc_unit {
00117   public:
00119     LLR_calc_unit();
00120 
00125     LLR_calc_unit(short int Dint1, short int Dint2, short int Dint3);
00126 
00143     void init_llr_tables(short int Dint1 = 12, short int Dint2 = 300, 
00144                          short int Dint3 = 7);
00145     // void init_llr_tables(const short int d1=10, const short int d2=300, const short int d3=5);
00146 
00148     inline QLLR to_qllr(const double &l);
00149 
00151     QLLRvec to_qllr(const vec &l);
00152 
00154     QLLRmat to_qllr(const mat &l);
00155   
00157     inline double to_double(const QLLR &l) const { return ( ((double) l) / ((double) (1<<Dint1))); };
00158 
00160     vec to_double(const QLLRvec &l);
00161 
00163     mat to_double(const QLLRmat &l);
00164 
00169     inline QLLR jaclog(QLLR a, QLLR b);
00170     // Note: a version of this function taking "double" values as input
00171     // is deliberately omitted, because this is rather slow. 
00172 
00179     inline QLLR Boxplus(QLLR a, QLLR b);
00180 
00184     QLLR logexp(QLLR x);
00185 
00187     ivec get_Dint();
00188   
00190     void operator=(const LLR_calc_unit&);
00191 
00193     friend std::ostream &operator<<(std::ostream &os, const LLR_calc_unit &lcu);
00194     
00195   private:
00197     ivec construct_logexp_table();
00198     
00200     ivec logexp_table;
00201     
00203     short int Dint1, Dint2, Dint3;
00204 
00205   };
00206 
00210   std::ostream &operator<<(std::ostream &os, const LLR_calc_unit &lcu);
00211 
00212   
00213 
00214   // ================ IMPLEMENTATIONS OF SOME LIKELIHOOD ALGEBRA FUNCTIONS =========== 
00215 
00216   inline QLLR LLR_calc_unit::to_qllr(const double &l)
00217     { 
00218       it_assert0(l<=to_double(QLLR_MAX),"LLR_calc_unit::to_qllr(): overflow");
00219       it_assert0(l>=-to_double(QLLR_MAX),"LLR_calc_unit::to_qllr(): overflow");
00220       return ( (int) std::floor(0.5+(1<<Dint1)*l) ); 
00221     };
00222   
00223   inline QLLR LLR_calc_unit::Boxplus(QLLR a, QLLR b)
00224     {
00225       /* These boundary checks are meant to completely eliminate the
00226          possibility of any numerical problem.  Perhaps they are not
00227          strictly necessary.  */
00228       if (a>=QLLR_MAX) { 
00229         return b;
00230       }
00231       if (b>=QLLR_MAX) {
00232         return a;
00233       }
00234       if (a<=-QLLR_MAX) {
00235         return (-b);
00236       }
00237       if (b<=-QLLR_MAX) {
00238         return (-a);
00239       }
00240 
00241       QLLR a_abs = (a>0 ? a : -a);
00242       QLLR b_abs = (b>0 ? b : -b);
00243       QLLR minabs = (a_abs>b_abs ? b_abs : a_abs);
00244       QLLR term1 = (a>0 ? (b>0 ? minabs : -minabs) : (b>0 ? -minabs : minabs));
00245       QLLR apb = a+b;
00246       QLLR term2 = logexp((apb>0 ? apb : -apb));
00247       QLLR amb = a-b;
00248       QLLR term3 = logexp((amb>0 ? amb : -amb));
00249 
00250       QLLR result = term1+term2-term3;
00251       if (result>=QLLR_MAX) {         
00252         return QLLR_MAX; 
00253       } else if (result<=-QLLR_MAX) {
00254         return (-QLLR_MAX);
00255       } else {
00256         return result;
00257       }
00258     }
00259 
00260   inline QLLR LLR_calc_unit::logexp(QLLR x)
00261     {
00262       it_assert0(x>=0,"LLR_calc_unit::logexp() is not defined for negative LLR values");
00263       int ind = x>>Dint3;
00264       if (ind>=Dint2) {
00265         return 0;
00266       }
00267 
00268       it_assert0(ind>=0,"LLR_calc_unit::logexp() internal error");
00269       it_assert0(ind<Dint2,"LLR_calc_unit::logexp() internal error");
00270 
00271       // With interpolation 
00272       // int delta=x-(ind<<Dint3);
00273       // return ((delta*logexp_table(ind+1) + ((1<<Dint3)-delta)*logexp_table(ind)) >> Dint3);
00274 
00275       // Without interpolation
00276       return logexp_table(ind);
00277     }
00278 
00279   inline QLLR LLR_calc_unit::jaclog(QLLR a, QLLR b )
00280     {
00281       QLLR x,maxab;
00282   
00283       if (a>b) {
00284         maxab = a;
00285         x=a-b;
00286       } else {
00287         maxab = b;
00288         x=b-a;
00289       }
00290  
00291       if (maxab>=QLLR_MAX) {
00292         return QLLR_MAX; 
00293       } else {
00294         return (maxab + logexp(x));
00295       }
00296     };
00297 
00298 }
00299 
00300 #endif
SourceForge Logo

Generated on Wed Apr 18 11:23:32 2007 for IT++ by Doxygen 1.5.2