~ubuntu-branches/ubuntu/precise/libbpp-phyl/precise

« back to all changes in this revision

Viewing changes to src/Bpp/Phyl/OptimizationTools.h

  • Committer: Bazaar Package Importer
  • Author(s): Julien Dutheil
  • Date: 2011-06-09 11:00:00 UTC
  • Revision ID: james.westby@ubuntu.com-20110609110000-yvx78svv6w7xxgph
Tags: upstream-2.0.2
ImportĀ upstreamĀ versionĀ 2.0.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
//
 
2
// File: OptimizationTools.h
 
3
// Created by: Julien Dutheil
 
4
// Created on: Sun Dec 14 09:43:32 2003
 
5
//
 
6
 
 
7
/*
 
8
  Copyright or Ā© or Copr. CNRS, (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
#ifndef _OPTIMIZATIONTOOLS_H_
 
41
#define _OPTIMIZATIONTOOLS_H_
 
42
 
 
43
#include "Likelihood/ClockTreeLikelihood.h"
 
44
#include "Likelihood/NNIHomogeneousTreeLikelihood.h"
 
45
#include "Likelihood/ClockTreeLikelihood.h"
 
46
#include "NNITopologySearch.h"
 
47
#include "Parsimony/DRTreeParsimonyScore.h"
 
48
#include "TreeTemplate.h"
 
49
#include "Distance/DistanceEstimation.h"
 
50
#include "Distance/AgglomerativeDistanceMethod.h"
 
51
 
 
52
#include <Bpp/Io/OutputStream.h>
 
53
#include <Bpp/App/ApplicationTools.h>
 
54
#include <Bpp/Numeric/Function/SimpleNewtonMultiDimensions.h>
 
55
 
 
56
namespace bpp
 
57
{
 
58
 
 
59
  
 
60
  /**
 
61
   * @brief Debugging wrapper: output data to a log file if a 0 likelihood is obtained during optimization.
 
62
   */
 
63
  class NaNWatcher: public DerivableSecondOrderWrapper
 
64
  {
 
65
 
 
66
  public:
 
67
    NaNWatcher(TreeLikelihood* tl) :
 
68
      DerivableSecondOrderWrapper(tl) {}
 
69
 
 
70
    NaNWatcher* clone() const { return new NaNWatcher(*this); }
 
71
 
 
72
    /**
 
73
     * @brief This function is redefined to output an error message if there is a 0 likelihood.
 
74
     *
 
75
     * @return The value of the likelihood function;
 
76
     */
 
77
    double getValue() const throw (Exception);
 
78
 
 
79
  };
 
80
 
 
81
 
 
82
 
 
83
  /**
 
84
   * @brief Listener used internally by the optimizeTreeNNI method.
 
85
   */
 
86
  class NNITopologyListener:
 
87
    public virtual TopologyListener
 
88
  {
 
89
  private:
 
90
    NNITopologySearch* topoSearch_;
 
91
    ParameterList parameters_;
 
92
    double tolerance_;
 
93
    OutputStream* messenger_;
 
94
    OutputStream* profiler_;
 
95
    unsigned int verbose_;
 
96
    unsigned int optimizeCounter_;
 
97
    unsigned int optimizeNumerical_;
 
98
    std::string optMethod_;
 
99
    unsigned int nStep_;
 
100
    bool reparametrization_;
 
101
 
 
102
  public:
 
103
    /**
 
104
     * @brief Build a new NNITopologyListener object.
 
105
     *
 
106
     * This listener listens to a NNITopologySearch object, and optimizes numerical parameters every *n* topological movements.
 
107
     * Optimization is performed using the optimizeNumericalParameters method (see there documentation for more details).
 
108
     *
 
109
     * @param ts         The NNITopologySearch object attached to this listener.
 
110
     * @param parameters The list of parameters to optimize. Use tl->getIndependentParameters() in order to estimate all parameters.
 
111
     * @param tolerance  Tolerance to use during optimizaton.
 
112
     * @param messenger  Where to output messages.
 
113
     * @param profiler   Where to output optimization steps.
 
114
     * @param verbose    Verbose level during optimization.
 
115
     * @param optMethod  Optimization method to use.
 
116
     * @param nStep      The number of optimization steps to perform.
 
117
     * @param reparametrization Tell if parameters should be transformed in order to remove constraints.
 
118
     *                          This can improve optimization, but is a bit slower.
 
119
     */
 
120
    NNITopologyListener(
 
121
                        NNITopologySearch* ts,
 
122
                        const ParameterList& parameters,
 
123
                        double tolerance,
 
124
                        OutputStream* messenger,
 
125
                        OutputStream* profiler,
 
126
                        unsigned int verbose,
 
127
                        const std::string& optMethod,
 
128
                        unsigned int nStep,
 
129
                        bool reparametrization) :
 
130
      topoSearch_(ts), parameters_(parameters), tolerance_(tolerance),
 
131
      messenger_(messenger), profiler_(profiler),
 
132
      verbose_(verbose),
 
133
      optimizeCounter_(0), optimizeNumerical_(1),
 
134
      optMethod_(optMethod), nStep_(nStep), reparametrization_(reparametrization) {}
 
135
 
 
136
    NNITopologyListener(const NNITopologyListener& tl) :
 
137
      topoSearch_(tl.topoSearch_),
 
138
      parameters_(tl.parameters_),
 
139
      tolerance_(tl.tolerance_),
 
140
      messenger_(tl.messenger_),
 
141
      profiler_(tl.profiler_),
 
142
      verbose_(tl.verbose_),
 
143
      optimizeCounter_(tl.optimizeCounter_),
 
144
      optimizeNumerical_(tl.optimizeNumerical_),
 
145
      optMethod_(tl.optMethod_),
 
146
      nStep_(tl.nStep_),
 
147
      reparametrization_(tl.reparametrization_)
 
148
    {}
 
149
 
 
150
    NNITopologyListener& operator=(const NNITopologyListener& tl)
 
151
    {
 
152
      topoSearch_        = tl.topoSearch_;
 
153
      parameters_        = tl.parameters_;
 
154
      tolerance_         = tl.tolerance_;
 
155
      messenger_         = tl.messenger_;
 
156
      profiler_          = tl.profiler_;
 
157
      verbose_           = tl.verbose_;
 
158
      optimizeCounter_   = tl.optimizeCounter_;
 
159
      optimizeNumerical_ = tl.optimizeNumerical_;
 
160
      optMethod_         = tl.optMethod_;
 
161
      nStep_             = tl.nStep_;
 
162
      reparametrization_ = tl.reparametrization_;
 
163
      return *this;
 
164
    }
 
165
 
 
166
    NNITopologyListener* clone() const { return new NNITopologyListener(*this); }
 
167
 
 
168
    virtual ~NNITopologyListener() {}
 
169
 
 
170
  public:
 
171
    void topologyChangeTested(const TopologyChangeEvent& event) {}
 
172
    void topologyChangeSuccessful(const TopologyChangeEvent& event);
 
173
    void setNumericalOptimizationCounter(unsigned int c) { optimizeNumerical_ = c; }
 
174
 
 
175
  };
 
176
 
 
177
  /**
 
178
   * @brief Listener used internally by the optimizeTreeNNI2 method.
 
179
   */
 
180
  class NNITopologyListener2:
 
181
    public TopologyListener
 
182
  {
 
183
  private:
 
184
    NNITopologySearch* topoSearch_;
 
185
    ParameterList parameters_;
 
186
    double tolerance_;
 
187
    OutputStream* messenger_;
 
188
    OutputStream* profiler_;
 
189
    unsigned int verbose_;
 
190
    unsigned int optimizeCounter_;
 
191
    unsigned int optimizeNumerical_;
 
192
    std::string optMethod_;
 
193
    bool reparametrization_;
 
194
 
 
195
  public:
 
196
    /**
 
197
     * @brief Build a new NNITopologyListener2 object.
 
198
     *
 
199
     * This listener listens to a NNITopologySearch object, and optimizes numerical parameters every *n* topological movements.
 
200
     * Optimization is performed using the optimizeNumericalParameters2 method (see there documentation for more details).
 
201
     *
 
202
     * @param ts         The NNITopologySearch object attached to this listener.
 
203
     * @param parameters The list of parameters to optimize. Use ts->getIndependentParameters() in order to estimate all parameters.
 
204
     * @param tolerance  Tolerance to use during optimizaton.
 
205
     * @param messenger  Where to output messages.
 
206
     * @param profiler   Where to output optimization steps.
 
207
     * @param verbose    Verbose level during optimization.
 
208
     * @param optMethod  Optimization method to use.
 
209
     * @param reparametrization Tell if parameters should be transformed in order to remove constraints.
 
210
     *                          This can improve optimization, but is a bit slower.
 
211
     */
 
212
    NNITopologyListener2(
 
213
                         NNITopologySearch* ts,
 
214
                         const ParameterList& parameters,
 
215
                         double tolerance,
 
216
                         OutputStream* messenger,
 
217
                         OutputStream* profiler,
 
218
                         unsigned int verbose,
 
219
                         const std::string& optMethod,
 
220
                         bool reparametrization) :
 
221
      topoSearch_(ts), parameters_(parameters), tolerance_(tolerance),
 
222
      messenger_(messenger), profiler_(profiler),
 
223
      verbose_(verbose),
 
224
      optimizeCounter_(0), optimizeNumerical_(1),
 
225
      optMethod_(optMethod), reparametrization_(reparametrization) {}
 
226
 
 
227
    NNITopologyListener2(const NNITopologyListener2& tl) :
 
228
      topoSearch_(tl.topoSearch_),
 
229
      parameters_(tl.parameters_),
 
230
      tolerance_(tl.tolerance_),
 
231
      messenger_(tl.messenger_),
 
232
      profiler_(tl.profiler_),
 
233
      verbose_(tl.verbose_),
 
234
      optimizeCounter_(tl.optimizeCounter_),
 
235
      optimizeNumerical_(tl.optimizeNumerical_),
 
236
      optMethod_(tl.optMethod_),
 
237
      reparametrization_(tl.reparametrization_)
 
238
    {}
 
239
 
 
240
    NNITopologyListener2& operator=(const NNITopologyListener2& tl)
 
241
    {
 
242
      topoSearch_        = tl.topoSearch_;
 
243
      parameters_        = tl.parameters_;
 
244
      tolerance_         = tl.tolerance_;
 
245
      messenger_         = tl.messenger_;
 
246
      profiler_          = tl.profiler_;
 
247
      verbose_           = tl.verbose_;
 
248
      optimizeCounter_   = tl.optimizeCounter_;
 
249
      optimizeNumerical_ = tl.optimizeNumerical_;
 
250
      optMethod_         = tl.optMethod_;
 
251
      reparametrization_ = tl.reparametrization_;
 
252
      return *this;
 
253
    }
 
254
 
 
255
    NNITopologyListener2* clone() const { return new NNITopologyListener2(*this); }
 
256
 
 
257
    virtual ~NNITopologyListener2() {}
 
258
 
 
259
  public:
 
260
    void topologyChangeTested(const TopologyChangeEvent& event) {}
 
261
    void topologyChangeSuccessful(const TopologyChangeEvent& event);
 
262
    void setNumericalOptimizationCounter(unsigned int c) { optimizeNumerical_ = c; }
 
263
 
 
264
  };
 
265
 
 
266
 
 
267
 
 
268
 
 
269
 
 
270
 
 
271
 
 
272
  /**
 
273
   * @brief Optimization methods for phylogenetic inference.
 
274
   *
 
275
   * This class provides optimization methods for phylogenetics.
 
276
   * Parameters of the optimization methods are set to work with TreeLikelihood
 
277
   * object. Some non trivial parameters are left to the user choice (tolerance, maximum
 
278
   * number of function evaluation, output streams).
 
279
   */
 
280
  class OptimizationTools
 
281
  {
 
282
  public:
 
283
    OptimizationTools();
 
284
    virtual ~OptimizationTools();
 
285
        
 
286
  public:
 
287
 
 
288
    static std::string OPTIMIZATION_GRADIENT;
 
289
    static std::string OPTIMIZATION_NEWTON;
 
290
    static std::string OPTIMIZATION_BRENT;
 
291
    static std::string OPTIMIZATION_BFGS;
 
292
                
 
293
    /**
 
294
     * @brief Optimize numerical parameters (branch length, substitution model & rate distribution) of a TreeLikelihood function.
 
295
     *
 
296
     * Uses Newton's method for branch length and Brent or BFGS one dimensional method for other parameters.
 
297
     *
 
298
     * A condition over function values is used as a stop condition for the algorithm.
 
299
     *
 
300
     * @see BrentOneDimension, BFGSMultiDimensions
 
301
     *
 
302
     * @param tl             A pointer toward the TreeLikelihood object to optimize.
 
303
     * @param parameters     The list of parameters to optimize. Use tl->getIndependentParameters() in order to estimate all parameters.
 
304
     * @param listener       A pointer toward an optimization listener, if needed.
 
305
     * @param nstep          The number of progressive steps to perform (see NewtonBrentMetaOptimizer). 1 means full precision from start.
 
306
     * @param tolerance      The tolerance to use in the algorithm.
 
307
     * @param tlEvalMax      The maximum number of function evaluations.
 
308
     * @param messageHandler The massage handler.
 
309
     * @param profiler       The profiler.
 
310
     * @param reparametrization Tell if parameters should be transformed in order to remove constraints.
 
311
     *                          This can improve optimization, but is a bit slower.
 
312
     * @param verbose        The verbose level.
 
313
     * @param optMethodDeriv Optimization type for derivable parameters (first or second order derivatives).
 
314
     * @see OPTIMIZATION_NEWTON, OPTIMIZATION_GRADIENT
 
315
     * @param optMethodModel Optimization type for model parameters (Brent or BFGS).
 
316
     * @see OPTIMIZATION_BRENT, OPTIMIZATION_BFGS     
 
317
     * @throw Exception any exception thrown by the Optimizer.
 
318
     */
 
319
    static unsigned int optimizeNumericalParameters(
 
320
                                                    DiscreteRatesAcrossSitesTreeLikelihood* tl,
 
321
                                                    const ParameterList& parameters,
 
322
                                                    OptimizationListener* listener = 0,
 
323
                                                    unsigned int nstep             = 1,
 
324
                                                    double tolerance               = 0.000001,
 
325
                                                    unsigned int tlEvalMax         = 1000000,
 
326
                                                    OutputStream* messageHandler   = ApplicationTools::message,
 
327
                                                    OutputStream* profiler         = ApplicationTools::message,
 
328
                                                    bool reparametrization         = false,
 
329
                                                    unsigned int verbose           = 1,
 
330
                                                    const std::string& optMethodDeriv      = OPTIMIZATION_NEWTON,
 
331
                                                    const std::string& optMethodModel      = OPTIMIZATION_BFGS)
 
332
      throw (Exception);
 
333
        
 
334
    /**
 
335
     * @brief Optimize numerical parameters (branch length, substitution model & rate distribution) of a TreeLikelihood function.
 
336
     *
 
337
     * Uses Newton's method for all parameters, branch length derivatives are computed analytically, derivatives for other parameters numerically.
 
338
     *
 
339
     * @see PseudoNewtonOptimizer
 
340
     *
 
341
     * @param tl             A pointer toward the TreeLikelihood object to optimize.
 
342
     * @param parameters     The list of parameters to optimize. Use tl->getIndependentParameters() in order to estimate all parameters.
 
343
     * @param listener       A pointer toward an optimization listener, if needed.
 
344
     * @param tolerance      The tolerance to use in the algorithm.
 
345
     * @param tlEvalMax      The maximum number of function evaluations.
 
346
     * @param messageHandler The massage handler.
 
347
     * @param profiler       The profiler.
 
348
     * @param reparametrization Tell if parameters should be transformed in order to remove constraints.
 
349
     *                          This can improve optimization, but is a bit slower.
 
350
     * @param verbose        The verbose level.
 
351
     * @param optMethodDeriv Optimization type for derivable parameters (first or second order derivatives).
 
352
     * @see OPTIMIZATION_NEWTON, OPTIMIZATION_GRADIENT
 
353
     * @throw Exception any exception thrown by the Optimizer.
 
354
     */
 
355
    static unsigned int optimizeNumericalParameters2(
 
356
                                                     DiscreteRatesAcrossSitesTreeLikelihood* tl,
 
357
                                                     const ParameterList& parameters,
 
358
                                                     OptimizationListener* listener = 0,
 
359
                                                     double tolerance               = 0.000001,
 
360
                                                     unsigned int tlEvalMax         = 1000000,
 
361
                                                     OutputStream* messageHandler   = ApplicationTools::message,
 
362
                                                     OutputStream* profiler         = ApplicationTools::message,
 
363
                                                     bool reparametrization         = false,
 
364
                                                     unsigned int verbose           = 1,
 
365
                                                     const std::string& optMethodDeriv  = OPTIMIZATION_NEWTON)
 
366
      throw (Exception);
 
367
 
 
368
    /**
 
369
     * @brief Optimize branch lengths parameters of a TreeLikelihood function.
 
370
     *
 
371
     * Uses Newton's method.
 
372
     *
 
373
     * A condition over function values is used as a stop condition for the algorithm.
 
374
     *
 
375
     * @see NewtonBrentMetaOptimizer
 
376
     *
 
377
     * @param tl             A pointer toward the TreeLikelihood object to optimize.
 
378
     * @param parameters     The list of parameters to optimize. The intersection of branch length parameters and the input set will be used. Use tl->getBranchLengthsParameters() in order to estimate all branch length parameters.
 
379
     * @param listener       A pointer toward an optimization listener, if needed.
 
380
     * @param tolerance      The tolerance to use in the algorithm.
 
381
     * @param tlEvalMax      The maximum number of function evaluations.
 
382
     * @param messageHandler The massage handler.
 
383
     * @param profiler       The profiler.
 
384
     * @param verbose        The verbose level.
 
385
     * @param optMethodDeriv Optimization type for derivable parameters (first or second order derivatives).
 
386
     * @see OPTIMIZATION_NEWTON, OPTIMIZATION_GRADIENT
 
387
     * @throw Exception any exception thrown by the Optimizer.
 
388
     */
 
389
    static unsigned int optimizeBranchLengthsParameters(
 
390
                                                        DiscreteRatesAcrossSitesTreeLikelihood * tl,
 
391
                                                        const ParameterList& parameters,
 
392
                                                        OptimizationListener* listener = 0,
 
393
                                                        double tolerance               = 0.000001,
 
394
                                                        unsigned int tlEvalMax         = 1000000,
 
395
                                                        OutputStream* messageHandler   = ApplicationTools::message,
 
396
                                                        OutputStream* profiler         = ApplicationTools::message,
 
397
                                                        unsigned int verbose           = 1,
 
398
                                                        const std::string& optMethodDeriv  = OPTIMIZATION_NEWTON)
 
399
      throw (Exception);
 
400
                
 
401
    /**
 
402
     * @brief Optimize numerical parameters assuming a global clock (branch heights, substitution model & rate distribution) of a ClockTreeLikelihood function.
 
403
     *
 
404
     * Uses Newton or conjugate gradient method for branch length and Brent's one dimensional method for other parameters
 
405
     * (NewtonBrentMetaOptimizer).
 
406
     * Derivatives are computed analytically.
 
407
     *
 
408
     * A condition over function values is used as a stop condition for the algorithm.
 
409
     *
 
410
     * @see NewtonBrentMetaOptimizer
 
411
     *
 
412
     * @param cl             A pointer toward the ClockTreeLikelihood object to optimize.
 
413
     * @param parameters     The list of parameters to optimize. Use cl->getIndependentParameters() in order to estimate all parameters.
 
414
     * @param listener       A pointer toward an optimization listener, if needed.
 
415
     * @param nstep          The number of progressive steps to perform (see NewtonBrentMetaOptimizer). 1 means full precision from start.
 
416
     * @param tolerance      The tolerance to use in the algorithm.
 
417
     * @param tlEvalMax      The maximum number of function evaluations.
 
418
     * @param messageHandler The massage handler.
 
419
     * @param profiler       The profiler.
 
420
     * @param verbose        The verbose level.
 
421
     * @param optMethodDeriv Optimization type for derivable parameters (first or second order derivatives).
 
422
     * @see OPTIMIZATION_NEWTON, OPTIMIZATION_GRADIENT
 
423
     * @throw Exception any exception thrown by the Optimizer.
 
424
     */
 
425
    static unsigned int optimizeNumericalParametersWithGlobalClock(
 
426
                                                                   DiscreteRatesAcrossSitesClockTreeLikelihood* cl,
 
427
                                                                   const ParameterList& parameters,
 
428
                                                                   OptimizationListener* listener = 0,
 
429
                                                                   unsigned int nstep = 1,
 
430
                                                                   double tolerance = 0.000001,
 
431
                                                                   unsigned int tlEvalMax = 1000000,
 
432
                                                                   OutputStream* messageHandler = ApplicationTools::message,
 
433
                                                                   OutputStream* profiler       = ApplicationTools::message,
 
434
                                                                   unsigned int verbose = 1,
 
435
                                                                   const std::string& optMethodDeriv = OPTIMIZATION_GRADIENT)
 
436
      throw (Exception);
 
437
 
 
438
    /**
 
439
     * @brief Optimize numerical parameters assuming a global clock (branch heights, substitution model & rate distribution) of a ClockTreeLikelihood function.
 
440
     *
 
441
     * Uses Newton or conjugate gradient method for all parameters, branch length derivatives are computed analytically, derivatives for other parameters numerically.
 
442
     *
 
443
     * @see PseudoNewtonOptimizer
 
444
     *
 
445
     * @param cl             A pointer toward the ClockTreeLikelihood object to optimize.
 
446
     * @param parameters     The list of parameters to optimize. Use cl->getIndependentParameters() in order to estimate all parameters.
 
447
     * @param listener       A pointer toward an optimization listener, if needed.
 
448
     * @param tolerance      The tolerance to use in the algorithm.
 
449
     * @param tlEvalMax      The maximum number of function evaluations.
 
450
     * @param messageHandler The massage handler.
 
451
     * @param profiler       The profiler.
 
452
     * @param verbose        The verbose level.
 
453
     * @param optMethodDeriv Optimization type for derivable parameters (first or second order derivatives).
 
454
     * @see OPTIMIZATION_NEWTON, OPTIMIZATION_GRADIENT
 
455
     * @throw Exception any exception thrown by the Optimizer.
 
456
     */
 
457
    static unsigned int optimizeNumericalParametersWithGlobalClock2(
 
458
                                                                    DiscreteRatesAcrossSitesClockTreeLikelihood* cl,
 
459
                                                                    const ParameterList& parameters,
 
460
                                                                    OptimizationListener* listener = 0,
 
461
                                                                    double tolerance = 0.000001,
 
462
                                                                    unsigned int tlEvalMax = 1000000,
 
463
                                                                    OutputStream* messageHandler = ApplicationTools::message,
 
464
                                                                    OutputStream* profiler       = ApplicationTools::message,
 
465
                                                                    unsigned int verbose = 1,
 
466
                                                                    const std::string& optMethodDeriv = OPTIMIZATION_GRADIENT)
 
467
      throw (Exception);
 
468
 
 
469
 
 
470
  private:
 
471
                
 
472
    class ScaleFunction:
 
473
      public virtual Function,
 
474
      public ParametrizableAdapter
 
475
    {                           
 
476
    private:
 
477
      TreeLikelihood* tl_;
 
478
      mutable ParameterList brLen_, lambda_;
 
479
                                
 
480
    public:
 
481
      ScaleFunction(TreeLikelihood* tl);
 
482
                                
 
483
      ScaleFunction(const ScaleFunction& sf) :
 
484
        tl_(sf.tl_), brLen_(sf.brLen_), lambda_(sf.lambda_)
 
485
      {}
 
486
 
 
487
      ScaleFunction& operator=(const ScaleFunction& sf)
 
488
      {
 
489
        tl_     = sf.tl_;
 
490
        brLen_  = sf.brLen_;
 
491
        lambda_ = sf.lambda_;
 
492
        return *this;
 
493
      }
 
494
 
 
495
      virtual ~ScaleFunction();
 
496
 
 
497
#ifndef NO_VIRTUAL_COV
 
498
      ScaleFunction*
 
499
#else
 
500
      Clonable*
 
501
#endif
 
502
      clone() const { return new ScaleFunction(*this); }
 
503
                                
 
504
    public:
 
505
      void setParameters(const ParameterList& lambda) throw (ParameterNotFoundException, ConstraintException);
 
506
      double getValue() const throw (ParameterException);
 
507
      const ParameterList& getParameters() const throw (Exception) { return lambda_; }
 
508
      const Parameter& getParameter(const std::string& name) const throw (ParameterNotFoundException)
 
509
      {
 
510
        if(name == "lambda") return lambda_[0];
 
511
        else throw ParameterNotFoundException("ScaleFunction::getParameter.", name);
 
512
      }
 
513
      double getParameterValue(const std::string& name) const throw (ParameterNotFoundException)
 
514
      {
 
515
        return lambda_.getParameter(name).getValue();
 
516
      }
 
517
      unsigned int getNumberOfParameters() const { return 1; }
 
518
      unsigned int getNumberOfIndependentParameters() const { return 1; }
 
519
    };
 
520
        
 
521
  public:
 
522
 
 
523
    /**
 
524
     * @brief Optimize the scale of a TreeLikelihood.
 
525
     *
 
526
     * This method only works on branch lengths parameters.
 
527
     * It multiply all branch length by a factor 'x' which is optimized
 
528
     * using Brent's algorithm in one dimension.
 
529
     * This method may be usefull for scaling a tree whose branch lengths
 
530
     * come from the Neighbor-Joining algorithm for instance.
 
531
     *
 
532
     * Practically, and contrarily to what one may expect, this does not
 
533
     * speed up the optimization!
 
534
     *
 
535
     * A condition over parameters is used as a stop condition for the algorithm.
 
536
     *
 
537
     * @param tl             A pointer toward the TreeLikelihood object to optimize.
 
538
     * @param tolerance      The tolerance to use in the algorithm.
 
539
     * @param tlEvalMax      The maximum number of function evaluations.
 
540
     * @param messageHandler The massage handler.
 
541
     * @param profiler       The profiler.
 
542
     * @param verbose        The verbose level.
 
543
     * @throw Exception any exception thrown by the optimizer.
 
544
     */
 
545
    static unsigned int optimizeTreeScale(
 
546
                                          TreeLikelihood* tl,
 
547
                                          double tolerance = 0.000001,
 
548
                                          int tlEvalMax = 1000000,
 
549
                                          OutputStream* messageHandler = ApplicationTools::message,
 
550
                                          OutputStream* profiler       = ApplicationTools::message,
 
551
                                          unsigned int verbose = 1)
 
552
      throw (Exception);
 
553
 
 
554
    /**
 
555
     * @brief Optimize all parameters from a TreeLikelihood object, including tree topology using Nearest Neighbor Interchanges.
 
556
     *
 
557
     * This function takes as input a TreeLikelihood object implementing the NNISearchable interface.
 
558
     *
 
559
     * Details:
 
560
     * A NNITopologySearch object is instanciated and is associated an additional TopologyListener.
 
561
     * This listener is used to re-estimate numerical parameters after one or several topology change.
 
562
     * By default, the PHYML option is used for the NNITopologySearch object, and numerical parameters are re-estimated
 
563
     * every 4 NNI runs (as in the phyml software).
 
564
     *
 
565
     * The optimizeNumericalParameters method is used for estimating numerical parameters.
 
566
     * The tolerance passed to this function is specified as input parameters.
 
567
     * They are generally very high to avoid local optima.
 
568
     *
 
569
     * @param tl                A pointer toward the TreeLikelihood object to optimize.
 
570
     * @param parameters        The list of parameters to optimize. Use tl->getIndependentParameters() in order to estimate all parameters.
 
571
     * @param optimizeNumFirst  Tell if we must optimize numerical parameters before searching topology.
 
572
     * @param tolBefore         The tolerance to use when estimating numerical parameters before topology search (if optimizeNumFirst is set to 'true').
 
573
     * @param tolDuring         The tolerance to use when estimating numerical parameters during the topology search.
 
574
     * @param tlEvalMax         The maximum number of function evaluations.
 
575
     * @param numStep           Number of NNI rounds before re-estimating numerical parameters.
 
576
     * @param messageHandler    The massage handler.
 
577
     * @param profiler          The profiler.
 
578
     * @param reparametrization Tell if parameters should be transformed in order to remove constraints.
 
579
     *                          This can improve optimization, but is a bit slower.
 
580
     * @param verbose           The verbose level.
 
581
     * @param optMethod         Option passed to optimizeNumericalParameters.
 
582
     * @param nStep             Option passed to optimizeNumericalParameters.
 
583
     * @param nniMethod         NNI algorithm to use.
 
584
     * @return A pointer toward the final likelihood object.
 
585
     * This pointer may be the same as passed in argument (tl), but in some cases the algorithm
 
586
     * clone this object. We may change this bahavior in the future...
 
587
     * You hence should write something like
 
588
     * @code
 
589
     * tl = OptimizationTools::optimizeTreeNNI(tl, ...);
 
590
     * @endcode
 
591
     * @throw Exception any exception thrown by the optimizer.
 
592
     */
 
593
    static NNIHomogeneousTreeLikelihood* optimizeTreeNNI(
 
594
                                                         NNIHomogeneousTreeLikelihood* tl,
 
595
                                                         const ParameterList& parameters,
 
596
                                                         bool optimizeNumFirst        = true,
 
597
                                                         double tolBefore             = 100,
 
598
                                                         double tolDuring             = 100,
 
599
                                                         int tlEvalMax                = 1000000,
 
600
                                                         unsigned int numStep         = 1,
 
601
                                                         OutputStream* messageHandler = ApplicationTools::message,
 
602
                                                         OutputStream* profiler       = ApplicationTools::message,
 
603
                                                         bool reparametrization       = false,
 
604
                                                         unsigned int verbose         = 1,
 
605
                                                         const std::string& optMethod = OptimizationTools::OPTIMIZATION_NEWTON,
 
606
                                                         unsigned int nStep           = 1,
 
607
                                                         const std::string& nniMethod = NNITopologySearch::PHYML)
 
608
      throw (Exception);
 
609
 
 
610
    /**
 
611
     * @brief Optimize all parameters from a TreeLikelihood object, including tree topology using Nearest Neighbor Interchanges.
 
612
     *
 
613
     * This function takes as input a TreeLikelihood object implementing the NNISearchable interface.
 
614
     *
 
615
     * Details:
 
616
     * A NNITopologySearch object is instanciated and is associated an additional TopologyListener.
 
617
     * This listener is used to re-estimate numerical parameters after one or several topology change.
 
618
     * By default, the PHYML option is used for the NNITopologySearch object, and numerical parameters are re-estimated
 
619
     * every 4 NNI runs (as in the phyml software).
 
620
     *
 
621
     * The optimizeNumericalParameters2 method is used for estimating numerical parameters.
 
622
     * The tolerance passed to this function is specified as input parameters.
 
623
     * They are generally very high to avoid local optima.
 
624
     *
 
625
     * @param tl                A pointer toward the TreeLikelihood object to optimize.
 
626
     * @param parameters        The list of parameters to optimize. Use tl->getIndependentParameters() in order to estimate all parameters.
 
627
     * @param optimizeNumFirst  Tell if we must optimize numerical parameters before searching topology.
 
628
     * @param tolBefore         The tolerance to use when estimating numerical parameters before topology search (if optimizeNumFirst is set to 'true').
 
629
     * @param tolDuring         The tolerance to use when estimating numerical parameters during the topology search.
 
630
     * @param tlEvalMax         The maximum number of function evaluations.
 
631
     * @param numStep           Number of NNI rounds before re-estimating numerical parameters.
 
632
     * @param messageHandler    The massage handler.
 
633
     * @param profiler          The profiler.
 
634
     * @param reparametrization Tell if parameters should be transformed in order to remove constraints.
 
635
     *                          This can improve optimization, but is a bit slower.
 
636
     * @param verbose           The verbose level.
 
637
     * @param optMethod         Option passed to optimizeNumericalParameters2.
 
638
     * @param nniMethod         NNI algorithm to use.
 
639
     * @return A pointer toward the final likelihood object.
 
640
     * This pointer may be the same as passed in argument (tl), but in some cases the algorithm
 
641
     * clone this object. We may change this bahavior in the future...
 
642
     * You hence should write something like
 
643
     * @code
 
644
     * tl = OptimizationTools::optimizeTreeNNI2(tl, ...);
 
645
     * @endcode
 
646
     * @throw Exception any exception thrown by the optimizer.
 
647
     */
 
648
    static NNIHomogeneousTreeLikelihood* optimizeTreeNNI2(
 
649
                                                          NNIHomogeneousTreeLikelihood* tl,
 
650
                                                          const ParameterList& parameters,
 
651
                                                          bool optimizeNumFirst        = true,
 
652
                                                          double tolBefore             = 100,
 
653
                                                          double tolDuring             = 100,
 
654
                                                          int tlEvalMax                = 1000000,
 
655
                                                          unsigned int numStep         = 1,
 
656
                                                          OutputStream* messageHandler = ApplicationTools::message,
 
657
                                                          OutputStream* profiler       = ApplicationTools::message,
 
658
                                                          bool reparametrization       = false,
 
659
                                                          unsigned int verbose         = 1,
 
660
                                                          const std::string& optMethod = OptimizationTools::OPTIMIZATION_NEWTON,
 
661
                                                          const std::string& nniMethod = NNITopologySearch::PHYML)
 
662
      throw (Exception);
 
663
 
 
664
    /**
 
665
     * @brief Optimize tree topology from a DRTreeParsimonyScore using Nearest Neighbor Interchanges.
 
666
     *
 
667
     * @param tp               A pointer toward the DRTreeParsimonyScore object to optimize.
 
668
     * @param verbose          The verbose level.
 
669
     * @return A pointer toward the final parsimony score object.
 
670
     * This pointer may be the same as passed in argument (tl), but in some cases the algorithm
 
671
     * clone this object. We may change this bahavior in the future...
 
672
     * You hence should write something like
 
673
     * @code
 
674
     * tp = OptimizationTools::optimizeTreeNNI(tp, ...);
 
675
     * @endcode
 
676
     */
 
677
    
 
678
    static DRTreeParsimonyScore* optimizeTreeNNI(
 
679
                                                 DRTreeParsimonyScore* tp,
 
680
                                                 unsigned int verbose = 1);
 
681
 
 
682
    /**
 
683
     * @brief Build a tree using a distance method.
 
684
     *
 
685
     * This method estimate a distance matrix using a DistanceEstimation object, and then builds the phylogenetic tree using a AgglomerativeDistanceMethod object.
 
686
     * The main issue here is to estimate non-branch lengths parameters, as substitution model and rate distribution parameters.
 
687
     * Three options are provideed here:
 
688
     * - DISTANCEMETHOD_INIT (default) keep parameters to there initial value,
 
689
     * - DISTANCEMETHOD_PAIRWISE estimated parameters in a pairwise manner, which is standard but not good at all...
 
690
     * - DISTANCEMETHOD_ITERATIONS uses Ninio et al's iterative algorithm, which uses Maximum Likelihood to estimate these parameters, and then update the distance matrix.
 
691
     * Ninio M, Privman E, Pupko T, Friedman N. 
 
692
     * Phylogeny reconstruction: increasing the accuracy of pairwise distance estimation using Bayesian inference of evolutionary rates.
 
693
     * Bioinformatics. 2007 Jan 15;23(2):e136-41.
 
694
     *
 
695
     * @param estimationMethod The distance estimation object to use.
 
696
     * @param reconstructionMethod The tree reconstruction object to use.
 
697
     * @param parametersToIgnore A list of parameters to ignore while optimizing parameters.
 
698
     * @param optimizeBrLen Tell if branch lengths should be optimized together with other parameters. This may lead to more accurate parameter estimation, but is slower.
 
699
     * @param rooted Tell if rooted trees must be produced.
 
700
     * @param param String describing the type of optimization to use.
 
701
     * @param tolerance Threshold on likelihood for stopping the iterative procedure. Used only with param=DISTANCEMETHOD_ITERATIONS.
 
702
     * @param tlEvalMax Maximum number of likelihood computations to perform when optimizing parameters. Used only with param=DISTANCEMETHOD_ITERATIONS.
 
703
     * @param profiler Output stream used by optimizer. Used only with param=DISTANCEMETHOD_ITERATIONS.
 
704
     * @param messenger Output stream used by optimizer. Used only with param=DISTANCEMETHOD_ITERATIONS.
 
705
     * @param verbose Verbose level.
 
706
     */
 
707
    static TreeTemplate<Node>* buildDistanceTree(
 
708
                                                 DistanceEstimation& estimationMethod,
 
709
                                                 AgglomerativeDistanceMethod& reconstructionMethod,
 
710
                                                 const ParameterList& parametersToIgnore,
 
711
                                                 bool optimizeBrLen = false,
 
712
                                                 bool rooted = false,
 
713
                                                 const std::string& param = DISTANCEMETHOD_INIT,
 
714
                                                 double tolerance = 0.000001,
 
715
                                                 unsigned int tlEvalMax = 1000000,
 
716
                                                 OutputStream* profiler = 0,
 
717
                                                 OutputStream* messenger = 0,
 
718
                                                 unsigned int verbose = 0) throw (Exception);
 
719
 
 
720
  public:
 
721
    static std::string DISTANCEMETHOD_INIT;
 
722
    static std::string DISTANCEMETHOD_PAIRWISE;
 
723
    static std::string DISTANCEMETHOD_ITERATIONS;
 
724
  };
 
725
 
 
726
 
 
727
} //end of namespace bpp.
 
728
 
 
729
#endif  //_OPTIMIZATIONTOOLS_H_
 
730