Commit 6377bc04 authored by omid's avatar omid
Browse files

merging master-eval up to WriteResults.h

parent 4cc63271
......@@ -33,8 +33,8 @@ void Forest::loadTreeAndGenes(
CodonFrequencies::CodonFrequencyModel aCodonFrequencyModel) {
// Collect global data that refers to the tree and that should not be
// duplicated on each tree of the forest
aTree.collectGlobalTreeData(mNodeNames, mBranchLengths,
&mMarkedInternalBranch);
aTree.collectGlobalTreeData(mNodeNames, mBranchLengths, mInternalBranches,
&mMarkedInternalBranch, mMarkedBranches);
// Number of branches of one tree
mNumBranches = aTree.getNumBranches();
......@@ -133,8 +133,9 @@ void Forest::loadTreeAndGenes(
map_internal_to_branchID.begin());
const std::map<unsigned int, unsigned int>::const_iterator end(
map_internal_to_branchID.end());
for (; im != end; ++im)
for (; im != end; ++im) {
mTableInternalToBranchID[im->first] = im->second;
}
// Save the new to original site number map
mSitesMappingToOriginal = aGenes.getSitesMappingToOriginal();
......@@ -271,69 +272,85 @@ void Forest::postLoad(void) {
#endif
bool Forest::getBranchRange(const CmdLine &aCmdLine, size_t &aBranchStart,
size_t &aBranchEnd) const {
const size_t num_branches = getNumInternalBranches();
const size_t marked_branch = getMarkedInternalBranch();
size_t &aBranchEnd, std::set<int> &aFgSet,
std::set<int> &aIbSet) const {
const size_t num_branches = getNumBranches(); /* getNumInternalBranches(); */
// const size_t num_internal_branches = getNumInternalBranches();
aFgSet = getMarkedBranches();
aIbSet = getInternalBranches();
// Check if the request make sense
if (num_branches == 0) {
throw FastCodeMLFatal("No internal branches present. Quitting.");
throw FastCodeMLFatal("No branches present. Quitting.");
}
// By default do all branches
bool do_all = true;
// Adjust the number of branches to compute
if (aCmdLine.mBranchFromFile) {
// Branch from file, verify if valid
if (marked_branch >= num_branches) {
if (aCmdLine.mVerboseLevel >= VERBOSE_INFO_OUTPUT)
std::cout << std::endl
<< "Invalid branch marked in tree file. Ignoring"
<< std::endl;
aBranchStart = 0;
aBranchEnd = num_branches - 1;
} else {
aBranchStart = marked_branch;
aBranchEnd = marked_branch;
do_all = false;
}
} else if (aCmdLine.mBranchStart < UINT_MAX &&
aCmdLine.mBranchStart >= num_branches) {
// Invalid start value, ignoring, do all branches
if (aCmdLine.mVerboseLevel >= VERBOSE_INFO_OUTPUT)
std::cout << std::endl
<< "Invalid branch requested. Ignoring" << std::endl;
aBranchStart = 0;
aBranchEnd = num_branches - 1;
} else if (aCmdLine.mBranchStart < UINT_MAX &&
aCmdLine.mBranchEnd == UINT_MAX) {
// Only start branch set. Do from it to the end.
aBranchStart = static_cast<size_t>(aCmdLine.mBranchStart);
aBranchEnd = num_branches - 1;
if (aBranchStart > 0)
do_all = false;
} else if (aCmdLine.mBranchStart < UINT_MAX &&
aCmdLine.mBranchEnd < UINT_MAX) {
// Both start and end branch (already tested start <= end)
aBranchStart = static_cast<size_t>(aCmdLine.mBranchStart);
if (aCmdLine.mBranchEnd >= num_branches) {
if (aCmdLine.mVerboseLevel >= VERBOSE_INFO_OUTPUT)
std::cout << std::endl
<< "Invalid end branch requested. Ignoring" << std::endl;
aBranchEnd = num_branches - 1;
if (aBranchStart > 0)
do_all = false;
} else {
aBranchEnd = static_cast<size_t>(aCmdLine.mBranchEnd);
if (aBranchStart > 0 && aBranchEnd < num_branches - 1)
do_all = false;
}
} else {
// No limit set, do all branches
aBranchStart = 0;
aBranchEnd = num_branches - 1;
}
/*if(aCmdLine.mBranchFromFile)
{
// Branch from file, verify if valid
if(marked_branch >= num_branches)
{
if(aCmdLine.mVerboseLevel >= VERBOSE_INFO_OUTPUT) std::cout << std::endl <<
"Invalid branch marked in tree file. Ignoring" << std::endl;
aBranchStart = 0;
aBranchEnd = num_branches-1;
}
else
{
aBranchStart = marked_branch;
aBranchEnd = marked_branch;
do_all = false;
}
}*/
/*else if(aCmdLine.mBranchStart < UINT_MAX && aCmdLine.mBranchStart >=
num_branches)
{
// Invalid start value, ignoring, do all branches
if(aCmdLine.mVerboseLevel >= VERBOSE_INFO_OUTPUT) std::cout << std::endl <<
"Invalid branch requested. Ignoring" << std::endl;
aBranchStart = 0;
aBranchEnd = num_branches-1;
}
else if(aCmdLine.mBranchStart < UINT_MAX && aCmdLine.mBranchEnd == UINT_MAX)
{
// Only start branch set. Do from it to the end.
aBranchStart = static_cast<size_t>(aCmdLine.mBranchStart);
aBranchEnd = num_branches-1;
if(aBranchStart > 0) do_all = false;
}
else if(aCmdLine.mBranchStart < UINT_MAX && aCmdLine.mBranchEnd < UINT_MAX)
{
// Both start and end branch (already tested start <= end)
aBranchStart = static_cast<size_t>(aCmdLine.mBranchStart);
if(aCmdLine.mBranchEnd >= num_branches)
{
if(aCmdLine.mVerboseLevel >= VERBOSE_INFO_OUTPUT) std::cout << std::endl <<
"Invalid end branch requested. Ignoring" << std::endl;
aBranchEnd = num_branches-1;
if(aBranchStart > 0) do_all = false;
}
else
{
aBranchEnd = static_cast<size_t>(aCmdLine.mBranchEnd);
if(aBranchStart > 0 && aBranchEnd < num_branches-1) do_all = false;
}
}*/
// if (aCmdLine.mBranchAll)
//{
// No limit set, do all branches
aBranchStart = 0;
aBranchEnd = num_branches - 1;
//}
// else
//{
// default, do all internal branches
// aBranchStart = 0;
// aBranchEnd = num_internal_branches-1;
//}
return do_all;
}
......@@ -414,13 +431,13 @@ void Forest::cleanReductionWorkingData(ForestNode *aNode) {
std::ostream &operator<<(std::ostream &aOut, const Forest &aForest) {
// General forest statistics
aOut << std::endl;
aOut << "Num branches: " << std::setw(7) << aForest.mNumBranches
<< std::endl;
aOut << "Internal branches: " << std::setw(7) << aForest.mNumInternalBranches
aOut << "Num branches: " << std::setw(7)
<< aForest.mNumBranches << std::endl;
aOut << "Internal branches: " << std::setw(7)
<< aForest.mNumInternalBranches << std::endl;
aOut << "Unique sites: " << std::setw(7) << aForest.mNumSites
<< std::endl;
aOut << "Unique sites: " << std::setw(7) << aForest.mNumSites
<< std::endl;
aOut << "Total branches: " << std::setw(7)
aOut << "Total branches: " << std::setw(7)
<< aForest.mNumBranches * aForest.mNumSites << std::endl;
// Count total branches on the reduced forest
......@@ -602,6 +619,12 @@ void Forest::mapInternalToBranchIdWalker(
mapInternalToBranchIdWalker(m, aMapInternalToBranchID);
}
// for(std::map<unsigned int, unsigned int >::const_iterator it =
// aMapInternalToBranchID.begin(); it != aMapInternalToBranchID.end(); ++it)
//{
// std::cout << "MAP " << it->first << " " << it->second << "\n";
//}
}
#ifndef NEW_LIKELIHOOD
......@@ -867,6 +890,52 @@ void Forest::computeLikelihoods(const ProbabilityMatrixSet &aSet,
}
}
void Forest::computeLikelihoods(const mfgProbabilityMatrixSet &aSet,
CacheAlignedDoubleVector &aLikelihoods,
const ListDependencies &aDependencies) {
// To speedup the inner OpenMP parallel loop, this is precomputed
// so in the call to computeLikelihoodsWalkerTC &mRoots[site] becomes
// tmp_roots+site
const ForestNode *tmp_roots = &mRoots[0];
double *likelihoods = &aLikelihoods[0];
ListDependencies::const_iterator ivs(aDependencies.begin());
const ListDependencies::const_iterator end(aDependencies.end());
for (; ivs != end; ++ivs) {
// Things that do not change in the parallel loop
const int len = static_cast<int>(ivs->size());
const unsigned int *tmp_ivs = &(*ivs)[0];
#ifdef _MSC_VER
#pragma omp parallel for default(none) shared(aSet, len, tmp_ivs, tmp_roots, \
likelihoods) schedule(static)
#else
#pragma omp parallel for default(shared)
#endif
for (int i = 0; i < len; ++i) {
#ifndef _MSC_VER
#pragma omp task untied
#endif
{
// Compute likelihood array at the root of one tree (the access order is
// the fastest)
const unsigned int tmp = tmp_ivs[i];
const unsigned int site = TreeAndSetsDependencies::getSiteNum(tmp);
const unsigned int set_idx = TreeAndSetsDependencies::getSetNum(tmp);
const double *g =
computeLikelihoodsWalkerTC(tmp_roots + site, aSet, set_idx);
#ifdef USE_CPV_SCALING
likelihoods[set_idx * mNumSites + site] = dot(mCodonFreq, g) * g[N];
#else
likelihoods[set_idx * mNumSites + site] = dot(mCodonFreq, g);
#endif
}
}
}
}
double *Forest::computeLikelihoodsWalkerTC(const ForestNode *aNode,
const ProbabilityMatrixSet &aSet,
unsigned int aSetIdx) {
......@@ -997,6 +1066,136 @@ double *Forest::computeLikelihoodsWalkerTC(const ForestNode *aNode,
return anode_prob;
}
double *Forest::computeLikelihoodsWalkerTC(const ForestNode *aNode,
const mfgProbabilityMatrixSet &aSet,
unsigned int aSetIdx) {
double *anode_prob = aNode->mProb[aSetIdx];
const unsigned int nc = aNode->mChildrenCount;
// Shortcut (on the leaves return immediately the probability vector)
if (nc == 0)
return anode_prob;
bool first = true;
for (unsigned int idx = 0; idx < nc; ++idx) {
// Copy to local var to avoid aliasing
double *anode_other_tree_prob = aNode->mOtherTreeProb[idx];
// If the node is in the same tree recurse and eventually save the value,
// else use the value
if (aNode->isSameTree(idx)) {
const ForestNode *m = aNode->mChildrenList[idx];
const unsigned int branch_id = m->mBranchId;
const int leaf_codon = m->mLeafCodon;
if (leaf_codon >= 0) {
if (first) {
aSet.doTransitionAtLeaf(aSetIdx, branch_id, leaf_codon, anode_prob);
#ifdef USE_CPV_SCALING
anode_prob[N] = normalizeVector(anode_prob);
if (anode_other_tree_prob)
memcpy(anode_other_tree_prob + VECTOR_SLOT * aSetIdx, anode_prob,
(N + 1) * sizeof(double));
#else
if (anode_other_tree_prob)
memcpy(anode_other_tree_prob + VECTOR_SLOT * aSetIdx, anode_prob,
N * sizeof(double));
#endif
first = false;
} else {
#ifdef USE_CPV_SCALING
double ALIGN64 temp[N + 1];
#else
double ALIGN64 temp[N];
#endif
double *x = anode_other_tree_prob
? anode_other_tree_prob + VECTOR_SLOT * aSetIdx
: temp;
aSet.doTransitionAtLeaf(aSetIdx, branch_id, leaf_codon, x);
#ifdef USE_CPV_SCALING
x[N] = normalizeVector(x);
anode_prob[N] *= x[N];
#endif
elementWiseMult(anode_prob, x);
}
} else {
if (first) {
double *g = computeLikelihoodsWalkerTC(m, aSet, aSetIdx);
aSet.doTransition(aSetIdx, branch_id, g, anode_prob);
#ifdef USE_CPV_SCALING
anode_prob[N] = normalizeVector(anode_prob) * g[N];
if (anode_other_tree_prob)
memcpy(anode_other_tree_prob + VECTOR_SLOT * aSetIdx, anode_prob,
(N + 1) * sizeof(double));
#else
if (anode_other_tree_prob)
memcpy(anode_other_tree_prob + VECTOR_SLOT * aSetIdx, anode_prob,
N * sizeof(double));
#endif
first = false;
} else {
#ifdef USE_CPV_SCALING
double ALIGN64 temp[N + 1];
#else
double ALIGN64 temp[N];
#endif
double *x = anode_other_tree_prob
? anode_other_tree_prob + VECTOR_SLOT * aSetIdx
: temp;
double *g = computeLikelihoodsWalkerTC(m, aSet, aSetIdx);
aSet.doTransition(aSetIdx, branch_id, g, x);
#ifdef USE_CPV_SCALING
x[N] = normalizeVector(x) * g[N];
anode_prob[N] *= x[N];
#endif
elementWiseMult(anode_prob, x);
}
}
} else {
if (first) {
#ifdef USE_CPV_SCALING
memcpy(anode_prob, anode_other_tree_prob + VECTOR_SLOT * aSetIdx,
(N + 1) * sizeof(double));
#else
memcpy(anode_prob, anode_other_tree_prob + VECTOR_SLOT * aSetIdx,
N * sizeof(double));
#endif
first = false;
} else {
#ifdef USE_CPV_SCALING
anode_prob[N] *= anode_other_tree_prob[VECTOR_SLOT * aSetIdx + N];
#endif
elementWiseMult(anode_prob,
anode_other_tree_prob + VECTOR_SLOT * aSetIdx);
}
}
}
#ifdef USE_LAPACK
switch (nc) {
case 1:
elementWiseMult(anode_prob, mInvCodonFreq);
break;
case 2:
elementWiseMult(anode_prob, mInv2CodonFreq);
break;
case 3:
elementWiseMult(anode_prob, mInv2CodonFreq);
elementWiseMult(anode_prob, mInvCodonFreq);
break;
case 4:
elementWiseMult(anode_prob, mInv2CodonFreq);
elementWiseMult(anode_prob, mInv2CodonFreq);
break;
default:
for (unsigned int idx = 0; idx < nc; ++idx)
elementWiseMult(anode_prob, mInvCodonFreq);
break;
}
#endif
return anode_prob;
}
#endif
#ifdef USE_DAG
......
......@@ -87,13 +87,16 @@ public:
/// program
/// @param[out] aBranchStart The first branch to be marked as foreground
/// @param[out] aBranchEnd The last branch to be marked as foreground
/// @param[out] aFgSet list of id of foreground branches
/// @param[out] aIbSet list of id of internal branches
///
/// @return True if all branches are selected
///
/// @exception FastCodeMLFatal Invalid range from command line
///
bool getBranchRange(const CmdLine &aCmdLine, size_t &aBranchStart,
size_t &aBranchEnd) const;
size_t &aBranchEnd, std::set<int> &aFgSet,
std::set<int> &aIbSet) const;
/// Reduce common subtrees on the whole forest.
///
......@@ -128,6 +131,10 @@ public:
void computeLikelihoods(const ProbabilityMatrixSet &aSet,
CacheAlignedDoubleVector &aLikelihoods,
const ListDependencies &aDependencies);
void computeLikelihoods(const mfgProbabilityMatrixSet &aSet,
CacheAlignedDoubleVector &aLikelihoods,
const ListDependencies &aDependencies);
#endif
#ifdef NON_RECURSIVE_VISIT
......@@ -194,6 +201,20 @@ public:
///
size_t getMarkedInternalBranch(void) const { return mMarkedInternalBranch; }
/// Get all marked branches
///
/// @return The internal branch index of the branches marked in the tree file.
/// UINT_MAX otherwise.
///
std::set<int> getMarkedBranches(void) const { return mMarkedBranches; }
/// Get all internal branches
///
/// @return The internal branch index of the internal branches in the tree
/// file. UINT_MAX otherwise.
///
std::set<int> getInternalBranches(void) const { return mInternalBranches; }
/// Get site multiplicity values.
///
/// @return Reference to the array of site multiplicities
......@@ -334,6 +355,10 @@ private:
double *computeLikelihoodsWalkerTC(const ForestNode *aNode,
const ProbabilityMatrixSet &aSet,
unsigned int aSetIdx);
double *computeLikelihoodsWalkerTC(const ForestNode *aNode,
const mfgProbabilityMatrixSet &aSet,
unsigned int aSetIdx);
#endif
#ifdef NON_RECURSIVE_VISIT
......@@ -378,22 +403,25 @@ private:
const double *mInv2CodonFreq; ///< Squared inverse of the codon frequencies
size_t mNumBranches; ///< Total number of branches of the original tree
std::vector<ForestNode> mRoots; ///< The roots of the forest's trees. Its
/// length is the number of valid sites
///length is the number of valid sites
std::vector<double> mSiteMultiplicity; ///< Multiplicity of the valid sites
size_t
mNumInternalBranches; ///< Total number of branches of the original tree
std::vector<unsigned int> mTableInternalToBranchID; ///< Map from internal
/// branch number to branch
/// number
///branch number to branch
///number
/// Here are global data that will be removed from the various (site) trees
std::vector<std::string> mNodeNames; ///< List of node names. Zero is the
/// root, then its first child and so on
std::vector<double> mBranchLengths; ///< List of branch lengths (read from
/// file or stored here to be exported in
/// the tree file)
size_t mMarkedInternalBranch; ///< Number of the internal branch as marked in
/// the tree file
///root, then its first child and so on
std::vector<double> mBranchLengths; ///< List of branch lengths (read from
///file or stored here to be exported in
///the tree file)
size_t mMarkedInternalBranch; ///< Number of the internal branch as marked in
///the tree file
std::set<int> mMarkedBranches; ///< id of the marked branches in the tree file
std::set<int>
mInternalBranches; ///< id of the internal branches in the tree file
#ifdef NEW_LIKELIHOOD
......
REQUIREMENTS to generate the executable:
REQUIREMENTS:
* Computer system:
Linux preferred, but sources are portable to other platforms
* C++ compiler, e.g. GCC 4
* CMake 2.8.0 (including ccmake) or later recommended, although compilation possible without
* Boost::Spirit 1.42.0 or later recommended, see http://boost-spirit.com/home/
* Reasonably new BLAS implementation (e.g. OpenBLAS, Goto2, ACML, MKL)
* Boost::Spirit library 1.42.0 or later recommended, see http://boost-spirit.com/home/
* Reasonably new BLAS library (e.g. OpenBLAS, Goto2, ACML, MKL)
packages from various Linux distributions can be used, but this deteriorates performance
recommended: OpenBLAS (http://xianyi.github.io/OpenBLAS/) or Intel MKL
* Reasonably new LAPACK library (e.g. original LAPACK or ACML, MKL)
packages from various Linux distributions can be used, but this deteriorates performance
* Nlopt library (see http://ab-initio.mit.edu/wiki/index.php/NLopt)
* MPI library (see http://www.mpich.org/static/downloads/)
INSTALLATION:
1-Prepring libraries:
INSTALLATION: How to generate the FastCodeML executable:
0)
* Generate BLAS, if necessary
Download OpenBLAS from https://github.com/xianyi/OpenBLAS/
make USE_THREAD=0 USE_OPENMP=1
make USE_OPENMP=1
make PREFIX=path/to/installation/blas install
* Generate LAPACK, if necessary
Download LAPACK from http://www.netlib.org/lapack/
Take the gfortran version of the makefile:
cp INSTALL/make.inc.gfortran ./make.inc
In the make.inc file set BLASLIB to the full path to the generated BLAS
e.g. BLASLIB = path/to/installation/blas/lib/libopenblas.a
Two libraries will be created: lapack.a and tmglib.a. We put them in a folder and merge them:
Run $make
Two libraries will be created: lapack.a and tmglib.a.
Put them inside a folder and merge them:
unpack in a new directory
ar x liblapack.a and ar x tmglib.a
merge into a new LAPACK library
ar x liblapack.a
ar x tmglib.a
merge all produced object files into a new LAPACK library
ar rcs liblapack.a *.o
copy the liblapack.a into the library folder
* Generate NLopt library, if necessary
Download NLopt from http://ab-initio.mit.edu/wiki/index.php/NLopt
Install it
./configure --prefix=path/to/installation/nlopt
make USE_THREAD=0
make
make install
* Generate BOOST library (>= 1.42.0), if necessary
Download BOOST from http://sourceforge.net/projects/boost/files/latest/download?source=files
Install it
./bootstrap.sh --prefix=path/to/installation/boost
./b2 install
* Generate MPI library, if necessary
Download MPI from http://www.mpich.org/static/downloads/
Install it
......@@ -44,24 +59,40 @@ INSTALLATION: How to generate the FastCodeML executable:
make
make install
1)
2-Generating executable:
* Open FastCodeML tarball
tar xvfz FastCodeML-X.Y.Z.tar.gz
cd FastCodeML-X.Y.Z/
* Edit CMakeLists.txt to enable/disable (ON/OFF) optimization libraries (in "Get the configuration switches" section)
switch USE_MPI and USE_OPENMP ON/OFF (other default settings should be ok)
* Set paths for libraries (change and execute SETPATHS) (use only what you require)
default:
USE_CPV_SCALING ON
USE_IDENTITY_MATRIX OFF
USE_LAPACK ON
USE_LIKELIHOOD_METHOD Original
USE_MKL_VML OFF
USE_MPI OFF
USE_OPENMP ON
USE_ORIGINAL_PROPORTIONS ON
* Set paths for required libraries (using export)
BLAS_LIB_DIR = path/to/installation/blas/lib
LAPACK_LIB_DIR = path/to/installation/blas/include
LAPACK_LIB_DIR = path/to/installation/lapack/lib
MKL_INCLUDE_DIR = path/to/installation/mkl/include
NLOPT_INCLUDE_DIR = path/to/installation/nlopt/include
NLOPT_LIB_DIR = path/to/installation/nlopt/lib
MATH_LIB_NAMES = "openblas;lapack;gfortranbegin;gfortran"
* Build executable
ccmake .
ccmake .
make
* An executable "fast" is created
Computer system:
* Linux preferred, but sources are portable to other platforms
......@@ -6,6 +6,9 @@
#undef __SSE2__
#endif
#include <iostream>
#include <stdio.h>
#ifdef __SSE2__
#include <emmintrin.h>
#endif
......@@ -18,6 +21,64 @@
#include <cmath>
#endif
// Codeml random number generator (setting seed, returning random number)
static unsigned int z_rndu = 1237;
static int w_rndu = 1237;
inline void SetSeedCodeml(int seed, int PrintSeed) {
int i;
FILE *frand, *fseed;
if (sizeof(unsigned int) != 4)
std::cout << "oh-oh, we are in trouble. int not 32-bit?";
if (seed <= 0) {
frand = fopen("/dev/urandom", "r");
if (frand) {
for (i = 0, seed = 0; i < sizeof(unsigned int); i++)
seed += (seed << 8) + getc(frand);
seed = 2 * seed + 1;
fclose(frand);
} else {
seed = 1234567891 * (int)time(NULL) + 1;
}
seed = abs(seed);
// if(PrintSeed) {
// fseed = fopen("SeedUsed", "w");
// if(fseed == NULL) std::cout << "can't open file SeedUsed.";
// fprin