00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00028 #ifndef SEQPP_MARKOV_H
00029 #define SEQPP_MARKOV_H
00030
00031 #include <seqpp/PhasedMarkov.h>
00032 #include <seqpp/arnoldi.h>
00033
00043 class Markov : public PhasedMarkov
00044 {
00045 public :
00046
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063 Markov( const char * ConfFile,
00064 bool calc_rank = false );
00066
00071 Markov( const SequenceSet & seqset,
00072 bool calc_rank = false,
00073 const string& prior_alpha_file = string() );
00074
00076
00081 Markov( const Sequence & seq,
00082 bool calc_rank = false,
00083 const string& prior_alpha_file = string() );
00084
00086 Markov( const Markov & );
00087
00089 Markov();
00090
00092
00099 Markov( short size, short order, bool alloc=true,
00100 const string& prior_alpha_file = string() ):PhasedMarkov(size,order,1,alloc, prior_alpha_file)
00101 {
00102 if (alloc){
00103 _Pi=_Pis[0];_container=_containers[0];
00104 }
00105 else{
00106 _Pi = NULL; _container=NULL;
00107 }
00108 _Mu=NULL;_PowPi=NULL;
00109 }
00110
00112
00117 Markov(const Markov &M1, const Markov &M2,const float p);
00118
00120
00138 Markov( const gsl_rng * r,
00139 short size, short order,
00140 bool calc_rank = false)
00141 : PhasedMarkov( r, size, order, 1, calc_rank ){
00142 _Pi=_Pis[0];_container=_containers[0];_Mu=NULL;_PowPi=NULL;
00143 };
00144
00146
00153 Markov( unsigned long * count,
00154 short size, short order,
00155 const string& prior_alpha_file = string(),
00156 bool calc_rank = false );
00157
00159 virtual Markov::~Markov();
00160
00161
00163
00169 template <class TSeq>
00170 void estimate( const TSeq & tseq,
00171 unsigned long beg, unsigned long end,
00172 bool calc_rank ) {
00173 this->PhasedMarkov::estimate( tseq, 1,0, beg, end, calc_rank );
00174 }
00175
00177
00182 void estimate( unsigned long * count, bool decal_required, bool calc_rank= false ){
00183 this->PhasedMarkov::estimate( &count, decal_required, calc_rank );
00184 }
00185
00187
00191 void estimate( const string& count_file, bool calc_rank= false )
00192 {
00193 this->PhasedMarkov::estimate( count_file, calc_rank );
00194 }
00195
00196
00197
00199 const double * markov_matrix() const{
00200 return PhasedMarkov::markov_matrix(0);
00201 }
00202
00204
00223 void draw_markov_matrix(const gsl_rng * r){
00224 PhasedMarkov::draw_markov_matrices(r);
00225 }
00226
00228 void free_markov_matrix(){
00229 PhasedMarkov::free_markov_matrices();
00230 }
00232
00233
00234
00235 void compute_stat_law( bool force=false ){
00236 PhasedMarkov::compute_stat_laws( force );
00237 _Mu = _Mus[0];
00238 }
00240 void free_stat_law(){
00241 PhasedMarkov::free_stat_laws();
00242 _Mu = NULL;
00243 }
00245 const double * stat_law() const{
00246 return PhasedMarkov::stat_law(0);
00247 }
00248
00250 virtual int compute_rank();
00251
00253 void compute_power();
00254
00256 int free_power();
00257
00259
00265
00267
00272
00273
00275
00281 double proba_step( long w1, long w2, int step ){
00282 compute_rank();
00283 compute_power();
00284 if (step<=_rank)
00285 return(_PowPi[step-1][w1][w2]);
00286 else
00287 return(_Mu[w2]);
00288 }
00289
00290
00292 bool isPi() const { return(_Pi != NULL);};
00294 bool isPow() const { return(_PowPi != NULL);};
00296 bool isMu() const { return(_Mu != NULL);};
00297
00299
00302 double & operator() (int index){
00303 return(_Pi[index]);
00304 }
00306
00309 double Mu(int index) const {
00310 return(_Mu[index]);
00311 }
00312
00313 protected :
00314
00316 double *_Pi;
00318 double *_container;
00319
00321 double *_Mu;
00322
00324 double ***_PowPi;
00325 };
00326 #endif
00327