1 //
2 // phylotreemixlen.cpp
3 // iqtree
4 //
5 // Created by Minh Bui on 24/08/15.
6 //
7 //
8
9 #include "phylotreemixlen.h"
10 #include "phylonodemixlen.h"
11 #include "model/modelfactorymixlen.h"
12 #include "model/modelmixture.h"
13 #include "model/ratefree.h"
14 #include "utils/MPIHelper.h"
15
16 #ifdef USE_CPPOPTLIB
17 #include "cppoptlib/solver/newtondescentsolver.h"
18 #include "cppoptlib/solver/lbfgsbsolver.h"
19 #endif
20
PhyloTreeMixlen()21 PhyloTreeMixlen::PhyloTreeMixlen() : IQTree()
22 #ifdef USE_CPPOPTLIB
23 , cppoptlib::BoundedProblem<double>()
24 #endif
25 {
26 mixlen = 1;
27 cur_mixture = -1;
28 // relative_treelen = NULL;
29 initializing_mixlen = false;
30 }
31
PhyloTreeMixlen(Alignment * aln,int mixlen)32 PhyloTreeMixlen::PhyloTreeMixlen(Alignment *aln, int mixlen) : IQTree(aln)
33 #ifdef USE_CPPOPTLIB
34 , cppoptlib::BoundedProblem<double>(mixlen)
35 #endif
36 {
37 // cout << "Initializing heterotachy mixture branch lengths" << endl;
38 cur_mixture = -1;
39 // relative_treelen = NULL;
40 initializing_mixlen = false;
41 setMixlen(mixlen);
42 }
43
~PhyloTreeMixlen()44 PhyloTreeMixlen::~PhyloTreeMixlen() {
45 // if (relative_treelen)
46 // aligned_free(relative_treelen);
47 }
48
startCheckpoint()49 void PhyloTreeMixlen::startCheckpoint() {
50 if (mixlen > 0)
51 checkpoint->startStruct("PhyloTreeMixlen" + convertIntToString(getMixlen()));
52 else
53 PhyloTree::startCheckpoint();
54 }
55
saveCheckpoint()56 void PhyloTreeMixlen::saveCheckpoint() {
57 if (mixlen > 0) {
58 startCheckpoint();
59 if (this->relative_treelen.size() > 0) {
60 ASSERT(mixlen == this->relative_treelen.size());
61 double relative_treelen[mixlen];
62 for (int i = 0; i < mixlen; i++)
63 relative_treelen[i] = this->relative_treelen[i];
64 CKP_ARRAY_SAVE(mixlen, relative_treelen);
65 }
66 endCheckpoint();
67 }
68 IQTree::saveCheckpoint();
69 }
70
71 /**
72 restore object from the checkpoint
73 */
restoreCheckpoint()74 void PhyloTreeMixlen::restoreCheckpoint() {
75 if (mixlen > 0) {
76 startCheckpoint();
77 double relative_treelen[mixlen];
78 if (CKP_ARRAY_RESTORE(mixlen, relative_treelen)) {
79 this->relative_treelen.resize(mixlen);
80 for (int i = 0; i < mixlen; i++)
81 this->relative_treelen[i] = relative_treelen[i];
82 }
83 endCheckpoint();
84 }
85 IQTree::restoreCheckpoint();
86 if (!root) {
87 // if not success, try to restore from PhyloTree
88 int orig_mixlen = mixlen;
89 mixlen = 0;
90 PhyloTree::restoreCheckpoint();
91 mixlen = orig_mixlen;
92 }
93 }
94
newNode(int node_id,const char * node_name)95 Node* PhyloTreeMixlen::newNode(int node_id, const char* node_name) {
96 return (Node*) (new PhyloNodeMixlen(node_id, node_name));
97 }
98
newNode(int node_id,int node_name)99 Node* PhyloTreeMixlen::newNode(int node_id, int node_name) {
100 return (Node*) (new PhyloNodeMixlen(node_id, node_name));
101 }
102
setMixlen(int mixlen)103 void PhyloTreeMixlen::setMixlen(int mixlen) {
104 this->mixlen = mixlen;
105 }
106
readTreeString(const string & tree_string)107 void PhyloTreeMixlen::readTreeString(const string &tree_string) {
108 IQTree::readTreeString(tree_string);
109 treeLengths(relative_treelen);
110 if (mixlen > 0 && relative_treelen[0] == 0.0)
111 relative_treelen.clear();
112 }
113
initializeModel(Params & params,string model_name,ModelsBlock * models_block)114 void PhyloTreeMixlen::initializeModel(Params ¶ms, string model_name, ModelsBlock *models_block) {
115 try {
116 if (!getModelFactory()) {
117 setModelFactory(new ModelFactoryMixlen(params, model_name, this, models_block));
118 }
119 } catch (string & str) {
120 outError(str);
121 }
122 IQTree::initializeModel(params, model_name, models_block);
123 }
124
treeLengths(DoubleVector & lenvec,Node * node,Node * dad)125 void PhyloTreeMixlen::treeLengths(DoubleVector &lenvec, Node *node, Node *dad) {
126 if (lenvec.empty())
127 lenvec.resize(mixlen, 0.0);
128 if (!node) node = root;
129 FOR_NEIGHBOR_IT(node, dad, it) {
130 treeLengths(lenvec, (*it)->node, node);
131 for (int i = 0; i < mixlen; i++)
132 lenvec[i] += (*it)->getLength(i);
133 }
134 }
135
136
initializeMixBranches(PhyloNode * node,PhyloNode * dad)137 void PhyloTreeMixlen::initializeMixBranches(PhyloNode *node, PhyloNode *dad) {
138 if (!node) {
139 node = (PhyloNode*)root;
140 // exit if already initialized
141 // if (!((PhyloNeighborMixlen*)root->neighbors[0])->lengths.empty())
142 // return;
143 }
144 int i;
145 FOR_NEIGHBOR_IT(node, dad, it) {
146 // assign length of left branch
147 PhyloNeighborMixlen *nei = (PhyloNeighborMixlen*)(*it);
148 PhyloNeighborMixlen *back_nei = (PhyloNeighborMixlen*)((*it)->node->findNeighbor(node));
149 if (nei->lengths.empty()) {
150 // no branch lengths, initialize with relative rates
151 ASSERT(nei->length >= 0);
152 nei->lengths.resize(mixlen, nei->length);
153 back_nei->lengths.resize(mixlen, back_nei->length);
154 for (i = 0; i < mixlen; i++) {
155 nei->lengths[i] = back_nei->lengths[i] = max(params->min_branch_length, nei->length * relative_treelen[i]);
156 }
157 } else if (nei->lengths.size() > mixlen) {
158 // too many lengths, cut down
159 nei->lengths.resize(mixlen);
160 back_nei->lengths.resize(mixlen);
161 } else {
162 // too few lengths, add more
163 int cur = nei->lengths.size();
164 nei->lengths.resize(mixlen, nei->length);
165 back_nei->lengths.resize(mixlen, back_nei->length);
166 double avglen = 0.0;
167 for (i = 0; i < cur; i++)
168 avglen += nei->lengths[i];
169 avglen /= cur;
170 for (i = cur; i < mixlen; i++) {
171 nei->lengths[i] = back_nei->lengths[i] = max(params->min_branch_length, avglen * relative_treelen[i]);
172 }
173 }
174
175 double mean_len = 0.0;
176 for (int i = 0; i < mixlen; i++)
177 mean_len += nei->lengths[i] * site_rate->getProp(i);
178 // mean_len /= mixlen;
179 nei->length = back_nei->length = mean_len;
180
181 // recursive call
182 initializeMixBranches((PhyloNode*)(*it)->node, node);
183 }
184 }
185
assignMeanMixBranches(Node * node,Node * dad)186 void PhyloTreeMixlen::assignMeanMixBranches(Node *node, Node *dad) {
187 if (!node) node = root;
188 FOR_NEIGHBOR_IT(node, dad, it) {
189 PhyloNeighborMixlen *nei = (PhyloNeighborMixlen*)(*it);
190 double mean_len = 0.0;
191 for (int i = 0; i < nei->lengths.size(); i++)
192 mean_len += nei->lengths[i] * site_rate->getProp(i);
193 // mean_len /= nei->lengths.size();
194 nei->length = mean_len;
195
196 nei = (PhyloNeighborMixlen*)(*it)->node->findNeighbor(node);
197 nei->length = mean_len;
198 assignMeanMixBranches((*it)->node, node);
199 }
200 }
201
initializeMixlen(double tolerance,bool write_info)202 void PhyloTreeMixlen::initializeMixlen(double tolerance, bool write_info) {
203 // initialize mixture branch lengths if empty
204
205 if (initializing_mixlen)
206 return;
207
208 int i;
209
210 initializing_mixlen = true;
211
212 if (relative_treelen.empty()) {
213 RateHeterogeneity *saved_rate = getRate();
214 bool saved_fused_mix_rate = model_factory->fused_mix_rate;
215
216 // create new rate model
217 // random alpha
218 // relative_rate = new RateGamma(mixlen, params->gamma_shape, params->gamma_median, this);
219
220 string param;
221
222 if (getRate()->getFixParams()) {
223 stringstream ss;
224 for (i = 0; i < mixlen; i++) {
225 if (i > 0) ss << ",";
226 ss << getRate()->getProp(i);
227 }
228 param = ss.str();
229 }
230
231 RateFree *relative_rate = new RateFree(mixlen, params->gamma_shape, param, false, params->optimize_alg, this);
232 relative_rate->setTree(this);
233
234 // setup new rate model
235 setRate(relative_rate);
236 model_factory->site_rate = relative_rate;
237 if (getModel()->isMixture()) {
238 // model_factory->fused_mix_rate = true;
239 setLikelihoodKernel(sse);
240 }
241
242 // optimize rate model
243 double tree_lh = relative_rate->optimizeParameters(tolerance);
244
245 // model_factory->optimizeParameters(params->fixed_branch_length, false, tolerance);
246
247 // optimizeModelParameters();
248
249 // 2016-07-22: BUGFIX should rescale rates
250 double mean_rate = relative_rate->rescaleRates();
251 if (fabs(mean_rate-1.0) > 1e-6 && params->fixed_branch_length != BRLEN_FIX) {
252 scaleLength(mean_rate);
253 }
254 if (write_info) {
255 cout << "Initial LogL: " << curScore << ", ";
256 // if (verbose_mode >= VB_MED)
257 relative_rate->writeInfo(cout);
258 }
259
260 // make the rates more distinct
261 if (mixlen > 1 && relative_rate->getRate(0) / relative_rate->getRate(mixlen-1) > 0.9) {
262 cout << "Making the rates more distinct..." << endl;
263 relative_rate->setRate(0, relative_rate->getRate(0)*0.95);
264 relative_rate->setRate(mixlen-1, relative_rate->getRate(mixlen-1)*1.05);
265 }
266
267 double treelen = treeLength();
268 relative_treelen.resize(mixlen);
269 // 2017-12-21: BUG: moved this out of write_info if
270 // otherwise, relative_treelen is not initialized
271 for (i = 0; i < mixlen; i++)
272 relative_treelen[i] = treelen * relative_rate->getRate(i);
273
274 if (write_info) {
275 cout << "relative_treelen:";
276 for (i = 0; i < mixlen; i++) {
277 cout << " " << relative_treelen[i];
278 }
279 cout << endl;
280 }
281
282
283 // restore rate model
284 setRate(saved_rate);
285 model_factory->site_rate = saved_rate;
286 model_factory->fused_mix_rate = saved_fused_mix_rate;
287 setLikelihoodKernel(sse);
288
289 // set the weights of heterotachy model
290 double pinvar = site_rate->getPInvar();
291 if (!site_rate->getFixParams())
292 for (i = 0; i < mixlen; i++)
293 site_rate->setProp(i, relative_rate->getProp(i)*(1.0-pinvar));
294
295 delete relative_rate;
296 clearAllPartialLH();
297 }
298
299 if (((PhyloNeighborMixlen*)root->neighbors[0])->lengths.size() != mixlen) {
300 // assign branch length from rate model
301 DoubleVector saved_treelen = relative_treelen;
302 DoubleVector lenvec;
303 treeLengths(lenvec);
304 for (i = 0; i < mixlen; i++) {
305 relative_treelen[i] = relative_treelen[i] / lenvec[i];
306 }
307 if (verbose_mode >= VB_MED) {
308 cout << "relative_ratio:";
309 for (i = 0; i < mixlen; i++)
310 cout << " " << relative_treelen[i];
311 cout << endl;
312 }
313 initializeMixBranches();
314 clearAllPartialLH();
315 relative_treelen = saved_treelen;
316 }
317
318 initializing_mixlen = false;
319
320 }
321
fixOneNegativeBranch(double branch_length,Neighbor * dad_branch,Node * dad)322 void PhyloTreeMixlen::fixOneNegativeBranch(double branch_length, Neighbor *dad_branch, Node *dad) {
323 PhyloTree::fixOneNegativeBranch(branch_length, dad_branch, dad);
324 /*
325 PhyloNeighborMixlen *br = (PhyloNeighborMixlen*)dad_branch;
326 if (br->lengths.empty())
327 return;
328 int i;
329 for (i = 0; i < br->lengths.size(); i++)
330 br->lengths[i] = branch_length * relative_rate->getRate(i);
331 br = (PhyloNeighborMixlen*)dad_branch->node->findNeighbor(dad);
332 for (i = 0; i < br->lengths.size(); i++)
333 br->lengths[i] = branch_length * relative_rate->getRate(i);
334 */
335 }
336
337
optimizeOneBranch(PhyloNode * node1,PhyloNode * node2,bool clearLH,int maxNRStep)338 void PhyloTreeMixlen::optimizeOneBranch(PhyloNode *node1, PhyloNode *node2, bool clearLH, int maxNRStep) {
339
340 if (initializing_mixlen)
341 return PhyloTree::optimizeOneBranch(node1, node2, clearLH, maxNRStep);
342
343 current_it = (PhyloNeighbor*) node1->findNeighbor(node2);
344 ASSERT(current_it);
345 current_it_back = (PhyloNeighbor*) node2->findNeighbor(node1);
346 ASSERT(current_it_back);
347
348 int i;
349
350 theta_computed = false;
351
352 #ifdef USE_CPPOPTLIB
353 if (params->optimize_alg_mixlen.find("cppopt") != string::npos) {
354 //----- using cppoptlib ------//
355
356 TVector lower_bound(mixlen), upper_bound(mixlen), variables(mixlen);
357
358 // variables.resize(mixlen);
359 for (i = 0; i < mixlen; i++) {
360 lower_bound[i] = params->min_branch_length;
361 variables[i] = current_it->getLength(i);
362 upper_bound[i] = params->max_branch_length;
363 }
364
365 setBoxConstraint(lower_bound, upper_bound);
366
367 cppoptlib::NewtonDescentSolver<PhyloTreeMixlen> solver;
368 // cppoptlib::LbfgsbSolver<PhyloTreeMixlen> solver;
369 solver.minimize(*this, variables);
370 for (i = 0; i < mixlen; i++) {
371 current_it->setLength(i, variables[i]);
372 current_it_back->setLength(i, variables[i]);
373 }
374 } else
375 #endif
376
377 if (params->optimize_alg_mixlen.find("newton") != string::npos) {
378
379 //----- Newton-Raphson -----//
380 double lower_bound[mixlen];
381 double upper_bound[mixlen];
382 double variables[mixlen];
383 for (i = 0; i < mixlen; i++) {
384 lower_bound[i] = params->min_branch_length;
385 variables[i] = current_it->getLength(i);
386 upper_bound[i] = params->max_branch_length;
387 }
388
389 double score = minimizeNewtonMulti(lower_bound, variables, upper_bound, params->min_branch_length, mixlen);
390 for (i = 0; i < mixlen; i++) {
391 current_it->setLength(i, variables[i]);
392 current_it_back->setLength(i, variables[i]);
393 }
394 } else if (params->optimize_alg_mixlen.find("BFGS") != string::npos) {
395
396 // BFGS method to simultaneously optimize all lengths per branch
397 // It is often better than the true Newton method (Numerical Recipes in C++, chap. 10.7)
398 double variables[mixlen+1];
399 double upper_bound[mixlen+1];
400 double lower_bound[mixlen+1];
401 bool bound_check[mixlen+1];
402 for (i = 0; i < mixlen; i++) {
403 lower_bound[i+1] = params->min_branch_length;
404 variables[i+1] = current_it->getLength(i);
405 upper_bound[i+1] = params->max_branch_length;
406 bound_check[i+1] = false;
407 }
408
409 double grad[mixlen+1], hessian[mixlen*mixlen];
410 computeFuncDervMulti(variables+1, grad, hessian);
411 double score;
412 if (params->optimize_alg.find("BFGS-B") != string::npos)
413 score = -L_BFGS_B(mixlen, variables+1, lower_bound+1, upper_bound+1, params->min_branch_length);
414 else
415 score = -minimizeMultiDimen(variables, mixlen, lower_bound, upper_bound, bound_check, params->min_branch_length, hessian);
416
417 for (i = 0; i < mixlen; i++) {
418 current_it->setLength(i, variables[i+1]);
419 current_it_back->setLength(i, variables[i+1]);
420 }
421 if (verbose_mode >= VB_DEBUG) {
422 cout << "Mixlen-LnL: " << score << endl;
423 }
424
425 } else {
426
427 if (!model_factory->fused_mix_rate && getModel()->isMixture())
428 outError("Please use option -optlen BFGS to disable EM algorithm");
429
430 // EM algorithm
431 size_t ptn, c;
432 size_t nptn = aln->getNPattern();
433 size_t nmix = site_rate->getNRate();
434 ASSERT(nmix == mixlen);
435
436 // first compute _pattern_lh_cat
437 double tree_lh = -DBL_MAX;
438
439 // 2 steps are empirically determined to be best!
440 for (int EM_step = 0; EM_step < 2; EM_step++) {
441
442 double new_tree_lh = computePatternLhCat(WSL_RATECAT);
443 if (new_tree_lh+1.0 < tree_lh)
444 cout << "WARN: at EM step " << EM_step << " new_tree_lh " << new_tree_lh << " worse than tree_lh " << tree_lh << endl;
445 if (new_tree_lh-params->min_branch_length < tree_lh)
446 break;
447 tree_lh = new_tree_lh;
448
449 // E-step
450 // decoupled weights (prop) from _pattern_lh_cat to obtain L_ci and compute pattern likelihood L_i
451 for (ptn = 0; ptn < nptn; ptn++) {
452 double *this_lk_cat = _pattern_lh_cat + ptn*nmix;
453 double lk_ptn = ptn_invar[ptn];
454 for (c = 0; c < nmix; c++) {
455 lk_ptn += this_lk_cat[c];
456 }
457 ASSERT(lk_ptn != 0.0);
458 lk_ptn = ptn_freq[ptn] / lk_ptn;
459 // transform _pattern_lh_cat into posterior probabilities of each category
460 for (c = 0; c < nmix; c++) {
461 this_lk_cat[c] *= lk_ptn;
462 }
463
464 }
465
466 double negative_lh;
467 double optx;
468 theta_computed = false;
469 computePtnFreq();
470
471 for (cur_mixture = 0; cur_mixture < mixlen; cur_mixture++) {
472
473 double *this_lk_cat = _pattern_lh_cat+cur_mixture;
474 for (ptn = 0; ptn < nptn; ptn++)
475 ptn_freq[ptn] = this_lk_cat[ptn*nmix];
476
477 double current_len = current_it->getLength(cur_mixture);
478 ASSERT(current_len >= 0.0);
479 // Newton-Raphson method
480 optx = minimizeNewton(params->min_branch_length, current_len, params->max_branch_length, params->min_branch_length, negative_lh, maxNRStep);
481
482 current_it->setLength(cur_mixture, optx);
483 current_it_back->setLength(cur_mixture, optx);
484
485 }
486
487 cur_mixture = -1;
488 // reset ptn_freq
489 ptn_freq_computed = false;
490 computePtnFreq();
491 } // for EM_step
492 }
493
494 if (clearLH) {
495 node1->clearReversePartialLh(node2);
496 node2->clearReversePartialLh(node1);
497 }
498
499 }
500
501 /**
502 return the number of dimensions
503 */
getNDim()504 int PhyloTreeMixlen::getNDim() {
505 return mixlen;
506 }
507
508
509 /**
510 the target function which needs to be optimized
511 @param x the input vector x
512 @return the function value at x
513 */
targetFunk(double x[])514 double PhyloTreeMixlen::targetFunk(double x[]) {
515 int i;
516 for (i = 0; i < mixlen; i++) {
517 current_it->setLength(i, x[i+1]);
518 current_it_back->setLength(i, x[i+1]);
519 }
520
521 if (theta_computed)
522 return -computeLikelihoodFromBuffer();
523 else
524 return -computeLikelihoodBranch(current_it, (PhyloNode*)current_it_back->node);
525 }
526
derivativeFunk(double x[],double dfx[])527 double PhyloTreeMixlen::derivativeFunk(double x[], double dfx[]) {
528 int i;
529 // cout.precision(10);
530 // cout << "x: ";
531 for (i = 0; i < mixlen; i++) {
532 ASSERT(!std::isnan(x[i+1]));
533 current_it->setLength(i, x[i+1]);
534 current_it_back->setLength(i, x[i+1]);
535 // cout << " " << x[i+1];
536 }
537 // cout << endl;
538 double df[mixlen+1], ddf[mixlen*mixlen];
539 computeLikelihoodDerv(current_it, (PhyloNode*)current_it_back->node, df, ddf);
540 for (i = 0; i < mixlen; i++)
541 df[i] = -df[i];
542 memcpy(dfx+1, df, sizeof(double)*mixlen);
543 return -df[mixlen];
544 }
545
computeFuncDervMulti(double * value,double * df,double * ddf)546 void PhyloTreeMixlen::computeFuncDervMulti(double *value, double *df, double *ddf) {
547 int i;
548 for (i = 0; i < mixlen; i++) {
549 current_it->setLength(i, value[i]);
550 current_it_back->setLength(i, value[i]);
551 }
552 computeLikelihoodDerv(current_it, (PhyloNode*)current_it_back->node, df, ddf);
553
554 // last element of df is the tree log-ikelihood
555 for (i = 0; i <= mixlen; i++) {
556 df[i] = -df[i];
557 }
558 int mixlen2 = mixlen * mixlen;
559 for (i = 0; i < mixlen2; i++)
560 ddf[i] = -ddf[i];
561 }
562
optimizeAllBranches(int my_iterations,double tolerance,int maxNRStep)563 double PhyloTreeMixlen::optimizeAllBranches(int my_iterations, double tolerance, int maxNRStep) {
564
565 initializeMixlen(tolerance, false);
566 clearAllPartialLH();
567
568 double tree_lh = PhyloTree::optimizeAllBranches(my_iterations, tolerance, maxNRStep);
569
570 if (!initializing_mixlen)
571 assignMeanMixBranches();
572
573 return tree_lh;
574 }
575
optimizeNNI(bool speedNNI)576 pair<int, int> PhyloTreeMixlen::optimizeNNI(bool speedNNI) {
577 int i, j;
578
579 DoubleVector meanlenvec;
580 treeLengths(meanlenvec);
581 // compute mean branch length
582 for (j = 0; j < mixlen; j++)
583 meanlenvec[j] /= (branchNum);
584
585 // scan over all branches and fix short/long branches
586 NodeVector nodes1, nodes2;
587 getBranches(nodes1, nodes2);
588 int num_fixed = 0;
589 for (i = 0; i < nodes1.size(); i++) {
590 PhyloNeighborMixlen* nei = (PhyloNeighborMixlen*)nodes1[i]->findNeighbor(nodes2[i]);
591 PhyloNeighborMixlen* nei_back = (PhyloNeighborMixlen*)nodes2[i]->findNeighbor(nodes1[i]);
592 for (j = 0; j < mixlen; j++)
593 // if (nei->lengths[j] < params->min_branch_length*2.0 || nei->lengths[j] > params->max_branch_length*0.9) {
594 if (nei->lengths[j] > params->max_branch_length*0.9) {
595 // if too long or too short branch, assign with mean branch length
596 nei->lengths[j] = nei_back->lengths[j] = meanlenvec[j];
597 num_fixed++;
598 }
599 }
600 if (num_fixed > 0)
601 optimizeBranches(num_fixed);
602
603 return IQTree::optimizeNNI(speedNNI);
604 }
605
printBranchLength(ostream & out,int brtype,bool print_slash,Neighbor * length_nei)606 void PhyloTreeMixlen::printBranchLength(ostream &out, int brtype, bool print_slash, Neighbor *length_nei) {
607 if (((PhyloNeighborMixlen*)length_nei)->lengths.empty())
608 return PhyloTree::printBranchLength(out, brtype, print_slash, length_nei);
609
610 if ((brtype & (WT_BR_LEN+WT_BR_SCALE)) == 0)
611 return;
612
613 PhyloNeighborMixlen *nei = (PhyloNeighborMixlen*) length_nei;
614
615 if (cur_mixture == -1) {
616 // print mixture branch lengths
617 out << "[";
618 for (int i = 0; i < mixlen; i++) {
619 if (i > 0) out << BRANCH_LENGTH_SEPARATOR;
620 double length = nei->lengths[i];
621 if (brtype & WT_BR_SCALE) length *= len_scale;
622 if (brtype & WT_BR_LEN_ROUNDING) length = round(length);
623 if (brtype & WT_BR_LEN) {
624 if (brtype & WT_BR_LEN_FIXED_WIDTH)
625 out << fixed << length;
626 else
627 out << length;
628 } else if (brtype & WT_BR_CLADE && length_nei->node->name != ROOT_NAME) {
629 out << length;
630 }
631 }
632 out << "]";
633 }
634
635 if (brtype & WT_BR_LEN)
636 out << ":";
637 else if ((brtype & WT_BR_CLADE) && print_slash && length_nei->node->name != ROOT_NAME)
638 out << "/";
639
640 double length = nei->length;
641
642 if (cur_mixture >= 0) {
643 // print branch length of a mixture component only!
644 length = nei->lengths[cur_mixture];
645 }
646
647 if (brtype & WT_BR_SCALE) length *= len_scale;
648 if (brtype & WT_BR_LEN_ROUNDING) length = round(length);
649 if (brtype & WT_BR_LEN) {
650 if (brtype & WT_BR_LEN_FIXED_WIDTH)
651 out << fixed << length;
652 else
653 out << length;
654 } else if (brtype & WT_BR_CLADE && length_nei->node->name != ROOT_NAME) {
655 out << length;
656 }
657 }
658
printResultTree(string suffix)659 void PhyloTreeMixlen::printResultTree(string suffix) {
660 if (MPIHelper::getInstance().isWorker()) {
661 return;
662 }
663 if (params->suppress_output_flags & OUT_TREEFILE)
664 return;
665 setRootNode(params->root);
666 string tree_file_name = params->out_prefix;
667 tree_file_name += ".treefile";
668 if (suffix.compare("") != 0) {
669 tree_file_name += "." + suffix;
670 }
671
672 try {
673 ofstream out;
674 out.exceptions(ios::failbit | ios::badbit);
675 out.open(tree_file_name.c_str());
676 cur_mixture = -1;
677 printTree(out, WT_BR_LEN | WT_BR_LEN_FIXED_WIDTH | WT_SORT_TAXA | WT_NEWLINE);
678 for (cur_mixture = 0; cur_mixture < mixlen; cur_mixture++) {
679 //out << "[Heterotachy class " << cur_mixture+1 << "]" << endl;
680 printTree(out, WT_BR_LEN | WT_BR_LEN_FIXED_WIDTH | WT_SORT_TAXA | WT_NEWLINE);
681 }
682 cur_mixture = -1;
683 out.close();
684 } catch (ios::failure) {
685 outError(ERR_WRITE_OUTPUT, tree_file_name);
686 }
687
688 if (verbose_mode >= VB_MED)
689 cout << "Best tree printed to " << tree_file_name << endl;
690 }
691
692
693 /*************** Using cppoptlib for branch length optimization ***********/
694
695 #ifdef USE_CPPOPTLIB
value(const TVector & x)696 double PhyloTreeMixlen::value(const TVector &x) {
697 double xx[mixlen+1];
698 for (int i = 0; i < mixlen; i++)
699 xx[i+1] = x(i);
700 return targetFunk(xx);
701 }
702
gradient(const TVector & x,TVector & grad)703 void PhyloTreeMixlen::gradient(const TVector &x, TVector &grad) {
704 int i;
705 double xx[mixlen];
706 for (i = 0; i < mixlen; i++)
707 xx[i] = x(i);
708 double df[mixlen+1], ddf[mixlen*mixlen];
709 computeFuncDervMulti(xx, df, ddf);
710 for (i = 0; i < mixlen; i++)
711 grad(i) = df[i];
712 }
713
hessian(const TVector & x,THessian & hessian)714 void PhyloTreeMixlen::hessian(const TVector &x, THessian &hessian) {
715 int i, j;
716 double xx[mixlen];
717 for (i = 0; i < mixlen; i++)
718 xx[i] = x(i);
719 int mixlen2 = mixlen*mixlen;
720 double df[mixlen+1], ddf[mixlen2];
721 computeFuncDervMulti(xx, df, ddf);
722
723 for (i = 0; i < mixlen; i++)
724 for (j = 0; j < mixlen; j++)
725 hessian(i, j) = ddf[i*mixlen+j];
726 }
727 #endif
728
729 // defining log-likelihood derivative function for EM algorithm
computeFuncDerv(double value,double & df,double & ddf)730 void PhyloTreeMixlen::computeFuncDerv(double value, double &df, double &ddf) {
731
732 if (cur_mixture < 0)
733 return PhyloTree::computeFuncDerv(value, df, ddf);
734
735 current_it->setLength(cur_mixture, value);
736 current_it_back->setLength(cur_mixture, value);
737
738 (this->*computeLikelihoodDervMixlenPointer)(current_it, (PhyloNode*) current_it_back->node, df, ddf);
739
740 df = -df;
741 ddf = -ddf;
742 return;
743
744
745 PhyloNeighbor* dad_branch = current_it;
746 PhyloNode *dad = (PhyloNode*) current_it_back->node;
747
748 PhyloNode *node = (PhyloNode*) dad_branch->node;
749 PhyloNeighbor *node_branch = (PhyloNeighbor*) node->findNeighbor(dad);
750 if (!central_partial_lh)
751 initializeAllPartialLh();
752 if (node->isLeaf()) {
753 PhyloNode *tmp_node = dad;
754 dad = node;
755 node = tmp_node;
756 PhyloNeighbor *tmp_nei = dad_branch;
757 dad_branch = node_branch;
758 node_branch = tmp_nei;
759 }
760
761 ASSERT((dad_branch->partial_lh_computed & 1) || node->isLeaf());
762 ASSERT((node_branch->partial_lh_computed & 1) || dad->isLeaf());
763 // if ((dad_branch->partial_lh_computed & 1) == 0)
764 // computePartialLikelihood(dad_branch, dad);
765 // if ((node_branch->partial_lh_computed & 1) == 0)
766 // computePartialLikelihood(node_branch, node);
767 size_t nstates = aln->num_states;
768 size_t ncat = site_rate->getNRate();
769 size_t nmixture = model->getNMixtures();
770
771 size_t block = ncat * nstates * nmixture;
772 size_t statemix = nstates * nmixture;
773 size_t statecat = nstates * ncat;
774 size_t ptn; // for big data size > 4GB memory required
775 size_t c, i, m = cur_mixture;
776 size_t orig_nptn = aln->size();
777 size_t nptn = aln->size()+model_factory->unobserved_ptns.size();
778 size_t maxptn = get_safe_upper_limit(nptn);
779 double *eval = model->getEigenvalues();
780 ASSERT(eval);
781
782 ASSERT(theta_all);
783 if (!theta_computed) {
784 // precompute theta for fast branch length optimization
785
786 if (dad->isLeaf()) {
787 // special treatment for TIP-INTERNAL NODE case
788 #ifdef _OPENMP
789 #pragma omp parallel for private(ptn, i, m)
790 #endif
791 for (ptn = 0; ptn < nptn; ptn++) {
792 double *partial_lh_dad = dad_branch->partial_lh + ptn*block;
793 double *theta = theta_all + ptn*block;
794
795 // TODO: check with vectorclass!
796 double *lh_tip = tip_partial_lh +
797 ((int)((ptn < orig_nptn) ? (aln->at(ptn))[dad->id] : model_factory->unobserved_ptns[ptn-orig_nptn][dad->id]))*statemix;
798 for (m = 0; m < nmixture; m++) {
799 for (i = 0; i < statecat; i++) {
800 theta[m*statecat+i] = lh_tip[m*nstates + i%nstates] * partial_lh_dad[m*statecat+i];
801 }
802 }
803
804 }
805 // ascertainment bias correction
806 } else {
807 // both dad and node are internal nodes
808 double *partial_lh_node = node_branch->partial_lh;
809 double *partial_lh_dad = dad_branch->partial_lh;
810
811 size_t all_entries = nptn*block;
812 #ifdef _OPENMP
813 #pragma omp parallel for
814 #endif
815 for (i = 0; i < all_entries; i++) {
816 theta_all[i] = partial_lh_node[i] * partial_lh_dad[i];
817 }
818 }
819 if (nptn < maxptn) {
820 // copy dummy values
821 for (ptn = nptn; ptn < maxptn; ptn++)
822 memcpy(&theta_all[ptn*block], &theta_all[(ptn-1)*block], block*sizeof(double));
823 }
824 theta_computed = true;
825 }
826
827 double *val0 = new double[statecat];
828 double *val1 = new double[statecat];
829 double *val2 = new double[statecat];
830 for (c = 0; c < ncat; c++) {
831 double prop = site_rate->getProp(c);
832 for (i = 0; i < nstates; i++) {
833 double cof = eval[cur_mixture*nstates+i]*site_rate->getRate(c);
834 // length for heterotachy model
835 double val = exp(cof*dad_branch->getLength(cur_mixture)) * prop * model->getMixtureWeight(cur_mixture);
836 double val1_ = cof*val;
837 val0[(c)*nstates+i] = val;
838 val1[(c)*nstates+i] = val1_;
839 val2[(c)*nstates+i] = cof*val1_;
840 }
841 }
842
843
844 double my_df = 0.0, my_ddf = 0.0, prob_const = 0.0, df_const = 0.0, ddf_const = 0.0;
845
846 #ifdef _OPENMP
847 #pragma omp parallel for reduction(+: my_df, my_ddf, prob_const, df_const, ddf_const) private(ptn, i)
848 #endif
849 for (ptn = 0; ptn < nptn; ptn++) {
850 double lh_ptn = ptn_invar[ptn], df_ptn = 0.0, ddf_ptn = 0.0;
851 double *theta = theta_all + ptn*block + cur_mixture*statecat;
852 for (i = 0; i < statecat; i++) {
853 lh_ptn += val0[i] * theta[i];
854 df_ptn += val1[i] * theta[i];
855 ddf_ptn += val2[i] * theta[i];
856 }
857
858 // assert(lh_ptn > 0.0);
859 lh_ptn = fabs(lh_ptn);
860
861 if (ptn < orig_nptn) {
862 double df_frac = df_ptn / lh_ptn;
863 double ddf_frac = ddf_ptn / lh_ptn;
864 double freq = ptn_freq[ptn];
865 double tmp1 = df_frac * freq;
866 double tmp2 = ddf_frac * freq;
867 my_df += tmp1;
868 my_ddf += tmp2 - tmp1 * df_frac;
869 } else {
870 // ascertainment bias correction
871 prob_const += lh_ptn;
872 df_const += df_ptn;
873 ddf_const += ddf_ptn;
874 }
875 }
876 df = my_df;
877 ddf = my_ddf;
878 if (std::isnan(df) || std::isinf(df)) {
879 df = 0.0;
880 ddf = 0.0;
881 // outWarning("Numerical instability (some site-likelihood = 0)");
882 }
883
884 if (orig_nptn < nptn) {
885 // ascertainment bias correction
886 prob_const = 1.0 - prob_const;
887 double df_frac = df_const / prob_const;
888 double ddf_frac = ddf_const / prob_const;
889 int nsites = aln->getNSite();
890 df += nsites * df_frac;
891 ddf += nsites *(ddf_frac + df_frac*df_frac);
892 }
893
894
895 delete [] val2;
896 delete [] val1;
897 delete [] val0;
898
899
900 df = -df;
901 ddf = -ddf;
902
903 // return lh;
904 }
905
906