IT++ Logo

llr.h

Go to the documentation of this file.
00001 
00030 #ifndef LLR_H
00031 #define LLR_H
00032 
00033 #include <limits>
00034 #include <itpp/base/vec.h>
00035 #include <itpp/base/mat.h>
00036 #include <itpp/base/specmat.h>
00037 #include <itpp/base/matfunc.h>
00038 #include <limits>
00039 
00040 namespace itpp {
00041 
00045   typedef signed int QLLR;
00046 
00050   typedef Vec<QLLR> QLLRvec;
00051 
00055   typedef Mat<QLLR> QLLRmat;
00056 
00060   const QLLR QLLR_MAX = (std::numeric_limits<QLLR>::max() >> 4);
00061   // added some margin to make sure the sum of two LLR is still permissible
00062 
00112   class LLR_calc_unit {
00113   public:
00115     LLR_calc_unit();
00116 
00121     LLR_calc_unit(short int Dint1, short int Dint2, short int Dint3);
00122 
00150     void init_llr_tables(short int Dint1 = 12, short int Dint2 = 300,
00151        short int Dint3 = 7);
00152     // void init_llr_tables(const short int d1=10, const short int d2=300,
00153     //            const short int d3=5);
00154 
00156     inline QLLR to_qllr(const double &l) const;
00157 
00159     QLLRvec to_qllr(const vec &l) const;
00160 
00162     QLLRmat to_qllr(const mat &l) const;
00163 
00165     inline double to_double(const QLLR &l) const { return ( ((double) l) / ((double) (1<<Dint1))); };
00166 
00168     vec to_double(const QLLRvec &l) const;
00169 
00171     mat to_double(const QLLRmat &l) const;
00172 
00177     inline QLLR jaclog(QLLR a, QLLR b) const;
00178     // Note: a version of this function taking "double" values as input
00179     // is deliberately omitted, because this is rather slow.
00180 
00187     inline QLLR Boxplus(QLLR a, QLLR b) const;
00188 
00192     inline QLLR logexp(QLLR x) const;
00193 
00195     ivec get_Dint();
00196 
00198     void operator=(const LLR_calc_unit&);
00199 
00201     friend std::ostream &operator<<(std::ostream &os, const LLR_calc_unit &lcu);
00202 
00203   private:
00205     ivec construct_logexp_table();
00206 
00208     ivec logexp_table;
00209 
00211     short int Dint1, Dint2, Dint3;
00212 
00213   };
00214 
00219   std::ostream &operator<<(std::ostream &os, const LLR_calc_unit &lcu);
00220 
00221   // ================ IMPLEMENTATIONS OF SOME LIKELIHOOD ALGEBRA FUNCTIONS ===========
00222 
00223   inline QLLR LLR_calc_unit::to_qllr(const double &l) const
00224     {
00225       it_assert_debug(l<=to_double(QLLR_MAX),"LLR_calc_unit::to_qllr(): overflow");
00226       it_assert_debug(l>=-to_double(QLLR_MAX),"LLR_calc_unit::to_qllr(): overflow");
00227       return ( (int) std::floor(0.5+(1<<Dint1)*l) );
00228     };
00229 
00230   inline QLLR LLR_calc_unit::Boxplus(QLLR a, QLLR b) const
00231     {
00232       it_assert_debug(a<QLLR_MAX,"LLR_calc_unit::Boxplus(): owerflow");
00233       it_assert_debug(b<QLLR_MAX,"LLR_calc_unit::Boxplus(): owerflow");
00234       it_assert_debug(a>-QLLR_MAX,"LLR_calc_unit::Boxplus(): owerflow");
00235       it_assert_debug(b>-QLLR_MAX,"LLR_calc_unit::Boxplus(): owerflow");
00236 
00237       QLLR a_abs = (a>0 ? a : -a);
00238       QLLR b_abs = (b>0 ? b : -b);
00239       QLLR minabs = (a_abs>b_abs ? b_abs : a_abs);
00240       QLLR term1 = (a>0 ? (b>0 ? minabs : -minabs) : (b>0 ? -minabs : minabs));
00241 
00242       if (Dint2==0) {  // logmax approximation - avoid looking into empty table
00243   it_assert_debug(term1<QLLR_MAX,"LLR_calc_unit::Boxplus(): owerflow");
00244   it_assert_debug(term1>-QLLR_MAX,"LLR_calc_unit::Boxplus(): owerflow");
00245   return term1;
00246       }
00247 
00248       QLLR apb = a+b;
00249       QLLR term2 = logexp((apb>0 ? apb : -apb));
00250       QLLR amb = a-b;
00251       QLLR term3 = logexp((amb>0 ? amb : -amb));
00252 
00253       QLLR result = term1+term2-term3;
00254 
00255       it_assert_debug(result<QLLR_MAX,"LLR_calc_unit::Boxplus(): owerflow");
00256       it_assert_debug(result>-QLLR_MAX,"LLR_calc_unit::Boxplus(): owerflow");
00257       return result;
00258     }
00259 
00260   inline QLLR LLR_calc_unit::logexp(QLLR x) const
00261     {
00262       it_assert_debug(x>=0,"LLR_calc_unit::logexp() is not defined for negative LLR values");
00263       int ind = x>>Dint3;
00264       if (ind>=Dint2) {  // outside table
00265   return 0;
00266       }
00267 
00268       it_assert_debug(ind>=0,"LLR_calc_unit::logexp() internal error");
00269       it_assert_debug(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 ) const
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 Sun Dec 9 17:26:18 2007 for IT++ by Doxygen 1.5.4