1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
|
/*
* modelliemarkov.h
*
* Created on: 24/05/2016
* Author: Michael Woodhams
*/
#ifndef MODELLIEMARKOV_H_
#define MODELLIEMARKOV_H_
#include "modelmarkov.h"
class ModelLieMarkov: public ModelMarkov {
public:
ModelLieMarkov(string model_name, PhyloTree *tree, string model_params, StateFreqType freq_type, string freq_params);
virtual ~ModelLieMarkov();
/**
this function is served for model testing
@param model_name name of the model
@param freq_type state frequency type, can be FREQ_USER_DEFINED, FREQ_EQUAL, FREQ_EMPIRICAL, or FREQ_ESTIMATE
*/
virtual void init(const char *model_name, string model_params, StateFreqType freq, string freq_params);
/**
start structure for checkpointing
*/
virtual void startCheckpoint();
/**
save object into the checkpoint
*/
virtual void saveCheckpoint();
/**
restore object from the checkpoint
*/
virtual void restoreCheckpoint();
/**
write information
@param out output stream
*/
virtual void writeInfo(ostream &out);
static void getLieMarkovModelInfo(string model_name, string &name, string &full_name, int &model_num, int &symmetry, StateFreqType &def_freq);
static string getModelInfo(string model_name, string &full_name, StateFreqType &def_freq);
// DO NOT override this function, because
// BQM, 2017-05-02: getNDimFreq should return degree of freedom, which is not included in getNDim()
// That's why 0 is returned for FREQ_ESTIMATE, num_states-1 for FREQ_EMPIRICAL
// virtual int getNDimFreq();
// this is redundant, there is already the same function below
// bool isTimeReversible();
/**
@return TRUE if model is time-reversible, FALSE otherwise
*/
virtual bool isReversible();
static bool validModelName(string model_name);
void setBounds(double *lower_bound, double *upper_bound, bool *bound_check);
/**
decompose the rate matrix into eigenvalues and eigenvectors
*/
virtual void decomposeRateMatrix();
/** decompose rate matrix using closed formula derived by Michael Woodhams */
void decomposeRateMatrixClosedForm();
/** decompose rate matrix using Eigen library */
void decomposeRateMatrixEigen3lib();
/**
compute the transition probability matrix.
@param time time between two events
@param mixture (optional) class for mixture model
@param trans_matrix (OUT) the transition matrix between all pairs of states.
Assume trans_matrix has size of num_states * num_states.
*/
virtual void computeTransMatrix(double time, double *trans_matrix, int mixture = 0);
// overrides Optimization::restartParameters
bool restartParameters(double guess[], int ndim, double lower[], double upper[], bool bound_check[], int iteration);
protected:
/**
Model parameters - cached so we know when they change, and thus when
recalculations are needed.
*/
double *model_parameters;
double **basis;
int symmetry; // RY->0, WS->1, MK->2
int model_num; // 0->1.1, etc to 36->12.12
void setBasis();
virtual void setRates();
bool nondiagonalizable; // will be set true for nondiagonalizable rate matrices, then will use scaled squaring method for matrix exponentiation.
/**
this function is served for the multi-dimension optimization. It should pack the model parameters
into a vector that is index from 1 (NOTE: not from 0)
@param variables (OUT) vector of variables, indexed from 1
*/
virtual void setVariables(double *variables);
/**
this function is served for the multi-dimension optimization. It should assign the model parameters
from a vector of variables that is index from 1 (NOTE: not from 0)
@param variables vector of variables, indexed from 1
@return TRUE if parameters are changed, FALSE otherwise (2015-10-20)
*/
virtual bool getVariables(double *variables);
static void parseModelName(string model_name, int* model_num, int* symmetry);
/*
* Overrides ModelMarkov::getName().
* Avoids appending +FO to name, as this is implied by how LM models
* work.
* Minh: you might chose to remove this override, if you like "+FO"
* to be on LM model names.
*/
string getName();
/*
const static double ***BASES;
const static int *MODEL_PARAMS;
const static string *SYMMETRY;
const static string *MODEL_NAMES;
const static int NUM_RATES;
const static int NUM_LM_MODELS;
*/
bool validFreqType();
};
#endif /* MODELLIEMARKOV_H_ */
|
