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
|
/*
* modelmorphology.cpp
*
* Created on: Apr 15, 2014
* Author: minh
*/
#include "modelmorphology.h"
ModelMorphology::ModelMorphology(const char *model_name, string model_params, StateFreqType freq, string freq_params, PhyloTree *tree)
: ModelMarkov(tree)
{
init(model_name, model_params, freq, freq_params);
}
void ModelMorphology::init(const char *model_name, string model_params, StateFreqType freq, string freq_params)
{
name = model_name;
full_name = model_name;
if (name == "MK") {
// all were initialized
num_params = 0;
} else if (name == "ORDERED") {
int i, j, k = 0;
// only allow for substitution from state i to state i+1 and back.
for (i = 0; i < num_states-1; i++) {
rates[k++] = 1.0;
for (j = i+2; j < num_states; j++, k++)
rates[k] = 0.0;
}
num_params = 0;
} else if (name == "GTR") {
outWarning("GTR multistate model will estimate " + convertIntToString(getNumRateEntries()-1) + " substitution rates that might be overfitting!");
outWarning("Please only use GTR with very large data and always test for model fit!");
} else {
// if name does not match, read the user-defined model
readParameters(model_name);
num_params = 0;
freq = FREQ_USER_DEFINED;
}
ModelMarkov::init(freq);
}
void ModelMorphology::readRates(istream &in) throw(const char*, string) {
int nrates = getNumRateEntries();
int row = 1, col = 0;
// since states for protein is stored in lower-triangle, special treatment is needed
for (int i = 0; i < nrates; i++, col++) {
if (col == row) {
row++; col = 0;
}
// switch col and row
int id = col*(2*num_states-col-1)/2 + (row-col-1);
if (id >= nrates) {
cout << row << " " << col << endl;
}
assert(id < nrates && id >= 0); // make sure that the conversion is correct
if (!(in >> rates[id]))
throw name+string(": Rate entries could not be read");
if (rates[id] < 0.0)
throw "Negative rates found";
}
}
int ModelMorphology::getNDim() {
int ndim = num_params;
if (freq_type == FREQ_ESTIMATE)
ndim += num_states-1;
return ndim;
}
ModelMorphology::~ModelMorphology() {
}
void ModelMorphology::startCheckpoint() {
checkpoint->startStruct("ModelMorph");
}
void ModelMorphology::saveCheckpoint() {
startCheckpoint();
if (num_params > 0)
CKP_ARRAY_SAVE(getNumRateEntries(), rates);
endCheckpoint();
ModelMarkov::saveCheckpoint();
}
void ModelMorphology::restoreCheckpoint() {
ModelMarkov::restoreCheckpoint();
startCheckpoint();
if (num_params > 0)
CKP_ARRAY_RESTORE(getNumRateEntries(), rates);
endCheckpoint();
decomposeRateMatrix();
if (phylo_tree)
phylo_tree->clearAllPartialLH();
}
string ModelMorphology::getNameParams() {
if (num_params == 0) return name;
ostringstream retname;
retname << name << '{';
int nrates = getNumRateEntries();
for (int i = 0; i < nrates; i++) {
if (i>0) retname << ',';
retname << rates[i];
}
retname << '}';
getNameParamsFreq(retname);
return retname.str();
}
void ModelMorphology::writeParameters(ostream &out) {
int i;
if (freq_type == FREQ_ESTIMATE) {
for (i = 0; i < num_states; i++)
out << "\t" << state_freq[i];
}
if (num_params == 0) return;
int nrateout = getNumRateEntries() - 1;
for (i = 0; i < nrateout; i++)
out << "\t" << rates[i];
}
void ModelMorphology::writeInfo(ostream &out) {
if (num_params > 0) {
out << "Rate parameters:";
int nrate = getNumRateEntries();
for (int i = 0; i < nrate; i++)
out << " " << rates[i];
out << endl;
}
if (freq_type != FREQ_EQUAL) {
out << "State frequencies:";
for (int i = 0; i < num_states; i++)
out << " " << state_freq[i];
out << endl;
}
}
|