fast.cpp 25 KB
Newer Older
Iakov Davydov's avatar
Iakov Davydov committed
1
2
3
4
5
/// @mainpage FastCodeML
///
/// @section intro_sect Introduction
///
/// FastCodeML is a rewrite of CodeML based directly on the pseudocode document.
Iakov Davydov's avatar
Iakov Davydov committed
6
7
8
9
/// It incorporates various parallelization strategies to be able to exploit
/// modern HPC machines architecture.
/// For this reason there are various parts of the code that can be selected at
/// compile time or run time to experiment with various, possible solutions.
Iakov Davydov's avatar
Iakov Davydov committed
10
11
12
///
/// @section contacts_sect Contacts
///
Iakov Davydov's avatar
Iakov Davydov committed
13
14
/// Contact us if you want more information on the project, want to collaborate
/// or suggest new ideas.
Iakov Davydov's avatar
Iakov Davydov committed
15
///
Iakov Davydov's avatar
Iakov Davydov committed
16
17
18
19
///- Ing. <a href="mailto:mvalle@cscs.ch">Mario Valle</a> - Swiss National
/// Supercomputing Centre (CSCS) - Switzerland
///- The HP2C <a href="mailto:selectome@hp2c.ch">Selectome</a> Project Group -
/// Mainly based in University of Lausanne - Switzerland
Iakov Davydov's avatar
Iakov Davydov committed
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
///

#include <iostream>
#include <iomanip>
#include <limits>
#include "CmdLine.h"
#include "Newick.h"
#include "Phylip.h"
#include "BayesTest.h"
#include "Forest.h"
#include "Exceptions.h"
#include "BranchSiteModel.h"
#include "ParseParameters.h"
#include "VerbosityLevels.h"
#include "WriteResults.h"

#ifndef VTRACE
#ifdef _OPENMP
#include <omp.h>
#endif
#endif
#ifdef USE_MKL_VML
#include <mkl_vml.h>
#endif
#include "Timer.h"
#ifdef USE_MPI
#include "HighLevelCoordinator.h"
#endif

/// Main program for FastCodeML.
///
Iakov Davydov's avatar
Iakov Davydov committed
51
52
53
///     @author Mario Valle - Swiss National Supercomputing Centre (CSCS)
///     @date 2010-12-22 (initial version)
///     @version 1.1
Iakov Davydov's avatar
Iakov Davydov committed
54
55
56
57
///
///	@param[in] aRgc Number of command line parameters
/// @param[in] aRgv Command line parameters
///
Iakov Davydov's avatar
Iakov Davydov committed
58
const char *version = "1.2.0";
Iakov Davydov's avatar
Iakov Davydov committed
59

Iakov Davydov's avatar
Iakov Davydov committed
60
61
int main(int aRgc, char **aRgv) {
  try {
Iakov Davydov's avatar
Iakov Davydov committed
62
63

#ifdef USE_MKL_VML
Iakov Davydov's avatar
Iakov Davydov committed
64
65
    // If used, intitialize the MKL VML library
    vmlSetMode(VML_HA | VML_DOUBLE_CONSISTENT);
Iakov Davydov's avatar
Iakov Davydov committed
66
67
68
#endif

#ifdef USE_MPI
Iakov Davydov's avatar
Iakov Davydov committed
69
70
    // Start the high level parallel executor (based on MPI)
    HighLevelCoordinator hlc(&aRgc, &aRgv);
Iakov Davydov's avatar
Iakov Davydov committed
71
72
#endif

Iakov Davydov's avatar
Iakov Davydov committed
73
74
75
    // Parse the command line
    CmdLine cmd;
    cmd.parseCmdLine(aRgc, aRgv);
Iakov Davydov's avatar
Iakov Davydov committed
76

Iakov Davydov's avatar
Iakov Davydov committed
77
// Adjust and report the number of threads that will be used
Iakov Davydov's avatar
Iakov Davydov committed
78
#ifdef _OPENMP
Iakov Davydov's avatar
Iakov Davydov committed
79
80
81
82
83
84
85
86
87
88
    int num_threads = omp_get_max_threads();
    // std::cout<<"max num of thr: "<< num_threads <<std::endl;

    if ((cmd.mNumThreads >= 1) &&
        (cmd.mNumThreads <= (unsigned int)num_threads))
      num_threads = cmd.mNumThreads;
    // std::cout<<"num of thr: "<< num_threads <<std::endl;

    omp_set_num_threads(num_threads);
/*if (num_threads < 2)
Iakov Davydov's avatar
Iakov Davydov committed
89
    cmd.mForceSerial = true;
Iakov Davydov's avatar
Iakov Davydov committed
90
 else
Iakov Davydov's avatar
Iakov Davydov committed
91
    cmd.mForceSerial = false;*/
Iakov Davydov's avatar
Iakov Davydov committed
92
#else
Iakov Davydov's avatar
Iakov Davydov committed
93
94
95
    cmd.mNumThreads = 1;
    int num_threads = 1;
    cmd.mForceSerial = true;
Iakov Davydov's avatar
Iakov Davydov committed
96
97
98
#endif

/*#ifdef _OPENMP
Iakov Davydov's avatar
Iakov Davydov committed
99
100
101
102
103
104
105
        int num_threads = omp_get_max_threads();
        if(num_threads < 2 || cmd.mForceSerial)
        {
                cmd.mForceSerial = true;
                num_threads = 1;
                omp_set_num_threads(1);
        }
Iakov Davydov's avatar
Iakov Davydov committed
106
#else
Iakov Davydov's avatar
Iakov Davydov committed
107
108
        cmd.mForceSerial = true;
        int num_threads = 1;
Iakov Davydov's avatar
Iakov Davydov committed
109
110
111
#endif*/

#ifdef USE_MPI
Iakov Davydov's avatar
Iakov Davydov committed
112
113
114
    // Shutdown messages from all MPI processes except the master
    if (!hlc.isMaster())
      cmd.mVerboseLevel = VERBOSE_NONE;
Iakov Davydov's avatar
Iakov Davydov committed
115
116
#endif

Iakov Davydov's avatar
Iakov Davydov committed
117
118
    //    std::cout <<std::endl<<"------------------"<< std::endl<<"FastCodeML
    //    V"<<version<<std::endl<<"------------------"<<std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
119
120
121
122
123
124
125
126
    // Write out command line parameters (if not quiet i.e. if verbose level >
    // 0)
    if (cmd.mVerboseLevel >= VERBOSE_INFO_OUTPUT) {

      std::cout << "------------------------------------" << std::endl;
      std::cout << "FastCodeML V" << version << std::endl;
      std::cout << "------------------------------------" << std::endl;
      std::cout << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
127
128
      std::cout << "Tree file:      " << cmd.mTreeFile << std::endl;
      std::cout << "Gene file:      " << cmd.mGeneFile << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
129
130
131
      std::cout << "Verbose level:  " << cmd.mVerboseLevel << " ("
                << decodeVerboseLevel(cmd.mVerboseLevel) << ')' << std::endl;
      if (cmd.mSeed)
Iakov Davydov's avatar
Iakov Davydov committed
132
        std::cout << "Seed:           " << cmd.mSeed << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
133
      if (cmd.mBranchFromFile)
Iakov Davydov's avatar
Iakov Davydov committed
134
        std::cout << "Branch:         From tree file" << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
135
136
      else if (cmd.mBranchStart != UINT_MAX &&
               cmd.mBranchStart == cmd.mBranchEnd)
Iakov Davydov's avatar
Iakov Davydov committed
137
        std::cout << "Branch:         " << cmd.mBranchStart << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
138
      else if (cmd.mBranchStart != UINT_MAX && cmd.mBranchEnd == UINT_MAX)
Iakov Davydov's avatar
Iakov Davydov committed
139
        std::cout << "Branches:       " << cmd.mBranchStart << "-end"
Iakov Davydov's avatar
Iakov Davydov committed
140
141
                  << std::endl;
      else if (cmd.mBranchStart != UINT_MAX && cmd.mBranchEnd != UINT_MAX)
Iakov Davydov's avatar
Iakov Davydov committed
142
        std::cout << "Branches:       " << cmd.mBranchStart << '-'
Iakov Davydov's avatar
Iakov Davydov committed
143
144
                  << cmd.mBranchEnd << std::endl;
      if (!cmd.mStopIfNotLRT)
Iakov Davydov's avatar
Iakov Davydov committed
145
        std::cout << "H0 pre stop:    No" << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
146
      if (cmd.mIgnoreFreq)
Iakov Davydov's avatar
Iakov Davydov committed
147
        std::cout << "Codon freq.:    Ignore" << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
      if (cmd.mDoNotReduceForest)
        std::cout << "Reduce forest:  Do not reduce" << std::endl;
      else
        std::cout << "Reduce forest:  Aggressive" << std::endl;
      if (cmd.mInitH0fromH1)
        std::cout << "Starting val.:  From H1" << std::endl;
      else if (cmd.mInitFromParams && cmd.mBranchLengthsFromFile)
        std::cout << "Starting val.:  Times from tree file and params from "
                     "const (see below)" << std::endl;
      else if (cmd.mInitFromParams)
        std::cout << "Starting val.:  Params from const (see below)"
                  << std::endl;
      else if (cmd.mBranchLengthsFromFile)
        std::cout << "Starting val.:  Times from tree file" << std::endl;
      if (cmd.mNoMaximization)
Iakov Davydov's avatar
Iakov Davydov committed
163
        std::cout << "Maximization:   No" << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
164
      if (cmd.mTrace)
Iakov Davydov's avatar
Iakov Davydov committed
165
        std::cout << "Trace:          On" << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
166
      if (cmd.mCleanData)
Iakov Davydov's avatar
Iakov Davydov committed
167
        std::cout << "Clean data:     On" << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
168
      else
Iakov Davydov's avatar
Iakov Davydov committed
169
        std::cout << "Clean data:     Off" << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
170
      if (cmd.mGraphFile)
Iakov Davydov's avatar
Iakov Davydov committed
171
        std::cout << "Graph file:     " << cmd.mGraphFile << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
172
      if (cmd.mGraphFile && cmd.mExportComputedTimes != UINT_MAX)
Iakov Davydov's avatar
Iakov Davydov committed
173
        std::cout << "Graph times:    From H" << cmd.mExportComputedTimes
Iakov Davydov's avatar
Iakov Davydov committed
174
175
                  << std::endl;
      if (!cmd.mNoMaximization)
Iakov Davydov's avatar
Iakov Davydov committed
176
        std::cout << "Optimizer:      " << cmd.mOptimizationAlgo << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
177
178
179
      if (cmd.mMaxIterations != MAX_ITERATIONS)
        std::cout << "Max iterations: " << cmd.mMaxIterations << std::endl;
      if (cmd.mDeltaValueForGradient > 0.0)
Iakov Davydov's avatar
Iakov Davydov committed
180
        std::cout << "Delta value:    " << cmd.mDeltaValueForGradient
Iakov Davydov's avatar
Iakov Davydov committed
181
182
183
                  << std::endl;
      std::cout << "Relative error: " << cmd.mRelativeError << std::endl;
      if (cmd.mResultsFile)
Iakov Davydov's avatar
Iakov Davydov committed
184
        std::cout << "Results file:   " << cmd.mResultsFile << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
185
186
      if (cmd.mFixedBranchLength)
        std::cerr << "Branch lengths are fixed" << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
187
#ifdef _OPENMP
Iakov Davydov's avatar
Iakov Davydov committed
188
      if (num_threads > 1) {
Iakov Davydov's avatar
Iakov Davydov committed
189
190
        std::cout << "Num. threads:   " << num_threads << std::endl
                  << "Num. cores:     " << omp_get_num_procs() << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
191
      } else
Iakov Davydov's avatar
Iakov Davydov committed
192
#endif
Iakov Davydov's avatar
Iakov Davydov committed
193
      {
Iakov Davydov's avatar
Iakov Davydov committed
194
195
        std::cout << "Num. threads:   1 serial" << std::endl
                  << "Num. cores:     1" << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
196
      }
Iakov Davydov's avatar
Iakov Davydov committed
197
#ifdef USE_MPI
Iakov Davydov's avatar
Iakov Davydov committed
198
199
200
201
202
203
      if (hlc.numJobs() > 2)
        std::cout << "Num. MPI proc:  1 (master) + " << hlc.numJobs() - 1
                  << " (workers)" << std::endl;
      else
        std::cout << "Num. MPI proc:  Insufficient, single task execution"
                  << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
204
#endif
Iakov Davydov's avatar
Iakov Davydov committed
205
      std::cout << "Compiled with:  ";
Iakov Davydov's avatar
Iakov Davydov committed
206
#ifdef _OPENMP
Iakov Davydov's avatar
Iakov Davydov committed
207
      std::cout << "USE_OPENMP ";
Iakov Davydov's avatar
Iakov Davydov committed
208
209
#endif
#ifdef USE_MPI
Iakov Davydov's avatar
Iakov Davydov committed
210
      std::cout << "USE_MPI ";
Iakov Davydov's avatar
Iakov Davydov committed
211
212
#endif
#ifdef USE_CPV_SCALING
Iakov Davydov's avatar
Iakov Davydov committed
213
      std::cout << "USE_CPV_SCALING ";
Iakov Davydov's avatar
Iakov Davydov committed
214
215
#endif
#ifdef NEW_LIKELIHOOD
Iakov Davydov's avatar
Iakov Davydov committed
216
      std::cout << "NEW_LIKELIHOOD ";
Iakov Davydov's avatar
Iakov Davydov committed
217
218
#endif
#ifdef NON_RECURSIVE_VISIT
Iakov Davydov's avatar
Iakov Davydov committed
219
      std::cout << "NON_RECURSIVE_VISIT ";
Iakov Davydov's avatar
Iakov Davydov committed
220
221
#endif
#ifdef USE_DAG
Iakov Davydov's avatar
Iakov Davydov committed
222
      std::cout << "USE_DAG ";
Iakov Davydov's avatar
Iakov Davydov committed
223
224
#endif
#ifdef USE_ORIGINAL_PROPORTIONS
Iakov Davydov's avatar
Iakov Davydov committed
225
      std::cout << "USE_ORIGINAL_PROPORTIONS ";
Iakov Davydov's avatar
Iakov Davydov committed
226
227
#endif
#ifdef USE_LAPACK
Iakov Davydov's avatar
Iakov Davydov committed
228
      std::cout << "USE_LAPACK ";
Iakov Davydov's avatar
Iakov Davydov committed
229
230
#endif
#ifdef USE_MKL_VML
Iakov Davydov's avatar
Iakov Davydov committed
231
      std::cout << "USE_MKL_VML";
Iakov Davydov's avatar
Iakov Davydov committed
232
#endif
Iakov Davydov's avatar
Iakov Davydov committed
233
234
235
236
237
238
239
240
241
242
243
      std::cout << std::endl
                << std::endl;
      if (cmd.mInitFromParams) {
        std::cout << "Param initial values:" << std::endl
                  << std::endl
                  << ParseParameters::getInstance();
      }
    }

// Initialize the random number generator (0 means it is not set on the command
// line)
Iakov Davydov's avatar
Iakov Davydov committed
244
#ifdef USE_MPI
Iakov Davydov's avatar
Iakov Davydov committed
245
246
247
248
    // Insure that each MPI process starts with a different seed
    if (cmd.mSeed == 0)
      cmd.mSeed = static_cast<unsigned int>(time(NULL)) +
                  static_cast<unsigned int>(hlc.getRank()) * 1000;
Iakov Davydov's avatar
Iakov Davydov committed
249
#else
Iakov Davydov's avatar
Iakov Davydov committed
250
251
    if (cmd.mSeed == 0)
      cmd.mSeed = static_cast<unsigned int>(time(NULL));
Iakov Davydov's avatar
Iakov Davydov committed
252
#endif
Iakov Davydov's avatar
Iakov Davydov committed
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
    srand(cmd.mSeed);

    // Verify the optimizer algorithm selected on the command line
    if (!cmd.mNoMaximization)
      BranchSiteModel::verifyOptimizerAlgo(cmd.mOptimizationAlgo);

    // Start a timer (to measure serial part over parallel one)
    Timer timer;
    if (cmd.mVerboseLevel >= VERBOSE_INFO_OUTPUT)
      timer.start();

    // Create the forest
    Forest forest(cmd.mVerboseLevel);

    // Enclose file loading into a block so temporary structures could be
    // deleted when no more needed
    {
      // Load the multiple sequence alignment (MSA)
      Phylip msa(cmd.mVerboseLevel);
      msa.readFile(cmd.mGeneFile, cmd.mCleanData);

      // Load the phylogenetic tree
      Newick tree(cmd.mVerboseLevel);
      tree.readFile(cmd.mTreeFile);

      // Check coherence between the two files
      msa.checkNameCoherence(tree.getSpecies());

      // Check root
      tree.checkRootBranches();

      // If times from file then check for null branch lengths for any leaf
      if (cmd.mBranchLengthsFromFile) {
        int zero_on_leaf_cnt = 0;
        int zero_on_int_cnt = 0;
        tree.countNullBranchLengths(zero_on_leaf_cnt, zero_on_int_cnt);

        if (zero_on_leaf_cnt > 0 || zero_on_int_cnt > 0) {
          std::cout
              << "Found null or missing branch length in tree file. On leaves: "
              << zero_on_leaf_cnt
              << "  on internal branches: " << zero_on_int_cnt << std::endl;
        }

        if (zero_on_leaf_cnt > 0) {
          throw FastCodeMLFatal("Null or missing branch length in tree file");
        }
      }

      // Print the tree with the numbering of internal branches
      if (cmd.mVerboseLevel >= VERBOSE_INFO_OUTPUT)
        tree.printTreeAnnotated(std::cout);

      // Load the forest
      forest.loadTreeAndGenes(
          tree, msa, cmd.mIgnoreFreq ? CodonFrequencies::CODON_FREQ_MODEL_UNIF
                                     : CodonFrequencies::CODON_FREQ_MODEL_F3X4);
    }

    // Reduce the forest merging common subtrees. Add also more reduction, then
    // clean the no more useful data.
    if (!cmd.mDoNotReduceForest) {
      // bool sts = forest.reduceSubtrees(cmd.mNumReductionBlocks);
      forest.reduceSubtrees();
Iakov Davydov's avatar
Iakov Davydov committed
317
318

#ifndef NEW_LIKELIHOOD
Iakov Davydov's avatar
Iakov Davydov committed
319
      forest.addAggressiveReduction();
Iakov Davydov's avatar
Iakov Davydov committed
320
#endif
Iakov Davydov's avatar
Iakov Davydov committed
321
      forest.cleanReductionWorkingData();
Iakov Davydov's avatar
Iakov Davydov committed
322
#ifdef NEW_LIKELIHOOD
Iakov Davydov's avatar
Iakov Davydov committed
323
      forest.prepareNewReduction();
Iakov Davydov's avatar
Iakov Davydov committed
324
#endif
Iakov Davydov's avatar
Iakov Davydov committed
325
    }
Iakov Davydov's avatar
Iakov Davydov committed
326
#ifdef NEW_LIKELIHOOD
Iakov Davydov's avatar
Iakov Davydov committed
327
328
329
    else {
      forest.prepareNewReductionNoReuse();
    }
Iakov Davydov's avatar
Iakov Davydov committed
330
331
332
#endif

#ifdef NON_RECURSIVE_VISIT
Iakov Davydov's avatar
Iakov Davydov committed
333
334
    // Prepare the pointers to visit the trees without recursion
    forest.prepareNonRecursiveVisit();
Iakov Davydov's avatar
Iakov Davydov committed
335
336
#endif

Iakov Davydov's avatar
Iakov Davydov committed
337
338
// Subdivide the trees in groups based on dependencies
// forest.prepareDependencies(cmd.mForceSerial || cmd.mDoNotReduceForest);
Iakov Davydov's avatar
Iakov Davydov committed
339
340

#ifdef USE_DAG
Iakov Davydov's avatar
Iakov Davydov committed
341
342
    // Load the forest into a DAG
    forest.loadForestIntoDAG(Nt);
Iakov Davydov's avatar
Iakov Davydov committed
343
344
#endif

Iakov Davydov's avatar
Iakov Davydov committed
345
346
347
348
349
350
351
    // Get the time needed by data preprocessing
    if (cmd.mVerboseLevel >= VERBOSE_INFO_OUTPUT) {
      timer.stop();
      std::cout << std::endl
                << "TIMER (preprocessing) ncores: " << std::setw(2)
                << num_threads << " time: " << timer.get() << std::endl;
    }
Iakov Davydov's avatar
Iakov Davydov committed
352

Iakov Davydov's avatar
Iakov Davydov committed
353
354
355
    // Print few statistics
    if (cmd.mVerboseLevel >= VERBOSE_INFO_OUTPUT)
      std::cout << forest;
Iakov Davydov's avatar
Iakov Davydov committed
356
357

#ifdef USE_MPI
Iakov Davydov's avatar
Iakov Davydov committed
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
    // Distribute the work. If run under MPI then finish, else return to the
    // standard execution flow
    if (cmd.mVerboseLevel >= VERBOSE_INFO_OUTPUT)
      timer.start();
    bool has_run_under_MPI = hlc.startWork(forest, cmd);

    // If executed under MPI report the time spent, otherwise stop the timer so
    // it can be restarted around the serial execution
    if (has_run_under_MPI) {
      if (cmd.mVerboseLevel >= VERBOSE_INFO_OUTPUT) {
        timer.stop();
        std::cout << std::endl
                  << "TIMER (processing) ncores: " << std::setw(2)
                  << num_threads * (hlc.numJobs() - 1) + 1
                  << " time: " << timer.get() << std::endl;
      }
      return 0;
    } else {
      timer.stop();
    }
Iakov Davydov's avatar
Iakov Davydov committed
378
379
#endif

Iakov Davydov's avatar
Iakov Davydov committed
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
    // Initialize the output results file (if the argument is null, no file is
    // created)
    WriteResults output_results(cmd.mResultsFile);

    // Compute the range of branches to mark as foreground
    size_t branch_start, branch_end;
    forest.getBranchRange(cmd, branch_start, branch_end);

    // Start timing parallel part
    if (cmd.mVerboseLevel >= VERBOSE_INFO_OUTPUT)
      timer.start();

    // Initialize the models
    BranchSiteModelNullHyp h0(forest, cmd);
    BranchSiteModelAltHyp h1(forest, cmd);

    // Initialize the test
    BayesTest beb(forest, cmd.mVerboseLevel, cmd.mDoNotReduceForest);

    // For all requested internal branches
    for (size_t fg_branch = branch_start; fg_branch <= branch_end;
         ++fg_branch) {
      if (cmd.mVerboseLevel >= VERBOSE_ONLY_RESULTS)
        std::cout << std::endl
                  << "Doing branch " << fg_branch << std::endl;

      // Compute the alternate model maximum loglikelihood
      double lnl1 = 0.;
      if (cmd.mComputeHypothesis != 0) {
        if (cmd.mInitFromParams)
          h1.initFromParams();
        if (cmd.mBranchLengthsFromFile)
          h1.initFromTree();

        lnl1 = h1(fg_branch);

        // Save the value for formatted output
        output_results.saveLnL(fg_branch, lnl1, 1);
      }

      // Compute the null model maximum loglikelihood
      double lnl0 = 0.;
      if (cmd.mComputeHypothesis != 1) {
        if (cmd.mInitH0fromH1)
          h0.initFromResult(h1.getVariables());
        else {
          if (cmd.mInitFromParams)
            h0.initFromParams();
          if (cmd.mBranchLengthsFromFile)
            h0.initFromTree();
        }

        lnl0 = h0(fg_branch, cmd.mStopIfNotLRT && cmd.mComputeHypothesis != 0,
                  lnl1 - THRESHOLD_FOR_LRT);

        // Save the value for formatted output (only if has not be forced to
        // stop)
        if (lnl0 < DBL_MAX)
          output_results.saveLnL(fg_branch, lnl0, 0);
      }

      if (cmd.mVerboseLevel >= VERBOSE_ONLY_RESULTS) {
        std::cout << std::endl;
        if (cmd.mComputeHypothesis != 1) {
          std::cout << "LnL0: ";
          if (lnl0 == std::numeric_limits<double>::infinity())
            std::cout << "**Invalid result**";
          else if (lnl0 < DBL_MAX)
            std::cout << std::setprecision(15) << std::fixed << lnl0;
          else
            std::cout << "(Doesn't pass LRT, skipping)";
Iakov Davydov's avatar
Iakov Davydov committed
451
          std::cout << " Function calls: " << h0.getNumEvaluations() << "   ";
Iakov Davydov's avatar
Iakov Davydov committed
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
          std::cout << std::endl
                    << std::endl;
          if (lnl0 != std::numeric_limits<double>::infinity()) {
            std::string s0 = h0.printFinalVars(std::cout);
            // std::cout<<"EDW0: "<< s0 <<std::endl;
            output_results.saveParameters(fg_branch, s0, 0);
          }
          std::cout << std::endl;
        }
        if (cmd.mComputeHypothesis != 0) {
          std::cout << "LnL1: ";
          if (lnl1 == std::numeric_limits<double>::infinity())
            std::cout << "**Invalid result**";
          else
            std::cout << std::setprecision(15) << std::fixed << lnl1;
          std::cout << " Function calls: " << h1.getNumEvaluations();
          std::cout << std::endl
                    << std::endl;
          if (lnl1 != std::numeric_limits<double>::infinity()) {
            std::string s1 = h1.printFinalVars(std::cout);
            // std::cout<<"EDW1: "<< s1 <<std::endl;
            output_results.saveParameters(fg_branch, s1, 1);
          }
          std::cout << std::endl;
        }
        if (cmd.mComputeHypothesis > 1) {
          if (lnl0 == std::numeric_limits<double>::infinity() ||
              lnl1 == std::numeric_limits<double>::infinity())
            std::cout << "LRT: **Invalid result**";
          else if (lnl0 < DBL_MAX)
            std::cout << "LRT: " << std::setprecision(15) << std::fixed
                      << lnl1 - lnl0
Iakov Davydov's avatar
Iakov Davydov committed
484
485
                      << "  (threshold: " << std::setprecision(15) << std::fixed
                      << THRESHOLD_FOR_LRT << ')';
Iakov Davydov's avatar
Iakov Davydov committed
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
          else
            std::cout << "LRT: < " << std::setprecision(15) << std::fixed
                      << THRESHOLD_FOR_LRT;
          std::cout << std::endl;
        }
      }

      // If requested set the time in the forest and export to a graph
      // visualization tool
      if (cmd.mGraphFile) {
        switch (cmd.mExportComputedTimes) {
        case 0:
          h0.saveComputedTimes();
          break;

        case 1:
          h1.saveComputedTimes();
          break;

        default:
          break;
        }

        // Use the forest export class
        ForestExport fe(forest);
        fe.exportForest(cmd.mGraphFile, fg_branch);
      }

      // If the two hypothesis are computed, H0 has not been stopped and the run
      // passes the LRT, then compute the BEB
      if (cmd.mComputeHypothesis > 1 && lnl0 < DBL_MAX &&
          BranchSiteModel::performLRT(lnl0, lnl1)) {
        // Get the scale values from the latest optimized h1.
        std::vector<double> scales(2);
        h1.getScales(scales);

        // Run the BEB test
        beb.computeBEB(h1.getVariables(), fg_branch, scales);

        // Output the sites under positive selection (if any)
        if (cmd.mVerboseLevel >= VERBOSE_ONLY_RESULTS)
          beb.printPositiveSelSites(fg_branch);

        // Get the sites under positive selection for printing in the results
        // file (if defined)
        if (output_results.isWriteResultsEnabled()) {
          std::vector<unsigned int> positive_sel_sites;
          std::vector<double> positive_sel_sites_prob;
          beb.extractPositiveSelSites(positive_sel_sites,
                                      positive_sel_sites_prob);
          output_results.savePositiveSelSites(fg_branch, positive_sel_sites,
                                              positive_sel_sites_prob);
        }
      }
    }

    // Get the time needed by the parallel part
    if (cmd.mVerboseLevel >= VERBOSE_INFO_OUTPUT) {
      timer.stop();
      std::cout << std::endl
                << "TIMER (processing) ncores: " << std::setw(2) << num_threads
                << " time: " << timer.get() << std::endl;
    }

    // Output the results
    output_results.outputResults();

    ////////////////////////////////////////////////////////////////////
    // Catch all exceptions
  } catch (const FastCodeMLSuccess &) {
    return 0;
  } catch (const FastCodeMLFatal &e) {
    // If a message associated (i.e. no empty string), display it
    if (e.what()[0])
      std::cout << std::endl
                << e.what() << std::endl
                << std::endl;
    return 1;
  } catch (const FastCodeMLMemoryError &e) {
    std::cout << std::endl
              << e.what() << std::endl
              << std::endl;
    return 1;
  } catch (const std::bad_alloc &e) {
    std::cout << std::endl
              << e.what() << std::endl
              << std::endl;
    return 1;
  } catch (...) {
    std::cout << std::endl
              << "Default exception caught." << std::endl
              << std::endl;
    return 1;
  }

  return 0;
Iakov Davydov's avatar
Iakov Davydov committed
582
583
584
585
586
587
588
589
590
591
592
593
}

/// @page cppstd_page C++ Coding Standard
/// Here are collected few rules for coding this project.
///
/// @section cnames_sect Class names
/// Class names are CamelCase with first letter uppercase.
///
/// Ex: %PhyloTree
///
/// @section cmeth_sect Class methods
/// Class methods names are CamelCase with the first letter lowercase.
Iakov Davydov's avatar
Iakov Davydov committed
594
595
/// Only very common and specific names should be all lowercase, like read,
/// clean, size.
Iakov Davydov's avatar
Iakov Davydov committed
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
///
/// Ex: testFillQ
///
/// @section cdatamemb_sect Class data members
/// Class member variables names start with 'm' followed by CamelCase name.
///
/// Ex: mFgBranch
///
/// @section carg_sect Function arguments
/// Function arguments names start with 'a' followed by CamelCase name.
///
/// Ex: aFgBranch
///
/// @section const_sect Constants and enumeration
/// Constants and enumerations are all uppercase with words separated by '_'.
Iakov Davydov's avatar
Iakov Davydov committed
611
612
/// The first letters specify the kind of constant (like: STS_ for status, OPT_
/// for option value).
Iakov Davydov's avatar
Iakov Davydov committed
613
614
615
616
617
618
619
620
621
622
623
///
/// Ex: STS_CANT_OPEN
///
/// @section stack_sect Temporary variables
/// All the other variables are all lower case with parts separated by '_'.
///
/// Ex: branch_list
///
/// @section misc_sect Miscellaneous rules
/// In case of error main should return 1.
///
Iakov Davydov's avatar
Iakov Davydov committed
624
625
/// Array sizes and corresponding indexes should be size_t. The remaining
/// counters should be unsigned int.
Iakov Davydov's avatar
Iakov Davydov committed
626
627
628
629
630
631
632
633
///
/// The null pointer should be written as NULL, not 0 to make clear its purpose.
///

/**
@page cmd_page Command Line Switches
Here is a quick list of the valid command line switches for FastCodeML.

Iakov Davydov's avatar
Iakov Davydov committed
634
635
The input `tree_file` is in %Newick format with the file containing only one
tree. The `alignment_file` instead is in %Phylip format.
Iakov Davydov's avatar
Iakov Davydov committed
636
637
638
639

@verbatim

Usage:
Iakov Davydov's avatar
Iakov Davydov committed
640
    FastCodeML [options] tree_file alignment_file
Iakov Davydov's avatar
Iakov Davydov committed
641

Iakov Davydov's avatar
Iakov Davydov committed
642
643
644
-d  --debug  -v  --verbose (required argument)
        Verbosity level (0: none; 1: results only; 2: normal info; 3: MPI trace;
4: more debug) (default: 1)
Iakov Davydov's avatar
Iakov Davydov committed
645

Iakov Davydov's avatar
Iakov Davydov committed
646
647
-q  --quiet (no argument)
        No messages except results
Iakov Davydov's avatar
Iakov Davydov committed
648

Iakov Davydov's avatar
Iakov Davydov committed
649
650
-?  -h  --help (no argument)
        This help
Iakov Davydov's avatar
Iakov Davydov committed
651

Iakov Davydov's avatar
Iakov Davydov committed
652
653
-s  --seed (required argument)
        Random number generator seed (0 < seed < 1000000000)
Iakov Davydov's avatar
Iakov Davydov committed
654

Iakov Davydov's avatar
Iakov Davydov committed
655
656
-b  --branch (required argument)
        Do only this branch as foreground branch (count from 0)
Iakov Davydov's avatar
Iakov Davydov committed
657

Iakov Davydov's avatar
Iakov Davydov committed
658
659
660
-bs  --branch-start (required argument)
        Start computing from this branch as foreground one (count from 0)
(default: first one)
Iakov Davydov's avatar
Iakov Davydov committed
661

Iakov Davydov's avatar
Iakov Davydov committed
662
663
664
-be  --branch-end (required argument)
        End computing at this branch as foreground one (count from 0) (default:
last one)
Iakov Davydov's avatar
Iakov Davydov committed
665

Iakov Davydov's avatar
Iakov Davydov committed
666
667
-i  --ignore-freq (no argument)
        Ignore computed codon frequency and set all of them to 1/61
Iakov Davydov's avatar
Iakov Davydov committed
668

Iakov Davydov's avatar
Iakov Davydov committed
669
670
671
-e  --export (required argument)
        Export forest in GML format (if %03d or @03d is present, one is created
for each fg branch)
Iakov Davydov's avatar
Iakov Davydov committed
672

Iakov Davydov's avatar
Iakov Davydov committed
673
674
-nr  --no-reduce (no argument)
        Do not reduce forest by merging common subtrees
Iakov Davydov's avatar
Iakov Davydov committed
675

Iakov Davydov's avatar
Iakov Davydov committed
676
677
-l  --lengths-from-file  --times-from-file (no argument)
        Initial branch lengths from tree file
Iakov Davydov's avatar
Iakov Davydov committed
678

Iakov Davydov's avatar
Iakov Davydov committed
679
680
-o  --initial-step (no argument)
        Only the initial step is performed (no maximization)
Iakov Davydov's avatar
Iakov Davydov committed
681

Iakov Davydov's avatar
Iakov Davydov committed
682
683
684
-c  --export-comp-times (required argument)
        Export the computed times from H0 if 0, H1 if 1, otherwise the one read
in the phylo tree
Iakov Davydov's avatar
Iakov Davydov committed
685

Iakov Davydov's avatar
Iakov Davydov committed
686
687
-r  --trace (no argument)
        Trace the maximization run
Iakov Davydov's avatar
Iakov Davydov committed
688

Iakov Davydov's avatar
Iakov Davydov committed
689
690
-nt  --number-of-threads (required argument)
        Number of threads (1 for non parallel execution)
Iakov Davydov's avatar
Iakov Davydov committed
691

Iakov Davydov's avatar
Iakov Davydov committed
692
693
-bf  --branch-from-file (no argument)
        Do only the branch marked in the file as foreground branch
Iakov Davydov's avatar
Iakov Davydov committed
694

Iakov Davydov's avatar
Iakov Davydov committed
695
696
-hy  --only-hyp (required argument)
        Compute only H0 if 0, H1 if 1
Iakov Davydov's avatar
Iakov Davydov committed
697

Iakov Davydov's avatar
Iakov Davydov committed
698
699
-i1  --init-from-h1 (no argument)
        Start H0 optimization from H1 results
Iakov Davydov's avatar
Iakov Davydov committed
700

Iakov Davydov's avatar
Iakov Davydov committed
701
702
703
-m  --maximizer (required argument)
        Optimizer algorithm (0:LBFGS, 1:VAR1, 2:VAR2, 3:SLSQP, 11:BOBYQA,
22:FromCodeML, 99:MLSL_LDS) (default: 22)
Iakov Davydov's avatar
Iakov Davydov committed
704

Iakov Davydov's avatar
Iakov Davydov committed
705
706
-sd  --small-diff (required argument)
        Delta used in gradient computation (default: 1.49e-8)
Iakov Davydov's avatar
Iakov Davydov committed
707

Iakov Davydov's avatar
Iakov Davydov committed
708
709
710
-p  --init-param (required argument)
        Pass initialization parameter in the form: P=value (P: w0, k, p0, p1,
w2)
Iakov Davydov's avatar
Iakov Davydov committed
711

Iakov Davydov's avatar
Iakov Davydov committed
712
713
-ic  --init-default (no argument)
        Start from default parameter values and times from tree file
Iakov Davydov's avatar
Iakov Davydov committed
714

Iakov Davydov's avatar
Iakov Davydov committed
715
716
-x  --extra-debug (required argument)
        Extra debug parameter (zero disables it)
Iakov Davydov's avatar
Iakov Davydov committed
717

Iakov Davydov's avatar
Iakov Davydov committed
718
719
-re  --relative-error (required argument)
        Relative error where to stop maximization (default: 1e-3)
Iakov Davydov's avatar
Iakov Davydov committed
720

Iakov Davydov's avatar
Iakov Davydov committed
721
722
-ou  --output (required argument)
        Write results formatted to this file
Iakov Davydov's avatar
Iakov Davydov committed
723

Iakov Davydov's avatar
Iakov Davydov committed
724
725
-cl  --clean-data (no argument)
        Remove ambiguous or missing sites from the MSA (default: no)
Iakov Davydov's avatar
Iakov Davydov committed
726

Iakov Davydov's avatar
Iakov Davydov committed
727
728
-ps  --no-pre-stop (no argument)
        Don't stop H0 maximization even if it cannot satisfy LRT (default: stop)
Iakov Davydov's avatar
Iakov Davydov committed
729

Iakov Davydov's avatar
Iakov Davydov committed
730
731
-mi  --max-iterations (required argument)
        Maximum number of iterations for the maximizer (default: 10000)
Iakov Davydov's avatar
Iakov Davydov committed
732

Iakov Davydov's avatar
Iakov Davydov committed
733
734
-bl  --branch-lengths-fixed (no argument)
        The length of the brances is fixed
Iakov Davydov's avatar
Iakov Davydov committed
735
736
737
738
739

@endverbatim
*/

/// @page vampir_page Using Vampir for profiling
Iakov Davydov's avatar
Iakov Davydov committed
740
741
/// On Linux we use VampirTrace to collect profile data and Vampir to display
/// the results (http://www.vampir.eu/).
Iakov Davydov's avatar
Iakov Davydov committed
742
743
744
745
746
747
748
749
///
/// Before invoking CMAKE define CXX=vtCC
///
/// Define CMAKE_BUILD_TYPE as: RelWithDebInfo
///
/// Run CMAKE and configure.
///
/// Then define CMAKE_CXX_FLAGS as: -vt:mt -vt:noopari
Iakov Davydov's avatar
Iakov Davydov committed
750
751
/// If you want to trace also the OpenMP calls then change it to: -vt:mt
/// -vt:preprocess -DVTRACE
Iakov Davydov's avatar
Iakov Davydov committed
752
753
754
755
756
///
/// Then proceed as usual to build the executable.
///
/// Before running the executable, define the following environment variables:
///
Iakov Davydov's avatar
Iakov Davydov committed
757
758
759
760
761
///     export VT_BUFFER_SIZE=512M
///     export VT_MAX_FLUSHES=0
///     export VT_SYNCH_FLUSH=yes
///     export VT_GPUTRACE=no
///     export VT_UNIFY=no
Iakov Davydov's avatar
Iakov Davydov committed
762
///
Iakov Davydov's avatar
Iakov Davydov committed
763
764
/// Due to a VampirTrace bug, at the end of the execution, run the vtunify
/// executable by itself.
Iakov Davydov's avatar
Iakov Davydov committed
765
///
Iakov Davydov's avatar
Iakov Davydov committed
766
767
/// Now you can analyze the results by running vampir on the *.otf file
/// generated.
Iakov Davydov's avatar
Iakov Davydov committed
768
///