fast.cpp 23.8 KB
Newer Older
1
2
3
/// @mainpage FastCodeML
///
/// @section intro_sect Introduction
4
///
5
/// FastCodeML is a rewrite of CodeML based directly on the pseudocode document.
mvalle's avatar
mvalle committed
6
7
/// 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.
8
9
10
11
///
/// @section contacts_sect Contacts
///
/// Contact us if you want more information on the project, want to collaborate or suggest new ideas.
12
///
13
14
///- 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
15
16
17
18
///

#include <iostream>
#include <iomanip>
19
#include <limits>
20
21
22
23
24
25
26
#include "CmdLine.h"
#include "Newick.h"
#include "Phylip.h"
#include "BayesTest.h"
#include "Forest.h"
#include "Exceptions.h"
#include "BranchSiteModel.h"
27
#include "ParseParameters.h"
28
#include "VerbosityLevels.h"
29
#include "WriteResults.h"
30

mvalle's avatar
mvalle committed
31
#ifndef VTRACE
32
33
34
#ifdef _OPENMP
#include <omp.h>
#endif
mvalle's avatar
mvalle committed
35
#endif
36
#ifdef USE_MKL_VML
mvalle's avatar
mvalle committed
37
#include <mkl_vml.h>
38
39
40
41
42
43
44
45
46
47
#endif
#include "Timer.h"
#ifdef USE_MPI
#include "HighLevelCoordinator.h"
#endif

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

55

56
int main(int aRgc, char **aRgv)
57
58
{
	try
59
	{
60

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

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

	// Parse the command line
	CmdLine cmd;
73
	cmd.parseCmdLine(aRgc, aRgv);
74
75
76

	// Adjust and report the number of threads that will be used
#ifdef _OPENMP
77
78
79
80
81
82
83
84
85
	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)
86
        cmd.mForceSerial = true;
87
88
89
     else
        cmd.mForceSerial = false;*/
#else
90
	cmd.mNumThreads=1;
91
	int num_threads = 1;
92
	cmd.mForceSerial = true;
93
94
95
96
97
98
99
100
101
102
103
104
105
106
#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*/
107
108
109

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

//    std::cout <<std::endl<<"------------------"<< std::endl<<"FastCodeML V"<<version<<std::endl<<"------------------"<<std::endl;
114
	// Write out command line parameters (if not quiet i.e. if verbose level > 0)
115
	if(cmd.mVerboseLevel >= VERBOSE_INFO_OUTPUT)
116
	{
117
118
119
120
121
	
		std::cout << "------------------------------------" << std::endl;
		std::cout << "FastCodeML V"<<version << std::endl;
		std::cout << "------------------------------------"<<std::endl;
		std::cout << std::endl;
122
123
124
125
126
													std::cout << "Tree file:      " << cmd.mTreeFile << std::endl;
													std::cout << "Gene file:      " << cmd.mGeneFile << std::endl;
													std::cout << "Verbose level:  " << cmd.mVerboseLevel << " (" << decodeVerboseLevel(cmd.mVerboseLevel) << ')' << std::endl;
		if(cmd.mSeed)								std::cout << "Seed:           " << cmd.mSeed << std::endl;
		if(cmd.mBranchFromFile)						std::cout << "Branch:         From tree file" << std::endl;
mvalle's avatar
mvalle committed
127
128
129
130
		else if(cmd.mBranchStart != UINT_MAX && cmd.mBranchStart == cmd.mBranchEnd)
			                                        std::cout << "Branch:         " << cmd.mBranchStart << std::endl;
		else if(cmd.mBranchStart != UINT_MAX && cmd.mBranchEnd == UINT_MAX)
			                                        std::cout << "Branches:       " << cmd.mBranchStart << "-end" << std::endl;
131
132
		else if(cmd.mBranchStart != UINT_MAX && cmd.mBranchEnd != UINT_MAX)
													std::cout << "Branches:       " << cmd.mBranchStart << '-' << cmd.mBranchEnd << std::endl;
133
		if(!cmd.mStopIfNotLRT)						std::cout << "H0 pre stop:    No" << std::endl;
134
135
136
137
		if(cmd.mIgnoreFreq)							std::cout << "Codon freq.:    Ignore" << std::endl;
		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;
138
139
140
		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;
141
		else if(cmd.mBranchLengthsFromFile)			std::cout << "Starting val.:  Times from tree file" << std::endl;
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;
mvalle's avatar
mvalle committed
147
		if(cmd.mGraphFile && cmd.mExportComputedTimes != UINT_MAX)
148
149
													std::cout << "Graph times:    From H" << cmd.mExportComputedTimes << std::endl;
		if(!cmd.mNoMaximization)					std::cout << "Optimizer:      " << cmd.mOptimizationAlgo << std::endl;
150
		if(cmd.mMaxIterations != MAX_ITERATIONS)	std::cout << "Max iterations: " << cmd.mMaxIterations << std::endl;
151
152
153
		if(cmd.mDeltaValueForGradient > 0.0)		std::cout << "Delta value:    " << cmd.mDeltaValueForGradient << std::endl;
													std::cout << "Relative error: " << cmd.mRelativeError << std::endl;
		if(cmd.mResultsFile)						std::cout << "Results file:   " << cmd.mResultsFile << std::endl;
154
        if(cmd.mNumThreads)                         std::cout << "Number of threads: " << cmd.mNumThreads << std::endl;
155
        if(cmd.mFixedBranchLength)                        std::cerr << "Branch lengths are fixed" << std::endl;
156
157
158
#ifdef _OPENMP
		if(num_threads > 1)
		{
159
													std::cout << "Num. threads:   " << num_threads << std::endl
mvalle's avatar
mvalle committed
160
		                                                      << "Num. cores:     " << omp_get_num_procs() << std::endl;
161
162
		}
		else
mvalle's avatar
mvalle committed
163
#endif
164
		{
165
													std::cout << "Num. threads:   1 serial" << std::endl
mvalle's avatar
mvalle committed
166
		                                                      << "Num. cores:     1"  << std::endl;
167
168
		}
#ifdef USE_MPI
169
170
		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;
171
#endif
172
													std::cout << "Compiled with:  ";
173
#ifdef _OPENMP
174
													std::cout << "USE_OPENMP ";
175
176
#endif
#ifdef USE_MPI
177
													std::cout << "USE_MPI ";
178
#endif
mvalle's avatar
mvalle committed
179
180
181
#ifdef USE_CPV_SCALING
													std::cout << "USE_CPV_SCALING ";
#endif
182
#ifdef NEW_LIKELIHOOD
183
													std::cout << "NEW_LIKELIHOOD ";
184
185
#endif
#ifdef NON_RECURSIVE_VISIT
186
													std::cout << "NON_RECURSIVE_VISIT ";
187
#endif
mvalle's avatar
mvalle committed
188
#ifdef USE_DAG
189
													std::cout << "USE_DAG ";
mvalle's avatar
mvalle committed
190
#endif
191
#ifdef USE_ORIGINAL_PROPORTIONS
192
													std::cout << "USE_ORIGINAL_PROPORTIONS ";
193
194
#endif
#ifdef USE_LAPACK
195
													std::cout << "USE_LAPACK ";
196
197
#endif
#ifdef USE_MKL_VML
198
													std::cout << "USE_MKL_VML";
199
#endif
200
													std::cout << std::endl << std::endl;
201
													if(cmd.mInitFromParams)
202
													{
203
														std::cout << "Param initial values:" << std::endl << std::endl
204
																  << ParseParameters::getInstance();
205
													}
206
207
208
	}

	// Initialize the random number generator (0 means it is not set on the command line)
209
210
211
212
#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
213
	if(cmd.mSeed == 0) cmd.mSeed = static_cast<unsigned int>(time(NULL));
214
#endif
215
216
	srand(cmd.mSeed);

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

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

224
	// Create the forest
225
	Forest forest(cmd.mVerboseLevel);
226
227
228

	// Enclose file loading into a block so temporary structures could be deleted when no more needed
	{
mvalle's avatar
mvalle committed
229
230
	// Load the multiple sequence alignment (MSA)
	Phylip msa(cmd.mVerboseLevel);
mvalle's avatar
mvalle committed
231
	msa.readFile(cmd.mGeneFile, cmd.mCleanData);
232
233
234

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

	// Check coherence between the two files
mvalle's avatar
mvalle committed
238
	msa.checkNameCoherence(tree.getSpecies());
239

240
	// Check root
241
242
	tree.checkRootBranches();

mvalle's avatar
mvalle committed
243
	// If times from file then check for null branch lengths for any leaf
244
245
246
247
248
249
250
251
	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)
		{
252
			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;
253
254
255
256
257
258
259
		}

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

	// Print the tree with the numbering of internal branches
	if(cmd.mVerboseLevel >= VERBOSE_INFO_OUTPUT) tree.printTreeAnnotated(std::cout);
mvalle's avatar
mvalle committed
263

264
	// Load the forest
mvalle's avatar
mvalle committed
265
	forest.loadTreeAndGenes(tree, msa, cmd.mIgnoreFreq ? CodonFrequencies::CODON_FREQ_MODEL_UNIF : CodonFrequencies::CODON_FREQ_MODEL_F3X4);
266
	}
267
268
269
270

	// Reduce the forest merging common subtrees. Add also more reduction, then clean the no more useful data.
	if(!cmd.mDoNotReduceForest)
	{
271
272
		//bool sts = forest.reduceSubtrees(cmd.mNumReductionBlocks);
		forest.reduceSubtrees();
273
274

#ifndef NEW_LIKELIHOOD
275
		forest.addAggressiveReduction();
276
#endif
277
		forest.cleanReductionWorkingData();
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
#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
295
	//forest.prepareDependencies(cmd.mForceSerial || cmd.mDoNotReduceForest);
296

mvalle's avatar
mvalle committed
297
298
299
300
301
#ifdef USE_DAG
	// Load the forest into a DAG
	forest.loadForestIntoDAG(Nt);
#endif

302
	// Get the time needed by data preprocessing
303
	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;}
304

mvalle's avatar
mvalle committed
305
	// Print few statistics
306
	if(cmd.mVerboseLevel >= VERBOSE_INFO_OUTPUT) std::cout << forest;
mvalle's avatar
mvalle committed
307

308
309
#ifdef USE_MPI
	// Distribute the work. If run under MPI then finish, else return to the standard execution flow
310
	if(cmd.mVerboseLevel >= VERBOSE_INFO_OUTPUT) timer.start();
mvalle's avatar
mvalle committed
311
	bool has_run_under_MPI = hlc.startWork(forest, cmd);
312
313

	// If executed under MPI report the time spent, otherwise stop the timer so it can be restarted around the serial execution
mvalle's avatar
mvalle committed
314
	if(has_run_under_MPI)
315
	{
316
		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;}
317
318
319
320
321
322
323
324
		return 0;
	}
	else
	{
		timer.stop();
	}
#endif

325
326
	// Initialize the output results file (if the argument is null, no file is created)
	WriteResults output_results(cmd.mResultsFile);
327

mvalle's avatar
mvalle committed
328
	// Compute the range of branches to mark as foreground
329
	size_t branch_start, branch_end;
mvalle's avatar
mvalle committed
330
	forest.getBranchRange(cmd, branch_start, branch_end);
331
332

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

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

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

342
	// For all requested internal branches
mvalle's avatar
mvalle committed
343
	for(size_t fg_branch=branch_start; fg_branch <= branch_end; ++fg_branch)
344
	{
345
		if(cmd.mVerboseLevel >= VERBOSE_ONLY_RESULTS) std::cout << std::endl << "Doing branch " << fg_branch << std::endl;
346
347

		// Compute the alternate model maximum loglikelihood
mvalle's avatar
mvalle committed
348
		double lnl1 = 0.;
349
350
		if(cmd.mComputeHypothesis != 0)
		{
351
352
			if(cmd.mInitFromParams)			h1.initFromParams();
			if(cmd.mBranchLengthsFromFile)	h1.initFromTree();
353
354

			lnl1 = h1(fg_branch);
355
356
357

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

mvalle's avatar
mvalle committed
360
361
362
363
		// Compute the null model maximum loglikelihood
		double lnl0 = 0.;
		if(cmd.mComputeHypothesis != 1)
		{
364
			if(cmd.mInitH0fromH1)				h0.initFromResult(h1.getVariables());
365
366
367
368
369
			else
			{
				if(cmd.mInitFromParams)			h0.initFromParams();
				if(cmd.mBranchLengthsFromFile)	h0.initFromTree();
			}
mvalle's avatar
mvalle committed
370

371
			lnl0 = h0(fg_branch, cmd.mStopIfNotLRT && cmd.mComputeHypothesis != 0, lnl1-THRESHOLD_FOR_LRT);
mvalle's avatar
mvalle committed
372
373
374
375
376

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

377
		if(cmd.mVerboseLevel >= VERBOSE_ONLY_RESULTS)
378
		{
379
			std::cout << std::endl;
380
381
			if(cmd.mComputeHypothesis != 1)
			{
382
				std::cout << "LnL0: ";
383
384
385
				if(lnl0 == std::numeric_limits<double>::infinity())
					std::cout << "**Invalid result**";
				else if(lnl0 < DBL_MAX)
386
					std::cout << std::setprecision(15) << std::fixed << lnl0;
mvalle's avatar
mvalle committed
387
				else
388
389
390
					std::cout << "(Doesn't pass LRT, skipping)";
				std::cout << " Function calls: " << h0.getNumEvaluations() << "   ";
				std::cout << std::endl << std::endl;
391
392
393
394
395
				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);
396
				}
397
				std::cout << std::endl;
398
399
400
			}
			if(cmd.mComputeHypothesis != 0)
			{
401
				std::cout << "LnL1: ";
402
403
404
405
				if(lnl1 == std::numeric_limits<double>::infinity())
					std::cout << "**Invalid result**";
				else
					std::cout << std::setprecision(15) << std::fixed << lnl1;
406
407
				std::cout << " Function calls: " << h1.getNumEvaluations();
				std::cout << std::endl << std::endl;
408
409
410
411
412
				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);
413
				}
414
				std::cout << std::endl;
415
			}
mvalle's avatar
mvalle committed
416
417
			if(cmd.mComputeHypothesis > 1)
			{
418
419
420
				if(lnl0 == std::numeric_limits<double>::infinity() || lnl1 == std::numeric_limits<double>::infinity())
					std::cout << "LRT: **Invalid result**";
				else if(lnl0 < DBL_MAX)
mvalle's avatar
mvalle committed
421
422
423
424
425
					std::cout << "LRT: " << std::setprecision(15) << std::fixed << lnl1 - lnl0 << "  (threshold: " << std::setprecision(15) << std::fixed << THRESHOLD_FOR_LRT << ')';
				else
					std::cout << "LRT: < " << std::setprecision(15) << std::fixed << THRESHOLD_FOR_LRT;
				std::cout << std::endl;
			}
426
427
428
429
430
431
432
433
434
435
436
437
438
439
		}

		// 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;
mvalle's avatar
mvalle committed
440
441
442

			default:
				break;
443
444
445
446
447
448
449
			}

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

450
		// If the two hypothesis are computed, H0 has not been stopped and the run passes the LRT, then compute the BEB
mvalle's avatar
mvalle committed
451
		if(cmd.mComputeHypothesis > 1 && lnl0 < DBL_MAX && BranchSiteModel::performLRT(lnl0, lnl1))
452
		{
mvalle's avatar
mvalle committed
453
454
455
456
457
			// Get the scale values from the latest optimized h1.
			std::vector<double> scales(2);
			h1.getScales(scales);

			// Run the BEB test
458
			beb.computeBEB(h1.getVariables(), fg_branch, scales);
mvalle's avatar
mvalle committed
459
460

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

463
464
465
466
467
468
469
470
			// 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);
			}
471
472
473
474
		}
	}

	// Get the time needed by the parallel part
475
	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;}
476

477
478
479
	// Output the results
	output_results.outputResults();

480
481
482
	////////////////////////////////////////////////////////////////////
	// Catch all exceptions
	}
mvalle's avatar
mvalle committed
483
	catch(const FastCodeMLSuccess&)
484
485
486
	{
		return 0;
	}
mvalle's avatar
mvalle committed
487
	catch(const FastCodeMLFatal& e)
488
	{
489
490
		// 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;
491
492
		return 1;
	}
mvalle's avatar
mvalle committed
493
494
	catch(const FastCodeMLMemoryError& e)
	{
495
		std::cout << std::endl << e.what() << std::endl << std::endl;
mvalle's avatar
mvalle committed
496
497
498
499
		return 1;
	}
	catch(const std::bad_alloc& e)
	{
500
		std::cout << std::endl << e.what() << std::endl << std::endl;
mvalle's avatar
mvalle committed
501
502
		return 1;
	}
503
504
	catch(...)
	{
505
		std::cout << std::endl << "Default exception caught." << std::endl << std::endl;
506
507
508
509
510
511
512
513
514
515
516
517
		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.
///
518
/// Ex: %PhyloTree
519
///
520
/// @section cmeth_sect Class methods
521
/// Class methods names are CamelCase with the first letter lowercase.
522
/// Only very common and specific names should be all lowercase, like read, clean, size.
523
524
525
///
/// Ex: testFillQ
///
526
/// @section cdatamemb_sect Class data members
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
/// 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.
///
550
/// Array sizes and corresponding indexes should be size_t. The remaining counters should be unsigned int.
551
///
mvalle's avatar
mvalle committed
552
/// The null pointer should be written as NULL, not 0 to make clear its purpose.
553
///
mvalle's avatar
mvalle committed
554
555
556
557
558

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

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

mvalle's avatar
mvalle committed
561
@verbatim
562

563
Usage:
mvalle's avatar
mvalle committed
564
    FastCodeML [options] tree_file alignment_file
565

566
-d  --debug  -v  --verbose (required argument)
567
        Verbosity level (0: none; 1: results only; 2: normal info; 3: MPI trace; 4: more debug) (default: 1)
568

569
570
-q  --quiet (no argument)
        No messages except results
571

572
573
-?  -h  --help (no argument)
        This help
574

575
576
-s  --seed (required argument)
        Random number generator seed (0 < seed < 1000000000)
577

578
-b  --branch (required argument)
mvalle's avatar
mvalle committed
579
580
581
582
583
584
585
        Do only this branch as foreground branch (count from 0)

-bs  --branch-start (required argument)
        Start computing from this branch as foreground one (count from 0) (default: first one)

-be  --branch-end (required argument)
        End computing at this branch as foreground one (count from 0) (default: last one)
586

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

590
591
-e  --export (required argument)
        Export forest in GML format (if %03d or @03d is present, one is created for each fg branch)
592

593
594
-nr  --no-reduce (no argument)
        Do not reduce forest by merging common subtrees
595

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

599
600
-o  --initial-step (no argument)
        Only the initial step is performed (no maximization)
601

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
604

605
606
-r  --trace (no argument)
        Trace the maximization run
607

608
609
-nt  --number-of-threads (required argument)
        Number of threads (1 for non parallel execution)
610

611
612
-bf  --branch-from-file (no argument)
        Do only the branch marked in the file as foreground branch
613

614
615
-hy  --only-hyp (required argument)
        Compute only H0 if 0, H1 if 1
616

617
618
-i1  --init-from-h1 (no argument)
        Start H0 optimization from H1 results
619

620
-m  --maximizer (required argument)
621
        Optimizer algorithm (0:LBFGS, 1:VAR1, 2:VAR2, 3:SLSQP, 11:BOBYQA, 22:FromCodeML, 99:MLSL_LDS) (default: 22)
622

623
-sd  --small-diff (required argument)
624
        Delta used in gradient computation (default: 1.49e-8)
625

626
627
-p  --init-param (required argument)
        Pass initialization parameter in the form: P=value (P: w0, k, p0, p1, w2)
628

629
630
-ic  --init-default (no argument)
        Start from default parameter values and times from tree file
631

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

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

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

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

644
-ps  --no-pre-stop (no argument)
mvalle's avatar
mvalle committed
645
        Don't stop H0 maximization even if it cannot satisfy LRT (default: stop)
646

647
-mi  --max-iterations (required argument)
648
649
650
        Maximum number of iterations for the maximizer (default: 10000)

-bl  --branch-lengths-fixed (no argument)
651
        The length of the brances is fixed
652

653
654
@endverbatim
*/
mvalle's avatar
mvalle committed
655

mvalle's avatar
mvalle committed
656
/// @page vampir_page Using Vampir for profiling
657
/// On Linux we use VampirTrace to collect profile data and Vampir to display the results (http://www.vampir.eu/).
mvalle's avatar
mvalle committed
658
659
660
661
662
663
664
///
/// Before invoking CMAKE define CXX=vtCC
///
/// Define CMAKE_BUILD_TYPE as: RelWithDebInfo
///
/// Run CMAKE and configure.
///
665
/// Then define CMAKE_CXX_FLAGS as: -vt:mt -vt:noopari
666
/// If you want to trace also the OpenMP calls then change it to: -vt:mt -vt:preprocess -DVTRACE
mvalle's avatar
mvalle committed
667
668
669
///
/// Then proceed as usual to build the executable.
///
670
/// Before running the executable, define the following environment variables:
mvalle's avatar
mvalle committed
671
///
672
///     export VT_BUFFER_SIZE=512M
mvalle's avatar
mvalle committed
673
///     export VT_MAX_FLUSHES=0
674
675
///     export VT_SYNCH_FLUSH=yes
///     export VT_GPUTRACE=no
mvalle's avatar
mvalle committed
676
677
///     export VT_UNIFY=no
///
678
/// Due to a VampirTrace bug, at the end of the execution, run the vtunify executable by itself.
mvalle's avatar
mvalle committed
679
680
681
///
/// Now you can analyze the results by running vampir on the *.otf file generated.
///