modelbase.h

 
#pragma once
 
#include <cassert>
#include <istream>
#include <string>
#include "../alignment/alignment.h"
#include "../utils/matrix.h"
#include "../utils/tools.h"
 
class ModelsBlock;  // do not pull in external headers!
 
namespace cmaple {
class SeqRegions;
 
class ModelBase {
 private:
  void initEqualStateFreqs();
 
  cmaple::SeqRegion::SeqType getSeqType();
 
  void readStateFreq(std::istream& in);
 
  void updateMutMatbyMutCount();
 
 protected:
  // NHANLT: we can change to use unique_ptr(s) instead of normal pointers
  ModelsBlock* model_block = nullptr;
 
  void init();
 
  ModelsBlock* readModelsDefinition(const char* builtin_models);
 
  bool readParametersString(std::string& model_str);
 
  void readRates(std::istream& in, const bool is_reversible);
 
  void normalizeQMatrix();
 
  virtual void extractRootFreqs(const Alignment* aln);
 
  void initPointers();
 
  template <cmaple::StateType num_states>
  void updateMutationMat();
 
  template <cmaple::StateType num_states>
  bool updateMutationMatEmpiricalTemplate();
 public:
  enum SubModel {
    JC,
    GTR,
    UNREST,
 
    GTR20,
    NONREV,
    LG,
    WAG,
    JTT,
    Q_PFAM,
    Q_BIRD,
    Q_MAMMAL,
    Q_INSECT,
    Q_PLANT,
    Q_YEAST,
    JTTDCMUT,
    DCMUT,
    VT,
    PMB,
    BLOSUM62,
    DAYHOFF,
    MTREV,
    MTART,
    MTZOA,
    MTMET,
    MTVER,
    MTINV,
    MTMAM,
    FLAVI,
    HIVB,
    HIVW,
    FLU,
    RTREV,
    CPREV,
    NQ_PFAM,
    NQ_BIRD,
    NQ_MAMMAL,
    NQ_INSECT,
    NQ_PLANT,
    NQ_YEAST,
 
    DEFAULT,  
    UNKNOWN,  
  };
 
  ModelBase() = delete;
 
  ModelBase(const SubModel sub_model, const cmaple::StateType num_states);
 
   virtual ~ModelBase();
 
  const static std::map<std::string, SubModel> dna_models_mapping;
 
  const static std::map<std::string, SubModel> aa_models_mapping;
 
  // NHANLT: we can change to use unique_ptr(s) instead of normal pointers in
  // the following
  SubModel sub_model;
 
  const cmaple::StateType num_states_ = 0;
 
  cmaple::RealNumType* pseu_mutation_count = nullptr;
 
  cmaple::RealNumType* root_freqs = nullptr;
 
  cmaple::RealNumType* mutation_mat = nullptr;
 
  cmaple::RealNumType* root_log_freqs = nullptr;
 
  cmaple::RealNumType* diagonal_mut_mat = nullptr;
 
  cmaple::RealNumType* transposed_mut_mat = nullptr;
 
  cmaple::RealNumType* inverse_root_freqs = nullptr;
 
  cmaple::RealNumType* freqi_freqj_qij = nullptr;
 
  cmaple::RealNumType* freq_j_transposed_ij = nullptr;
 
  cmaple::StateType* row_index = nullptr;
    
  cmaple::RealNumType normalized_factor = 1.0;
 
  bool fixed_params = false;
 
  std::string getModelName() const;
 
  std::string exportRootFrequenciesStr();
 
  std::string exportQMatrixStr();
 
  void extractRefInfo(const Alignment* aln);
 
  virtual void initMutationMat(){}
 
  virtual bool updateMutationMatEmpirical() {
    return false;
  }
 
  virtual void updatePesudoCount(const Alignment* aln,
                                 const SeqRegions& regions1,
                                 const SeqRegions& regions2);
 
  static cmaple::SeqRegion::SeqType detectSeqType(const SubModel sub_model);
 
  static cmaple::ModelBase::SubModel parseModel(const std::string& model_name);
};
 
template <StateType num_states>
void cmaple::ModelBase::updateMutationMat() {
  // update Mutation matrix regarding the pseudo muation count
  updateMutMatbyMutCount();
 
  // compute the total rate regarding the root freqs
  RealNumType total_rate = 0;
  total_rate -= dotProduct<num_states>(root_freqs, diagonal_mut_mat);
    
  // record sum as the normalized factor
  normalized_factor = total_rate;
 
  // inverse total_rate
  total_rate = 1.0 / total_rate;
 
  // normalize the mutation_mat
  RealNumType* mutation_mat_row = mutation_mat;
  RealNumType* freqi_freqj_qij_row = freqi_freqj_qij;
  for (StateType i = 0; i < num_states_; ++i, mutation_mat_row += num_states_,
                 freqi_freqj_qij_row += num_states_) {
    for (StateType j = 0; j < num_states_; ++j) {
      mutation_mat_row[j] *= total_rate;
      // mutation_mat_row[j] /= total_rate;
 
      // update freqi_freqj_qij
      if (i != j) {
        freqi_freqj_qij_row[j] =
            root_freqs[i] * inverse_root_freqs[j] * mutation_mat_row[j];
            // root_freqs[i] / root_freqs[j] * mutation_mat_row[j];
      } else {
        freqi_freqj_qij_row[j] = mutation_mat_row[j];
      }
 
      // update the transposed mutation matrix
      transposed_mut_mat[row_index[j] + i] = mutation_mat_row[j];
    }
 
    // update diagonal
    diagonal_mut_mat[i] = mutation_mat_row[i];
  }
 
  // pre-compute matrix to speedup
  RealNumType* transposed_mut_mat_row = transposed_mut_mat;
  RealNumType* freq_j_transposed_ij_row = freq_j_transposed_ij;
 
  for (StateType i = 0; i < num_states_; ++i,
                 transposed_mut_mat_row += num_states_,
                 freq_j_transposed_ij_row += num_states_) {
    setVecByProduct<num_states>(freq_j_transposed_ij_row, root_freqs,
                                transposed_mut_mat_row);
  }
}
 
template <StateType num_states>
auto cmaple::ModelBase::updateMutationMatEmpiricalTemplate()
    -> bool {
  bool update = false;
 
  if (!fixed_params) {
    // clone the current mutation matrix
    RealNumType* tmp_mut_mat = new RealNumType[row_index[num_states_]];
    memcpy(tmp_mut_mat, mutation_mat,
           row_index[num_states_] * sizeof(RealNumType));
 
    // update the mutation matrix regarding the pseu_mutation_count
    updateMutationMat<num_states>();
 
    // set update = true if the mutation matrix changes more than a threshold
    const RealNumType change_thresh = 1e-3;
    RealNumType sum_change = 0;
    RealNumType* tmp_mut_mat_ptr = tmp_mut_mat;
    RealNumType* mutation_mat_ptr = mutation_mat;
    for (StateType i = 0; i < num_states_;
         ++i, tmp_mut_mat_ptr += num_states_, mutation_mat_ptr += num_states_) {
      for (StateType j = 0; j < num_states_; ++j) {
        sum_change += fabs(tmp_mut_mat_ptr[j] - mutation_mat_ptr[j]);
      }
 
      sum_change -= fabs(tmp_mut_mat_ptr[i] - mutation_mat_ptr[i]);
    }
 
    update = sum_change > change_thresh;
 
    // delete tmp_diagonal_mutation_mat
    delete[] tmp_mut_mat;
  }
 
  // return update
  return update;
}
}  // namespace cmaple