merging master-eval up to CodeMLOptimizer.h

BUGS

-1. Parameter optimizers other than "-m 22" are untested and might not converge properly.
-2. -nt 1 yields results inconsistent with the test case when specified without MPI
+1. Look at the Gitlab page / issues section (https://gitlab.isb-sib.ch/phylo/fastcodeml/issues)

BayesTest.h

@@ -184,4 +184,150 @@ private:
       mBEBset; ///< Probability matrix set to be used for likelihood computation
 };
 
+class MfgBayesTest {
+public:
+  /// Constructor.
+  ///
+  /// @param[in] aForest The forest
+  /// @param[in] aVerbose The verbosity level
+  /// @param[in] aNoReduction If true the dependencies computed are for no
+  /// reduced forests
+  ///
+  explicit MfgBayesTest(Forest &aForest, unsigned int aVerbose = 0,
+                        bool aNoReduction = true)
+      : mForest(aForest), mNumSites(mForest.getNumSites()),
+        mSiteClassProb(BEB_DIMS * mNumSites), mVerbose(aVerbose),
+        mPriors(mNumSites * BEB_NUM_CAT), mDependencies(mForest, aVerbose),
+        mBEBset(mForest.getNumBranches()) {
+    // Create the dependency list for forest likelihood computation
+    mDependencies.computeDependencies(1, aNoReduction);
+    if (mVerbose >= VERBOSE_ONLY_RESULTS)
+      mDependencies.print("TEST FOR BEB (before optimization)");
+    mDependencies.optimizeDependencies();
+    if (mVerbose >= VERBOSE_ONLY_RESULTS)
+      mDependencies.print("TEST FOR BEB");
+  }
+
+  /// Destructor.
+  ///
+  ~MfgBayesTest() {}
+
+  /// Bayes Empirical Bayes (BEB) test.
+  ///
+  /// @param[in] aVars The variables optimized at the end of the H1 run.
+  /// @param[in] aFgBranchSet The foreground branch set under test.
+  /// @param[in] aScales The two scales ([0] bg; [1] fg) to rescale the branch
+  /// lengths. They are computed in H1.
+  ///
+  void computeBEB(const std::vector<double> &aVars, std::set<int> aFgBranchSet,
+                  const std::vector<double> &aScales);
+
+  /// Print the sites under positive selection.
+  ///
+  /// @param[in] aFgBranchSet Set of the branches marked as foreground branch
+  ///
+  void printPositiveSelSites(std::set<int> aFgBranchSet) const;
+
+  /// Extract the sites under positive selection and the corresponding
+  /// probabilities.
+  ///
+  /// @param[out] aPositiveSelSites Vector of sites under positive selection
+  /// @param[out] aPositiveSelSitesProb Corresponding probabilities
+  ///
+  void
+  extractPositiveSelSites(std::vector<unsigned int> &aPositiveSelSites,
+                          std::vector<double> &aPositiveSelSitesProb) const;
+
+private:
+  /// This sets up the grid (mPara[][]) according to the priors.
+  /// It calculates the probability of data at each site given w: f(f_h|w).
+  /// This is calculated using the branch model (NSsites = 0 model = 2), with
+  /// BayesEB=2 used to force the use of the correct scale factors in
+  /// GetPMatBranch().
+  ///
+  ///@verbatim
+  /// Order of site classes for iw or f(x_h|w):
+  ///						back   fore		num.sets
+  /// Branchsite A (121 sets)
+  ///	  site class 0:		 w0		w0		  10
+  ///	  site class 1:		 w1=1	w1=1	   1
+  ///	  site class 2a:	 w0		w2		 100
+  ///	  site class 2b:	 w1=1	w2		  10
+  ///@endverbatim
+  ///
+  /// @param[in] aVars  The variables optimized at the end of the H1 run.
+  /// @param[in] aSiteMultiplicity  The site multiplicity vector.
+  /// @param[in] aFgBranchSet  The foreground branches under test.
+  /// @param[in] aScales	The two scales ([0] bg; [1] fg) to rescale the branch
+  /// lengths. They are computed in H1.
+  ///
+  /// @return The computed scale.
+  ///
+  double getGridParams(const std::vector<double> &aVars,
+                       const std::vector<double> &aSiteMultiplicity,
+                       std::set<int> aFgBranchSet,
+                       const std::vector<double> &aScales);
+
+  /// This gives the indices (ix, iy) and the coordinates (aProbX, aProbY,
+  /// 1-aProbX-aProbY) for
+  /// the aTriangleIdx-th triangle, with aTriangleIdx from 0, 1, ...,
+  /// BEB_N1D*BEB_N1D-1.
+  ///
+  /// The ternary graph (0-1 on each axis) is partitioned into BEB_N1D*BEB_N1D
+  /// equal-sized triangles.
+  /// In the first row (ix=0), there is one triangle (iy=0);
+  /// In the second row (ix=1), there are 3 triangles (iy=0,1,2);
+  /// In the i-th row (ix=i), there are 2*i+1 triangles (iy=0,1,...,2*i).
+  ///
+  /// aProbX rises when ix goes up, but aProbY decreases when iy increases.
+  /// (aProbX, aProbY) is the
+  /// centroid in the ij-th small triangle. aProbX and aProbY each takes on
+  /// 2*BEB_N1D-1 possible values.
+  ///
+  /// @param[out] aProbX The p0 value on the X axis of the triangular grid.
+  /// @param[out] aProbY The p1 value on the Y axis of the triangular grid.
+  /// @param[in] aTriangleIdx The index inside the triangular grid.
+  ///
+  void getIndexTernary(double *aProbX, double *aProbY,
+                       unsigned int aTriangleIdx);
+
+private:
+  /// Disabled assignment operator to avoid warnings on Windows.
+  ///
+  /// @fn BayesTest& operator=(const BayesTest& aObj)
+  ///
+  /// @param[in] aObj The object to be assigned
+  ///
+  /// @return The object receiving the assignment
+  ///
+  MfgBayesTest &operator=(const MfgBayesTest &);
+
+private:
+  const static unsigned int BEB_N1D = 10; ///< Number of intervals for w0 and w2
+  const static unsigned int BEB_DIMS =
+      4; ///< Number of codon classes (0, 1, 2a, 2b)
+  const static unsigned int BEB_NUM_CAT =
+      BEB_N1D + 1 + BEB_N1D * BEB_N1D + BEB_N1D; ///< Total number of categories
+                                                 ///for w0 and w2 (it is
+                                                 ///com.ncatG in codeml.c)
+  const static unsigned int BEB_NGRID =
+      Pow<BEB_N1D, BEB_DIMS>::value; ///< Number of points in the grid used to
+                                     ///evaluate the integral. It is
+                                     ///BEB_N1D^BEB_DIMS
+
+private:
+  Forest &mForest;                    ///< The forest.
+  size_t mNumSites;                   ///< Number of sites.
+  std::vector<double> mSiteClassProb; ///< Probability of a site to pertain to a
+                                      ///given class (one row per class (4
+                                      ///classes), one column per site).
+  unsigned int mVerbose;              ///< If greater than zero prints more info
+  std::vector<double>
+      mPriors; ///< Computed priors (each points to a list, one for each site)
+  TreeAndSetsDependencies
+      mDependencies; ///< Dependency list for likelihood computation
+  mfgProbabilityMatrixSetBEB
+      mBEBset; ///< Probability matrix set to be used for likelihood computation
+};
+
 #endif

BranchSiteModel.h

@@ -117,6 +117,25 @@ public:
   double maximizeLikelihood(size_t aFgBranch, bool aStopIfBigger,
                             double aThreshold);
 
+  /// Compute the maximum likelihood for the given forest (multiple foregrounds)
+  ///
+  /// @param[in] aFgBranch The number of the internal branch to be marked as
+  /// foreground
+  /// @param[in] aStopIfBigger If true stop computation as soon as lnl is over
+  /// aThreshold
+  /// @param[in] aThreshold The threshold at which the maximization should be
+  /// stopped
+  ///
+  /// @return The maximum Likelihood value
+  ///
+  /// @exception FastCodeMLFatal Exception in Ming2 computation
+  /// @exception FastCodeMLFatal Invalid optimization algorithm identifier on
+  /// the command line
+  /// @exception FastCodeMLFatal Exception in computation
+  ///
+  double maximizeLikelihood(std::set<int> aFgBranchSet, bool aStopIfBigger,
+                            double aThreshold);
+
   /// Compute the likelihood for the given forest and the given set of
   /// parameters when computing gradient.
   ///
@@ -169,7 +188,16 @@ public:
   ///
   /// @return The optimized variables
   ///
-  const std::vector<double> &getVariables(void) const { return mVar; }
+  /// const std::vector<double>* getVariables(void) const {return &mVar;} //
+  /// should return refference ?
+
+  /// Get variable values
+  ///
+  /// @return The optimized variables
+  ///
+  const std::vector<double> &getVariables(void) const {
+    return mVar;
+  } // should return refference ?
 
   /// Get the number of function evaluation.
   ///
@@ -214,6 +242,10 @@ public:
     aScales[1] = mFgScale;
   }
 
+  /// Get the array of branch lengths.
+  ///
+  std::vector<double> getBranchLengths(void) const { return mBranches; }
+
   /// Verify if the code entered on command line is valid.
   ///
   /// @param[in] aOptimizationAlgo The code set on the command line.
@@ -242,11 +274,44 @@ protected:
   ///
   void getProportions(double aV0, double aV1, double *aProportions) const {
 #ifdef USE_ORIGINAL_PROPORTIONS
+
+    /*   // codeml code //
+     int f_and_x(double x[], double f[], int n, int fromf, int LastItem)
+     {
+     This transforms between x and f.  x and f can be identical.
+     If (fromf), f->x
+     else        x->f.
+     The iterative variable x[] and frequency f[0,1,n-2] are related as:
+     freq[k] = exp(x[k])/(1+SUM(exp(x[k]))), k=0,1,...,n-2,
+     x[] and freq[] may be the same vector.
+     The last element (f[n-1] or x[n-1]=1) is updated only if(LastItem).
+
+     int i;
+     double tot;
+
+     if (fromf) {   f => x
+     if((tot=1-sum(f,n-1))<1e-80) error2("f[n-1]==1, not dealt with.");
+     tot = 1/tot;
+     for(i=0; i<n-1; i++)  x[i] = log(f[i]*tot);
+     if(LastItem) x[n-1] = 0;
+     }
+     else {         x => f
+     for(i=0,tot=1; i<n-1; i++)  tot  += (f[i]=exp(x[i]));
+     for(i=0; i<n-1; i++)        f[i] /= tot;
+     if(LastItem) f[n-1] = 1/tot;
+     }
+     return(0);
+     }*/
+
     aProportions[0] = exp(aV0);
     aProportions[1] = exp(aV1);
     double tot = aProportions[0] + aProportions[1] + 1;
     aProportions[0] /= tot;
     aProportions[1] /= tot;
+
+    // std::cout << "p[0]" << aProportions[0]<< "p[1]" << aProportions[1] <<
+    // std::endl;
+
     tot = aProportions[0] + aProportions[1];
 
     aProportions[2] = (1. - tot) * aProportions[0] / tot;
@@ -265,6 +330,12 @@ protected:
   ///
   void initVariables(void);
 
+  /// Initialize variables to be optimized in the same way as codeml.
+  /// It uses mInitType to know what has been already initialized by
+  /// initFromTree() or initFromResult()
+  ///
+  void initVariablesCodeml(void);
+
 private:
   /// Set upper and lower limits for the maximization domain
   ///
@@ -592,4 +663,107 @@ private:
   double mScaleQ1;     ///< Scale value for Q1
 };
 
+/// Alternate Hypothesis test (multiple fg grounds).
+///
+///		@author Mario Valle - Swiss National Supercomputing Centre
+///(CSCS)
+///		@date 2010-12-23 (initial version)
+///		@version 1.1
+///
+///
+class MfgBranchSiteModelAltHyp : public BranchSiteModel {
+public:
+  /// Constructor.
+  ///
+  /// @param[in] aForest The forest for which the maximum likelihood should be
+  /// computed
+  /// @param[in] aCmdLine The command line parameters
+  ///
+  MfgBranchSiteModelAltHyp(Forest &aForest, const CmdLine &aCmdLine)
+      : BranchSiteModel(
+            aForest, aForest.getNumBranches(), aForest.getNumSites(),
+            aCmdLine.mSeed, 5, aCmdLine.mNoMaximization, aCmdLine.mTrace,
+            aCmdLine.mOptimizationAlgo, aCmdLine.mDeltaValueForGradient,
+            aCmdLine.mRelativeError,
+            aCmdLine.mForceSerial || aCmdLine.mDoNotReduceForest,
+            aCmdLine.mVerboseLevel, aCmdLine.mExtraDebug,
+            aCmdLine.mMaxIterations, aCmdLine.mFixedBranchLength),
+        mfgmSet(aForest.getNumBranches()),
+        mfgmSetForGradient(aForest.getNumBranches()), mPrevK(DBL_MAX),
+        mPrevOmega0(DBL_MAX), mPrevOmega2(DBL_MAX) {
+    // Initialize the dependency set
+    mDependencies.computeDependencies(4, mNoParallel);
+    mDependencies.print("TEST FOR H1 (before optimization)");
+    mDependencies.optimizeDependencies();
+    mDependencies.print("TEST FOR H1");
+  }
+
+  /// Compute the alternative hypothesis log likelihood.
+  ///
+  /// @param[in] aFgBranch The identifier for the branch marked as foreground
+  /// branch
+  ///
+  /// @return The log likelihood under the alternative hypothesis
+  ///
+  double operator()(std::set<int> aFgBranchSet);
+
+  /// Compute the likelihood for the given forest and the given set of
+  /// parameters when computing gradient.
+  ///
+  /// @param[in] aVar The optimizer variables
+  /// @param[in] aTrace If set visualize the best result so far
+  /// @param[in] aGradientVar Used in gradient computation to avoid unneeded
+  /// computations. If set to UINT_MAX it is ignored
+  ///
+  /// @return The maximum Likelihood value
+  ///
+  double computeLikelihoodForGradient(const std::vector<double> &aVar,
+                                      bool aTrace, size_t aGradientVar);
+
+  /// Compute the likelihood for the given forest and the given set of
+  /// parameters.
+  ///
+  /// @param[in] aVar The optimizer variables
+  /// @param[in] aTrace If set visualize the best result so far
+  ///
+  /// @return The maximum Likelihood value
+  ///
+  double computeLikelihood(const std::vector<double> &aVar, bool aTrace);
+
+private:
+  /// Disabled assignment operator to avoid warnings on Windows.
+  ///
+  /// @fn BranchSiteModelAltHyp& operator=(const BranchSiteModelAltHyp& aObj)
+  ///
+  /// @param[in] aObj The object to be assigned
+  ///
+  /// @return The object receiving the assignment
+  ///
+  MfgBranchSiteModelAltHyp &
+  operator=(const MfgBranchSiteModelAltHyp & /*aObj*/);
+
+  /// Combine the sites' various codon classes likelihoods into one
+  /// log-likelihood value
+  ///
+  /// @return The tree log-likelihood value
+  ///
+  double combineSiteLikelihoods(void);
+
+private:
+  CheckpointableTransitionMatrix mQw0; ///< Q matrix for the omega0 case
+  TransitionMatrix mQw2;               ///< Q matrix for the omega2 case
+  CheckpointableTransitionMatrix mQ1;  ///< Q matrix for the omega1 == 1 case
+  mfgProbabilityMatrixSetH1 mfgmSet;   ///< Set of matrices used for the tree
+                                       ///visits (multiple foregrounds)
+  mfgProbabilityMatrixSetH1 mfgmSetForGradient; ///< Set of matrices used for
+                                                ///the tree visits (multiple
+                                                ///foregrounds)
+  double mPrevK;      ///< Previous k value used to compute matrices
+  double mPrevOmega0; ///< Previous w0 value used to compute matrices
+  double mPrevOmega2; ///< Previous w2 value used to compute matrices
+  double mScaleQw0;   ///< Scale value for Qw0
+  double mScaleQw2;   ///< Scale value for Qw2
+  double mScaleQ1;    ///< Scale value for Q1
+};
+
 #endif

CHANGELOG

+FastCodeML 1.3.0:
+- Support for trees with length 0 branches
+- Supporting rooted trees by unrooting them first
+- Support for terminal branches as foreground
+- Supporting multiple foreground branches
+- More consistent results with codeml by unifying initial points and distributions
+-- Detail
+- Better installations guide, standard output and command line options
+- Correcting bugs in converting external and internal proportions
+
 FastCodeML 1.2.0:
 - (Major) BG/Q compatibility (detail below)
 - (Minor) Removal of -np option (now controlled by -nt), cleanup of std::cout in main.
@@ -25,6 +35,7 @@ FastCodeML 1.2.0:
 - fast.cpp
 	> Removal of -np option.
 - SETPATHS - example explaining how to set envionment variables per machine.
+
 FastCodeML 1.1.0:
 - Few corrections to accept pipe character also in sequence names
 - New option to fixed branch lengths (-bl)

CmdLine.cpp

@@ -116,9 +116,9 @@ void CmdLine::parseCmdLine(int aCnt, char **aVal) {
     OPT_QUIET,
     OPT_HELP,
     OPT_SEED,
-    OPT_BRANCH,
-    OPT_BRANCH_START,
-    OPT_BRANCH_END,
+    OPT_BRANCH_ALL,
+    // OPT_BRANCH_START,
+    // OPT_BRANCH_END,
     OPT_IGNORE_FREQ,
     OPT_EXPORT,
     OPT_NOT_REDUCE,
@@ -146,69 +146,63 @@ void CmdLine::parseCmdLine(int aCnt, char **aVal) {
 
   // Then the definitions of each command line option
   CSimpleOpt::SOption parser_options[] = {
-      {OPT_VERBOSE, "-d", SO_REQ_SEP, "Verbosity level (0: none; 1: results "
+
+      // reporting parameters
+      {OPT_HELP, "-h", SO_NONE, "This help"},
+      //{ OPT_HELP,				"-h",
+      //SO_NONE,	"" },
+      {OPT_HELP, "--help", SO_NONE, ""},
+      //{ OPT_VERBOSE,			"-d",
+      //SO_REQ_SEP, "Verbosity level (0: none; 1: results only; 2: normal info;
+      //3: MPI trace; 4: more debug) (default: 1)" },
+      //{ OPT_VERBOSE,			"--debug",
+      //SO_REQ_SEP, "" },
+      {OPT_VERBOSE, "-v", SO_REQ_SEP, "Verbosity level (0: none; 1: results "
                                       "only; 2: normal info; 3: MPI trace; 4: "
                                       "more debug) (default: 1)"},
-      {OPT_VERBOSE, "--debug", SO_REQ_SEP, ""},
-      {OPT_VERBOSE, "-v", SO_REQ_SEP, ""},
       {OPT_VERBOSE, "--verbose", SO_REQ_SEP, ""},
       {OPT_QUIET, "-q", SO_NONE, "No messages except results"},
       {OPT_QUIET, "--quiet", SO_NONE, ""},
-      {OPT_HELP, "-?", SO_NONE, "This help"},
-      {OPT_HELP, "-h", SO_NONE, ""},
-      {OPT_HELP, "--help", SO_NONE, ""},
-      {OPT_SEED, "-s", SO_REQ_SEP,
-       "Random number generator seed (0 < seed < 1000000000)"},
-      {OPT_SEED, "--seed", SO_REQ_SEP, ""},
-      {OPT_BRANCH, "-b", SO_REQ_SEP,
-       "Do only this branch as foreground branch (count from 0)"},
-      {OPT_BRANCH, "--branch", SO_REQ_SEP, ""},
-      {OPT_BRANCH_START, "-bs", SO_REQ_SEP, "Start computing from this branch "
-                                            "as foreground one (count from 0) "
-                                            "(default: first one)"},
-      {OPT_BRANCH_START, "--branch-start", SO_REQ_SEP, ""},
-      {OPT_BRANCH_END, "-be", SO_REQ_SEP, "End computing at this branch as "
-                                          "foreground one (count from 0) "
-                                          "(default: last one)"},
-      {OPT_BRANCH_END, "--branch-end", SO_REQ_SEP, ""},
-      {OPT_IGNORE_FREQ, "-i", SO_NONE,
-       "Ignore computed codon frequency and set all of them to 1/61"},
-      {OPT_IGNORE_FREQ, "--ignore-freq", SO_NONE, ""},
       {OPT_EXPORT, "-e", SO_REQ_SEP, "Export forest in GML format (if %03d or "
                                      "@03d is present, one is created for each "
                                      "fg branch)"},
       {OPT_EXPORT, "--export", SO_REQ_SEP, ""},
-      {OPT_NOT_REDUCE, "-nr", SO_NONE,
-       "Do not reduce forest by merging common subtrees"},
-      {OPT_NOT_REDUCE, "--no-reduce", SO_NONE, ""},
-      {OPT_TIMES_FROM_FILE, "-l", SO_NONE,
-       "Initial branch lengths from tree file"},
-      {OPT_TIMES_FROM_FILE, "--lengths-from-file", SO_NONE, ""},
-      {OPT_TIMES_FROM_FILE, "--times-from-file", SO_NONE, ""},
-      {OPT_ONE_STEP, "-o", SO_NONE,
-       "Only the initial step is performed (no maximization)"},
-      {OPT_ONE_STEP, "--initial-step", SO_NONE, ""},
-      {OPT_COMP_TIMES, "-c", SO_REQ_SEP, "Export the computed times from H0 if "
-                                         "0, H1 if 1, otherwise the one read "
-                                         "in the phylo tree"},
-      {OPT_COMP_TIMES, "--export-comp-times", SO_REQ_SEP, ""},
       {OPT_TRACE, "-r", SO_NONE, "Trace the maximization run"},
       {OPT_TRACE, "--trace", SO_NONE, ""},
+      {OPT_EXTRA_DEBUG, "-x", SO_REQ_SEP,
+       "Extra debug parameter (zero disables it)"},
+      {OPT_EXTRA_DEBUG, "--extra-debug", SO_REQ_SEP, ""},
+      {OPT_OUT_RESULTS, "-ou", SO_REQ_SEP,
+       "Write results formatted to this file"},
+      {OPT_OUT_RESULTS, "--output", SO_REQ_SEP, ""},
+      {OPT_CLEAN_DATA, "-cl", SO_NONE,
+       "Remove ambiguous or missing sites from the MSA (default: no)"},
+      {OPT_CLEAN_DATA, "--clean-data", SO_NONE, ""},
+
+      // configuration parameters
+
+      {OPT_SEED, "-s", SO_REQ_SEP,
+       "Random number generator seed (0 < seed < 1000000000)"},
+      {OPT_SEED, "--seed", SO_REQ_SEP, ""},
       {OPT_NUM_THREADS, "-nt", SO_REQ_SEP,
        "Number of threads (1 for non parallel execution)"},
       {OPT_NUM_THREADS, "--number-of-threads", SO_REQ_SEP, ""},
-      //{ OPT_FORCE_SERIAL,		"-np",
-      // SO_NONE,	"Don't use parallel execution" },
-      //{ OPT_FORCE_SERIAL,		"--no-parallel",
-      // SO_NONE,	"" },
-      {OPT_BRANCH_FROM_FILE, "-bf", SO_NONE,
-       "Do only the branch marked in the file as foreground branch"},
-      {OPT_BRANCH_FROM_FILE, "--branch-from-file", SO_NONE, ""},
-      {OPT_ONE_HYP_ONLY, "-hy", SO_REQ_SEP, "Compute only H0 if 0, H1 if 1"},
-      {OPT_ONE_HYP_ONLY, "--only-hyp", SO_REQ_SEP, ""},
-      {OPT_INIT_H0_FROM_H1, "-i1", SO_NONE,
-       "Start H0 optimization from H1 results"},
-      {OPT_INIT_H0_FROM_H1, "--init-from-h1", SO_NONE, ""},
+      //	{ OPT_BRANCH_START,		"-bs",
+      //SO_REQ_SEP, "Start computing from this branch as foreground one (count
+      //from 0) (default: first one)" },
+      //	{ OPT_BRANCH_START,		"--branch-start",
+      //SO_REQ_SEP, "" },
+      //	{ OPT_BRANCH_END,		"-be",
+      //SO_REQ_SEP, "End computing at this branch as foreground one (count from
+      //0) (default: last one)" },
+      //	{ OPT_BRANCH_END,		"--branch-end",
+      //SO_REQ_SEP, "" },
+      {OPT_IGNORE_FREQ, "-i", SO_NONE,
+       "Ignore computed codon frequency and set all of them to 1/61"},
+      {OPT_IGNORE_FREQ, "--ignore-freq", SO_NONE, ""},
+      {OPT_NOT_REDUCE, "-nr", SO_NONE,
+       "Do not reduce forest by merging common subtrees"},
+      {OPT_NOT_REDUCE, "--no-reduce", SO_NONE, ""},
       {OPT_OPTIM_ALGO, "-m", SO_REQ_SEP,
        "Optimizer algorithm (0:LBFGS, 1:VAR1, 2:VAR2, 3:SLSQP, 11:BOBYQA, "
        "22:FromCodeML, 99:MLSL_LDS) (default: 22)"},
@@ -216,33 +210,58 @@ void CmdLine::parseCmdLine(int aCnt, char **aVal) {
       {OPT_DELTA_VAL, "-sd", SO_REQ_SEP,
        "Delta used in gradient computation (default: 1.49e-8)"},
       {OPT_DELTA_VAL, "--small-diff", SO_REQ_SEP, ""},
-      {OPT_INIT_PARAM, "-p", SO_REQ_SEP, "Pass initialization parameter in the "
-                                         "form: P=value (P: w0, k, p0, p1, "
-                                         "w2)"},
-      {OPT_INIT_PARAM, "--init-param", SO_REQ_SEP, ""},
-      {OPT_INIT_DEFAULT, "-ic", SO_NONE,
-       "Start from default parameter values and times from tree file"},
-      {OPT_INIT_DEFAULT, "--init-default", SO_NONE, ""},
-      {OPT_EXTRA_DEBUG, "-x", SO_REQ_SEP,
-       "Extra debug parameter (zero disables it)"},
-      {OPT_EXTRA_DEBUG, "--extra-debug", SO_REQ_SEP, ""},
       {OPT_REL_ERROR, "-re", SO_REQ_SEP,
        "Relative error where to stop maximization (default: 1e-3)"},
       {OPT_REL_ERROR, "--relative-error", SO_REQ_SEP, ""},
-      {OPT_OUT_RESULTS, "-ou", SO_REQ_SEP,
-       "Write results formatted to this file"},
-      {OPT_OUT_RESULTS, "--output", SO_REQ_SEP, ""},
-      {OPT_CLEAN_DATA, "-cl", SO_NONE,
-       "Remove ambiguous or missing sites from the MSA (default: no)"},
-      {OPT_CLEAN_DATA, "--clean-data", SO_NONE, ""},
-      {OPT_NO_PRE_STOP, "-ps", SO_NONE, "Don't stop H0 maximization even if it "
-                                        "cannot satisfy LRT (default: stop)"},
-      {OPT_NO_PRE_STOP, "--no-pre-stop", SO_NONE, ""},
       {OPT_MAX_ITER, "-mi", SO_REQ_SEP,
        "Maximum number of iterations for the maximizer (default: 10000)"},
       {OPT_MAX_ITER, "--max-iterations", SO_REQ_SEP, ""},
-      {OPT_BRANCH_LENGTH, "-bl", SO_NONE, "The length of the brances is fixed"},
+
+      //{ OPT_TIMES_FROM_FILE,	"--times-from-file",	SO_NONE,	"" },
+
+      //{ OPT_FORCE_SERIAL,		"-np",
+      //SO_NONE,	"Don't use parallel execution" },
+      //{ OPT_FORCE_SERIAL,		"--no-parallel",
+      //SO_NONE,	"" },
+      //{ OPT_BRANCH_FROM_FILE,	"-bf",
+      //SO_NONE,	"Do only the branch marked in the file as foreground
+      //branch" },
+      //{ OPT_BRANCH_FROM_FILE,	"--branch-from-file",	SO_NONE,	"" },
+      {OPT_ONE_HYP_ONLY, "-hy", SO_REQ_SEP, "Compute only H0 if 0, H1 if 1"},
+      {OPT_ONE_HYP_ONLY, "--only-hyp", SO_REQ_SEP, ""},
+      {OPT_INIT_H0_FROM_H1, "-i1", SO_NONE,
+       "Start H0 optimization from H1 results"},
+      {OPT_INIT_H0_FROM_H1, "--init-from-h1", SO_NONE, ""},
+      {OPT_ONE_STEP, "-o", SO_NONE,
+       "Only the initial step is performed (no maximization)"},
+      {OPT_ONE_STEP, "--initial-step", SO_NONE, ""},
+      {OPT_NO_PRE_STOP, "-ps", SO_NONE, "Don't stop H0 maximization even if it "
+                                        "cannot satisfy LRT (default: stop)"},
+      {OPT_NO_PRE_STOP, "--no-pre-stop", SO_NONE, ""},
+
+      // initialization parameters
+
+      {OPT_TIMES_FROM_FILE, "-l", SO_NONE,
+       "Initial branch lengths from tree file"},
+      {OPT_TIMES_FROM_FILE, "--blengths-from-file", SO_NONE, ""},
+      {OPT_BRANCH_LENGTH, "-bl", SO_NONE,
+       "The length of the branches is fixed"},
       {OPT_BRANCH_LENGTH, "--branch-lengths-fixed", SO_NONE, ""},
+      {OPT_BRANCH_ALL, "-ba", SO_NONE,
+       "Do for all branches as foreground branch (including leaves)"},
+      {OPT_BRANCH_ALL, "--branch-all", SO_NONE, ""},
+      {OPT_COMP_TIMES, "-c", SO_REQ_SEP, "Export the computed times from H0 if "
+                                         "0, H1 if 1, otherwise the one read "
+                                         "in the phylo tree"},
+      {OPT_COMP_TIMES, "--export-comp-times", SO_REQ_SEP, ""},
+      {OPT_INIT_PARAM, "-p", SO_REQ_SEP, "Pass initialization parameter in the "
+                                         "form: P=value (P: w0, k, p0, p1, w2) "
+                                         "-> e.g. -p w0=1.1 -p p0=2.1 ..."},
+      {OPT_INIT_PARAM, "--init-param", SO_REQ_SEP, ""},
+      {OPT_INIT_DEFAULT, "-ic", SO_NONE,
+       "Start from default parameter values and times from tree file"},
+      {OPT_INIT_DEFAULT, "--init-default", SO_NONE, ""},
+
       SO_END_OF_OPTIONS};
 
   // Setup the usage string
@@ -303,26 +322,27 @@ void CmdLine::parseCmdLine(int aCnt, char **aVal) {
       mSeed = static_cast<unsigned int>(tmpi);
       break;
 
-    case OPT_BRANCH:
-      tmpi = atoi(args.OptionArg());
-      if (tmpi < 0)
-        throw FastCodeMLFatal("Invalid branch value");
-      mBranchStart = mBranchEnd = static_cast<unsigned int>(tmpi);
+    case OPT_BRANCH_ALL:
+      mBranchAll = true;
+      //		tmpi = atoi(args.OptionArg());
+      //		if(tmpi < 0) throw FastCodeMLFatal("Invalid branch
+      //value");
+      //		mBranchStart = mBranchEnd = static_cast<unsigned
+      //int>(tmpi);
       break;
 
-    case OPT_BRANCH_START:
-      tmpi = atoi(args.OptionArg());
-      if (tmpi < 0)
-        throw FastCodeMLFatal("Invalid start branch value");
-      mBranchStart = static_cast<unsigned int>(tmpi);
-      break;
+    //	case OPT_BRANCH_START:
+    //		tmpi = atoi(args.OptionArg());