Commit 82d4e2bc authored by mvalle's avatar mvalle
Browse files

Print annotated tree with the internal branch number.

Changed name to the internal variable mTimesFromFile to mBranchLengthFromFile.


git-svn-id: https://svn.vital-it.ch/svn/hp2c/trunk/Codeml_Variants/Fastcodeml@5567 95c5a404-1f88-0410-a6b0-c3f062b6f34a
parent 3b7e20cb
......@@ -266,7 +266,7 @@ void CmdLine::parseCmdLine(int aCnt, char **aVal)
break;
case OPT_TIMES_FROM_FILE:
mTimesFromFile = true;
mBranchLengthsFromFile = true;
break;
case OPT_ONE_STEP:
......
......@@ -35,7 +35,7 @@ public:
mExtraDebug(0),
mIgnoreFreq(false),
mDoNotReduceForest(false),
mTimesFromFile(false),
mBranchLengthsFromFile(false),
mNoMaximization(false),
mTrace(false),
mForceSerial(false),
......@@ -78,7 +78,7 @@ public:
unsigned int mExtraDebug; ///< Extra debug parameter for development tests
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 mTimesFromFile; ///< The initial value of the branch lengths is taken from the phylo tree file
bool mBranchLengthsFromFile; ///< The initial value of the branch lengths is taken from the phylo tree file
bool mNoMaximization; ///< Only the first step of the likelihood maximization is taken
bool mTrace; ///< Trace the optimization steps
bool mForceSerial; ///< Disable all parallelism
......
......@@ -656,8 +656,8 @@ void HighLevelCoordinator::doWorker(Forest& aForest, const CmdLine& aCmdLine)
case JOB_H0:
{
// Initialize maximizer
if(aCmdLine.mInitFromParams) h0.initFromTreeAndParams();
else if(aCmdLine.mTimesFromFile) h0.initFromTree();
if(aCmdLine.mInitFromParams) h0.initFromTreeAndParams();
else 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.;
......@@ -676,8 +676,8 @@ void HighLevelCoordinator::doWorker(Forest& aForest, const CmdLine& aCmdLine)
case JOB_H1:
{
// Initialize maximizer
if(aCmdLine.mInitFromParams) h1.initFromTreeAndParams();
else if(aCmdLine.mTimesFromFile) h1.initFromTree();
if(aCmdLine.mInitFromParams) h1.initFromTreeAndParams();
else if(aCmdLine.mBranchLengthsFromFile) h1.initFromTree();
// Compute H1
double lnl = h1(static_cast<size_t>(job[1]));
......
......@@ -169,8 +169,6 @@ void Newick::readFile(const char *aFilename)
{
p1 = str.find_first_of('(');
if(p1 != std::string::npos) break;
// p1 = str.find_first_not_of(" \t\r\n");
// if(p1 != std::string::npos && str[p1] == '(') break;
}
in.close();
if(p1 == std::string::npos)
......@@ -276,3 +274,44 @@ void Newick::printTreeUnformatted(std::ostream& aOut, TreeNode *aNode) const
}
}
int Newick::printTreeAnnotated(std::ostream& aOut, TreeNode *aNode, int aInternalBranch) const
{
TreeNode *m;
unsigned int idx;
int internal_branch_idx = aInternalBranch;
// Special case for the root
if(!aNode)
{
aOut << std::endl << "Annotated Newick tree (*N mark the internal branch N)" << std::endl;
aOut << '(';
for(idx=0; (m = mTreeRoot.getChild(idx)) != NULL; ++idx)
{
if(idx > 0) aOut << ',';
internal_branch_idx = printTreeAnnotated(aOut, m, internal_branch_idx);
}
aOut << ')';
mTreeRoot.printNode();
aOut << ';' << std::endl;
}
else if(aNode->isLeaf())
{
aNode->printNode();
}
else
{
internal_branch_idx = aInternalBranch+1;
aOut << '(';
for(idx=0; (m = aNode->getChild(idx)) != NULL; ++idx)
{
if(idx > 0) aOut << ',';
internal_branch_idx = printTreeAnnotated(aOut, m, internal_branch_idx);
}
aOut << ')';
aNode->printNode();
aOut << '*' << aInternalBranch;
}
return internal_branch_idx;
}
......@@ -50,6 +50,16 @@ public:
///
virtual void printTreeUnformatted(std::ostream& aOut, TreeNode *aNode=NULL) const;
/// Print the phylogenetic tree completed with all the info loaded in the same format as read in and annotated with the internal branch number.
///
/// @param[in] aOut Output stream
/// @param[in] aNode The node from which to start. If null starts from the root.
/// @param[in] aInternalBranch Internal branch identifier to annotate the current branch.
///
/// @return The new internal branch id
///
virtual int printTreeAnnotated(std::ostream& aOut, TreeNode *aNode=NULL, int aInternalBranch=0) const;
private:
/// Load a phylo tree definition from a Newick formatted string.
......
......@@ -46,6 +46,16 @@ public:
///
virtual void printTreeUnformatted(std::ostream& aOut, TreeNode *aNode=NULL) const =0;
/// Print the phylogenetic tree completed with all the info loaded in the same format as read in and annotated with the internal branch number.
///
/// @param[in] aOut Output stream
/// @param[in] aNode The node from which to start. If null starts from the root.
/// @param[in] aInternalBranch Internal branch identifier to annotate the current branch.
///
/// @return The new internal branch id
///
virtual int printTreeAnnotated(std::ostream& aOut, TreeNode *aNode=NULL, int aInternalBranch=0) const =0;
/// Return the list of species.
///
/// @return Reference to a vector of species names as read from the leaves of the tree
......
......@@ -101,7 +101,7 @@ int main(int aRgc, char **aRgv)
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.mTimesFromFile) std::cout << "Starting val.: Times from tree file" << 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;
......@@ -188,7 +188,10 @@ int main(int aRgc, char **aRgv)
msa.checkNameCoherence(tree.getSpecies());
// If times from file then check for null branch lengths for any leaf
if(cmd.mTimesFromFile) tree.checkNullBranchLengths();
if(cmd.mBranchLengthsFromFile) tree.checkNullBranchLengths();
// 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);
......@@ -308,8 +311,8 @@ int main(int aRgc, char **aRgv)
double lnl1 = 0.;
if(cmd.mComputeHypothesis != 0)
{
if(cmd.mInitFromParams) h1.initFromTreeAndParams();
else if(cmd.mTimesFromFile) h1.initFromTree();
if(cmd.mInitFromParams) h1.initFromTreeAndParams();
else if(cmd.mBranchLengthsFromFile) h1.initFromTree();
lnl1 = h1(fg_branch);
......@@ -321,9 +324,9 @@ int main(int aRgc, char **aRgv)
double lnl0 = 0.;
if(cmd.mComputeHypothesis != 1)
{
if(cmd.mInitH0fromH1) h0.initFromResult(h1.getVariables());
else if(cmd.mInitFromParams) h0.initFromTreeAndParams();
else if(cmd.mTimesFromFile) h0.initFromTree();
if(cmd.mInitH0fromH1) h0.initFromResult(h1.getVariables());
else if(cmd.mInitFromParams) h0.initFromTreeAndParams();
else 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