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
144
145
146
147
148
149
|
// $Id: ssrvDistanceSeqs2Tree.cpp 962 2006-11-07 15:13:34Z privmane $
#include "ssrvDistanceSeqs2Tree.h"
//#include "bestAlphaAndNu.h"
#include "bestParamUSSRV.h"
#include "someUtil.h"
#include <float.h>
tree ssrvDistanceSeqs2Tree::seqs2TreeIterative(const sequenceContainer &sc, MDOUBLE initAlpha, MDOUBLE initNu, const Vdouble *weights, const tree* constraintTreePtr) {
_constraintTreePtr=constraintTreePtr;
_alpha = initAlpha;
_newNu = _nu = initNu;
_weights = weights;
return seqs2TreeIterativeInternal(sc, true);
}
tree ssrvDistanceSeqs2Tree::seqs2TreeIterative(const sequenceContainer &sc, const Vdouble *weights, const tree* constraintTreePtr) {
_constraintTreePtr=constraintTreePtr;
_weights = weights;
return seqs2TreeIterativeInternal(sc, false);
}
tree ssrvDistanceSeqs2Tree::seqs2TreeIterative(const sequenceContainer &sc, const tree &initTree, const Vdouble *weights, const tree* constraintTreePtr) {
_constraintTreePtr=constraintTreePtr;
_weights = weights;
return seqs2TreeIterativeInternalInitTreeGiven(sc, initTree);
}
tree ssrvDistanceSeqs2Tree::seqs2TreeIterative(const sequenceContainer &sc, const tree &initTree, MDOUBLE initAlpha, const Vdouble *weights, const tree* constraintTreePtr) {
_alpha = initAlpha;
_weights = weights;
_constraintTreePtr=constraintTreePtr;
return seqs2TreeIterativeInternalInitTreeGiven(sc, false, initTree, initAlpha);
}
tree ssrvDistanceSeqs2Tree::seqs2TreeIterative(const sequenceContainer &sc, const tree &initTree, MDOUBLE initAlpha, MDOUBLE initNu, const Vdouble *weights, const tree* constraintTreePtr) {
_alpha = initAlpha;
_newNu = _nu = initNu;
_weights = weights;
_constraintTreePtr=constraintTreePtr;
return seqs2TreeIterativeInternalInitTreeGiven(sc, true, initTree, initAlpha);
}
// NOTE! This version is a NON-ITERATIVE version that uses the side info supplied by the user
tree ssrvDistanceSeqs2Tree::seqs2Tree(const sequenceContainer &sc, MDOUBLE alpha, MDOUBLE nu, const Vdouble *weights, const tree* constraintTreePtr) {
_weights = weights;
_alpha = alpha;
_newNu = _nu = nu;
_constraintTreePtr=constraintTreePtr;
seqs2TreeOneIterationInternal(sc, true);
return _newTree;
}
tree ssrvDistanceSeqs2Tree::seqs2TreeBootstrap(const sequenceContainer &sc, const MDOUBLE alpha, MDOUBLE nu, const Vdouble *weights, const tree* constraintTreePtr) {
_weights = weights;
_alpha = alpha;
_newNu = _nu = nu;
return static_cast<iterativeDistanceSeqs2Tree *>(this)->seqs2TreeBootstrap(sc, weights, constraintTreePtr);
}
// NOTE! This version calls ITERATIVE seqs2Tree because side info is not given by the user, so we have to generate and optimize it
tree ssrvDistanceSeqs2Tree::seqs2Tree(const sequenceContainer &sc, const Vdouble *weights, const tree* constraintTreePtr) {
return seqs2TreeIterative(sc,weights,constraintTreePtr);
}
MDOUBLE ssrvDistanceSeqs2Tree::optimizeSideInfo(const sequenceContainer &sc, tree &et)
{
if (!dynamic_cast<tamura92*>(
static_cast<replacementModelSSRV*>(_spPtr->getPijAccelerator()->getReplacementModel())
->getBaseRM()
)
) {
bestParamSSRV optimizer(true,true,false,true); // optimize alpha, nu, NOT tamura92 params, and bbl
optimizer(et,sc,*static_cast<stochasticProcessSSRV*>(_spPtr),_weights,
15,15,0.5,_epsilonLikelihoodImprovement4alphaOptimiz,_epsilonLikelihoodImprovement,
_epsilonLikelihoodImprovement4BBL,_maxIterationsBBL,5);
_newAlpha=optimizer.getBestAlpha();
_newNu=optimizer.getBestNu();
return(optimizer.getBestL());
} else {
bestParamSSRV optimizer(true,true,true,true); // optimize alpha, nu, tamura92 params, and bbl
optimizer(et,sc,*static_cast<stochasticProcessSSRV*>(_spPtr),_weights,
15,15,0.5,_epsilonLikelihoodImprovement4alphaOptimiz,_epsilonLikelihoodImprovement,
_epsilonLikelihoodImprovement4BBL,_maxIterationsBBL,5);
_newAlpha=optimizer.getBestAlpha();
_newNu=optimizer.getBestNu();
return(optimizer.getBestL());
}
}
MDOUBLE ssrvDistanceSeqs2Tree::calcSideInfoGivenTreeAndAlpha(const sequenceContainer &sc, const tree &et, MDOUBLE alpha)
{
_newAlpha = alpha;
(static_cast<gammaDistribution*>(_spPtr->distr()))->setAlpha(alpha);
// optimize only nu (and tamura92 params, if relevant)
if (!dynamic_cast<tamura92*>(
static_cast<replacementModelSSRV*>(_spPtr->getPijAccelerator()->getReplacementModel())
->getBaseRM()
)
) {
bestParamSSRV optimizer(false,true,false,false);
optimizer(et,sc,*(static_cast<stochasticProcessSSRV*>(_spPtr)),_weights,
15,15,_epsilonLikelihoodImprovement4alphaOptimiz,_epsilonLikelihoodImprovement,
_epsilonLikelihoodImprovement4BBL,_maxIterationsBBL,5);
_newNu=optimizer.getBestNu();
return(optimizer.getBestL());
} else {
bestParamSSRV optimizer(false,true,true,false);
optimizer(et,sc,*(static_cast<stochasticProcessSSRV*>(_spPtr)),_weights,
15,15,_epsilonLikelihoodImprovement4alphaOptimiz,_epsilonLikelihoodImprovement,
_epsilonLikelihoodImprovement4BBL,_maxIterationsBBL,5);
_newNu=optimizer.getBestNu();
return(optimizer.getBestL());
}
}
void ssrvDistanceSeqs2Tree::acceptSideInfo()
{
_alpha = _newAlpha;
_nu = _newNu;
}
void ssrvDistanceSeqs2Tree::utilizeSideInfo()
{
// set new alpha value in the sp that is used in _distM
LOG(10,<<"# utilizing alpha "<<_alpha<<" and nu "<<_nu<<endl);
(static_cast<gammaDistribution*>(_spPtr->distr()))->setAlpha(_alpha);
(static_cast<stochasticProcessSSRV*>(_spPtr))->setRateOfRate(_nu);
}
void ssrvDistanceSeqs2Tree::printSideInfo(ostream& out) const
{
out<<"Alpha: "<< _alpha <<" Nu: "<< _nu <<endl;
}
// non virtual
void ssrvDistanceSeqs2Tree::setSideInfo(const MDOUBLE alpha, MDOUBLE nu)
{
_alpha = alpha;
_nu = nu;
}
ssrvDistanceSeqs2Tree::alphaAndNu ssrvDistanceSeqs2Tree::getSideInfo() const
{
return alphaAndNu(_alpha, _nu);
}
|
