Commit 64f47277 authored by mvalle's avatar mvalle
Browse files

Added exception for LRT not satisfied.

Added maximization stop after given number of iterations.
Added new cmd line options to define maximum number of iterations.


git-svn-id: https://svn.vital-it.ch/svn/hp2c/trunk/Codeml_Variants/Fastcodeml@6593 95c5a404-1f88-0410-a6b0-c3f062b6f34a
parent 9bd7b64a
......@@ -960,7 +960,8 @@ public:
/// @param[in] aStopIfBigger If true stop computation as soon as value is over aThreshold
/// @param[in] aThreshold The threshold at which the maximization should be stopped
///
MaximizerFunction(BranchSiteModel* aModel, bool aTrace, const std::vector<double>& aUpper, double aDeltaForGradient, size_t aNumMatrixParams, bool aStopIfBigger, double aThreshold)
MaximizerFunction(BranchSiteModel* aModel, bool aTrace, const std::vector<double>& aUpper, double aDeltaForGradient, size_t aNumMatrixParams,
bool aStopIfBigger, double aThreshold)
: mModel(aModel), mTrace(aTrace), mUpper(aUpper), mDeltaForGradient(aDeltaForGradient),
mTotalNumVariables(aUpper.size()), mNumBranchLengths(aUpper.size() - aNumMatrixParams),
mStopIfBigger(aStopIfBigger), mThreshold(aThreshold)
......@@ -1133,7 +1134,7 @@ double BranchSiteModel::maximizeLikelihood(size_t aFgBranch, bool aStopIfBigger,
try
{
// Create the optimizer (instead of mRelativeError is used the fixed value from CodeML)
Ming2 optim(this, mTrace, mVerbose, mLowerBound, mUpperBound, mDeltaForGradient, 1e-8, aStopIfBigger, aThreshold);
Ming2 optim(this, mTrace, mVerbose, mLowerBound, mUpperBound, mDeltaForGradient, 1e-8, aStopIfBigger, aThreshold, mMaxIterations);
// Do the maximization
double maxl = optim.minimizeFunction(mVar);
......@@ -1146,7 +1147,7 @@ double BranchSiteModel::maximizeLikelihood(size_t aFgBranch, bool aStopIfBigger,
}
return maxl;
}
catch(FastCodeMLSuccess&)
catch(FastCodeMLEarlyStopLRT&)
{
if(mTrace) std::cout << "Optimization stopped because LRT not satisfied" << std::endl;
return DBL_MAX;
......@@ -1217,6 +1218,9 @@ double BranchSiteModel::maximizeLikelihood(size_t aFgBranch, bool aStopIfBigger,
opt->set_max_objective(MaximizerFunction::wrapFunction, &compute);
// If the user has set a maximum number of iterations set it
if(mMaxIterations != MAX_ITERATIONS) opt->set_maxeval(mMaxIterations);
nlopt::result result = opt->optimize(mVar, maxl);
// Print the final optimum value
......
......@@ -38,6 +38,7 @@ protected:
/// @param[in] aNoParallel If true no parallel execution support is setup
/// @param[in] aVerbose The verbosity level
/// @param[in] aExtraDebug Extra parameter for testing during development
/// @param[in] aMaxIterations Maximum number of iterations for the maximization
///
BranchSiteModel(Forest& aForest,
size_t aNumBranches,
......@@ -51,7 +52,8 @@ protected:
double aRelativeError,
bool aNoParallel,
unsigned int aVerbose,
unsigned int aExtraDebug)
unsigned int aExtraDebug,
unsigned int aMaxIterations)
: mForest(aForest),
mVar(aNumBranches+aNumVariables),
mFgScale(0.0),
......@@ -68,6 +70,7 @@ protected:
mExtraDebug(aExtraDebug),
mVerbose(aVerbose),
mNumEvaluations(0),
mMaxIterations(aMaxIterations),
mDependencies(aForest, aVerbose),
mNoParallel(aNoParallel),
mSeed(aSeed),
......@@ -333,6 +336,7 @@ protected:
unsigned int mExtraDebug; ///< Parameter for extra development testing
unsigned int mVerbose; ///< Parameter for extra development testing
unsigned int mNumEvaluations; ///< Counter of the likelihood function evaluations
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
......@@ -362,7 +366,8 @@ public:
: BranchSiteModel(aForest, aForest.getNumBranches(), aForest.getNumSites(),
aCmdLine.mSeed, 4, aCmdLine.mNoMaximization, aCmdLine.mTrace,
aCmdLine.mOptimizationAlgo, aCmdLine.mDeltaValueForGradient,
aCmdLine.mRelativeError, aCmdLine.mForceSerial || aCmdLine.mDoNotReduceForest, aCmdLine.mVerboseLevel, aCmdLine.mExtraDebug),
aCmdLine.mRelativeError, aCmdLine.mForceSerial || aCmdLine.mDoNotReduceForest,
aCmdLine.mVerboseLevel, aCmdLine.mExtraDebug, aCmdLine.mMaxIterations),
mSet(aForest.getNumBranches()), mSetForGradient(aForest.getNumBranches()),
mPrevK(DBL_MAX), mPrevOmega0(DBL_MAX)
{
......@@ -452,7 +457,8 @@ public:
: 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.mRelativeError, aCmdLine.mForceSerial || aCmdLine.mDoNotReduceForest,
aCmdLine.mVerboseLevel, aCmdLine.mExtraDebug, aCmdLine.mMaxIterations),
mSet(aForest.getNumBranches()), mSetForGradient(aForest.getNumBranches()),
mPrevK(DBL_MAX), mPrevOmega0(DBL_MAX), mPrevOmega2(DBL_MAX)
{
......
......@@ -136,7 +136,8 @@ void CmdLine::parseCmdLine(int aCnt, char **aVal)
OPT_REL_ERROR,
OPT_OUT_RESULTS,
OPT_CLEAN_DATA,
OPT_NO_PRE_STOP
OPT_NO_PRE_STOP,
OPT_MAX_ITER
};
// Then the definitions of each command line option
......@@ -199,6 +200,8 @@ void CmdLine::parseCmdLine(int aCnt, char **aVal)
{ 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, "" },
SO_END_OF_OPTIONS
};
......@@ -356,6 +359,12 @@ void CmdLine::parseCmdLine(int aCnt, char **aVal)
case OPT_NO_PRE_STOP:
mStopIfNotLRT = false;
break;
case OPT_MAX_ITER:
tmpi = atoi(args.OptionArg());
if(tmpi <= 0) throw FastCodeMLFatal("Invalid max number of iterations. Should be > 0");
mMaxIterations = static_cast<unsigned int>(tmpi);
break;
}
}
......
......@@ -5,6 +5,11 @@
#include <climits>
#include "VerbosityLevels.h"
/// The default maximum number of optimization steps.
///
static const unsigned int MAX_ITERATIONS=10000;
/// Parse the command line flags.
///
/// @author Mario Valle - Swiss National Supercomputing Centre (CSCS)
......@@ -34,6 +39,7 @@ public:
mComputeHypothesis(UINT_MAX),
mOptimizationAlgo(0),
mExtraDebug(0),
mMaxIterations(MAX_ITERATIONS),
mIgnoreFreq(false),
mDoNotReduceForest(false),
mBranchLengthsFromFile(false),
......@@ -78,6 +84,7 @@ public:
unsigned int mComputeHypothesis; ///< If set to 0 compute only H0, if set to 1 compute only H1, otherwise compute both
unsigned int mOptimizationAlgo; ///< Select the optimization algorithm to use
unsigned int mExtraDebug; ///< Extra debug parameter for development tests
unsigned int mMaxIterations; ///< Maximum number of iterations for the maximization
bool mIgnoreFreq; ///< Ignore the computed codon frequencies and set them all to 1/61
bool mDoNotReduceForest; ///< If true do not reduce the forest merging common subtrees
bool mBranchLengthsFromFile; ///< The initial value of the branch lengths is taken from the phylo tree file
......
......@@ -58,10 +58,6 @@ static inline double distance(const double* RESTRICT x, const double* RESTRICT y
return sqrt(t);
}
/// The number of maximal optimization steps (was variable 'maxround')
///
static const int MAX_ITERATIONS=10000;
double Ming2::minimizeFunction(std::vector<double>& aVars)
{
......@@ -176,10 +172,10 @@ int Ming2::ming2(FILE *fout, double *f, double x[], const double xl[], const dou
identityMatrix(H, nfree);
for(iround = 0; iround < MAX_ITERATIONS; ++iround)
for(iround = 0; iround < mMaxIterations; ++iround)
{
// Check if the optimization can be stopped in advance due to LRT non satisfied
if(mStopIfBigger && *f < mThreshold) throw FastCodeMLSuccess();
if(mStopIfBigger && *f < mThreshold) throw FastCodeMLEarlyStopLRT();
if(fout)
{
......@@ -232,7 +228,7 @@ int Ming2::ming2(FILE *fout, double *f, double x[], const double xl[], const dou
{
if(mAlwaysCenter)
{
iround = MAX_ITERATIONS;
iround = mMaxIterations;
break;
}
else
......@@ -383,7 +379,7 @@ int Ming2::ming2(FILE *fout, double *f, double x[], const double xl[], const dou
v += y[i] * s[i];
}
if (fabs(v) < small)
if(fabs(v) < small)
{
identityMatrix(H, nfree);
fail = 1;
......@@ -393,22 +389,23 @@ int Ming2::ming2(FILE *fout, double *f, double x[], const double xl[], const dou
for(i=0; i < nfree; ++i)
for(j=0; j < nfree; ++j)
H[i*nfree + j] += ((1 + w / v) * s[i] * s[j] - z[i] * s[j] - s[i] * z[j]) / v;
} /* for (iround,MAX_ITERATIONS) */
} // end of for(iround, mMaxIterations)
/* try to remove this after updating LineSearch2() */
//*f = (*fun) (x, n); ++mNumFunCall;
*f = -mModel->computeLikelihood(x, n, mTraceFun);
if (mNoisy > 2)
if(mNoisy > 2)
{
printf("\n");
}
if (iround == MAX_ITERATIONS)
if(iround == mMaxIterations)
{
if (fout)
{
fprintf(fout, "\ncheck convergence!\n");
fprintf(fout, "\ncheck convergence! (max number of iterations reached)\n");
}
return -1;
......
......@@ -26,8 +26,10 @@ public:
/// @param[in] aRelativeError Relative error to stop computation
/// @param[in] aStopIfBigger If true stop computation as soon as value is over aThreshold
/// @param[in] aThreshold The threshold at which the maximization should be stopped
/// @param[in] aMaxIterations Maximum number of iterations for the maximization
///
Ming2(BranchSiteModel* aModel, bool aTrace, unsigned int aVerbose, const std::vector<double>& aLowerBound, const std::vector<double>& aUpperBound, double aDeltaForGradient, double aRelativeError, bool aStopIfBigger, double aThreshold) :
Ming2(BranchSiteModel* aModel, bool aTrace, unsigned int aVerbose, const std::vector<double>& aLowerBound,
const std::vector<double>& aUpperBound, double aDeltaForGradient, double aRelativeError, bool aStopIfBigger, double aThreshold, int aMaxIterations) :
mModel(aModel),
mTrace(aTrace),
mTraceFun(false),
......@@ -38,6 +40,7 @@ public:
mVerbose(aVerbose),
mStopIfBigger(aStopIfBigger),
mThreshold(-aThreshold),
mMaxIterations(aMaxIterations),
mAlwaysCenter(false),
mNoisy((aTrace && aVerbose > 0) ? 9 : 0) {}
......@@ -65,7 +68,7 @@ private:
///
/// @return Optimization status (-1 check convergence; 0 success; 1 fail)
///
/// @exception FastCodeMLSuccess If the optimization is stop in advance because LRT not satisfied
/// @exception FastCodeMLEarlyStopLRT If the optimization has been stopped in advance because LRT is not satisfied
///
int ming2(FILE *fout, double *f, double x[], const double xl[], const double xu[], double space[], int ispace[], double rel_error, int n);
......@@ -150,6 +153,7 @@ private:
unsigned int mVerbose; ///< The verbose flag from the BranchSiteModel class
bool mStopIfBigger; ///< When true stop if lnL is bigger than mThreshold
double mThreshold; ///< Threshold for the early stop of optimization if LRT non satisfied (the value is stored with sign changed)
int mMaxIterations; ///< Maximum number of iterations for the maximization
private:
/// The following variables are from the original code
......
......@@ -34,7 +34,6 @@ public:
/// Early successful termination exception.
/// It does not signal a fatal error, but something like having printed the help.
/// Signals also the early termination of a maximization due to LRT no more satisfied.
///
class FastCodeMLSuccess : public std::exception
{
......@@ -45,6 +44,20 @@ public:
{}
};
/// Early termination due to LRT not satisfied exception.
/// It does not signal a fatal error, but that the optimization stopped because LRT cannot be satisfied.
///
class FastCodeMLEarlyStopLRT : public std::exception
{
public:
/// Early successful termination exception.
///
FastCodeMLEarlyStopLRT() : exception()
{}
};
/// Memory error in FastCodeML.
/// The message explains the reason in detail.
///
......
......@@ -119,6 +119,7 @@ int main(int aRgc, char **aRgv)
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;
......@@ -604,6 +605,9 @@ Usage:
-ps --no-pre-stop (no argument)
Don't stop H0 maximization even if it cannot satisfy LRT (default: stop)
-mi --max-iterations (required argument)
Maximum number of iterations for the maximizer (default: 10000)
@endverbatim
*/
......
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