Version 4.2.0

Markov.h

Go to the documentation of this file.

/* seqpp/Markov.h
 *
 * Copyright (C) 2003 Laboratoire Statistique & Génome
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#ifndef SEQPP_MARKOV_H
#define SEQPP_MARKOV_H

#include <seqpp/PhasedMarkov.h>
#include <seqpp/arnoldi.h>

class Markov : public PhasedMarkov
{
 public :
 
  /*
    \param markov_file description of an input model in a file (generated by print method) 
    \param calc_rank calculus of the convergence rank if true
    Exple of file:
    \code
    # 4 <- Order of the phased Markov chain<br>
    # 3 <- Phase<br>
    # 4 <- Alphabet size<br>
    # 10 steps <- Convergence to the stationnary distribution (!!Optionnal!!)<br>
    # Transition matrix:<br>
    ...<br>
    # Stationnary Probability: (!!Optionnal!!)<<br>
    ...
    \endcode
  */
  Markov( const char * ConfFile, 
          bool calc_rank = false );


  Markov( const SequenceSet & seqset,
          bool calc_rank = false,
          const string& prior_alpha_file = string() );


  Markov( const Sequence & seq,
          bool calc_rank = false,
          const string& prior_alpha_file = string() );

  Markov( const Markov & );
  
  Markov();


  Markov( short size, short order, bool alloc=true,
          const string& prior_alpha_file = string() ):PhasedMarkov(size,order,1,alloc, prior_alpha_file) 
  {//cout<<"Markov constr. 6"<<endl;
    if (alloc){
      _Pi=_Pis[0];_container=_containers[0];
    }
    else{
      _Pi = NULL; _container=NULL;
    }
    _Mu=NULL;_PowPi=NULL;
  }


  Markov(const Markov &M1, const Markov &M2,const float p);


  Markov( const gsl_rng * r,
          short size, short order,
          bool calc_rank = false)
    : PhasedMarkov( r, size, order, 1, calc_rank ){
    _Pi=_Pis[0];_container=_containers[0];_Mu=NULL;_PowPi=NULL;
  };


  Markov( unsigned long * count,
          short size, short order,
          const string& prior_alpha_file = string(),
          bool calc_rank = false );

  virtual ~Markov();

  //---------------------------------------

  template <class TSeq>
  void estimate( const TSeq & tseq,
                 unsigned long beg, unsigned long end,
                 bool calc_rank ) {
    this->PhasedMarkov::estimate( tseq, 1,0, beg, end, calc_rank );
  }


  void estimate( unsigned long * count, bool decal_required, bool calc_rank= false ){
    this->PhasedMarkov::estimate( &count, decal_required, calc_rank );
  }


  void estimate( const string& count_file, bool calc_rank= false )
  {
    this->PhasedMarkov::estimate( count_file, calc_rank );
  }


  //---------------------------------------
  const double * markov_matrix() const{
    return PhasedMarkov::markov_matrix(0);
  } 


  void draw_markov_matrix(const  gsl_rng * r){
    PhasedMarkov::draw_markov_matrices(r);
  }

  void free_markov_matrix(){
    PhasedMarkov::free_markov_matrices();
  }
  
  /*
    \param force "true" to force the calculation, even if the law already exists. Default => "false". 
  */ 
  void compute_stat_law( bool force=false ){
    PhasedMarkov::compute_stat_laws( force );
    _Mu = _Mus[0];
  }  
  void free_stat_law(){
    PhasedMarkov::free_stat_laws();
    _Mu = NULL;
  }
  const double * stat_law() const{
    return PhasedMarkov::stat_law(0);
  }
  //--------------------------------------- 
  virtual int compute_rank(); 
  
  void compute_power();
   
  int free_power();


  // double proba_step( const string & w1, const string & w2, int step ) const;

  //double proba_step( const vector<short> & w1, const vector<short> & w2, int step ) const;


  double proba_step( long w1, long w2, int step ){
    compute_rank();
    compute_power();
      if (step<=_rank)
        return(_PowPi[step-1][w1][w2]);
      else 
        return(_Mu[w2]);
  }
   
  //---------------------------------------
  bool isPi() const { return(_Pi != NULL);};
  bool isPow() const { return(_PowPi != NULL);};
  bool isMu()  const { return(_Mu != NULL);};


  double & operator() (int index){ 
    return(_Pi[index]); 
  }

  double Mu(int index) const { 
    return(_Mu[index]); 
  }

protected :
  
  double *_Pi;
  double *_container;
  
  double *_Mu; 
  
  double ***_PowPi;
}; 
#endif/*SEQPP_MARKOV_H*/

Download seq++ 4.2.0

Download previous versions

Statistique & Genome Home

Contributors : M.Baudry, P.Y.Bourguignon, M.Hoebeke, V.Miele, P.Nicolas, G.Nuel, H.Richard, D.Robelin

---