1 /*************************************************************************** 2 * Copyright (C) 2009 by BUI Quang Minh * 3 * minh.bui@univie.ac.at * 4 * * 5 * This program is free software; you can redistribute it and/or modify * 6 * it under the terms of the GNU General Public License as published by * 7 * the Free Software Foundation; either version 2 of the License, or * 8 * (at your option) any later version. * 9 * * 10 * This program is distributed in the hope that it will be useful, * 11 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * 13 * GNU General Public License for more details. * 14 * * 15 * You should have received a copy of the GNU General Public License * 16 * along with this program; if not, write to the * 17 * Free Software Foundation, Inc., * 18 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * 19 ***************************************************************************/ 20 #ifndef MODELMARKOV_H 21 #define MODELMARKOV_H 22 23 #include "tree/phylotree.h" 24 #include "modelsubst.h" 25 #include "utils/optimization.h" 26 #include "alignment/alignment.h" 27 #include "utils/eigendecomposition.h" 28 #include <complex> 29 30 const double MIN_RATE = 1e-4; 31 const double TOL_RATE = 1e-4; 32 const double MAX_RATE = 100; 33 34 string freqTypeString(StateFreqType freq_type, SeqType seq_type, bool full_str); 35 36 /** 37 General Markov model of substitution (reversible or non-reversible) 38 This works for all kind of data 39 40 @author BUI Quang Minh <minh.bui@univie.ac.at> 41 */ 42 class ModelMarkov : public ModelSubst, public EigenDecomposition 43 { 44 45 friend class ModelSet; 46 friend class ModelMixture; 47 friend class ModelPoMo; 48 friend class PartitionModel; 49 friend class PartitionModelPlen; 50 51 public: 52 /** 53 constructor 54 @param tree associated tree for the model 55 @param reversible TRUE (default) for reversible model, FALSE for non-reversible 56 @param adapt_tree TRUE (default) to convert rooted<->unrooted tree 57 */ 58 ModelMarkov(PhyloTree *tree, bool reversible = true, bool adapt_tree = true); 59 60 /** 61 @return TRUE if model is time-reversible, FALSE otherwise 62 */ isReversible()63 virtual bool isReversible() { return is_reversible; }; 64 65 /** 66 set the reversibility of the model 67 @param reversible TRUE to make model reversible, FALSE otherwise 68 @param adapt_tree TRUE (default) to convert between rooted and unrooted tree 69 */ 70 virtual void setReversible(bool reversible, bool adapt_tree = true); 71 72 73 /** 74 init the model and decompose the rate matrix. This function should always be called 75 after creating the class. Otherwise it will not work properly. 76 @param freq_type state frequency type, can be FREQ_USER_DEFINED, FREQ_EQUAL, FREQ_EMPIRICAL, or FREQ_ESTIMATE 77 */ 78 void init(StateFreqType freq_type); 79 80 /** 81 initializes ModelSubst::freq_type array according to freq_type 82 (can be FREQ_USER_DEFINED, FREQ_EQUAL, FREQ_EMPIRICAL, or FREQ_ESTIMATE) 83 */ 84 void init_state_freq(StateFreqType freq_type); 85 86 /** 87 this function is served for ModelDNA and ModelProtein 88 @param model_name name of the model 89 @param freq_type state frequency type, can be FREQ_USER_DEFINED, FREQ_EQUAL, FREQ_EMPIRICAL, or FREQ_ESTIMATE 90 */ init(const char * model_name,string model_params,StateFreqType freq,string freq_params)91 virtual void init(const char *model_name, string model_params, StateFreqType freq, string freq_params) {} 92 93 /** 94 destructor 95 */ 96 virtual ~ModelMarkov(); 97 98 /** 99 start structure for checkpointing 100 */ 101 virtual void startCheckpoint(); 102 103 /** 104 save object into the checkpoint 105 */ 106 virtual void saveCheckpoint(); 107 108 /** 109 restore object from the checkpoint 110 */ 111 virtual void restoreCheckpoint(); 112 113 /** 114 * @return model name 115 */ 116 virtual string getName(); 117 118 /** 119 * @return model name with parameters in form of e.g. GTR{a,b,c,d,e,f} 120 */ 121 virtual string getNameParams(); 122 123 /** 124 internal function: return string for frequency 125 @param retname output stream 126 */ 127 void getNameParamsFreq(ostream &retname); 128 129 /** 130 @return the number of rate entries, equal to the number of non-diagonal elements 131 of the rate matrix (since model is NOT reversible) 132 */ 133 virtual int getNumRateEntries(); 134 135 /** 136 set the associated tree 137 @param tree the associated tree 138 */ 139 void setTree(PhyloTree *tree); 140 141 /** 142 Read the upper-triangle rate matrix from an input stream. 143 It will throw error messages if failed 144 @param in input stream 145 */ 146 virtual void readRates(istream &in) throw(const char*, string); 147 148 /** 149 Read the rate parameters from a comma-separated string 150 It will throw error messages if failed 151 @param in input stream 152 */ 153 virtual void readRates(string str) throw(const char*); 154 155 /** 156 Read state frequencies from an input stream. 157 It will throw error messages if failed 158 @param in input stream 159 */ 160 virtual void readStateFreq(istream &in) throw(const char*); 161 162 /** 163 Read state frequencies from comma-separated string 164 It will throw error messages if failed 165 @param str input string 166 */ 167 virtual void readStateFreq(string str) throw(const char*); 168 169 /** 170 read model parameters from a file 171 @param file_name file containing rate matrix and state frequencies 172 */ 173 void readParameters(const char *file_name, bool adapt_tree = true); 174 175 /** 176 read model parameters from a string 177 @param model_str string containing rate matrix and state frequencies 178 */ 179 void readParametersString(string &model_str, bool adapt_tree = true); 180 181 /** 182 compute the transition probability matrix. 183 @param time time between two events 184 @param mixture (optional) class for mixture model 185 @param trans_matrix (OUT) the transition matrix between all pairs of states. 186 Assume trans_matrix has size of num_states * num_states. 187 */ 188 virtual void computeTransMatrix(double time, double *trans_matrix, int mixture = 0); 189 190 /** 191 compute the transition probability matrix for non-reversible model 192 @param time time between two events 193 @param mixture (optional) class for mixture model 194 @param trans_matrix (OUT) the transition matrix between all pairs of states. 195 Assume trans_matrix has size of num_states * num_states. 196 */ 197 virtual void computeTransMatrixNonrev(double time, double *trans_matrix, int mixture = 0); 198 199 /** 200 compute the transition probability between two states 201 @param time time between two events 202 @param state1 first state 203 @param state2 second state 204 */ 205 virtual double computeTrans(double time, int state1, int state2); 206 207 /** 208 compute the transition probability between two states 209 @param time time between two events 210 @param state1 first state 211 @param state2 second state 212 @param derv1 (OUT) 1st derivative 213 @param derv2 (OUT) 2nd derivative 214 */ 215 virtual double computeTrans(double time, int state1, int state2, double &derv1, double &derv2); 216 217 /** 218 Get the rate matrix. 219 @param rate_mat (OUT) upper-triagle rate matrix. Assume rate_mat has size of num_states*(num_states-1)/2 220 */ 221 virtual void getRateMatrix(double *rate_mat); 222 223 /** 224 Set the rate matrix. 225 @param rate_mat upper-triagle rate matrix. Assume rate_mat has size of num_states*(num_states-1)/2 226 */ 227 virtual void setRateMatrix(double *rate_mat); 228 229 /** 230 Set the full rate matrix of size num_states*num_states 231 @param rate_mat full rate matrix 232 @param freq state frequency 233 */ 234 virtual void setFullRateMatrix(double *rate_mat, double *freq); 235 236 /** 237 compute the state frequency vector 238 @param mixture (optional) class for mixture model 239 @param state_freq (OUT) state frequency vector. Assume state_freq has size of num_states 240 */ 241 virtual void getStateFrequency(double *state_freq, int mixture = 0); 242 243 /** 244 set the state frequency vector 245 @param state_freq (IN) state frequency vector. Assume state_freq has size of num_states 246 */ 247 virtual void setStateFrequency(double *state_freq); 248 249 /** 250 set the state frequency vector 251 @param state_freq (IN) state frequency vector. Assume state_freq has size of num_states 252 */ 253 virtual void adaptStateFrequency(double *state_freq); 254 255 /** 256 * compute Q matrix 257 * @param q_mat (OUT) Q matrix, assuming of size num_states * num_states 258 */ 259 virtual void getQMatrix(double *q_mat); 260 261 /** 262 rescale the state frequencies 263 @param sum_one TRUE to make frequencies sum to 1, FALSE to make last entry equal to 1 264 */ 265 void scaleStateFreq(bool sum_one); 266 267 /** 268 get frequency type 269 @return frequency type 270 */ getFreqType()271 virtual StateFreqType getFreqType() { return freq_type; } 272 273 274 /** 275 compute the transition probability matrix.and the derivative 1 and 2 276 @param time time between two events 277 @param mixture (optional) class for mixture model 278 @param trans_matrix (OUT) the transition matrix between all pairs of states. 279 Assume trans_matrix has size of num_states * num_states. 280 @param trans_derv1 (OUT) the 1st derivative matrix between all pairs of states. 281 @param trans_derv2 (OUT) the 2nd derivative matrix between all pairs of states. 282 */ 283 virtual void computeTransDerv(double time, double *trans_matrix, 284 double *trans_derv1, double *trans_derv2, int mixture = 0); 285 286 /** 287 @return the number of dimensions 288 */ 289 virtual int getNDim(); 290 291 /** 292 @return the number of dimensions corresponding to state frequencies, which is 293 not counted in getNDim(). This serves e.g. for computing AIC, BIC score 294 */ 295 virtual int getNDimFreq(); 296 297 298 /** 299 the target function which needs to be optimized 300 @param x the input vector x 301 @return the function value at x 302 */ 303 virtual double targetFunk(double x[]); 304 305 /** 306 * setup the bounds for joint optimization with BFGS 307 */ 308 virtual void setBounds(double *lower_bound, double *upper_bound, bool *bound_check); 309 310 /** 311 optimize model parameters 312 @return the best likelihood 313 */ 314 virtual double optimizeParameters(double gradient_epsilon); 315 316 /** 317 * @return TRUE if parameters are at the boundary that may cause numerical unstability 318 */ 319 virtual bool isUnstableParameters(); 320 321 // A simple helper function that prints the rates in a nice way and can be 322 // reused by children. The title is necessary, because, e.g., for PoMo, the 323 // rates are mutation rates and not substitution rates, and also 324 // exchangeabilities may be reported. 325 void report_rates(ostream &out, string title, double *r); 326 327 // Report the stationary frequencies state_freq or custom_state_freq (if 328 // given) to output stream out. 329 void report_state_freqs(ostream &out, double *custom_state_freq=NULL); 330 331 /** 332 write information 333 @param out output stream 334 */ 335 virtual void writeInfo(ostream &out); 336 337 /** 338 write parameters, used with modeltest 339 @param out output stream 340 */ writeParameters(ostream & out)341 virtual void writeParameters(ostream &out){} 342 343 /** decompose rate matrix for non-reversible models */ 344 virtual void decomposeRateMatrixNonrev(); 345 346 /** old version of decompose rate matrix for reversible models */ 347 void decomposeRateMatrixRev(); 348 349 /** 350 decompose the rate matrix into eigenvalues and eigenvectors 351 */ 352 virtual void decomposeRateMatrix(); 353 354 // double *getEigenCoeff() const; 355 356 virtual double *getEigenvalues() const; 357 358 virtual double *getEigenvectors() const; 359 virtual double *getInverseEigenvectors() const; 360 361 // void setEigenCoeff(double *eigenCoeff); 362 363 void setEigenvalues(double *eigenvalues); 364 365 void setEigenvectors(double *eigenvectors); 366 367 void setInverseEigenvectors(double *inv_eigenvectors); 368 369 /** 370 * compute the memory size for the model, can be large for site-specific models 371 * @return memory size required in bytes 372 */ getMemoryRequired()373 virtual uint64_t getMemoryRequired() { 374 return ModelSubst::getMemoryRequired() + sizeof(double)*num_states*num_states*3; 375 } 376 377 /** default TRUE: store only upper half of the rate matrix */ 378 bool half_matrix; 379 380 /****************************************************/ 381 /* NON-REVERSIBLE STUFFS */ 382 /****************************************************/ 383 384 /** 385 * Return a model of type given by model_name. (Will be some subclass of ModelMarkov.) 386 */ 387 static ModelMarkov* getModelByName(string model_name, PhyloTree *tree, string model_params, StateFreqType freq_type, string freq_params); 388 389 /** 390 * true if model_name is the name of some known non-reversible model 391 */ 392 static bool validModelName(string model_name); 393 394 // Mon Jul 3 14:47:08 BST 2017; added by Dominik. I had problems with mixture 395 // models together with PoMo and rate heterogeneity. E.g., a model 396 // "MIX{HKY+P+N9+G2,GTR+P+N9+G2}" leads to segmentation faults because the 397 // `ModelPoMoMixture` reports a /wrong/ number of states (i.e., it reports 52 398 // instead of 104). Consequently, the `initMem()` function of ModelMixture, 399 // messes up the `eigenvalues`, etc., variables of the `ModelPoMoMixture`s. I 400 // circumvent this, by adding this virtual function; for normal models, it 401 // just returns `num_states`, however, for mixture models, it returns 402 // `num_states*nmixtures`. 403 virtual int get_num_states_total(); 404 405 // Mon Jul 3 15:53:00 BST 2017; added by Dominik. Same problem as with 406 // `get_num_states_total()`. The pointers to the eigenvalues and eigenvectors 407 // need to be updated recursively, if the model is a mixture model. For a 408 // normal Markov model, only the standard pointers are set. This was done in 409 // `ModelMixture::initMem()` before. 410 virtual void update_eigen_pointers(double *eval, double *evec, double *inv_evec); 411 412 413 /** 414 set num_params variable 415 */ setNParams(int num_params)416 virtual void setNParams(int num_params) { 417 this->num_params = num_params; 418 } 419 420 /** 421 get num_params variable 422 */ getNParams()423 virtual int getNParams() { 424 return num_params; 425 } 426 427 protected: 428 429 /** 430 this function is served for the multi-dimension optimization. It should pack the model parameters 431 into a vector that is index from 1 (NOTE: not from 0) 432 @param variables (OUT) vector of variables, indexed from 1 433 */ 434 virtual void setVariables(double *variables); 435 436 /** 437 this function is served for the multi-dimension optimization. It should assign the model parameters 438 from a vector of variables that is index from 1 (NOTE: not from 0) 439 @param variables vector of variables, indexed from 1 440 @return TRUE if parameters are changed, FALSE otherwise (2015-10-20) 441 */ 442 virtual bool getVariables(double *variables); 443 444 445 /** 446 * Called from getVariables to update the rate matrix for the new 447 * model parameters. 448 */ 449 virtual void setRates(); 450 451 /** 452 free all allocated memory 453 */ 454 virtual void freeMem(); 455 456 /** TRUE if model is reversible */ 457 bool is_reversible; 458 459 /** 460 phylogenetic tree associated 461 */ 462 PhyloTree *phylo_tree; 463 464 /** 465 rates between pairs of states of the unit rate matrix Q. 466 In order A-C, A-G, A-T, C-G, C-T (rate G-T = 1 always) 467 */ 468 double *rates; 469 470 /** 471 the number of free rate parameters 472 */ 473 int num_params; 474 475 /** 476 eigenvalues of the rate matrix Q 477 */ 478 double *eigenvalues; 479 480 /** 481 eigenvectors of the rate matrix Q 482 */ 483 double *eigenvectors; 484 485 /** 486 inverse eigenvectors of the rate matrix Q 487 */ 488 double *inv_eigenvectors; 489 490 /** 491 coefficient cache, served for fast computation of the P(t) matrix 492 */ 493 // double *eigen_coeff; 494 495 /** state with highest frequency, used when optimizing state frequencies +FO */ 496 int highest_freq_state; 497 498 /****************************************************/ 499 /* NON-REVERSIBLE STUFFS */ 500 /****************************************************/ 501 502 /** 503 compute the transition probability matrix using (complex) eigenvalues 504 @param time time between two events 505 @param trans_matrix (OUT) the transition matrix between all pairs of states. 506 Assume trans_matrix has size of num_states * num_states. 507 */ 508 void computeTransMatrixEigen(double time, double *trans_matrix); 509 510 /** 511 unrestricted Q matrix. Note that Q is normalized to 1 and has row sums of 0. 512 no state frequencies are involved here since Q is a general matrix. 513 */ 514 double *rate_matrix; 515 516 /** imaginary part of eigenvalues */ 517 double *eigenvalues_imag; 518 519 /** 520 complex eigenvalues and eigenvectors, pointing to the same pointer 521 to the previous double *eigenvalues and double *eigenvectors 522 */ 523 std::complex<double> *ceval, *cevec, *cinv_evec; 524 525 /** will be set true for nondiagonalizable rate matrices, 526 then will use scaled squaring method for matrix exponentiation. 527 */ 528 bool nondiagonalizable; 529 530 }; 531 532 #endif 533