00001 #ifndef SEQPP_MTDCORE
00002 #define SEQPP_MTDCORE
00003
00004 #include <vector>
00005 #include <seqpp/Markov.h>
00006
00007 class mtd_core
00008 {
00009 public:
00010 mtd_core( short alphabet_size, short mtd_order, short mkv_order );
00011 ~mtd_core();
00012
00013
00014 void estimate( unsigned long * count, bool decal_required,
00015 short nbseed, int nbiter_max, double eps, bool log );
00016
00021 double E_step();
00022
00023 void M_step();
00024
00025 void mtd_to_matrix( double * Pi_markov );
00026
00027 int nb_params() const{
00028 return ((_npi)*_pow[_mkv_order]*(_alphabet_size-1) + _npi-1);
00029 }
00030
00031
00032 private:
00034 short _mtd_order;
00036 short _alphabet_size;
00037
00039 int _nbiter_max;
00040 double _eps;
00041 short _nbseed;
00042
00044 short _mkv_order;
00045
00047 short _npi;
00048 Markov** _pi;
00049 vector< double > _phi;
00051 void init_param( Markov** &pi, vector< double > &phi, const gsl_rng * r );
00053 void free_param( Markov** &pi, vector< double > &phi ){
00054 if (pi != NULL){
00055 for (short i=0; i<_npi; i++)
00056 delete pi[i];
00057 delete[] pi;
00058 pi = NULL;
00059 }
00060 }
00061
00062
00064 vector<unsigned long> _count;
00066 long _count_sum;
00067 void init_count( unsigned long * count, bool decal_required );
00068
00069
00070 vector<long> _pow;
00071 void init_pow(){
00072 _pow.push_back(1);
00073 for (short i=1; i<=(_mtd_order+1); i++)
00074 _pow.push_back( _pow[i-1]*_alphabet_size );
00075 }
00076
00077
00078 double ** _pcond;
00079 void init_pcond(){
00080 _pcond = new double* [_npi];
00081 for (short itpos=0; itpos<_npi; itpos++)
00082 _pcond[itpos] = new double[ _pow[_mtd_order+1] ];
00083 }
00084 void free_pcond(){
00085 for (short itpos=0; itpos<_npi; itpos++)
00086 delete[] _pcond[itpos];
00087 delete[] _pcond;
00088 }
00089
00090
00091 vector<long> _suiv;
00092 void init_suiv(){
00093 for (long i=0; i<(_pow[_mkv_order]-1);i++)
00094 _suiv.push_back( i+1 );
00095 _suiv.push_back( 0 );
00096 }
00097 void free_suiv(){
00098 _suiv.clear();
00099 }
00100
00101 vector<long> _predict_path;
00102 double ** _sum;
00103 void init_pile(){
00104 _predict_path.assign(_npi, 0);
00105 if (_sum == NULL){
00106 _sum = new double*[_npi];
00107 for (short itpos=0; itpos<_npi; itpos++){
00108 _sum[itpos] = new double[_alphabet_size];
00109 }
00110 }
00111 }
00112 void free_pile(){
00113 _predict_path.clear();
00114 for (short i=0; i<_npi; i++)
00115 delete[] _sum[i];
00116 delete[] _sum;
00117 _sum = NULL;
00118 }
00119
00120
00121 double numerator_M_step( short itpos, long predict, short i);
00122 };
00123 #endif //SEQPP_MTDCORE