bpp-phyl  2.4.0
OptimizationTools.cpp
Go to the documentation of this file.
1 //
2 // File: OptimizationTools.cpp
3 // Created by: Julien Dutheil
4 // Created on: Sun Dec 14 09:43:32 2003
5 //
6 
7 /*
8  Copyright or © or Copr. Bio++ Development Team, (November 16, 2004)
9 
10  This software is a computer program whose purpose is to provide classes
11  for phylogenetic data analysis.
12 
13  This software is governed by the CeCILL license under French law and
14  abiding by the rules of distribution of free software. You can use,
15  modify and/ or redistribute the software under the terms of the CeCILL
16  license as circulated by CEA, CNRS and INRIA at the following URL
17  "http://www.cecill.info".
18 
19  As a counterpart to the access to the source code and rights to copy,
20  modify and redistribute granted by the license, users are provided only
21  with a limited warranty and the software's author, the holder of the
22  economic rights, and the successive licensors have only limited
23  liability.
24 
25  In this respect, the user's attention is drawn to the risks associated
26  with loading, using, modifying and/or developing or reproducing the
27  software by the user in light of its specific status of free software,
28  that may mean that it is complicated to manipulate, and that also
29  therefore means that it is reserved for developers and experienced
30  professionals having in-depth computer knowledge. Users are therefore
31  encouraged to load and test the software's suitability as regards their
32  requirements in conditions enabling the security of their systems and/or
33  data to be ensured and, more generally, to use and operate it in the
34  same conditions as regards security.
35 
36  The fact that you are presently reading this means that you have had
37  knowledge of the CeCILL license and that you accept its terms.
38  */
39 
40 #include "OptimizationTools.h"
43 #include "NNISearchable.h"
44 #include "NNITopologySearch.h"
45 #include "Io/Newick.h"
46 
55 
56 // From bpp-seq:
57 #include <Bpp/Seq/Io/Fasta.h>
58 
59 using namespace bpp;
60 using namespace std;
61 
62 /******************************************************************************/
63 
65 
67 
68 /******************************************************************************/
69 
70 std::string OptimizationTools::OPTIMIZATION_NEWTON = "newton";
71 std::string OptimizationTools::OPTIMIZATION_GRADIENT = "gradient";
72 std::string OptimizationTools::OPTIMIZATION_BRENT = "Brent";
73 std::string OptimizationTools::OPTIMIZATION_BFGS = "BFGS";
74 
75 /******************************************************************************/
76 
78  tl_(tl),
79  brLen_(),
80  lambda_()
81 {
82  // We work only on the branch lengths:
84  if (brLen_.hasParameter("RootPosition"))
85  brLen_.deleteParameter("RootPosition");
86  lambda_.addParameter(Parameter("scale factor", 0));
87 }
88 
90 
92 {
93  if (lambda.size() != 1)
94  throw Exception("OptimizationTools::ScaleFunction::f(). This is a one parameter function!");
96 }
97 
99 {
100  // Scale the tree:
101  ParameterList brLen = brLen_;
102  double s = exp(lambda_[0].getValue());
103  for (unsigned int i = 0; i < brLen.size(); i++)
104  {
105  try
106  {
107  brLen[i].setValue(brLen[i].getValue() * s);
108  }
109  catch (ConstraintException& cex)
110  {
111  // Do nothing. Branch value is already at bound...
112  }
113  }
114  return tl_->f(brLen);
115 }
116 
117 /******************************************************************************/
118 
120  TreeLikelihood* tl,
121  double tolerance,
122  unsigned int tlEvalMax,
123  OutputStream* messageHandler,
124  OutputStream* profiler,
125  unsigned int verbose)
126 {
127  ScaleFunction sf(tl);
128  BrentOneDimension bod(&sf);
129  bod.setMessageHandler(messageHandler);
130  bod.setProfiler(profiler);
131  ParameterList singleParameter;
132  singleParameter.addParameter(Parameter("scale factor", 0));
133  bod.setInitialInterval(-0.5, 0.5);
134  bod.init(singleParameter);
135  ParametersStopCondition PS(&bod, tolerance);
136  bod.setStopCondition(PS);
137  bod.setMaximumNumberOfEvaluations(tlEvalMax);
138  bod.optimize();
140  if (verbose > 0)
141  ApplicationTools::displayResult("Tree scaled by", exp(sf.getParameters()[0].getValue()));
142  return bod.getNumberOfEvaluations();
143 }
144 
145 /******************************************************************************/
146 
149  const ParameterList& parameters,
150  OptimizationListener* listener,
151  unsigned int nstep,
152  double tolerance,
153  unsigned int tlEvalMax,
154  OutputStream* messageHandler,
155  OutputStream* profiler,
156  bool reparametrization,
157  unsigned int verbose,
158  const std::string& optMethodDeriv,
159  const std::string& optMethodModel)
160 {
161  DerivableSecondOrder* f = tl;
162  ParameterList pl = parameters;
163 
164  // Shall we reparametrize the function to remove constraints?
165  unique_ptr<DerivableSecondOrder> frep;
166  if (reparametrization)
167  {
168  frep.reset(new ReparametrizationDerivableSecondOrderWrapper(f, parameters));
169  f = frep.get();
170 
171  // Reset parameters to remove constraints:
172  pl = f->getParameters().subList(parameters.getParameterNames());
173  }
174 
175  // ///////////////
176  // Build optimizer:
177 
178  // Branch lengths
179 
181  MetaOptimizer* poptimizer = 0;
183 
184  if (optMethodDeriv == OPTIMIZATION_GRADIENT)
186  else if (optMethodDeriv == OPTIMIZATION_NEWTON)
188  else if (optMethodDeriv == OPTIMIZATION_BFGS)
190  else
191  throw Exception("OptimizationTools::optimizeNumericalParameters. Unknown optimization method: " + optMethodDeriv);
192 
193  // Other parameters
194 
195  if (optMethodModel == OPTIMIZATION_BRENT)
196  {
198  desc->addOptimizer("Substitution model parameter", new SimpleMultiDimensions(f), plsm.getParameterNames(), 0, MetaOptimizerInfos::IT_TYPE_STEP);
199 
200 
202  desc->addOptimizer("Rate distribution parameter", new SimpleMultiDimensions(f), plrd.getParameterNames(), 0, MetaOptimizerInfos::IT_TYPE_STEP);
203  poptimizer = new MetaOptimizer(f, desc, nstep);
204  }
205  else if (optMethodModel == OPTIMIZATION_BFGS)
206  {
207  vector<string> vNameDer;
208 
210  vNameDer = plsm.getParameterNames();
211 
213 
214  vector<string> vNameDer2 = plrd.getParameterNames();
215 
216  vNameDer.insert(vNameDer.begin(), vNameDer2.begin(), vNameDer2.end());
217  fnum->setParametersToDerivate(vNameDer);
218 
219  desc->addOptimizer("Rate & model distribution parameters", new BfgsMultiDimensions(fnum), vNameDer, 1, MetaOptimizerInfos::IT_TYPE_FULL);
220  poptimizer = new MetaOptimizer(fnum, desc, nstep);
221  }
222  else
223  throw Exception("OptimizationTools::optimizeNumericalParameters. Unknown optimization method: " + optMethodModel);
224 
225  poptimizer->setVerbose(verbose);
226  poptimizer->setProfiler(profiler);
227  poptimizer->setMessageHandler(messageHandler);
228  poptimizer->setMaximumNumberOfEvaluations(tlEvalMax);
229  poptimizer->getStopCondition()->setTolerance(tolerance);
230 
231  // Optimize TreeLikelihood function:
233  NaNListener* nanListener = new NaNListener(poptimizer, tl);
234  poptimizer->addOptimizationListener(nanListener);
235  if (listener)
236  poptimizer->addOptimizationListener(listener);
237  poptimizer->init(pl);
238  poptimizer->optimize();
239 
240  if (verbose > 0)
242 
243  // We're done.
244  unsigned int nb = poptimizer->getNumberOfEvaluations();
245  delete poptimizer;
246  return nb;
247 }
248 
249 /******************************************************************************/
250 
253  const ParameterList& parameters,
254  OptimizationListener* listener,
255  double tolerance,
256  unsigned int tlEvalMax,
257  OutputStream* messageHandler,
258  OutputStream* profiler,
259  bool reparametrization,
260  bool useClock,
261  unsigned int verbose,
262  const std::string& optMethodDeriv)
263 {
264  DerivableSecondOrder* f = tl;
265  ParameterList pl = parameters;
266  // Shall we use a molecular clock constraint on branch lengths?
267  unique_ptr<GlobalClockTreeLikelihoodFunctionWrapper> fclock;
268  if (useClock)
269  {
270  fclock.reset(new GlobalClockTreeLikelihoodFunctionWrapper(tl));
271  f = fclock.get();
272  if (verbose > 0)
273  ApplicationTools::displayResult("Log-likelihood after adding clock", -tl->getLogLikelihood());
274 
275  // Reset parameters to use new branch lengths. WARNING! 'old' branch parameters do not exist anymore and have been replaced by heights
276  pl = fclock->getParameters().getCommonParametersWith(parameters);
277  pl.addParameters(fclock->getHeightParameters());
278  }
279  // Shall we reparametrize the function to remove constraints?
280  unique_ptr<DerivableSecondOrder> frep;
281  if (reparametrization)
282  {
283  frep.reset(new ReparametrizationDerivableSecondOrderWrapper(f, pl));
284  f = frep.get();
285 
286  // Reset parameters to remove constraints:
287  pl = f->getParameters().subList(pl.getParameterNames());
288  }
289 
290  unique_ptr<AbstractNumericalDerivative> fnum;
291  // Build optimizer:
292  unique_ptr<Optimizer> optimizer;
293  if (optMethodDeriv == OPTIMIZATION_GRADIENT)
294  {
295  fnum.reset(new TwoPointsNumericalDerivative(f));
296  fnum->setInterval(0.0000001);
297  optimizer.reset(new ConjugateGradientMultiDimensions(fnum.get()));
298  }
299  else if (optMethodDeriv == OPTIMIZATION_NEWTON)
300  {
301  fnum.reset(new ThreePointsNumericalDerivative(f));
302  fnum->setInterval(0.0001);
303  optimizer.reset(new PseudoNewtonOptimizer(fnum.get()));
304  }
305  else if (optMethodDeriv == OPTIMIZATION_BFGS)
306  {
307  fnum.reset(new TwoPointsNumericalDerivative(f));
308  fnum->setInterval(0.0001);
309  optimizer.reset(new BfgsMultiDimensions(fnum.get()));
310  }
311  else
312  throw Exception("OptimizationTools::optimizeNumericalParameters2. Unknown optimization method: " + optMethodDeriv);
313 
314  // Numerical derivatives:
316  if (useClock)
317  tmp.addParameters(fclock->getHeightParameters());
319  optimizer->setVerbose(verbose);
320  optimizer->setProfiler(profiler);
321  optimizer->setMessageHandler(messageHandler);
322  optimizer->setMaximumNumberOfEvaluations(tlEvalMax);
323  optimizer->getStopCondition()->setTolerance(tolerance);
324 
325  // Optimize TreeLikelihood function:
327  NaNListener* nanListener = new NaNListener(optimizer.get(), tl);
328  optimizer->addOptimizationListener(nanListener);
329  if (listener)
330  optimizer->addOptimizationListener(listener);
331  optimizer->init(pl);
332  optimizer->optimize();
333 
334  if (verbose > 0)
336 
337  // We're done.
338  return optimizer->getNumberOfEvaluations();
339 }
340 
341 /******************************************************************************/
342 
345  const ParameterList& parameters,
346  OptimizationListener* listener,
347  double tolerance,
348  unsigned int tlEvalMax,
349  OutputStream* messageHandler,
350  OutputStream* profiler,
351  unsigned int verbose,
352  const std::string& optMethodDeriv)
353 {
354  // Build optimizer:
355  Optimizer* optimizer = 0;
356  if (optMethodDeriv == OPTIMIZATION_GRADIENT)
357  {
358  tl->enableFirstOrderDerivatives(true);
359  tl->enableSecondOrderDerivatives(false);
360  optimizer = new ConjugateGradientMultiDimensions(tl);
361  }
362  else if (optMethodDeriv == OPTIMIZATION_NEWTON)
363  {
364  tl->enableFirstOrderDerivatives(true);
366  optimizer = new PseudoNewtonOptimizer(tl);
367  }
368  else if (optMethodDeriv == OPTIMIZATION_BFGS)
369  {
370  tl->enableFirstOrderDerivatives(true);
371  tl->enableSecondOrderDerivatives(false);
372  optimizer = new BfgsMultiDimensions(tl);
373  }
374  else
375  throw Exception("OptimizationTools::optimizeBranchLengthsParameters. Unknown optimization method: " + optMethodDeriv);
376  optimizer->setVerbose(verbose);
377  optimizer->setProfiler(profiler);
378  optimizer->setMessageHandler(messageHandler);
379  optimizer->setMaximumNumberOfEvaluations(tlEvalMax);
380  optimizer->getStopCondition()->setTolerance(tolerance);
381 
382  // Optimize TreeLikelihood function:
385  if (listener)
386  optimizer->addOptimizationListener(listener);
387  optimizer->init(pl);
388  optimizer->optimize();
389  if (verbose > 0)
391 
392  // We're done.
393  unsigned int n = optimizer->getNumberOfEvaluations();
394  delete optimizer;
395  return n;
396 }
397 
398 /******************************************************************************/
399 
402  const ParameterList& parameters,
403  OptimizationListener* listener,
404  unsigned int nstep,
405  double tolerance,
406  unsigned int tlEvalMax,
407  OutputStream* messageHandler,
408  OutputStream* profiler,
409  unsigned int verbose,
410  const std::string& optMethodDeriv)
411 {
413 
414  // Build optimizer:
416  if (optMethodDeriv == OPTIMIZATION_GRADIENT)
417  {
418  fun = new TwoPointsNumericalDerivative(cl);
419  fun->setInterval(0.0000001);
421  }
422  else if (optMethodDeriv == OPTIMIZATION_NEWTON)
423  {
424  fun = new ThreePointsNumericalDerivative(cl);
425  fun->setInterval(0.0001);
427  }
428  else
429  throw Exception("OptimizationTools::optimizeNumericalParametersWithGlobalClock. Unknown optimization method: " + optMethodDeriv);
430 
431  // Numerical derivatives:
434 
436  if (plsm.size() < 10)
437  desc->addOptimizer("Substitution model parameter", new SimpleMultiDimensions(cl), plsm.getParameterNames(), 0, MetaOptimizerInfos::IT_TYPE_STEP);
438  else
439  desc->addOptimizer("Substitution model parameters", new DownhillSimplexMethod(cl), plsm.getParameterNames(), 0, MetaOptimizerInfos::IT_TYPE_FULL);
440 
442  if (plrd.size() < 10)
443  desc->addOptimizer("Rate distribution parameter", new SimpleMultiDimensions(cl), plrd.getParameterNames(), 0, MetaOptimizerInfos::IT_TYPE_STEP);
444  else
445  desc->addOptimizer("Rate dsitribution parameters", new DownhillSimplexMethod(cl), plrd.getParameterNames(), 0, MetaOptimizerInfos::IT_TYPE_FULL);
446 
447  MetaOptimizer optimizer(fun, desc, nstep);
448  optimizer.setVerbose(verbose);
449  optimizer.setProfiler(profiler);
450  optimizer.setMessageHandler(messageHandler);
451  optimizer.setMaximumNumberOfEvaluations(tlEvalMax);
452  optimizer.getStopCondition()->setTolerance(tolerance);
453 
454  // Optimize TreeLikelihood function:
456  if (listener)
457  optimizer.addOptimizationListener(listener);
458  optimizer.init(parameters);
459  optimizer.optimize();
460  if (verbose > 0)
462 
463  // We're done.
464  return optimizer.getNumberOfEvaluations();
465 }
466 
467 /******************************************************************************/
468 
471  const ParameterList& parameters,
472  OptimizationListener* listener,
473  double tolerance,
474  unsigned int tlEvalMax,
475  OutputStream* messageHandler,
476  OutputStream* profiler,
477  unsigned int verbose,
478  const std::string& optMethodDeriv)
479 {
481 
482  // Build optimizer:
483  Optimizer* optimizer = 0;
484  if (optMethodDeriv == OPTIMIZATION_GRADIENT)
485  {
486  fun = new TwoPointsNumericalDerivative(cl);
487  fun->setInterval(0.0000001);
488  optimizer = new ConjugateGradientMultiDimensions(fun);
489  }
490  else if (optMethodDeriv == OPTIMIZATION_NEWTON)
491  {
492  fun = new ThreePointsNumericalDerivative(cl);
493  fun->setInterval(0.0001);
494  optimizer = new PseudoNewtonOptimizer(fun);
495  }
496  else
497  throw Exception("OptimizationTools::optimizeBranchLengthsParameters. Unknown optimization method: " + optMethodDeriv);
498 
499  // Numerical derivatives:
500  ParameterList tmp = parameters.getCommonParametersWith(cl->getParameters());
502 
503  optimizer->setVerbose(verbose);
504  optimizer->setProfiler(profiler);
505  optimizer->setMessageHandler(messageHandler);
506  optimizer->setMaximumNumberOfEvaluations(tlEvalMax);
507  optimizer->getStopCondition()->setTolerance(tolerance);
508 
509  // Optimize TreeLikelihood function:
511  if (listener)
512  optimizer->addOptimizationListener(listener);
513  optimizer->init(parameters);
514  optimizer->optimize();
515  if (verbose > 0)
517 
518  // We're done.
519  unsigned int n = optimizer->getNumberOfEvaluations();
520  delete optimizer;
521 
522  // We're done.
523  return n;
524 }
525 
526 /******************************************************************************/
527 
529 {
530  optimizeCounter_++;
531  if (optimizeCounter_ == optimizeNumerical_)
532  {
533  DiscreteRatesAcrossSitesTreeLikelihood* likelihood = dynamic_cast<DiscreteRatesAcrossSitesTreeLikelihood*>(topoSearch_->getSearchableObject());
535  OptimizationTools::optimizeNumericalParameters(likelihood, parameters_, 0, nStep_, tolerance_, 1000000, messenger_, profiler_, reparametrization_, verbose_, optMethod_);
536  optimizeCounter_ = 0;
537  }
538 }
539 
540 /******************************************************************************/
541 
543 {
544  optimizeCounter_++;
545  if (optimizeCounter_ == optimizeNumerical_)
546  {
547  DiscreteRatesAcrossSitesTreeLikelihood* likelihood = dynamic_cast<DiscreteRatesAcrossSitesTreeLikelihood*>(topoSearch_->getSearchableObject());
549  OptimizationTools::optimizeNumericalParameters2(likelihood, parameters_, 0, tolerance_, 1000000, messenger_, profiler_, reparametrization_, false, verbose_, optMethod_);
550  optimizeCounter_ = 0;
551  }
552 }
553 
554 // ******************************************************************************/
555 
558  const ParameterList& parameters,
559  bool optimizeNumFirst,
560  double tolBefore,
561  double tolDuring,
562  unsigned int tlEvalMax,
563  unsigned int numStep,
564  OutputStream* messageHandler,
565  OutputStream* profiler,
566  bool reparametrization,
567  unsigned int verbose,
568  const std::string& optMethodDeriv,
569  unsigned int nStep,
570  const std::string& nniMethod)
571 {
572  // Roughly optimize parameter
573  if (optimizeNumFirst)
574  {
575  OptimizationTools::optimizeNumericalParameters(tl, parameters, NULL, nStep, tolBefore, 1000000, messageHandler, profiler, reparametrization, verbose, optMethodDeriv);
576  }
577  // Begin topo search:
578  NNITopologySearch topoSearch(*tl, nniMethod, verbose > 2 ? verbose - 2 : 0);
579  NNITopologyListener* topoListener = new NNITopologyListener(&topoSearch, parameters, tolDuring, messageHandler, profiler, verbose, optMethodDeriv, nStep, reparametrization);
580  topoListener->setNumericalOptimizationCounter(numStep);
581  topoSearch.addTopologyListener(topoListener);
582  topoSearch.search();
583  return dynamic_cast<NNIHomogeneousTreeLikelihood*>(topoSearch.getSearchableObject());
584 }
585 
586 /******************************************************************************/
587 
590  const ParameterList& parameters,
591  bool optimizeNumFirst,
592  double tolBefore,
593  double tolDuring,
594  unsigned int tlEvalMax,
595  unsigned int numStep,
596  OutputStream* messageHandler,
597  OutputStream* profiler,
598  bool reparametrization,
599  unsigned int verbose,
600  const std::string& optMethodDeriv,
601  const std::string& nniMethod)
602 {
603  // Roughly optimize parameter
604  if (optimizeNumFirst)
605  {
606  OptimizationTools::optimizeNumericalParameters2(tl, parameters, NULL, tolBefore, 1000000, messageHandler, profiler, reparametrization, false, verbose, optMethodDeriv);
607  }
608  // Begin topo search:
609  NNITopologySearch topoSearch(*tl, nniMethod, verbose > 2 ? verbose - 2 : 0);
610  NNITopologyListener2* topoListener = new NNITopologyListener2(&topoSearch, parameters, tolDuring, messageHandler, profiler, verbose, optMethodDeriv, reparametrization);
611  topoListener->setNumericalOptimizationCounter(numStep);
612  topoSearch.addTopologyListener(topoListener);
613  topoSearch.search();
614  return dynamic_cast<NNIHomogeneousTreeLikelihood*>(topoSearch.getSearchableObject());
615 }
616 
617 /******************************************************************************/
618 
621  unsigned int verbose)
622 {
623  NNISearchable* topo = dynamic_cast<NNISearchable*>(tp);
624  NNITopologySearch topoSearch(*topo, NNITopologySearch::PHYML, verbose);
625  topoSearch.search();
626  return dynamic_cast<DRTreeParsimonyScore*>(topoSearch.getSearchableObject());
627 }
628 
629 /******************************************************************************/
630 
631 std::string OptimizationTools::DISTANCEMETHOD_INIT = "init";
632 std::string OptimizationTools::DISTANCEMETHOD_PAIRWISE = "pairwise";
633 std::string OptimizationTools::DISTANCEMETHOD_ITERATIONS = "iterations";
634 
635 /******************************************************************************/
636 
638  DistanceEstimation& estimationMethod,
639  const ParameterList& parametersToIgnore,
640  const std::string& param,
641  unsigned int verbose)
642 {
643  if (param != DISTANCEMETHOD_PAIRWISE && param != DISTANCEMETHOD_INIT)
644  throw Exception("OptimizationTools::estimateDistanceMatrix. Invalid option param=" + param + ".");
645  estimationMethod.resetAdditionalParameters();
646  estimationMethod.setVerbose(verbose);
647  if (param == DISTANCEMETHOD_PAIRWISE)
648  {
649  ParameterList tmp = estimationMethod.getModel().getIndependentParameters();
651  tmp.deleteParameters(parametersToIgnore.getParameterNames());
652  estimationMethod.setAdditionalParameters(tmp);
653  }
654  // Compute matrice:
655  if (verbose > 0)
656  ApplicationTools::displayTask("Estimating distance matrix", true);
657  estimationMethod.computeMatrix();
658  unique_ptr<DistanceMatrix> matrix(estimationMethod.getMatrix());
659  if (verbose > 0)
661 
662  return matrix.release();
663 }
664 
665 /******************************************************************************/
666 
668  DistanceEstimation& estimationMethod,
669  AgglomerativeDistanceMethod& reconstructionMethod,
670  const ParameterList& parametersToIgnore,
671  bool optimizeBrLen,
672  const std::string& param,
673  double tolerance,
674  unsigned int tlEvalMax,
675  OutputStream* profiler,
676  OutputStream* messenger,
677  unsigned int verbose)
678 {
679  estimationMethod.resetAdditionalParameters();
680  estimationMethod.setVerbose(verbose);
681  if (param == DISTANCEMETHOD_PAIRWISE)
682  {
683  ParameterList tmp = estimationMethod.getModel().getIndependentParameters();
685  tmp.deleteParameters(parametersToIgnore.getParameterNames());
686  estimationMethod.setAdditionalParameters(tmp);
687  }
688  TreeTemplate<Node>* tree = NULL;
689  TreeTemplate<Node>* previousTree = NULL;
690  bool test = true;
691  while (test)
692  {
693  // Compute matrice:
694  if (verbose > 0)
695  ApplicationTools::displayTask("Estimating distance matrix", true);
696  estimationMethod.computeMatrix();
697  DistanceMatrix* matrix = estimationMethod.getMatrix();
698  if (verbose > 0)
700 
701  // Compute tree:
702  if (matrix->size() == 2) {
703  //Special case, there is only one possible tree:
704  Node* n1 = new Node(0);
705  Node* n2 = new Node(1, matrix->getName(0));
706  n2->setDistanceToFather((*matrix)(0,0) / 2.);
707  Node* n3 = new Node(2, matrix->getName(1));
708  n3->setDistanceToFather((*matrix)(0,0) / 2.);
709  n1->addSon(n2);
710  n1->addSon(n3);
711  tree = new TreeTemplate<Node>(n1);
712  break;
713  }
714  if (verbose > 0)
715  ApplicationTools::displayTask("Building tree");
716  reconstructionMethod.setDistanceMatrix(*matrix);
717  reconstructionMethod.computeTree();
718  previousTree = tree;
719  delete matrix;
720  tree = dynamic_cast<TreeTemplate<Node>*>(reconstructionMethod.getTree());
721  if (verbose > 0)
723  if (previousTree && verbose > 0)
724  {
725  int rf = TreeTools::robinsonFouldsDistance(*previousTree, *tree, false);
726  ApplicationTools::displayResult("Topo. distance with previous iteration", TextTools::toString(rf));
727  test = (rf == 0);
728  delete previousTree;
729  }
730  if (param != DISTANCEMETHOD_ITERATIONS)
731  break; // Ends here.
732 
733  // Now, re-estimate parameters:
734  unique_ptr<TransitionModel> model(estimationMethod.getModel().clone());
735  unique_ptr<DiscreteDistribution> rdist(estimationMethod.getRateDistribution().clone());
737  *estimationMethod.getData(),
738  model.get(),
739  rdist.get(),
740  true, verbose > 1);
741  tl.initialize();
742  ParameterList parameters = tl.getParameters();
743  if (!optimizeBrLen)
744  {
745  vector<string> vs = tl.getBranchLengthsParameters().getParameterNames();
746  parameters.deleteParameters(vs);
747  }
748  parameters.deleteParameters(parametersToIgnore.getParameterNames());
749  optimizeNumericalParameters(&tl, parameters, NULL, 0, tolerance, tlEvalMax, messenger, profiler, verbose > 0 ? verbose - 1 : 0);
750  if (verbose > 0)
751  {
752  ParameterList tmp = tl.getSubstitutionModelParameters();
753  for (unsigned int i = 0; i < tmp.size(); i++)
754  {
755  ApplicationTools::displayResult(tmp[i].getName(), TextTools::toString(tmp[i].getValue()));
756  }
757  tmp = tl.getRateDistributionParameters();
758  for (unsigned int i = 0; i < tmp.size(); i++)
759  {
760  ApplicationTools::displayResult(tmp[i].getName(), TextTools::toString(tmp[i].getValue()));
761  }
762  }
763  }
764  return tree;
765 }
766 
767 /******************************************************************************/
768 
static unsigned int optimizeNumericalParameters2(DiscreteRatesAcrossSitesTreeLikelihood *tl, const ParameterList &parameters, OptimizationListener *listener=0, double tolerance=0.000001, unsigned int tlEvalMax=1000000, OutputStream *messageHandler=ApplicationTools::message.get(), OutputStream *profiler=ApplicationTools::message.get(), bool reparametrization=false, bool useClock=false, unsigned int verbose=1, const std::string &optMethodDeriv=OPTIMIZATION_NEWTON)
Optimize numerical parameters (branch length, substitution model & rate distribution) of a TreeLikeli...
static unsigned int optimizeTreeScale(TreeLikelihood *tl, double tolerance=0.000001, unsigned int tlEvalMax=1000000, OutputStream *messageHandler=ApplicationTools::message.get(), OutputStream *profiler=ApplicationTools::message.get(), unsigned int verbose=1)
Optimize the scale of a TreeLikelihood.
virtual void addSon(size_t pos, Node *node)
Definition: Node.h:407
void setNumericalOptimizationCounter(unsigned int c)
static void displayTaskDone()
virtual std::vector< std::string > getParameterNames() const
virtual void setMessageHandler(OutputStream *mh)=0
static unsigned int optimizeNumericalParametersWithGlobalClock(DiscreteRatesAcrossSitesClockTreeLikelihood *cl, const ParameterList &parameters, OptimizationListener *listener=0, unsigned int nstep=1, double tolerance=0.000001, unsigned int tlEvalMax=1000000, OutputStream *messageHandler=ApplicationTools::message.get(), OutputStream *profiler=ApplicationTools::message.get(), unsigned int verbose=1, const std::string &optMethodDeriv=OPTIMIZATION_GRADIENT)
Optimize numerical parameters assuming a global clock (branch heights, substitution model & rate dist...
static DistanceMatrix * estimateDistanceMatrix(DistanceEstimation &estimationMethod, const ParameterList &parametersToIgnore, const std::string &param=DISTANCEMETHOD_INIT, unsigned int verbose=0)
Estimate a distance matrix using maximum likelihood.
static std::string DISTANCEMETHOD_PAIRWISE
virtual void setDistanceMatrix(const DistanceMatrix &matrix)=0
Set the distance matrix to use.
Listener used internally by the optimizeTreeNNI2 method.
void search()
Performs the search.
TransitionModel * clone() const =0
virtual double f(const ParameterList &parameters)
virtual ParameterList getBranchLengthsParameters() const =0
Get the branch lengths parameters.
OptimizationStopCondition * getStopCondition()
static unsigned int optimizeNumericalParameters(DiscreteRatesAcrossSitesTreeLikelihood *tl, const ParameterList &parameters, OptimizationListener *listener=0, unsigned int nstep=1, double tolerance=0.000001, unsigned int tlEvalMax=1000000, OutputStream *messageHandler=ApplicationTools::message.get(), OutputStream *profiler=ApplicationTools::message.get(), bool reparametrization=false, unsigned int verbose=1, const std::string &optMethodDeriv=OPTIMIZATION_NEWTON, const std::string &optMethodModel=OPTIMIZATION_BRENT)
Optimize numerical parameters (branch length, substitution model & rate distribution) of a TreeLikeli...
DistanceMatrix * getMatrix() const
Get the distance matrix.
virtual bool matchParametersValues(const ParameterList &params, std::vector< size_t > *updatedParameters=0)
virtual void enableSecondOrderDerivatives(bool yn)=0
virtual void deleteParameters(const std::vector< std::string > &names, bool mustExist=true)
virtual double optimize()=0
void setParameters(const ParameterList &lambda)
The TreeLikelihood interface.
virtual void addOptimizer(const std::string &name, Optimizer *optimizer, const std::vector< std::string > &params, unsigned short derivatives=0, const std::string &type=IT_TYPE_STEP)
static std::string IT_TYPE_STEP
A listener which capture NaN function values and throw an exception in case this happens.
Double recursive implementation of interface TreeParsimonyScore.
static void displayResult(const std::string &text, const T &result)
virtual const ParameterList & getParameters() const =0
static void displayTask(const std::string &text, bool eof=false)
static std::string DISTANCEMETHOD_INIT
STL namespace.
Interface for Nearest Neighbor Interchanges algorithms.
Definition: NNISearchable.h:96
The phylogenetic tree class.
static TreeTemplate< Node > * buildDistanceTree(DistanceEstimation &estimationMethod, AgglomerativeDistanceMethod &reconstructionMethod, const ParameterList &parametersToIgnore, bool optimizeBrLen=false, const std::string &param=DISTANCEMETHOD_INIT, double tolerance=0.000001, unsigned int tlEvalMax=1000000, OutputStream *profiler=0, OutputStream *messenger=0, unsigned int verbose=0)
Build a tree using a distance method.
static std::string OPTIMIZATION_GRADIENT
static unsigned int optimizeNumericalParametersWithGlobalClock2(DiscreteRatesAcrossSitesClockTreeLikelihood *cl, const ParameterList &parameters, OptimizationListener *listener=0, double tolerance=0.000001, unsigned int tlEvalMax=1000000, OutputStream *messageHandler=ApplicationTools::message.get(), OutputStream *profiler=ApplicationTools::message.get(), unsigned int verbose=1, const std::string &optMethodDeriv=OPTIMIZATION_GRADIENT)
Optimize numerical parameters assuming a global clock (branch heights, substitution model & rate dist...
unsigned int getNumberOfEvaluations() const
virtual void deleteParameter(const std::string &name)
This class adds support for NNI topology estimation to the DRHomogeneousTreeLikelihood class...
std::size_t size() const
static unsigned int optimizeBranchLengthsParameters(DiscreteRatesAcrossSitesTreeLikelihood *tl, const ParameterList &parameters, OptimizationListener *listener=0, double tolerance=0.000001, unsigned int tlEvalMax=1000000, OutputStream *messageHandler=ApplicationTools::message.get(), OutputStream *profiler=ApplicationTools::message.get(), unsigned int verbose=1, const std::string &optMethodDeriv=OPTIMIZATION_NEWTON)
Optimize branch lengths parameters of a TreeLikelihood function.
static std::string DISTANCEMETHOD_ITERATIONS
void setVerbose(size_t verbose)
virtual void setTolerance(double tolerance)=0
static std::string IT_TYPE_FULL
const DiscreteDistribution & getRateDistribution() const
void setMaximumNumberOfEvaluations(unsigned int max)
virtual void addParameters(const ParameterList &params)
virtual void enableFirstOrderDerivatives(bool yn)=0
Interface for agglomerative distance methods.
const TransitionModel & getModel() const
const std::string & getName(std::size_t i) const
virtual ParameterList getRateDistributionParameters() const =0
Get the parameters associated to the rate distirbution.
virtual void addOptimizationListener(OptimizationListener *listener)=0
const SiteContainer * getData() const
void addOptimizationListener(OptimizationListener *listener)
virtual bool hasParameter(const std::string &name) const
void addTopologyListener(TopologyListener *listener)
Add a listener to the list.
virtual ParameterList getNonDerivableParameters() const =0
All non derivable parameters.
Interface for rate across sites (RAS) implementation.
This class implements the likelihood computation for a tree using the double-recursive algorithm...
static std::string OPTIMIZATION_BRENT
void setNumericalOptimizationCounter(unsigned int c)
Interface for likelihood computation with a global clock and rate across sites variation.
void computeMatrix()
Perform the distance computation.
This Optimizer implements Newton&#39;s algorithm for finding a minimum of a function. This is in fact a m...
virtual const ParameterList & getIndependentParameters() const =0
static void displayMessage(const std::string &text)
virtual unsigned int getNumberOfEvaluations() const =0
Class for notifying new toplogy change events.
NNISearchable * getSearchableObject()
size_t size() const
void topologyChangeSuccessful(const TopologyChangeEvent &event)
Tell that a topology change is definitive.
The phylogenetic node class.
Definition: Node.h:90
virtual void setVerbose(unsigned int v)=0
virtual void setMaximumNumberOfEvaluations(unsigned int max)=0
virtual ParameterList getCommonParametersWith(const ParameterList &params) const
static const std::string PHYML
virtual void addParameter(const Parameter &param)
void setParametersToDerivate(const std::vector< std::string > &variables)
void setProfiler(OutputStream *profiler)
void setVerbose(unsigned int v)
static std::string CONSTRAINTS_AUTO
virtual void setDistanceToFather(double distance)
Set or update the distance toward the father node.
Definition: Node.h:299
virtual ParameterList subList(const std::vector< std::string > &names) const
void setAdditionalParameters(const ParameterList &parameters)
Specify a list of parameters to be estimated.
void setStopCondition(const OptimizationStopCondition &stopCondition)
virtual double getLogLikelihood() const =0
Get the logarithm of the likelihood for the whole dataset.
static NNIHomogeneousTreeLikelihood * optimizeTreeNNI(NNIHomogeneousTreeLikelihood *tl, const ParameterList &parameters, bool optimizeNumFirst=true, double tolBefore=100, double tolDuring=100, unsigned int tlEvalMax=1000000, unsigned int numStep=1, OutputStream *messageHandler=ApplicationTools::message.get(), OutputStream *profiler=ApplicationTools::message.get(), bool reparametrization=false, unsigned int verbose=1, const std::string &optMethod=OptimizationTools::OPTIMIZATION_NEWTON, unsigned int nStep=1, const std::string &nniMethod=NNITopologySearch::PHYML)
Optimize all parameters from a TreeLikelihood object, including tree topology using Nearest Neighbor ...
void setConstraintPolicy(const std::string &constraintPolicy)
void setMessageHandler(OutputStream *mh)
virtual void setProfiler(OutputStream *profiler)=0
void topologyChangeSuccessful(const TopologyChangeEvent &event)
Tell that a topology change is definitive.
void setInitialInterval(double inf, double sup)
Estimate a distance matrix from sequence data, according to a given model.
std::string toString(T t)
void resetAdditionalParameters()
Reset all additional parameters.
virtual void init(const ParameterList &params)=0
static std::string OPTIMIZATION_BFGS
static NNIHomogeneousTreeLikelihood * optimizeTreeNNI2(NNIHomogeneousTreeLikelihood *tl, const ParameterList &parameters, bool optimizeNumFirst=true, double tolBefore=100, double tolDuring=100, unsigned int tlEvalMax=1000000, unsigned int numStep=1, OutputStream *messageHandler=ApplicationTools::message.get(), OutputStream *profiler=ApplicationTools::message.get(), bool reparametrization=false, unsigned int verbose=1, const std::string &optMethod=OptimizationTools::OPTIMIZATION_NEWTON, const std::string &nniMethod=NNITopologySearch::PHYML)
Optimize all parameters from a TreeLikelihood object, including tree topology using Nearest Neighbor ...
virtual Tree * getTree() const =0
virtual void computeTree()=0
Perform the clustering.
DiscreteDistribution * clone() const =0
Listener used internally by the optimizeTreeNNI method.
const ParameterList & getParameters() const
virtual ParameterList getSubstitutionModelParameters() const =0
Get the parameters associated to substitution model(s).
static std::string OPTIMIZATION_NEWTON
virtual void setParametersValues(const ParameterList &params)
NNI topology search method.
virtual void setConstraintPolicy(const std::string &constraintPolicy)=0
static int robinsonFouldsDistance(const Tree &tr1, const Tree &tr2, bool checkNames=true, int *missing_in_tr2=NULL, int *missing_in_tr1=NULL)
Calculates the Robinson-Foulds topological distance between two trees.
Definition: TreeTools.cpp:887
virtual OptimizationStopCondition * getStopCondition()=0
void init(const ParameterList &params)