fast.cpp 30.8 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
///- Ing. <a href="mailto:mvalle@cscs.ch">Mario Valle</a> - Swiss National
17
///Supercomputing Centre (CSCS) - Switzerland
Iakov Davydov's avatar
Iakov Davydov committed
18
///- The HP2C <a href="mailto:selectome@hp2c.ch">Selectome</a> Project Group -
19
///Mainly based in University of Lausanne - Switzerland
Iakov Davydov's avatar
Iakov Davydov committed
20
21
22
23
24
///

#include <iostream>
#include <iomanip>
#include <limits>
omid's avatar
omid committed
25
#include <string>
Iakov Davydov's avatar
Iakov Davydov committed
26
27
28
29
30
31
32
33
34
35
#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"
omid's avatar
omid committed
36
#include "MathSupport.h"
Iakov Davydov's avatar
Iakov Davydov committed
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

#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.
///
53
54
55
56
///		@author Mario Valle - Swiss National Supercomputing Centre
///(CSCS)
///		@date 2010-12-22 (initial version)
///		@version 1.1
Iakov Davydov's avatar
Iakov Davydov committed
57
58
59
60
///
///	@param[in] aRgc Number of command line parameters
/// @param[in] aRgv Command line parameters
///
omid's avatar
omid committed
61
const char *version = "1.3.0";
Iakov Davydov's avatar
Iakov Davydov committed
62

Iakov Davydov's avatar
Iakov Davydov committed
63
int main(int aRgc, char **aRgv) {
omid's avatar
omid committed
64

65
	Timer timer_app;
omid's avatar
omid committed
66

67
	timer_app.start();
omid's avatar
omid committed
68

69
	try {
Iakov Davydov's avatar
Iakov Davydov committed
70
71

#ifdef USE_MKL_VML
72
73
		// If used, intitialize the MKL VML library
		vmlSetMode(VML_HA | VML_DOUBLE_CONSISTENT);
Iakov Davydov's avatar
Iakov Davydov committed
74
75
76
#endif

#ifdef USE_MPI
77
78
		// Start the high level parallel executor (based on MPI)
		HighLevelCoordinator hlc(&aRgc, &aRgv);
Iakov Davydov's avatar
Iakov Davydov committed
79
80
#endif

81
82
83
		// Parse the command line
		CmdLine cmd;
		cmd.parseCmdLine(aRgc, aRgv);
Iakov Davydov's avatar
Iakov Davydov committed
84

Iakov Davydov's avatar
Iakov Davydov committed
85
// Adjust and report the number of threads that will be used
Iakov Davydov's avatar
Iakov Davydov committed
86
#ifdef _OPENMP
87
		int num_threads = omp_get_max_threads();
omid's avatar
omid committed
88
		//std::cout<<"max num of thr: "<< num_threads <<std::endl;
89
90
91
92

		if ((cmd.mNumThreads >= 1) &&
				(cmd.mNumThreads <= (unsigned int)num_threads))
		num_threads = cmd.mNumThreads;
omid's avatar
omid committed
93
94
		//else
		//num_threads = num_threads/2 + 1;// to prevent possible false sharing when it is set to maximum
95
96
97
		// std::cout<<"num of thr: "<< num_threads <<std::endl;

		omp_set_num_threads(num_threads);
omid's avatar
omid committed
98
		//std::cout<<"current num of thr: "<< num_threads <<", command line num of thr: "<<cmd.mNumThreads<<std::endl;
99
100
101
102
		/*if (num_threads < 2)
		 cmd.mForceSerial = true;
		 else
		 cmd.mForceSerial = false;*/
Iakov Davydov's avatar
Iakov Davydov committed
103
#else
104
105
106
		cmd.mNumThreads = 1;
		int num_threads = 1;
		cmd.mForceSerial = true;
Iakov Davydov's avatar
Iakov Davydov committed
107
108
#endif

109
110
111
112
113
114
115
116
117
118
119
120
		/*#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*/
Iakov Davydov's avatar
Iakov Davydov committed
121
122

#ifdef USE_MPI
123
124
125
		// Shutdown messages from all MPI processes except the master
		if (!hlc.isMaster())
		cmd.mVerboseLevel = VERBOSE_NONE;
Iakov Davydov's avatar
Iakov Davydov committed
126
127
#endif

128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
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
		//	  std::cout <<std::endl<<"------------------"<< std::endl<<"FastCodeML
		//V"<<version<<std::endl<<"------------------"<<std::endl;
		// 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;
			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;
			else if (cmd.mBranchAll)
				std::cout << "FG Branches:	  All (internals + leaves) "
						<< std::endl;
			// 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;
			// else if(cmd.mBranchStart != UINT_MAX && cmd.mBranchEnd != UINT_MAX)
			//											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;
			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)
				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;
			if (cmd.mGraphFile && cmd.mExportComputedTimes != UINT_MAX)
				std::cout << "Graph times:	  From H" << cmd.mExportComputedTimes
						<< std::endl;
			if (!cmd.mNoMaximization)
				std::cout << "Optimizer:	  " << cmd.mOptimizationAlgo
						<< std::endl;
			if (cmd.mMaxIterations != MAX_ITERATIONS)
				std::cout << "Max iterations: " << cmd.mMaxIterations
						<< std::endl;
			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;
			if (cmd.mNumThreads)
				std::cout << "Number of threads: " << cmd.mNumThreads
						<< std::endl;
			if (cmd.mFixedBranchLength)
				std::cout << "Branch lengths are fixed" << std::endl;
Iakov Davydov's avatar
Iakov Davydov committed
212
#ifdef _OPENMP
213
			if (num_threads > 1) {
214
215
				std::cout << "Current num. threads:	  " << num_threads << std::endl
				<< "Total num. cores:	  " << omp_get_num_procs() << std::endl;
216
			} else
Iakov Davydov's avatar
Iakov Davydov committed
217
#endif
218
219
220
221
			{
				std::cout << "Num. threads:	  1 serial" << std::endl
						<< "Num. cores:	  1" << std::endl;
			}
Iakov Davydov's avatar
Iakov Davydov committed
222
#ifdef USE_MPI
223
224
225
226
227
228
			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
229
#endif
230
			std::cout << "Compiled with:  ";
Iakov Davydov's avatar
Iakov Davydov committed
231
#ifdef _OPENMP
232
			std::cout << "USE_OPENMP ";
Iakov Davydov's avatar
Iakov Davydov committed
233
234
#endif
#ifdef USE_MPI
235
			std::cout << "USE_MPI ";
Iakov Davydov's avatar
Iakov Davydov committed
236
237
#endif
#ifdef USE_CPV_SCALING
238
			std::cout << "USE_CPV_SCALING ";
Iakov Davydov's avatar
Iakov Davydov committed
239
240
#endif
#ifdef NEW_LIKELIHOOD
241
			std::cout << "NEW_LIKELIHOOD ";
Iakov Davydov's avatar
Iakov Davydov committed
242
243
#endif
#ifdef NON_RECURSIVE_VISIT
244
			std::cout << "NON_RECURSIVE_VISIT ";
Iakov Davydov's avatar
Iakov Davydov committed
245
246
#endif
#ifdef USE_DAG
247
			std::cout << "USE_DAG ";
Iakov Davydov's avatar
Iakov Davydov committed
248
249
#endif
#ifdef USE_ORIGINAL_PROPORTIONS
250
			std::cout << "USE_ORIGINAL_PROPORTIONS ";
Iakov Davydov's avatar
Iakov Davydov committed
251
252
#endif
#ifdef USE_LAPACK
253
			std::cout << "USE_LAPACK ";
Iakov Davydov's avatar
Iakov Davydov committed
254
255
#endif
#ifdef USE_MKL_VML
256
			std::cout << "USE_MKL_VML";
Iakov Davydov's avatar
Iakov Davydov committed
257
#endif
258
259
260
261
262
263
			std::cout << std::endl << std::endl;
			if (cmd.mInitFromParams) {
				std::cout << "Param initial values:" << std::endl << std::endl
						<< ParseParameters::getInstance();
			}
		}
Iakov Davydov's avatar
Iakov Davydov committed
264
265
266

// Initialize the random number generator (0 means it is not set on the command
// line)
Iakov Davydov's avatar
Iakov Davydov committed
267
#ifdef USE_MPI
268
269
270
271
		// 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
272
#else
273
274
		if (cmd.mSeed == 0)
			cmd.mSeed = static_cast<unsigned int>(time(NULL));
Iakov Davydov's avatar
Iakov Davydov committed
275
#endif
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
		// srand(cmd.mSeed); // fastcodeml seed
		SetSeedCodeml(cmd.mSeed, 0); // codeml seed is 1

		// 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 and unrooting if tree is rooted
		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) {
				if (cmd.mVerboseLevel >= VERBOSE_ONLY_RESULTS) {
					std::cout
							<< "Found null or missing branch length in tree file: on "
							<< zero_on_leaf_cnt << " leave(s) and on "
							<< zero_on_int_cnt << " internal branch(es)."
							<< std::endl;
				}
			}
		}

		// 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
340
341

#ifndef NEW_LIKELIHOOD
342
			forest.addAggressiveReduction();
Iakov Davydov's avatar
Iakov Davydov committed
343
#endif
344
			forest.cleanReductionWorkingData();
Iakov Davydov's avatar
Iakov Davydov committed
345
#ifdef NEW_LIKELIHOOD
346
			forest.prepareNewReduction();
Iakov Davydov's avatar
Iakov Davydov committed
347
#endif
348
		}
Iakov Davydov's avatar
Iakov Davydov committed
349
#ifdef NEW_LIKELIHOOD
350
351
352
		else {
			forest.prepareNewReductionNoReuse();
		}
Iakov Davydov's avatar
Iakov Davydov committed
353
354
355
#endif

#ifdef NON_RECURSIVE_VISIT
356
357
		// Prepare the pointers to visit the trees without recursion
		forest.prepareNonRecursiveVisit();
Iakov Davydov's avatar
Iakov Davydov committed
358
359
#endif

Iakov Davydov's avatar
Iakov Davydov committed
360
361
// Subdivide the trees in groups based on dependencies
// forest.prepareDependencies(cmd.mForceSerial || cmd.mDoNotReduceForest);
Iakov Davydov's avatar
Iakov Davydov committed
362
363

#ifdef USE_DAG
364
365
		// Load the forest into a DAG
		forest.loadForestIntoDAG(Nt);
Iakov Davydov's avatar
Iakov Davydov committed
366
367
#endif

368
369
370
371
372
373
374
		// 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
375

376
377
378
		// Print few statistics
		if (cmd.mVerboseLevel >= VERBOSE_INFO_OUTPUT)
			std::cout << forest;
Iakov Davydov's avatar
Iakov Davydov committed
379
380

#ifdef USE_MPI
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
		// 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
401
402
#endif

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
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
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
		// Compute the range of branches to mark as foreground
		size_t branch_start, branch_end;

		std::set<int> fg_set;     // to save a list of fg branches from the
								  // getBranchRange function
		std::set<int> ib_set;    // to save a list of internal branches from the
								 // getBranchRange function
		std::vector<double> mVar; // to save optimization variables

		forest.getBranchRange(cmd, branch_start, branch_end, fg_set, ib_set); // fgset is added to save a list of fg branches

		// for (std::set<int>::iterator it=ib_set.begin(); it!=ib_set.end(); ++it)
		// std::cout << " " << *it << ",";

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

		double lnl0, lnl1 = 0.;

		if (!fg_set.empty()) // in case of marked fg branches (one or multiple fg)
		{
			if (cmd.mVerboseLevel >= VERBOSE_ONLY_RESULTS)
				std::cout << std::endl
						<< "Doing foreground branch(es) from tree file "
						<< std::endl;
			if (cmd.mVerboseLevel >= VERBOSE_ONLY_RESULTS)
				std::cout << "-------------------------------------------"
						<< std::endl;
			if (cmd.mVerboseLevel >= VERBOSE_ONLY_RESULTS) {
				std::cout << "Doing foreground branch(es) ";
				for (std::set<int>::iterator it = fg_set.begin();
						it != fg_set.end(); ++it)
					std::cout << " " << *it << " ";
				std::cout << std::endl;
			}
			// Initialize the models
			MfgBranchSiteModelNullHyp h0(forest, cmd);
			MfgBranchSiteModelAltHyp h1(forest, cmd);

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

			// 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_set);
				// h1.saveComputedTimes();

				// std::cout << "lnl1 = " << lnl1 << std::endl;
				// Save the value for formatted output
				// output_results.saveLnL(fg_set, 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_set,
						cmd.mStopIfNotLRT && cmd.mComputeHypothesis != 0,
						lnl1 - THRESHOLD_FOR_LRT);

				// std::cout << "lnl0 = " << lnl0 << std::endl;
				// 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)";
					std::cout << " Function calls: " << h0.getNumEvaluations()
							<< "	  ";
					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 << "	 (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;
				}
			}

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

			// 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)) {
				if (cmd.mVerboseLevel >= VERBOSE_ONLY_RESULTS)
					std::cout << std::endl
							<< "LRT is significant. Computing sites under positive "
									"selection ... " << std::endl;

				// Get the scale values from the latest optimized h1.
				std::vector<double> scales(2);
				h1.getScales(scales);

				// Run the BEB test
				// if(cmd.mVerboseLevel >= VERBOSE_ONLY_RESULTS) std::cout << std::endl
				// << "LRT is significant. Computing sites under positive selection ...
				// " << std::endl ;
				beb.computeBEB(h1.getVariables(), fg_set, scales);

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

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

				 if (cmd.mVerboseLevel >= VERBOSE_ONLY_RESULTS) {
				 std::cout << std::endl
				 << "Positively selected sites and their probabilities : ";
				 std::cout << std::endl;
				 for (std::vector<unsigned int>::iterator it =
				 positive_sel_sites.begin();
				 it != positive_sel_sites.end(); ++it)
				 std::cout << " " << *it << ",";
				 std::cout << std::endl;
				 for (std::vector<double>::iterator it =
				 positive_sel_sites_prob.begin();
				 it != positive_sel_sites_prob.end(); ++it)
				 std::cout << " " << *it << ",";
				 std::cout << std::endl;
				 }

				 output_results.savePositiveSelSites(fg_set, positive_sel_sites,
				 positive_sel_sites_prob);
				 }*/
			}

			// if branches are fixed

			if (cmd.mFixedBranchLength)

			{
				if (cmd.mVerboseLevel >= VERBOSE_ONLY_RESULTS) {
					std::cout << std::endl << "Final ";
					tree.printTreeAnnotated(std::cout, NULL, 0, true);
omid's avatar
omid committed
620
					std::cout << std::endl;
omid's avatar
omid committed
621
					tree.printTreeAnnotated(std::cout, NULL, 0, true, false);
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
				}
			}

			else

			{
				if (cmd.mVerboseLevel >= VERBOSE_ONLY_RESULTS) {
					if (cmd.mComputeHypothesis == 0)
						mVar = h0.getVariables();
					else
						mVar = h1.getVariables();

					if (cmd.mComputeHypothesis == 0)
						std::cout << std::endl << "H0 Final ";
					else
omid's avatar
omid committed
637
						std::cout << std::endl << "H1 Final ";
638
639
					tree.printTreeAnnotatedWithEstLens(std::cout, NULL, 0, true,
							&mVar);
omid's avatar
omid committed
640
					std::cout << std::endl;
641
642
					tree.printTreeAnnotatedWithEstLens(std::cout, NULL, 0, true,
							&mVar, false);
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
				}
				// tree.printTreeAnnotatedWithEstLens(std::cout, NULL, 0, true,
				// &h1.getVariables());
				// std :: cout << std::endl;
			}

			/*if(cmd.mInitFromParams)			h0.initFromParams();
			 if(cmd.mBranchLengthsFromFile)	h0.initFromTree();

			 lnl0 = h0(fg_set, cmd.mStopIfNotLRT && cmd.mComputeHypothesis != 0,
			 0-THRESHOLD_FOR_LRT);
			 std::cout << "lnl0 (multiple fg) = " << lnl0 << std::endl;

			 if(cmd.mInitFromParams)			h1.initFromParams();
			 if(cmd.mBranchLengthsFromFile)	h1.initFromTree();

			 lnl1 = h1(fg_set);
			 std::cout << "lnl1 (multiple fg) = " << lnl1 << std::endl;*/

			timer_app.stop();
			if (cmd.mVerboseLevel >= VERBOSE_ONLY_RESULTS) {
664
665
				std::cout << std::endl << "Time used: " << timer_app.get() / 60000
									<< "m:" << (timer_app.get() / 1000) % 60 << "s" << std::endl;
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
				std::cout << "Cores used: " << num_threads << std::endl;
			}
			return 0;
		}

		// Else for all requested internal branches

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

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

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

		if (cmd.mVerboseLevel >= VERBOSE_ONLY_RESULTS) {
			if (cmd.mBranchAll)
686
				std::cout << std::endl << "Doing for all branches (internal + leaf) as foreground"
687
688
						<< std::endl;
			else
689
				std::cout << std::endl << "Doing for all internal branches as foreground"
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
						<< std::endl;
			std::cout << "------------------------------------" << std::endl;
		}

		for (size_t fg_branch = branch_start; fg_branch <= branch_end;
				++fg_branch) {

			if (cmd.mBranchAll
					or (!cmd.mBranchAll
							and ib_set.find(fg_branch) != ib_set.end()))

					{

				if (cmd.mVerboseLevel >= VERBOSE_ONLY_RESULTS)
					std::cout << "Doing foreground 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);
					// h1.saveComputedTimes();
					// h1.mBranches

					// 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)";
						std::cout << " Function calls: "
								<< h0.getNumEvaluations() << "	  ";
						std::cout << std::endl << std::endl;
						if (lnl0 != std::numeric_limits<double>::infinity()) {
							std::string s0 = h0.printFinalVars(std::cout);
							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);
							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
									<< "	 (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;
					}
				}

				// 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)) {
					if (cmd.mVerboseLevel >= VERBOSE_ONLY_RESULTS)
						std::cout << std::endl
								<< "LRT is significant. Computing sites under positive "
										"selection ... " << std::endl;

					// Get the scale values from the latest optimized h1.
					std::vector<double> scales(2);
					h1.getScales(scales);

					// Run the BEB test
					// if(cmd.mVerboseLevel >= VERBOSE_ONLY_RESULTS) std::cout <<
					// std::endl << "LRT is significant. Computing sites under positive
					// selection ... " << std::endl ;
					// 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);
					}
				}

				if (cmd.mFixedBranchLength)

				{
					if (cmd.mVerboseLevel >= VERBOSE_ONLY_RESULTS) {
						std::cout << std::endl << "Final ";
						tree.printTreeAnnotated(std::cout, NULL, 0, true);
omid's avatar
omid committed
859
						std::cout << std::endl;
omid's avatar
omid committed
860
						tree.printTreeAnnotated(std::cout, NULL, 0, true, false);
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
						std::cout << std::endl;
					}
				}

				else {
					// std :: cout << std::endl;
					if (cmd.mVerboseLevel >= VERBOSE_ONLY_RESULTS) {
						if (cmd.mComputeHypothesis == 0)
							mVar = h0.getVariables();
						else
							mVar = h1.getVariables();

						if (cmd.mComputeHypothesis == 0)
							std::cout << std::endl << "H0 Final ";
						else
omid's avatar
omid committed
876
							std::cout << std::endl << "H1 Final ";
877
878
						tree.printTreeAnnotatedWithEstLens(std::cout, NULL, 0,
								true, &mVar);
omid's avatar
omid committed
879
						std::cout << std::endl;
880
881
						tree.printTreeAnnotatedWithEstLens(std::cout, NULL, 0,
														true, &mVar, false);
omid's avatar
omid committed
882
						std::cout << std::endl;
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
					}
				}
			}
		}

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

		timer_app.stop();
		if (cmd.mVerboseLevel >= VERBOSE_ONLY_RESULTS) {
			std::cout << std::endl << "Time used: " << timer_app.get() / 60000
899
					<< "m:" << (timer_app.get() / 1000) % 60 << "s" << std::endl;
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
			std::cout << "Cores used: " << num_threads << 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
927
928
929
930
931
932
933
934
935
936
937
938
}

/// @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
939
940
/// Only very common and specific names should be all lowercase, like read,
/// clean, size.
Iakov Davydov's avatar
Iakov Davydov committed
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
///
/// 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
956
957
/// The first letters specify the kind of constant (like: STS_ for status, OPT_
/// for option value).
Iakov Davydov's avatar
Iakov Davydov committed
958
959
960
961
962
963
964
965
966
967
968
///
/// 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
969
970
/// Array sizes and corresponding indexes should be size_t. The remaining
/// counters should be unsigned int.
Iakov Davydov's avatar
Iakov Davydov committed
971
972
973
974
///
/// The null pointer should be written as NULL, not 0 to make clear its purpose.
///
/// @page vampir_page Using Vampir for profiling
Iakov Davydov's avatar
Iakov Davydov committed
975
976
/// 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
977
978
979
980
981
982
983
984
///
/// 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
985
986
/// 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
987
988
989
990
991
///
/// Then proceed as usual to build the executable.
///
/// Before running the executable, define the following environment variables:
///
992
993
994
995
996
///		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
997
///
Iakov Davydov's avatar
Iakov Davydov committed
998
999
/// Due to a VampirTrace bug, at the end of the execution, run the vtunify
/// executable by itself.
Iakov Davydov's avatar
Iakov Davydov committed
1000
///