Commit fa0dfd72 authored by mvalle's avatar mvalle
Browse files

Changed variable initialization logic.

Implemented initialize h0 from h1 also for MPI.


git-svn-id: https://svn.vital-it.ch/svn/hp2c/trunk/Codeml_Variants/Fastcodeml@5831 95c5a404-1f88-0410-a6b0-c3f062b6f34a
parent c99c26c9
......@@ -165,16 +165,13 @@ void BranchSiteModel::initFromTree(void)
mForest.setTimesFromLengths(mVar);
// Ask for initialization completion
mInitType = INIT_TYPE_TIMES;
mInitStatus |= INIT_TIMES;
}
void BranchSiteModel::initFromTreeAndParams(void)
void BranchSiteModel::initFromParams(void)
{
// Initialize branch lengths from the phylo tree
mForest.setTimesFromLengths(mVar);
// Get the parameters
// Get the parameters (the default values or the values set on the command line)
ParseParameters* params = ParseParameters::getInstance();
// Initialization as in CodeML (seems)
......@@ -192,18 +189,10 @@ void BranchSiteModel::initFromTreeAndParams(void)
mVar[mNumTimes+0] = p0+p1; // p0+p1
mVar[mNumTimes+1] = p0/(p0+p1); // p0/(p0+p1)
#endif
if(mNumVariables == 5 && mInitType != INIT_TYPE_RES_5)
{
mVar[mNumTimes+4] = params->getParameter("w2"); // w2
if(mNumVariables == 5) mVar[mNumTimes+4] = params->getParameter("w2"); // w2
// Ask for initialization completion
mInitType = INIT_TYPE_RES_5;
}
else
{
// Ask for initialization completion
mInitType = INIT_TYPE_RES_4;
}
// The parameters have been initializated
mInitStatus |= INIT_PARAMS;
}
void BranchSiteModel::initFromResult(const std::vector<double>& aPreviousResult, unsigned int aValidLen)
......@@ -215,7 +204,7 @@ void BranchSiteModel::initFromResult(const std::vector<double>& aPreviousResult,
if(aValidLen > mNumTimes+mNumVariables) aValidLen = mNumTimes+mNumVariables;
else if(aValidLen < mNumTimes)
{
mInitType = INIT_TYPE_NONE;
mInitStatus = INIT_NONE;
return;
}
else if(aValidLen < mNumTimes+4) aValidLen = mNumTimes;
......@@ -225,9 +214,9 @@ void BranchSiteModel::initFromResult(const std::vector<double>& aPreviousResult,
mVar.resize(static_cast<size_t>(mNumTimes+mNumVariables));
// Ask for initialization completion
if(aValidLen == mNumTimes) mInitType = INIT_TYPE_TIMES;
else if(aValidLen == mNumTimes+4) mInitType = INIT_TYPE_RES_4;
else mInitType = INIT_TYPE_RES_5;
if(aValidLen == mNumTimes) mInitStatus = INIT_TIMES;
else if(aValidLen == mNumTimes+4) mInitStatus = INIT_TIMES|INIT_PARAMS_H0;
else mInitStatus = INIT_TIMES|INIT_PARAMS;
}
......@@ -236,33 +225,33 @@ void BranchSiteModel::initVariables(void)
unsigned int i;
// Initialize time
if(mInitType == INIT_TYPE_NONE)
if((mInitStatus & INIT_TIMES) != INIT_TIMES)
{
for(i=0; i < mNumTimes; ++i) mVar[i] = randFrom0to1()*.1+0.01; // T
for(i=0; i < mNumTimes; ++i) mVar[i] = 0.01 + 0.1 * randFrom0to1(); // T
}
// Initialize w0, k, v1, v2
if(mInitType == INIT_TYPE_TIMES || mInitType == INIT_TYPE_NONE)
if((mInitStatus & INIT_PARAMS_H0) != INIT_PARAMS_H0)
{
#ifdef USE_ORIGINAL_PROPORTIONS
mVar[mNumTimes+0] = 1.0 + 0.2 * randFrom0to1(); // x0 -> p0
mVar[mNumTimes+1] = 0.2*randFrom0to1(); // x1 -> p1
mVar[mNumTimes+0] = 1.0 + 0.2 * randFrom0to1(); // x0 -> p0
mVar[mNumTimes+1] = 0.0 + 0.2 * randFrom0to1(); // x1 -> p1
#else
mVar[mNumTimes+0] = randFrom0to1(); // p0+p1
mVar[mNumTimes+1] = randFrom0to1(); // p0/(p0+p1)
mVar[mNumTimes+0] = randFrom0to1(); // p0+p1
mVar[mNumTimes+1] = randFrom0to1(); // p0/(p0+p1)
#endif
mVar[mNumTimes+2] = randFrom0to1()*0.8 + 0.1; // w0
mVar[mNumTimes+3] = 2.0; // k
mVar[mNumTimes+2] = 0.05 + 0.8 * randFrom0to1(); // w0
mVar[mNumTimes+3] = 0.2 + 2.1 * randFrom0to1(); // k
}
// Initialize w2 if needed
if(mNumVariables == 5 && mInitType != INIT_TYPE_RES_5)
if(mNumVariables == 5 && (mInitStatus & INIT_PARAM_W2) != INIT_PARAM_W2)
{
mVar[mNumTimes+4] = 1.001 + 0.149 * randFrom0to1(); // w2
mVar[mNumTimes+4] = 1.001 + 0.149 * randFrom0to1(); // w2
}
// Re-initialize the next time
mInitType = INIT_TYPE_NONE;
mInitStatus = INIT_NONE;
// Check the initial values to be inside the domain (otherwise clamp them to the domain)
unsigned int nv = mNumTimes+mNumVariables;
......
......@@ -62,7 +62,7 @@ protected:
mOnlyInitialStep(aOnlyInitialStep),
mTrace(aTrace),
mOptAlgo(aOptAlgo),
mInitType(INIT_TYPE_NONE),
mInitStatus(INIT_NONE),
mNumTimes(static_cast<unsigned int>(aNumBranches)),
mNumVariables(aNumVariables),
mExtraDebug(aExtraDebug),
......@@ -261,27 +261,28 @@ private:
OPTIM_MLSL_LDS = 99 ///< A global optimizer
};
/// Valid values for the mInitType variable depicting from where the variables have been initialized.
/// Valid values for the mInitStatus variable depicting from where the variables have been initialized.
///
enum InitVarStatus
{
INIT_TYPE_NONE, ///< All variables to optimize should be initialized
INIT_TYPE_TIMES, ///< All variables to optimize should be initialized except times
INIT_TYPE_RES_4, ///< All variables to optimize should be initialized except times, w0, k, v1, v2
INIT_TYPE_RES_5 ///< All variables to optimize have been initialized
INIT_NONE=0, ///< All variables to optimize should be initialized
INIT_TIMES=1, ///< All variables to optimize should be initialized except times
INIT_PARAMS_H0=2,
INIT_PARAM_W2=4,
INIT_PARAMS=6 ///< All variables to optimize should be initialized except w0, w2, k, p0, p1
};
public:
/// Initialize the times from the input phylogenetic tree
/// Initialize the times from the input phylogenetic tree.
///
void initFromTree(void);
/// Initialize the times from the input phylogenetic tree and set the other values to hardcoded constants with P0=1 and P1=0.
/// Set the other parameters values to hardcoded constants.
/// This routine could be used in place of initFromTree()
///
/// @exception FastCodeMLFatal Invalid p0 and p1 values
///
void initFromTreeAndParams(void);
void initFromParams(void);
/// Initialize variables from a previous optimization result
///
......@@ -317,7 +318,7 @@ protected:
bool mOnlyInitialStep; ///< Only the initial step is executed, no optimization
bool mTrace; ///< Enable maximization tracing
unsigned int mOptAlgo; ///< Optimization algorithm to use
InitVarStatus mInitType; ///< From where the variables have been initialized
unsigned int mInitStatus; ///< Which variables have been initialized
unsigned int mNumTimes; ///< Number of branch lengths values
unsigned int mNumVariables; ///< The number of extra variables (4 for H0 and 5 for H1)
unsigned int mExtraDebug; ///< Parameter for extra development testing
......
......@@ -332,6 +332,7 @@ void CmdLine::parseCmdLine(int aCnt, char **aVal)
break;
case OPT_INIT_DEFAULT:
mBranchLengthsFromFile = true;
mInitFromParams = true;
break;
......
......@@ -103,12 +103,13 @@ struct HighLevelCoordinator::WorkTable
///
/// @param[out] aJob The job request: [0] is set to the kind of job (JOB_H0, JOB_H1, JOB_BEB, JOB_SHUTDOWN);
/// [1] to the fg branch number (or zero for JOB_SHUTDOWN);
/// [2] zero or the number of variables for a JOB_BEB request
/// [2] zero or the number of variables for a JOB_BEB or JOB_H0 requests
/// @param[in] aRank The current worker rank
/// @param[in] aInitFromH1 If true for a H0 request send back also all H1 variables
///
/// @return True if a job has been assigned, false if the job is JOB_SHUTDOWN
///
bool getNextJob(int* aJob, int aRank);
bool getNextJob(int* aJob, int aRank, bool aInitFromH1=false);
/// Mark the job assigned to the aRank worker as finished
///
......@@ -150,7 +151,7 @@ struct HighLevelCoordinator::WorkTable
};
bool HighLevelCoordinator::WorkTable::getNextJob(int* aJob, int aRank)
bool HighLevelCoordinator::WorkTable::getNextJob(int* aJob, int aRank, bool aInitFromH1)
{
// Assign all H1 jobs
for(size_t branch=0; branch < mNumInternalBranches; ++branch)
......@@ -175,7 +176,7 @@ bool HighLevelCoordinator::WorkTable::getNextJob(int* aJob, int aRank)
{
aJob[0] = JOB_H0;
aJob[1] = static_cast<int>(branch);
aJob[2] = 1; // Send the lnl of the corresponding H1 step
aJob[2] = (aInitFromH1) ? static_cast<int>(mResults[branch].mHxVariables[1].size())+1 : 1; // Send the lnl of the corresponding H1 step or mResults[branch].mHxVariables[1]
mJobStatus[idx] = JOB_ASSIGNED;
mWorkList[idx] = aRank;
return true;
......@@ -448,7 +449,7 @@ bool HighLevelCoordinator::startWork(Forest& aForest, const CmdLine& aCmdLine)
WriteResults output_results(aCmdLine.mResultsFile);
// In the master process initialize the master
doMaster(output_results);
doMaster(output_results, aCmdLine);
}
else
{
......@@ -461,7 +462,7 @@ bool HighLevelCoordinator::startWork(Forest& aForest, const CmdLine& aCmdLine)
}
void HighLevelCoordinator::doMaster(WriteResults& aOutputResults)
void HighLevelCoordinator::doMaster(WriteResults& aOutputResults, const CmdLine& aCmdLine)
{
// Push work to free workers
unsigned int num_workers = 0;
......@@ -570,7 +571,7 @@ void HighLevelCoordinator::doMaster(WriteResults& aOutputResults)
// Send work packet or shutdown request (job[0] is the step to be done, job[1] is the fg branch, job[2] the length of the additional data)
int job[3];
mWorkTable->getNextJob(job, worker);
mWorkTable->getNextJob(job, worker, aCmdLine.mInitH0fromH1);
MPI_Send(static_cast<void*>(job), 3, MPI_INTEGER, worker, MSG_NEW_JOB, MPI_COMM_WORLD);
// For BEB send the variables from H1; for H0 send the lnl value from corresponding H1
......@@ -585,8 +586,17 @@ void HighLevelCoordinator::doMaster(WriteResults& aOutputResults)
break;
case JOB_H0:
v = &(mWorkTable->mResults[job[1]].mLnl[1]);
MPI_Send(static_cast<void*>(v), 1, MPI_DOUBLE, worker, MSG_NEW_JOB, MPI_COMM_WORLD);
if(job[2] == 1)
{
v = &(mWorkTable->mResults[job[1]].mLnl[1]);
}
else
{
results_double.assign(mWorkTable->mResults[job[1]].mHxVariables[1].begin(), mWorkTable->mResults[job[1]].mHxVariables[1].end());
results_double.push_back(mWorkTable->mResults[job[1]].mLnl[1]);
v = &results_double[0];
}
MPI_Send(static_cast<void*>(v), job[2], MPI_DOUBLE, worker, MSG_NEW_JOB, MPI_COMM_WORLD);
break;
}
}
......@@ -774,11 +784,15 @@ void HighLevelCoordinator::doWorker(Forest& aForest, const CmdLine& aCmdLine)
case JOB_H0:
{
// Initialize maximizer
if(aCmdLine.mInitFromParams) h0.initFromTreeAndParams();
else if(aCmdLine.mBranchLengthsFromFile) h0.initFromTree();
if(aCmdLine.mInitH0fromH1 && job[2] > 1) h0.initFromResult(values_double, static_cast<unsigned int>(values_double.size())-1);
else
{
if(aCmdLine.mInitFromParams) h0.initFromParams();
if(aCmdLine.mBranchLengthsFromFile) h0.initFromTree();
}
// Get the lnl from the corresponding H1 step if any
double threshold = (job[2] > 0) ? values_double[0]-THRESHOLD_FOR_LRT : 0.;
double threshold = (job[2] > 0) ? values_double.back()-THRESHOLD_FOR_LRT : 0.;
// Compute H0
double lnl = h0(static_cast<size_t>(job[1]), aCmdLine.mStopIfNotLRT && job[2] > 0, threshold);
......@@ -794,8 +808,8 @@ void HighLevelCoordinator::doWorker(Forest& aForest, const CmdLine& aCmdLine)
case JOB_H1:
{
// Initialize maximizer
if(aCmdLine.mInitFromParams) h1.initFromTreeAndParams();
else if(aCmdLine.mBranchLengthsFromFile) h1.initFromTree();
if(aCmdLine.mInitFromParams) h1.initFromParams();
if(aCmdLine.mBranchLengthsFromFile) h1.initFromTree();
// Compute H1
double lnl = h1(static_cast<size_t>(job[1]));
......
......@@ -61,10 +61,11 @@ private:
/// The master coordination job
///
/// @param[in] aOutputResults To collect and output results to a results file
/// @param[in] aCmdLine The parameters from the command line of the main program
///
/// @exception FastCodeMLFatal Invalid job request found
///
void doMaster(WriteResults& aOutputResults);
void doMaster(WriteResults& aOutputResults, const CmdLine& aCmdLine);
/// The worker high level loop
///
......
......@@ -99,12 +99,15 @@ int main(int aRgc, char **aRgv)
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 std::cout << "Branches: " << cmd.mBranchStart << '-' << cmd.mBranchEnd << std::endl;
else if(cmd.mBranchStart != UINT_MAX && cmd.mBranchEnd != UINT_MAX)
std::cout << "Branches: " << cmd.mBranchStart << '-' << cmd.mBranchEnd << 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) std::cout << "Starting val.: Times from tree file and params from const (see below)" << 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;
......@@ -302,8 +305,8 @@ int main(int aRgc, char **aRgv)
double lnl1 = 0.;
if(cmd.mComputeHypothesis != 0)
{
if(cmd.mInitFromParams) h1.initFromTreeAndParams();
else if(cmd.mBranchLengthsFromFile) h1.initFromTree();
if(cmd.mInitFromParams) h1.initFromParams();
if(cmd.mBranchLengthsFromFile) h1.initFromTree();
lnl1 = h1(fg_branch);
......@@ -316,8 +319,11 @@ int main(int aRgc, char **aRgv)
if(cmd.mComputeHypothesis != 1)
{
if(cmd.mInitH0fromH1) h0.initFromResult(h1.getVariables());
else if(cmd.mInitFromParams) h0.initFromTreeAndParams();
else if(cmd.mBranchLengthsFromFile) h0.initFromTree();
else
{
if(cmd.mInitFromParams) h0.initFromParams();
if(cmd.mBranchLengthsFromFile) h0.initFromTree();
}
lnl0 = h0(fg_branch, cmd.mStopIfNotLRT && cmd.mComputeHypothesis != 0, lnl1-THRESHOLD_FOR_LRT);
......
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