Commit 3b60da8f authored by akalantz's avatar akalantz
Browse files

branch lenghts fixed attribute added (-bl)

git-svn-id: https://svn.vital-it.ch/svn/hp2c/trunk/Codeml_Variants/Fastcodeml@7149 95c5a404-1f88-0410-a6b0-c3f062b6f34a
parent 1d5901d3
This diff is collapsed.
......@@ -53,9 +53,10 @@ protected:
bool aNoParallel,
unsigned int aVerbose,
unsigned int aExtraDebug,
unsigned int aMaxIterations)
: mForest(aForest),
mVar(aNumBranches+aNumVariables),
unsigned int aMaxIterations,
bool aFixedBranchLength)
: mForest(aForest),
mVar(aFixedBranchLength ? aNumVariables : aNumBranches+aNumVariables), //mVar(aNumBranches+aNumVariables),
mFgScale(0.0),
mBgScale(0.0),
mMaxLnL(-DBL_MAX),
......@@ -74,9 +75,11 @@ protected:
mDependencies(aForest, aVerbose),
mNoParallel(aNoParallel),
mSeed(aSeed),
mRelativeError(aRelativeError)
mRelativeError(aRelativeError),
mFixedBranchLength(aFixedBranchLength),
mBranches(aNumBranches)
{
setLimits(aNumBranches, static_cast<size_t>(aNumVariables));
setLimits(aNumBranches, static_cast<size_t>(aNumVariables), aFixedBranchLength);
}
/// Destructor.
......@@ -177,7 +180,7 @@ public:
/// @return True if the test is passed
///
static bool performLRT(double aLnL0, double aLnL1) {return (aLnL1 - aLnL0) > THRESHOLD_FOR_LRT;}
/// Get site multiplicity values.
///
/// @return Reference to the array of site multiplicities
......@@ -247,8 +250,9 @@ private:
///
/// @param[in] aNumTimes Number of times (i.e.\ branch lengths)
/// @param[in] aNumVariables Number of other variables (4 for H0, 5 for H1)
/// @param[in] aFixedBranch
///
void setLimits(size_t aNumTimes, size_t aNumVariables);
void setLimits(size_t aNumTimes, size_t aNumVariables, bool aFixedBranchLength);
/// Generate a double random number between 0 and 1
///
......@@ -339,6 +343,8 @@ protected:
unsigned int mMaxIterations; ///< Maximum number of iterations for the maximization
TreeAndSetsDependencies mDependencies; ///< The dependency list between trees to use in this run
bool mNoParallel; ///< True if no preparation for multithreading should be done
bool mFixedBranchLength; ///< True if branch lengths are kept fixed (not optimised)
std::vector<double> mBranches; ///< Variable with the branch lengths
private:
unsigned int mSeed; ///< Random number generator seed to be used also by the optimizer
......@@ -367,7 +373,7 @@ public:
aCmdLine.mSeed, 4, aCmdLine.mNoMaximization, aCmdLine.mTrace,
aCmdLine.mOptimizationAlgo, aCmdLine.mDeltaValueForGradient,
aCmdLine.mRelativeError, aCmdLine.mForceSerial || aCmdLine.mDoNotReduceForest,
aCmdLine.mVerboseLevel, aCmdLine.mExtraDebug, aCmdLine.mMaxIterations),
aCmdLine.mVerboseLevel, aCmdLine.mExtraDebug, aCmdLine.mMaxIterations, aCmdLine.mFixedBranchLength),
mSet(aForest.getNumBranches()), mSetForGradient(aForest.getNumBranches()),
mPrevK(DBL_MAX), mPrevOmega0(DBL_MAX)
{
......@@ -458,7 +464,7 @@ public:
aCmdLine.mSeed, 5, aCmdLine.mNoMaximization, aCmdLine.mTrace,
aCmdLine.mOptimizationAlgo, aCmdLine.mDeltaValueForGradient,
aCmdLine.mRelativeError, aCmdLine.mForceSerial || aCmdLine.mDoNotReduceForest,
aCmdLine.mVerboseLevel, aCmdLine.mExtraDebug, aCmdLine.mMaxIterations),
aCmdLine.mVerboseLevel, aCmdLine.mExtraDebug, aCmdLine.mMaxIterations, aCmdLine.mFixedBranchLength),
mSet(aForest.getNumBranches()), mSetForGradient(aForest.getNumBranches()),
mPrevK(DBL_MAX), mPrevOmega0(DBL_MAX), mPrevOmega2(DBL_MAX)
{
......
......@@ -186,7 +186,7 @@ void CmdLine::parseCmdLine(int aCnt, char **aVal)
{ 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_OPTIM_ALGO, "-m", SO_REQ_SEP, "Optimizer algorithm (0:LBFGS, 1:VAR1, 2:VAR2, 3:SLSQP, 11:BOBYQA, 22:FromCodeML, 99:MLSL_LDS) (default: 0)" },
{ 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)" },
{ OPT_OPTIM_ALGO, "--maximizer", SO_REQ_SEP, "" },
{ OPT_DELTA_VAL, "-sd", SO_REQ_SEP, "Delta used in gradient computation (default: 1.49e-8)" },
{ OPT_DELTA_VAL, "--small-diff", SO_REQ_SEP, "" },
......@@ -205,7 +205,7 @@ void CmdLine::parseCmdLine(int aCnt, char **aVal)
{ 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_MAX_ITER, "--max-iterations", SO_REQ_SEP, "" },
{OPT_BRANCH_LENGTH, "-bl", SO_NONE, "The length of the brances is fixed (only for KCM and algorithms LBFGS and SLSQP)"},
{OPT_BRANCH_LENGTH, "--branch-lengths-fixed", SO_NONE, ""},
SO_END_OF_OPTIONS
......@@ -377,10 +377,11 @@ void CmdLine::parseCmdLine(int aCnt, char **aVal)
/* if (mNumThreads == 1)
mForceSerial = true;*/
if (mNumThreads <=0) throw FastCodeMLFatal("Invalid number of threads");
break;
break;
case OPT_BRANCH_LENGTH:
//mFixedBranchLength = true;
mFixedBranchLength = true;
mBranchLengthsFromFile = true;
break;
}
}
......
......@@ -60,7 +60,8 @@ public:
mInitFromParams(false),
mCleanData(false),
mStopIfNotLRT(true),
mCmdLineImpl(NULL)
mCmdLineImpl(NULL),
mFixedBranchLength(false)
{}
/// Destructor.
......@@ -106,6 +107,7 @@ public:
bool mCleanData; ///< Remove ambiguous or missing sites from the MSA (genes)
bool mStopIfNotLRT; ///< Stop H0 maximization when LRT cannot be satisfied
unsigned int mNumThreads; ///< Number of threads (if 1 the parallelization is disabled)
bool mFixedBranchLength; ///<fixed branch lengths
private:
struct CmdLineImpl;
......
......@@ -193,7 +193,7 @@ void PhyloTree::checkRootBranches(void) const
{
std::cout << std::endl << "Root has " << cnt_root_branches << " children of which " << cnt_root_leaves << " are leaves" << std::endl;
}
if (cnt_root_branches == 2)
{
std::cout << std::endl << "This is a rooted tree. Please check!" << std::endl;
......
......@@ -70,23 +70,23 @@ int main(int aRgc, char **aRgv)
// Adjust and report the number of threads that will be used
#ifdef _OPENMP
int num_threads = omp_get_max_threads();
//std::cout<<"max num of thr: "<< num_threads <<std::endl;
if((cmd.mNumThreads >=1)&&(cmd.mNumThreads <= 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)
cmd.mForceSerial = true;
else
int num_threads = omp_get_max_threads();
//std::cout<<"max num of thr: "<< num_threads <<std::endl;
if((cmd.mNumThreads >=1)&&(cmd.mNumThreads <= 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)
cmd.mForceSerial = true;
else
cmd.mForceSerial = false;*/
#else
#else
cmd.mNumThreads=1;
int num_threads = 1;
int num_threads = 1;
cmd.mForceSerial = true;
#endif
#endif
/*#ifdef _OPENMP
int num_threads = omp_get_max_threads();
......@@ -99,7 +99,7 @@ int main(int aRgc, char **aRgv)
#else
cmd.mForceSerial = true;
int num_threads = 1;
#endif*/
#endif*/
#ifdef USE_MPI
// Shutdown messages from all MPI processes except the master
......@@ -142,8 +142,8 @@ int main(int aRgc, char **aRgv)
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::cerr << "Branch lengths are fixed" << std::endl;
if(cmd.mNumThreads) std::cout << "Number of threads: " << cmd.mNumThreads << std::endl;
if(cmd.mFixedBranchLength) std::cerr << "Branch lengths are fixed" << std::endl;
#ifdef _OPENMP
if(num_threads > 1)
{
......@@ -379,11 +379,11 @@ int main(int aRgc, char **aRgv)
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);
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;
}
......@@ -396,11 +396,11 @@ int main(int aRgc, char **aRgv)
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);
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;
}
......@@ -596,8 +596,8 @@ Usage:
-r --trace (no argument)
Trace the maximization run
-nt --number-of-threads (required argument)
Number of threads (1 for non parallel execution)
-nt --number-of-threads (required argument)
Number of threads (1 for non parallel execution)
-bf --branch-from-file (no argument)
Do only the branch marked in the file as foreground branch
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment