fast.cpp 22.6 KB
Newer Older
Iakov Davydov's avatar
Iakov Davydov committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
/// @mainpage FastCodeML
///
/// @section intro_sect Introduction
///
/// FastCodeML is a rewrite of CodeML based directly on the pseudocode document.
/// 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.
///
/// @section contacts_sect Contacts
///
/// Contact us if you want more information on the project, want to collaborate or suggest new ideas.
///
///- 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
///

#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
46
47
48
///		@author Mario Valle - Swiss National Supercomputing Centre (CSCS)
///		@date 2010-12-22 (initial version)
///		@version 1.1
Iakov Davydov's avatar
Iakov Davydov committed
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
///
///	@param[in] aRgc Number of command line parameters
/// @param[in] aRgv Command line parameters
///
const char* version="1.2.0";


int main(int aRgc, char **aRgv)
{
	try
	{

#ifdef USE_MKL_VML
	// If used, intitialize the MKL VML library
	vmlSetMode(VML_HA|VML_DOUBLE_CONSISTENT);
#endif

#ifdef USE_MPI
	// Start the high level parallel executor (based on MPI)
	HighLevelCoordinator hlc(&aRgc, &aRgv);
#endif

	// Parse the command line
	CmdLine cmd;
	cmd.parseCmdLine(aRgc, aRgv);

	// Adjust and report the number of threads that will be used
#ifdef _OPENMP
	int num_threads = omp_get_max_threads();
Iakov Davydov's avatar
Iakov Davydov committed
78
	//std::cout<<"max num of thr: "<< num_threads <<std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
79

Iakov Davydov's avatar
Iakov Davydov committed
80
81
82
	 if((cmd.mNumThreads >=1)&&(cmd.mNumThreads <= (unsigned int)num_threads))
		num_threads = cmd.mNumThreads;
	// std::cout<<"num of thr: "<< num_threads <<std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
83

Iakov Davydov's avatar
Iakov Davydov committed
84
85
86
87
88
	 omp_set_num_threads(num_threads);
	/*if (num_threads < 2)
		cmd.mForceSerial = true;
	 else
		cmd.mForceSerial = false;*/
Iakov Davydov's avatar
Iakov Davydov committed
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
#else
	cmd.mNumThreads=1;
	int num_threads = 1;
	cmd.mForceSerial = true;
#endif

/*#ifdef _OPENMP
	int num_threads = omp_get_max_threads();
	if(num_threads < 2 || cmd.mForceSerial)
	{
		cmd.mForceSerial = true;
		num_threads = 1;
		omp_set_num_threads(1);
	}
#else
	cmd.mForceSerial = true;
	int num_threads = 1;
#endif*/

#ifdef USE_MPI
	// Shutdown messages from all MPI processes except the master
	if(!hlc.isMaster()) cmd.mVerboseLevel = VERBOSE_NONE;
#endif

Iakov Davydov's avatar
Iakov Davydov committed
113
//	  std::cout <<std::endl<<"------------------"<< std::endl<<"FastCodeML V"<<version<<std::endl<<"------------------"<<std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
114
115
116
	// Write out command line parameters (if not quiet i.e. if verbose level > 0)
	if(cmd.mVerboseLevel >= VERBOSE_INFO_OUTPUT)
	{
117

Iakov Davydov's avatar
Iakov Davydov committed
118
119
120
121
		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
122
123
													std::cout << "Tree file:	  " << cmd.mTreeFile << std::endl;
													std::cout << "Gene file:	  " << cmd.mGeneFile << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
124
													std::cout << "Verbose level:  " << cmd.mVerboseLevel << " (" << decodeVerboseLevel(cmd.mVerboseLevel) << ')' << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
125
126
		if(cmd.mSeed)								std::cout << "Seed:			  " << cmd.mSeed << std::endl;
		if(cmd.mBranchFromFile)						std::cout << "Branch:		  From tree file" << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
127
		else if(cmd.mBranchStart != UINT_MAX && cmd.mBranchStart == cmd.mBranchEnd)
Iakov Davydov's avatar
Iakov Davydov committed
128
													std::cout << "Branch:		  " << cmd.mBranchStart << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
129
		else if(cmd.mBranchStart != UINT_MAX && cmd.mBranchEnd == UINT_MAX)
Iakov Davydov's avatar
Iakov Davydov committed
130
													std::cout << "Branches:		  " << cmd.mBranchStart << "-end" << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
131
		else if(cmd.mBranchStart != UINT_MAX && cmd.mBranchEnd != UINT_MAX)
Iakov Davydov's avatar
Iakov Davydov committed
132
133
134
													std::cout << "Branches:		  " << cmd.mBranchStart << '-' << cmd.mBranchEnd << std::endl;
		if(!cmd.mStopIfNotLRT)						std::cout << "H0 pre stop:	  No" << std::endl;
		if(cmd.mIgnoreFreq)							std::cout << "Codon freq.:	  Ignore" << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
135
136
137
138
139
140
141
		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;
Iakov Davydov's avatar
Iakov Davydov committed
142
143
144
145
146
		if(cmd.mNoMaximization)						std::cout << "Maximization:	  No" << std::endl;
		if(cmd.mTrace)								std::cout << "Trace:		  On" << std::endl;
		if(cmd.mCleanData)							std::cout << "Clean data:	  On" << std::endl;
		else										std::cout << "Clean data:	  Off" << std::endl;
		if(cmd.mGraphFile)							std::cout << "Graph file:	  " << cmd.mGraphFile << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
147
		if(cmd.mGraphFile && cmd.mExportComputedTimes != UINT_MAX)
Iakov Davydov's avatar
Iakov Davydov committed
148
149
													std::cout << "Graph times:	  From H" << cmd.mExportComputedTimes << std::endl;
		if(!cmd.mNoMaximization)					std::cout << "Optimizer:	  " << cmd.mOptimizationAlgo << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
150
		if(cmd.mMaxIterations != MAX_ITERATIONS)	std::cout << "Max iterations: " << cmd.mMaxIterations << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
151
		if(cmd.mDeltaValueForGradient > 0.0)		std::cout << "Delta value:	  " << cmd.mDeltaValueForGradient << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
152
													std::cout << "Relative error: " << cmd.mRelativeError << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
153
154
155
		if(cmd.mResultsFile)						std::cout << "Results file:	  " << cmd.mResultsFile << std::endl;
		if(cmd.mNumThreads)							std::cout << "Number of threads: " << cmd.mNumThreads << std::endl;
		if(cmd.mFixedBranchLength)						  std::cerr << "Branch lengths are fixed" << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
156
157
158
#ifdef _OPENMP
		if(num_threads > 1)
		{
Iakov Davydov's avatar
Iakov Davydov committed
159
160
													std::cout << "Num. threads:	  " << num_threads << std::endl
															  << "Num. cores:	  " << omp_get_num_procs() << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
161
162
163
164
		}
		else
#endif
		{
Iakov Davydov's avatar
Iakov Davydov committed
165
166
													std::cout << "Num. threads:	  1 serial" << std::endl
															  << "Num. cores:	  1"  << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
		}
#ifdef USE_MPI
		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;
#endif
													std::cout << "Compiled with:  ";
#ifdef _OPENMP
													std::cout << "USE_OPENMP ";
#endif
#ifdef USE_MPI
													std::cout << "USE_MPI ";
#endif
#ifdef USE_CPV_SCALING
													std::cout << "USE_CPV_SCALING ";
#endif
#ifdef NEW_LIKELIHOOD
													std::cout << "NEW_LIKELIHOOD ";
#endif
#ifdef NON_RECURSIVE_VISIT
													std::cout << "NON_RECURSIVE_VISIT ";
#endif
#ifdef USE_DAG
													std::cout << "USE_DAG ";
#endif
#ifdef USE_ORIGINAL_PROPORTIONS
													std::cout << "USE_ORIGINAL_PROPORTIONS ";
#endif
#ifdef USE_LAPACK
													std::cout << "USE_LAPACK ";
#endif
#ifdef USE_MKL_VML
													std::cout << "USE_MKL_VML";
#endif
													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)
#ifdef USE_MPI
	// 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;
#else
	if(cmd.mSeed == 0) cmd.mSeed = static_cast<unsigned int>(time(NULL));
#endif
	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;
Iakov Davydov's avatar
Iakov Davydov committed
247
		int zero_on_int_cnt	 = 0;
Iakov Davydov's avatar
Iakov Davydov committed
248
249
250
251
252
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
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
		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();

#ifndef NEW_LIKELIHOOD
		forest.addAggressiveReduction();
#endif
		forest.cleanReductionWorkingData();
#ifdef NEW_LIKELIHOOD
		forest.prepareNewReduction();
#endif
	}
#ifdef NEW_LIKELIHOOD
	else
	{
		forest.prepareNewReductionNoReuse();
	}
#endif

#ifdef NON_RECURSIVE_VISIT
	// Prepare the pointers to visit the trees without recursion
	forest.prepareNonRecursiveVisit();
#endif

	// Subdivide the trees in groups based on dependencies
	//forest.prepareDependencies(cmd.mForceSerial || cmd.mDoNotReduceForest);

#ifdef USE_DAG
	// Load the forest into a DAG
	forest.loadForestIntoDAG(Nt);
#endif

	// 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;}

	// Print few statistics
	if(cmd.mVerboseLevel >= VERBOSE_INFO_OUTPUT) std::cout << forest;

#ifdef USE_MPI
	// 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();
	}
#endif

	// 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
389
				std::cout << " Function calls: " << h0.getNumEvaluations() << "	  ";
Iakov Davydov's avatar
Iakov Davydov committed
390
391
392
				std::cout << std::endl << std::endl;
				if(lnl0 != std::numeric_limits<double>::infinity())
				{
Iakov Davydov's avatar
Iakov Davydov committed
393
394
395
					std::string s0 = h0.printFinalVars(std::cout);
					//std::cout<<"EDW0: "<< s0 <<std::endl;
					output_results.saveParameters(fg_branch, s0, 0);
Iakov Davydov's avatar
Iakov Davydov committed
396
397
398
399
400
401
402
403
404
405
406
407
408
409
				}
				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())
				{
Iakov Davydov's avatar
Iakov Davydov committed
410
411
412
					std::string s1= h1.printFinalVars(std::cout);
					//std::cout<<"EDW1: "<< s1 <<std::endl;
					output_results.saveParameters(fg_branch, s1, 1);
Iakov Davydov's avatar
Iakov Davydov committed
413
414
415
416
417
418
419
420
				}
				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)
Iakov Davydov's avatar
Iakov Davydov committed
421
					std::cout << "LRT: " << std::setprecision(15) << std::fixed << lnl1 - lnl0 << "	 (threshold: " << std::setprecision(15) << std::fixed << THRESHOLD_FOR_LRT << ')';
Iakov Davydov's avatar
Iakov Davydov committed
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
451
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
484
485
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
				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;
}

/// @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.
/// Only very common and specific names should be all lowercase, like read, clean, size.
///
/// 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 '_'.
/// The first letters specify the kind of constant (like: STS_ for status, OPT_ for option value).
///
/// 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.
///
/// Array sizes and corresponding indexes should be size_t. The remaining counters should be unsigned int.
///
/// 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.

The input `tree_file` is in %Newick format with the file containing only one tree. The `alignment_file` instead is in %Phylip format.

@verbatim

Usage:
Iakov Davydov's avatar
Iakov Davydov committed
564
	FastCodeML [options] tree_file alignment_file
Iakov Davydov's avatar
Iakov Davydov committed
565

Iakov Davydov's avatar
Iakov Davydov committed
566
567
-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
568

Iakov Davydov's avatar
Iakov Davydov committed
569
570
-q	--quiet (no argument)
		No messages except results
Iakov Davydov's avatar
Iakov Davydov committed
571

Iakov Davydov's avatar
Iakov Davydov committed
572
573
-?	-h	--help (no argument)
		This help
Iakov Davydov's avatar
Iakov Davydov committed
574

Iakov Davydov's avatar
Iakov Davydov committed
575
576
-s	--seed (required argument)
		Random number generator seed (0 < seed < 1000000000)
Iakov Davydov's avatar
Iakov Davydov committed
577

Iakov Davydov's avatar
Iakov Davydov committed
578
579
-b	--branch (required argument)
		Do only this branch as foreground branch (count from 0)
Iakov Davydov's avatar
Iakov Davydov committed
580

Iakov Davydov's avatar
Iakov Davydov committed
581
582
-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
583

Iakov Davydov's avatar
Iakov Davydov committed
584
585
-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
586

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

Iakov Davydov's avatar
Iakov Davydov committed
590
591
-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
592

Iakov Davydov's avatar
Iakov Davydov committed
593
594
-nr	 --no-reduce (no argument)
		Do not reduce forest by merging common subtrees
Iakov Davydov's avatar
Iakov Davydov committed
595

Iakov Davydov's avatar
Iakov Davydov committed
596
597
-l	--lengths-from-file	 --times-from-file (no argument)
		Initial branch lengths from tree file
Iakov Davydov's avatar
Iakov Davydov committed
598

Iakov Davydov's avatar
Iakov Davydov committed
599
600
-o	--initial-step (no argument)
		Only the initial step is performed (no maximization)
Iakov Davydov's avatar
Iakov Davydov committed
601

Iakov Davydov's avatar
Iakov Davydov committed
602
603
-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
604

Iakov Davydov's avatar
Iakov Davydov committed
605
606
-r	--trace (no argument)
		Trace the maximization run
Iakov Davydov's avatar
Iakov Davydov committed
607

Iakov Davydov's avatar
Iakov Davydov committed
608
609
-nt	 --number-of-threads (required argument)
		Number of threads (1 for non parallel execution)
Iakov Davydov's avatar
Iakov Davydov committed
610

Iakov Davydov's avatar
Iakov Davydov committed
611
612
-bf	 --branch-from-file (no argument)
		Do only the branch marked in the file as foreground branch
Iakov Davydov's avatar
Iakov Davydov committed
613

Iakov Davydov's avatar
Iakov Davydov committed
614
615
-hy	 --only-hyp (required argument)
		Compute only H0 if 0, H1 if 1
Iakov Davydov's avatar
Iakov Davydov committed
616

Iakov Davydov's avatar
Iakov Davydov committed
617
618
-i1	 --init-from-h1 (no argument)
		Start H0 optimization from H1 results
Iakov Davydov's avatar
Iakov Davydov committed
619

Iakov Davydov's avatar
Iakov Davydov committed
620
621
-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
622

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

Iakov Davydov's avatar
Iakov Davydov committed
626
627
-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
628

Iakov Davydov's avatar
Iakov Davydov committed
629
630
-ic	 --init-default (no argument)
		Start from default parameter values and times from tree file
Iakov Davydov's avatar
Iakov Davydov committed
631

Iakov Davydov's avatar
Iakov Davydov committed
632
633
-x	--extra-debug (required argument)
		Extra debug parameter (zero disables it)
Iakov Davydov's avatar
Iakov Davydov committed
634

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

Iakov Davydov's avatar
Iakov Davydov committed
638
639
-ou	 --output (required argument)
		Write results formatted to this file
Iakov Davydov's avatar
Iakov Davydov committed
640

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

Iakov Davydov's avatar
Iakov Davydov committed
644
645
-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
646

Iakov Davydov's avatar
Iakov Davydov committed
647
648
-mi	 --max-iterations (required argument)
		Maximum number of iterations for the maximizer (default: 10000)
Iakov Davydov's avatar
Iakov Davydov committed
649

Iakov Davydov's avatar
Iakov Davydov committed
650
651
-bl	 --branch-lengths-fixed (no argument)
		The length of the brances is fixed
Iakov Davydov's avatar
Iakov Davydov committed
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671

@endverbatim
*/

/// @page vampir_page Using Vampir for profiling
/// On Linux we use VampirTrace to collect profile data and Vampir to display the results (http://www.vampir.eu/).
///
/// 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
/// If you want to trace also the OpenMP calls then change it to: -vt:mt -vt:preprocess -DVTRACE
///
/// Then proceed as usual to build the executable.
///
/// Before running the executable, define the following environment variables:
///
Iakov Davydov's avatar
Iakov Davydov committed
672
673
674
675
676
///		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
677
678
679
680
681
///
/// Due to a VampirTrace bug, at the end of the execution, run the vtunify executable by itself.
///
/// Now you can analyze the results by running vampir on the *.otf file generated.
///