Version 4.2.0

PhasedMarkov.h

Go to the documentation of this file.

/* seqpp/PhasedMarkov.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_PHASEDMARKOV_H
#define SEQPP_PHASEDMARKOV_H

using namespace std;
#include <seqpp/SequenceSet.h>
#include <seqpp/arnoldi.h>
#include <vector>
#include <string>

#include <gsl/gsl_sf_gamma.h>

class 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>
    # Phase n°: 0<br>
    # 10 steps <- Convergence to the stationnary distribution (!!Optionnal!!)<br>
    # Transition matrix:<br>
    ...<br>
    # Stationnary Probability: (!!Optionnal!!)<<br>
    ...
    \endcode
  */
  PhasedMarkov(  const string & markov_file,
                 bool calc_rank = false  );
 

  PhasedMarkov( const SequenceSet & seqset,
                short phase, short initial_phase = 0,
                bool calc_rank = false,
                const string& prior_alpha_file = string() );
  

  PhasedMarkov( const Sequence & seq,
                short phase, short initial_phase = 0,
                bool calc_rank = false,
                const string& prior_alpha_file = string() );

  PhasedMarkov( const PhasedMarkov & phm );

  PhasedMarkov();


  PhasedMarkov( short size, short order, short phase, bool alloc=true,
                const string& prior_alpha_file = string() );

  /*
    \param M1 first Markov chain object
    \param M2 second Markov chain object
    \param p weight of M1 in the resulting M(with 0<=p<=1)
  */
  PhasedMarkov(const PhasedMarkov &M1, const PhasedMarkov &M2, const float p);


  PhasedMarkov(const SequenceSet & seqset, 
               const vector<int> &Indseq, 
               short phase, short initial_phase = 0,
               bool calc_rank = false,
               const string& prior_alpha_file = string() ); 


  PhasedMarkov( const gsl_rng * r,
                short size, short order, 
                short phase, 
                bool calc_rank = false );


  PhasedMarkov( unsigned long * * count,
                short size, short order,
                short phase, short initial_phase = 0,
                bool calc_rank = false,
                const string& prior_alpha_file = string() );

  virtual ~PhasedMarkov();
  
  //---------------------------------------

  template <class TSeq>
    void estimate( const TSeq & tseq,
                   short phase, short initial_phase, 
                   unsigned long beg, unsigned long end,
                   bool calc_rank = false, bool count_again = true )
    {
      if (count_again)
        tseq.clear_count();  
      tseq.count_p_occurencies(_phase,initial_phase,beg,end);
      
      estimate( tseq.get_p_count(), true, calc_rank );
    } 



  void estimate( const string& count_file, bool calc_rank= false )
  {
    unsigned long ** count = new unsigned long*[_phase];
    for (short p=0; p<_phase; p++ )
      count[p] = new unsigned long[_nPi];    

    if (!load_count( count_file, count )){
      cerr<<"Bad count file"<<endl;
      exit(1);
    }
      
    estimate( count, false, calc_rank );

    for (short p=0; p<_phase; p++ )
      delete[] count[p];
    delete [] count;    
  }


  void estimate( unsigned long * * count, bool decal_required, bool calc_rank= false );


  //---------------------------------------
  const double ** markov_matrices( ) const{
    return const_cast< const double** >(_Pis);
  }
  const double * markov_matrix( short numphase ) const{
    if ( numphase<_phase )
      return _Pis[numphase];
    else 
      return NULL;
  }
  

  void draw_markov_matrices(const  gsl_rng * r);


  inline virtual void new_markov_matrices();
  inline virtual void free_markov_matrices();

  double total_variation( const PhasedMarkov & M );


  //---------------------------------------
  /*
    \param force "true" to force the calculation, even if the laws already exist. Default => "false". 
  */ 
  void compute_stat_laws( bool force = false ); 
  const double * stat_law( short numphase = 0 ) const{
    if ( numphase<_phase )
      return _Mus[numphase];
    else 
      return NULL;
  }
  inline void free_stat_laws();
  
  /*
    \param MuInit initial law, must have allocated
    \seqset set of sequences
  */
  void compute_init_law( double * MuInit, const SequenceSet & seqset ) const;


  //---------------------------------------
  virtual int compute_rank();   
  virtual long nb_parameters() const{
    return  _nb_param;
  }


  //---------------------------------------
  void link_to_translator( const Translator & trans ){
    _trans = &trans;
  }

  double proba_c( const string & word, 
                  Coder & coder, short numphase=0 ) const;

  double proba( const string & word, 
                Coder & coder, short numphase=0 ) const;



  double proba_c( const vector<short> & word, 
                  Coder & coder, short numphase=0 ) const;

  double proba( const vector<short> & word, 
                Coder & coder, short numphase=0 ) const;



  double proba_c( long word, int lw=-1, long jump=-1, short numphase=0 ) const;

  double proba( long word, int lw=-1, long jump=-1, short numphase=0 ) const;

  

  double proba_c(const long * seq, long tbeg, long tend, short numphase=0) const;


  double proba(const long * seq, long tbeg, long tend, short numphase=0) const;


  double proba_c(const Sequence& seq, long tbeg, long tend, short numphase=0) const{
    return proba_c(seq.seq(), tbeg, tend, numphase);
  }


  double proba(const Sequence& seq, long tbeg, long tend, short numphase=0) const{
    return  proba(seq.seq(), tbeg, tend, numphase);
  }
 
 
  //---------------------------------------

  double log_likelihood( const SequenceSet & seqset,
                         short initial_phase = 0, short numphase=-1 ) const; 


  double log_ratio_likelihood( const SequenceSet & seqset,
                               const PhasedMarkov &M,
                               short initial_phase1=0,                               
                               short initial_phase2=0 ) const ; 



  double log_likelihood( const Sequence & seq,
                         short initial_phase = 0, short numphase=-1 ) const; 


  double log_ratio_likelihood( const Sequence & seq,
                               const PhasedMarkov &M,
                               short initial_phase1=0,                               
                               short initial_phase2=0 ) const ;
  

  template <class TSeq>
    double BIC( const TSeq & tseq, short initial_phase=0 ) const{   
    //BIC = -2*loglikelihood + nbparam*log(length)
    return( -2*log_likelihood( tseq, initial_phase ) + _nb_param*log( (double)tseq.tell_length() ) );
  }

 

  template <class TSeq>
    double AIC( const TSeq & tseq, short initial_phase=0 ) const{
    //AIC = -2*loglikelihood + 2*nbparam
    return( -2*log_likelihood( tseq, initial_phase )  + 2*_nb_param );
  }


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

   template <class TSeq1, class TSeq2>
   double post_log_likelihood( const TSeq1& tseq_train, const TSeq2& tseq_eval,
                               bool force=false,
                               short initial_phase_train = 0, short initial_phase_eval = 0 )
   {
    double LogLik_te = 0., LogLik_t = 0.,sum1=0., sum2=0.;
    short p, u;
    long j, word;

    if (( tseq_train.tell_order() != _order )||( tseq_eval.tell_order() != _order )){
      cerr << "Order incompatibility beetween  sequence and model " <<endl;
      exit(-1);
    }

    if (_prior_alpha.size() == 0){
       cerr << "Prior input required for post_log_likelihood." <<endl;
      exit(-1);
    }
      
    // A FAIRE attention seq > order...

    tseq_train.count_p_occurencies( _phase, initial_phase_train );
    tseq_eval.count_p_occurencies(_phase, initial_phase_eval);
       
    unsigned long * countp_t, * countp_e;
    long jump = tseq_train.tell_jump();

   
    vector<double> lng1; 
    lng1.reserve(_size);
    for (p=0; p<_phase; p++){
      
      double sum = 0;
      for (u=0; u<_size; u++){
        lng1[u] = gsl_sf_lngamma(_prior_alpha[p][u]);
        sum += _prior_alpha[p][u];
      }
      double lng2 = gsl_sf_lngamma(sum);
      
      countp_t = (tseq_train.get_p_count())[p] + jump;
      countp_e = (tseq_eval.get_p_count())[p] + jump;

        for(j = 0; j < _nMu; j++){      
          sum1=0., sum2=0.;
          word = j*_size;
          for (u=0; u<_size; u++){
            sum1 = countp_t[word+u] + countp_e[word+u] + _prior_alpha[p][u];
            LogLik_te += gsl_sf_lngamma(sum1) - lng1[u];
            sum2 += sum1;
          }
          LogLik_te += lng2 - gsl_sf_lngamma(sum2);
        } 
    }
    
    if ((_postloglike>0)||(force==true)){
      for (p=0; p<_phase; p++){

        double sum = 0;
        for (u=0; u<_size; u++){
          lng1[u] = gsl_sf_lngamma(_prior_alpha[p][u]);
          sum += _prior_alpha[p][u];
        }
        double lng2 = gsl_sf_lngamma(sum);
        
        countp_t = (tseq_train.get_p_count())[p] + jump;
        for(j = 0; j < _nMu; j++){      
          sum1=0., sum2=0.;
          word = j*_size;
          for (u=0; u<_size; u++){
            sum1 = countp_t[word+u] + _prior_alpha[p][u];
            LogLik_t += gsl_sf_lngamma(sum1) - lng1[u];
            sum2 += sum1;
          }
          LogLik_t += lng2 - gsl_sf_lngamma(sum2);
        }       
      }
      _postloglike = LogLik_t;
    }
    
    return(LogLik_te - _postloglike);
   }

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

   inline void print( const string &FileOut );
  
  void print(ofstream &Out) const;

  //---------------------------------------
  int tell_size() const { return(_size);};
  int tell_rank() const { return(_rank);};
  int tell_order() const { return(_order);};
  int tell_phase() const { return(_phase); }

  //---------------------------------------
  int nMu() const { return(_nMu); }
  int nPi() const { return(_nPi); }


  double Pi(int index, int p=0) const { 
    return(_Pis[p][index]); 
  }

  double & operator() (int index, int p=0){ 
    return(_Pis[p][index]); 
  }


  double Mu(int index, int p=0) const { 
    return(_Mus[p][index]); 
  }  
  bool isPis() const { return(_Pis != NULL);};
  bool isMus()  const { return(_Mus != NULL);};


  inline short nextPhase(short p) const;
  inline short prevPhase(short p) const;
  //--------------------------------------- 

  bool Stochasticity();


  void file_to_count( const string& src_file, unsigned long ** dest_count);
   
 protected :
  
  short _phase;

  double **_Pis;
  double **_containers; 

  double **_Mus;   
  
  short _size; 
  
  short _order;
  
  long _nPi; 

  long _nMu; 
 
  long _nb_param;

  int _rank;

  long _jump;

  short *_nextPhase;

  short *_prevPhase;

  const Translator * _trans;

  double _postloglike;

  vector< vector<double> > _prior_alpha;


  inline void FreeNextPhase();
  inline void FreePrevPhase();

  bool isNextPhase()  const { return(_nextPhase != NULL);};
  bool isPrevPhase()  const { return(_prevPhase != NULL);}; 

  bool load_count( const string & count_file, unsigned long** data );

  bool load_prior_alpha( const string & prior_alpha_file );
};


inline short PhasedMarkov::nextPhase(short p) const
{
  return (_nextPhase[p]);
}

inline short PhasedMarkov::prevPhase(short p) const
{
  return (_prevPhase[p]);
}

inline void PhasedMarkov::new_markov_matrices()
{   
  if (!_Pis){
    _containers = new double*[_phase];
    _Pis = new double*[_phase];
    for (short p=0; p<_phase; p++){
      _containers[p] = new double[_nPi];
      _Pis[p] = _containers[p];
    }
  }
}

inline void PhasedMarkov::free_markov_matrices()
{
  if (_Pis){
    for (int p=0;p<_phase;p++)
      delete []_containers[p];
    delete []_containers;
    
    delete []_Pis;
  }
  _Pis = NULL;
  _containers = NULL;
}

inline void PhasedMarkov::free_stat_laws()
{
  if (_Mus)
    {
      for (int p=0;p<_phase;p++)
        delete []_Mus[p];
      delete []_Mus;
    }
  _Mus = NULL;
}

inline void PhasedMarkov::FreeNextPhase()
{
  if (isNextPhase())
    delete [] _nextPhase;
  _nextPhase = NULL;
}

inline void PhasedMarkov::FreePrevPhase()
{
  if (isPrevPhase())
    delete [] _prevPhase;
  _prevPhase = NULL;
}


inline void PhasedMarkov::print(const string & FileOut )
{
  ofstream Out(FileOut.c_str());
  print(Out);
}

#endif/*SEQPP_PHASEDMARKOV_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

---