Commit 22133b10 authored by mvalle's avatar mvalle
Browse files

Added guard against infinites in results

Added print of no pre stop switch
Fixed global scaling code


git-svn-id: https://svn.vital-it.ch/svn/hp2c/trunk/Codeml_Variants/Fastcodeml@5861 95c5a404-1f88-0410-a6b0-c3f062b6f34a
parent 177d77f5
......@@ -510,9 +510,6 @@ double BranchSiteModelNullHyp::combineSiteLikelihoods(void)
const std::vector<double>& mult = mForest.getSiteMultiplicity();
double lnl = 0;
double* likelihoods = &mLikelihoods[0];
#ifdef USE_GLOBAL_SCALING
double scale = 0;
#endif
for(size_t site=0; site < num_sites; ++site)
{
// The following computation is split to avoid negative values
......@@ -527,12 +524,9 @@ double BranchSiteModelNullHyp::combineSiteLikelihoods(void)
x = likelihoods[2*num_sites+site];
if(x > 0) p += p2a*x;
//x = (p > 0) ? log(p)-(mForest.getNumBranches()-mForest.getNumInternalBranches())*log(GLOBAL_SCALING_FACTOR) : mMaxLnL-100000;
x = (p > 0) ? log(p) : mMaxLnL-100000;
lnl += x*mult[site];
#ifdef USE_GLOBAL_SCALING
scale += mult[site]*(mForest.getNumBranches()-mForest.getNumInternalBranches());
#endif
if(mExtraDebug > 1)
{
std::cout << std::setw(4) << site << ' ';
......@@ -543,7 +537,7 @@ double BranchSiteModelNullHyp::combineSiteLikelihoods(void)
}
}
#ifdef USE_GLOBAL_SCALING
lnl -= scale*log(GLOBAL_SCALING_FACTOR);
lnl -= mGlobalScale;
#endif
return lnl;
......@@ -846,9 +840,6 @@ double BranchSiteModelAltHyp::combineSiteLikelihoods(void)
const std::vector<double>& mult = mForest.getSiteMultiplicity();
double lnl = 0;
double* likelihoods = &mLikelihoods[0];
#ifdef USE_GLOBAL_SCALING
double scale = 0;
#endif
for(size_t site=0; site < num_sites; ++site)
{
// The following computation is split to avoid negative values
......@@ -867,13 +858,9 @@ double BranchSiteModelAltHyp::combineSiteLikelihoods(void)
x = likelihoods[3*num_sites+site];
if(x > 0) p += p2b*x;
//x = (p > 0) ? log(p)-(mForest.getNumBranches()-mForest.getNumInternalBranches())*log(GLOBAL_SCALING_FACTOR) : mMaxLnL-100000;
//lnl += x*mult[site];
x = (p > 0) ? log(p) : mMaxLnL-100000;
lnl += x*mult[site];
#ifdef USE_GLOBAL_SCALING
scale += mult[site]*(mForest.getNumBranches()-mForest.getNumInternalBranches());
#endif
if(mExtraDebug > 1)
{
std::cout << std::setw(4) << site << ' ';
......@@ -885,7 +872,7 @@ double BranchSiteModelAltHyp::combineSiteLikelihoods(void)
}
}
#ifdef USE_GLOBAL_SCALING
lnl -= scale*log(GLOBAL_SCALING_FACTOR);
lnl -= mGlobalScale;
#endif
return lnl;
......@@ -957,7 +944,6 @@ public:
private:
#if 1
/// Compute the function gradient
///
/// @param[in] aPointValue The value of the function at aVars
......@@ -1016,7 +1002,8 @@ private:
vars_working_copy[i] = aVars[i];
}
}
#else
#if 0
/// Compute the function gradient (original method)
///
/// @param[in] aPointValue The value of the function at aVars
......
......@@ -6,8 +6,8 @@
#include <cstdlib>
#include <cfloat>
#include "TransitionMatrix.h"
#include "ProbabilityMatrixSet.h"
#include "Forest.h"
#include "ProbabilityMatrixSet.h"
#include "CmdLine.h"
#include "TreeAndSetsDependencies.h"
......@@ -74,6 +74,14 @@ protected:
mRelativeError(aRelativeError)
{
setLimits(aNumBranches, static_cast<size_t>(aNumVariables));
#ifdef USE_GLOBAL_SCALING
const size_t num_sites = mForest.getNumSites();
const std::vector<double>& mult = mForest.getSiteMultiplicity();
double total_mult = 0.;
for(size_t site=0; site < num_sites; ++site) total_mult += mult[site];
mGlobalScale = (mForest.getNumBranches()-mForest.getNumInternalBranches())*log(GLOBAL_SCALING_FACTOR)*total_mult;
#endif
}
/// Destructor.
......@@ -314,6 +322,9 @@ protected:
double mBgScale; ///< The computed background branches scale
double mMaxLnL; ///< Maximum value of LnL found during optimization
double mDeltaForGradient; ///< Value used to change the variables to compute gradient
#ifdef USE_GLOBAL_SCALING
double mGlobalScale;
#endif
CacheAlignedDoubleVector mLikelihoods; ///< Computed likelihoods at the root of all trees. Defined here to make it aligned.
bool mOnlyInitialStep; ///< Only the initial step is executed, no optimization
bool mTrace; ///< Enable maximization tracing
......
......@@ -23,11 +23,6 @@
#include "TreeAndSetsDependencies.h"
#include "CmdLine.h"
// Global scaling factor (uncomment the define only if the scaling factor is different from one)
//#define USE_GLOBAL_SCALING
#ifdef USE_GLOBAL_SCALING
static const double GLOBAL_SCALING_FACTOR = 1.0;
#endif
/// The phylogenetic tree's forest.
/// This class encapsulates the forest of phylogenetic tree that will be used for computing the tree's maximum likelihood
......
......@@ -12,6 +12,7 @@
#include <iostream>
#include <iomanip>
#include <limits>
#ifndef VTRACE
#ifdef _OPENMP
#include <omp.h>
......@@ -648,7 +649,11 @@ void HighLevelCoordinator::doMaster(WriteResults& aOutputResults, const CmdLine&
if(mWorkTable->mJobStatus[branch*JOBS_PER_BRANCH+JOB_H1] == JOB_SKIP) continue;
std::cout << "Branch: " << std::fixed << std::setw(3) << branch;
if(mWorkTable->mResults[branch].mLnl[0] == DBL_MAX)
if(mWorkTable->mResults[branch].mLnl[0] == std::numeric_limits<double>::infinity())
{
std::cout << " Lnl H0: " << std::setw(24) << "Inf";
}
else if(mWorkTable->mResults[branch].mLnl[0] == DBL_MAX)
{
std::cout << " Lnl H0: " << std::setw(24) << "NA";
}
......@@ -656,16 +661,26 @@ void HighLevelCoordinator::doMaster(WriteResults& aOutputResults, const CmdLine&
{
std::cout << " Lnl H0: " << std::setw(24) << std::setprecision(15) << mWorkTable->mResults[branch].mLnl[0];
}
std::cout << " Lnl H1: " << std::setw(24) << std::setprecision(15) << mWorkTable->mResults[branch].mLnl[1];
if(mWorkTable->mResults[branch].mLnl[0] < DBL_MAX)
if(mWorkTable->mResults[branch].mLnl[1] == std::numeric_limits<double>::infinity())
{
std::cout << " Lnl H1: " << std::setw(24) << "Inf";
}
else
{
std::cout << " Lnl H1: " << std::setw(24) << std::setprecision(15) << mWorkTable->mResults[branch].mLnl[1];
}
if(mWorkTable->mResults[branch].mLnl[0] == std::numeric_limits<double>::infinity() || mWorkTable->mResults[branch].mLnl[1] == std::numeric_limits<double>::infinity())
{
std::cout << " LRT: " << std::setw(24) << "*Invalid*";
}
else if(mWorkTable->mResults[branch].mLnl[0] < DBL_MAX)
std::cout << " LRT: " << std::setw(24) << std::setprecision(15) << std::fixed << mWorkTable->mResults[branch].mLnl[1] - mWorkTable->mResults[branch].mLnl[0] << " (threshold: " << std::setprecision(15) << std::fixed << THRESHOLD_FOR_LRT << ')';
else
std::cout << " LRT: < " << std::setprecision(15) << std::fixed << THRESHOLD_FOR_LRT;
std::cout << std::endl;
std::cout << std::endl;
if(mWorkTable->mResults[branch].mLnl[0] != DBL_MAX) mWorkTable->printVariables(branch, 0, std::cout);
mWorkTable->printVariables(branch, 1, std::cout);
if(mWorkTable->mResults[branch].mLnl[0] != std::numeric_limits<double>::infinity() && mWorkTable->mResults[branch].mLnl[0] != DBL_MAX) mWorkTable->printVariables(branch, 0, std::cout);
if(mWorkTable->mResults[branch].mLnl[1] != std::numeric_limits<double>::infinity()) mWorkTable->printVariables(branch, 1, std::cout);
std::cout << std::endl;
}
......
......@@ -17,6 +17,12 @@
#include <mkl_vml_functions.h>
#endif
// Global scaling factor (uncomment the define only if the scaling factor is different from one)
#define USE_GLOBAL_SCALING
#ifdef USE_GLOBAL_SCALING
static const double GLOBAL_SCALING_FACTOR = 4.25;
#endif
/// Set of probability matrices for all branches of a tree.
///
/// @author Mario Valle - Swiss National Supercomputing Centre (CSCS)
......@@ -126,6 +132,10 @@ public:
break;
}
#ifdef USE_GLOBAL_SCALING
for(int i=0; i < N; ++i) aGout[i] *= GLOBAL_SCALING_FACTOR;
#endif
// The element wise multiplication has been moved up one level
//#if !defined(BUNDLE_ELEMENT_WISE_MULT)
// elementWiseMult(aGout, mInvCodonFreq);
......@@ -133,7 +143,11 @@ public:
#else
for(int r=0; r < N; ++r)
{
#ifdef USE_GLOBAL_SCALING
aGout[r] = mMatrices[aSetIdx*mNumMatrices+aBranch][r*N+aCodon]*GLOBAL_SCALING_FACTOR;
#else
aGout[r] = mMatrices[aSetIdx*mNumMatrices+aBranch][r*N+aCodon];
#endif
}
#endif
}
......
......@@ -16,6 +16,7 @@
#include <iostream>
#include <iomanip>
#include <limits>
#include "CmdLine.h"
#include "Newick.h"
#include "Phylip.h"
......@@ -101,6 +102,7 @@ int main(int aRgc, char **aRgv)
std::cout << "Branches: " << cmd.mBranchStart << "-end" << std::endl;
else if(cmd.mBranchStart != UINT_MAX && cmd.mBranchEnd != UINT_MAX)
std::cout << "Branches: " << cmd.mBranchStart << '-' << cmd.mBranchEnd << std::endl;
if(!cmd.mStopIfNotLRT) std::cout << "H0 pre stop: No" << 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;
......@@ -211,7 +213,7 @@ int main(int aRgc, char **aRgv)
if(zero_on_leaf_cnt > 0 || zero_on_int_cnt > 0)
{
std::cout << "Found Null or missing branch length in tree file. On leaves: " << zero_on_leaf_cnt << " on internal branches: " << zero_on_int_cnt << std::endl;
std::cout << "Found null or missing branch length in tree file. On leaves: " << zero_on_leaf_cnt << " on internal branches: " << zero_on_int_cnt << std::endl;
}
if(zero_on_leaf_cnt > 0)
......@@ -342,27 +344,34 @@ int main(int aRgc, char **aRgv)
if(cmd.mComputeHypothesis != 1)
{
std::cout << "LnL0: ";
if(lnl0 < DBL_MAX)
if(lnl0 == std::numeric_limits<double>::infinity())
std::cout << "**Invalid result**";
else if(lnl0 < DBL_MAX)
std::cout << std::setprecision(15) << std::fixed << lnl0;
else
std::cout << "(Doesn't pass LRT, skipping)";
std::cout << " Function calls: " << h0.getNumEvaluations() << " ";
std::cout << std::endl << std::endl;
h0.printFinalVars(std::cout);
if(lnl0 != std::numeric_limits<double>::infinity()) h0.printFinalVars(std::cout);
std::cout << std::endl;
}
if(cmd.mComputeHypothesis != 0)
{
std::cout << "LnL1: ";
std::cout << std::setprecision(15) << std::fixed << lnl1;
if(lnl1 == std::numeric_limits<double>::infinity())
std::cout << "**Invalid result**";
else
std::cout << std::setprecision(15) << std::fixed << lnl1;
std::cout << " Function calls: " << h1.getNumEvaluations();
std::cout << std::endl << std::endl;
h1.printFinalVars(std::cout);
if(lnl1 != std::numeric_limits<double>::infinity()) h1.printFinalVars(std::cout);
std::cout << std::endl;
}
if(cmd.mComputeHypothesis > 1)
{
if(lnl0 < DBL_MAX)
if(lnl0 == std::numeric_limits<double>::infinity() || lnl1 == std::numeric_limits<double>::infinity())
std::cout << "LRT: **Invalid result**";
else if(lnl0 < DBL_MAX)
std::cout << "LRT: " << std::setprecision(15) << std::fixed << lnl1 - lnl0 << " (threshold: " << std::setprecision(15) << std::fixed << THRESHOLD_FOR_LRT << ')';
else
std::cout << "LRT: < " << std::setprecision(15) << std::fixed << 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