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
|
/***************************************************************************
* Copyright (C) 2009 by BUI Quang Minh *
* minh.bui@univie.ac.at *
* *
* 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., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#include "rateinvar.h"
RateInvar::RateInvar(double p_invar_sites, PhyloTree *tree)
: RateHeterogeneity()
{
if (tree) {
if (tree->aln->frac_const_sites == 0.0)
p_invar = 0.0;
else
p_invar = max(tree->aln->frac_const_sites/2.0, MIN_PINVAR);
// p_invar = MIN_PINVAR;
} else
p_invar = MIN_PINVAR;
fix_p_invar = false;
optimize_p_invar = true;
phylo_tree = tree;
name = "+I";
full_name = "Invar";
if (p_invar_sites >= 0) {
p_invar = p_invar_sites;
fix_p_invar = true;
}
}
void RateInvar::saveCheckpoint() {
checkpoint->startStruct("RateInvar");
CKP_SAVE(p_invar);
// CKP_SAVE(fix_p_invar);
// CKP_SAVE(optimize_p_invar);
checkpoint->endStruct();
RateHeterogeneity::saveCheckpoint();
}
void RateInvar::restoreCheckpoint() {
RateHeterogeneity::restoreCheckpoint();
checkpoint->startStruct("RateInvar");
CKP_RESTORE(p_invar);
// CKP_RESTORE(fix_p_invar);
// CKP_RESTORE(optimize_p_invar);
checkpoint->endStruct();
}
string RateInvar::getNameParams() {
ostringstream str;
str << "+I{" << p_invar << '}';
return str.str();
}
double RateInvar::computeFunction(double p_invar_value) {
p_invar = p_invar_value;
// fix bug: computeTip... will update ptn_invar vector
// phylo_tree->computePtnInvar();
phylo_tree->clearAllPartialLH();
return -phylo_tree->computeLikelihood();
}
double RateInvar::targetFunk(double x[]) {
getVariables(x);
// fix bug: computeTip... will update ptn_invar vector
phylo_tree->computePtnInvar();
return -phylo_tree->computeLikelihood();
}
void RateInvar::setBounds(double *lower_bound, double *upper_bound, bool *bound_check) {
if (getNDim() == 0) return;
lower_bound[1] = MIN_PINVAR;
upper_bound[1] = phylo_tree->aln->frac_const_sites;
bound_check[1] = false;
}
double RateInvar::optimizeParameters(double gradient_epsilon) {
if (phylo_tree->aln->frac_const_sites == 0.0)
return -computeFunction(0.0);
if (fix_p_invar || !optimize_p_invar)
return -computeFunction(p_invar);
if (verbose_mode >= VB_MAX)
cout << "Optimizing proportion of invariable sites..." << endl;
double negative_lh;
double ferror;
p_invar = minimizeOneDimen(MIN_PINVAR, p_invar, min(phylo_tree->aln->frac_const_sites, 1.0-MIN_PINVAR), max(gradient_epsilon, TOL_PINVAR), &negative_lh, &ferror);
//p_invar = minimizeOneDimen(MIN_PINVAR, p_invar, 1.0 - MIN_PINVAR, TOL_PINVAR, &negative_lh, &ferror);
// phylo_tree->clearAllPartialLH();
// phylo_tree->computePtnInvar();
// return -negative_lh;
return -computeFunction(p_invar);
}
void RateInvar::writeInfo(ostream &out) {
out << "Proportion of invariable sites: " << p_invar << endl;
}
void RateInvar::writeParameters(ostream &out) {
out << "\t" << p_invar;
}
void RateInvar::setVariables(double *variables) {
if (RateInvar::getNDim() == 0) return;
variables[1] = p_invar;
}
bool RateInvar::getVariables(double *variables) {
if (RateInvar::getNDim() == 0) return false;
bool changed = (p_invar != variables[1]);
p_invar = variables[1];
return changed;
}
|
