Skip to content
GitLab
Commit 654e04d9 authored by Iakov Davydov's avatar Iakov Davydov
Browse files

format all the code

parent e493e3f0
Expand all Hide whitespace changes
Inline Side-by-side
AlignedAllocator.cpp
View file @ 654e04d9
// The following headers are required for all allocators.
#include <cstddef> // Required for size_t and ptrdiff_t and NULL
#include <cstddef> // Required for size_t and ptrdiff_t and NULL
#include <stdexcept> // Required for std::length_error
// The following headers contain stuff that AlignedAllocator uses.
#include <cstdlib> // For malloc() and free()
#include <cstdlib> // For malloc() and free()
// For XMT the following function should be supplied till the bug is fixed by Cray
// For XMT the following function should be supplied till the bug is fixed by
// Cray
#ifdef __MTA__
//extern "C" int posix_memalign(void **memptr, size_t alignment, size_t size);
static int posix_memalign(void **memptr, size_t alignment, size_t size)
{
*memptr = malloc(size);
return 0;
// extern "C" int posix_memalign(void **memptr, size_t alignment, size_t size);
static int posix_memalign(void **memptr, size_t alignment, size_t size) {
*memptr = malloc(size);
return 0;
}
#endif
......@@ -23,8 +23,7 @@ static int posix_memalign(void **memptr, size_t alignment, size_t size)
#endif
// Alignment must be power of 2 (1,2,4,8,16...)
void* alignedMalloc(size_t aSize, size_t aAlignment)
{
void *alignedMalloc(size_t aSize, size_t aAlignment) {
#ifdef _MSC_VER
#if 0
--aAlignment;
......@@ -35,34 +34,34 @@ void* alignedMalloc(size_t aSize, size_t aAlignment)
reinterpret_cast<uintptr_t*>(o)[-1] = r;
return reinterpret_cast<void*>(o);
#endif
return _aligned_malloc(aSize, aAlignment);
return _aligned_malloc(aSize, aAlignment);
#else
void* ptr = NULL;
if(posix_memalign(&ptr, aAlignment, aSize)) return NULL;
return ptr;
void *ptr = NULL;
if (posix_memalign(&ptr, aAlignment, aSize))
return NULL;
return ptr;
#endif
}
void alignedFree(void* aPtr)
{
if(!aPtr) return;
void alignedFree(void *aPtr) {
if (!aPtr)
return;
#ifdef _MSC_VER
#if 0
free(reinterpret_cast<void*>(reinterpret_cast<uintptr_t*>(aPtr)[-1]));
#endif
_aligned_free(aPtr);
_aligned_free(aPtr);
#else
free(aPtr);
free(aPtr);
#endif
}
#if 0
// The following headers contain stuff that main() uses.
#include <iostream> // For std::cout
#include <ostream> // For std::endl
#include <vector> // For std::vector
#include <iostream> // For std::cout
#include <ostream> // For std::endl
#include <vector> // For std::vector
#include "AlignedAllocator.h"
int main()
......@@ -96,4 +95,3 @@ int main()
cout << endl << "Destroying l:" << endl;
}
#endif
AlignedAllocator.h
View file @ 654e04d9
......@@ -3,7 +3,7 @@
#define ALIGNED_ALLOCATOR_H
// The following headers are required for all allocators.
#include <cstddef> // Required for size_t and ptrdiff_t and NULL
#include <cstddef> // Required for size_t and ptrdiff_t and NULL
//#include <new> // Required for placement new and std::bad_alloc
#include <stdexcept> // Required for std::length_error
......@@ -11,180 +11,166 @@
/// Aligned allocator definition.
/// It will be used to obtain a vector aligned to a given power of 2.
/// Example allocation aligned to 64: std::vector<double, AlignedAllocator<double, 64> > aligned_vector;
/// Example allocation aligned to 64: std::vector<double,
/// AlignedAllocator<double, 64> > aligned_vector;
///
/// @author Mario Valle - Swiss National Supercomputing Centre (CSCS)
/// @date 2010-12-22 (initial version)
/// @version 1.1
///
template <typename T, size_t A> class AlignedAllocator
{
template <typename T, size_t A> class AlignedAllocator {
public:
// The following will be the same for virtually all allocators.
typedef T * pointer;
typedef const T * const_pointer;
typedef T& reference;
typedef const T& const_reference;
typedef T value_type;
typedef size_t size_type;
typedef ptrdiff_t difference_type;
T * address(T& r) const
{
return &r;
}
const T * address(const T& s) const
{
return &s;
}
size_t max_size() const
{
// The following has been carefully written to be independent of
// the definition of size_t and to avoid signed/unsigned warnings.
return (static_cast<size_t>(0) - static_cast<size_t>(1)) / sizeof(T);
}
/// Internal definition for AlignedAllocator.
/// The following must be the same for all allocators.
///
template <typename U> struct rebind
{
typedef AlignedAllocator<U, A> other;
};
bool operator!=(const AlignedAllocator& other) const
{
return !(*this == other);
}
void construct(T * const p, const T& t) const
{
void * const pv = static_cast<void *>(p);
new (pv) T(t);
}
void destroy(T * const p) const; // Defined below.
/// Returns true if and only if storage allocated from *this
/// can be deallocated from other, and vice versa.
/// Always returns true for stateless allocators.
///
/// @param[in] other The other allocator to be compared
///
/// @return Always true, this is a stateless allocator.
///
bool operator==(const AlignedAllocator& other) const
{
return true;
}
/// Default constructor. It should be empty for stateless allocators.
///
AlignedAllocator() { }
/// Default copy constructor. It should be empty for stateless allocators.
///
AlignedAllocator(const AlignedAllocator&) { }
/// Default rebinding constructor. It should be empty for stateless allocators.
///
template <typename U> AlignedAllocator(const AlignedAllocator<U, A>&) { }
/// Default destructor. It should be empty for stateless allocators.
///
~AlignedAllocator() { }
/// The main allocator routine.
/// The following will be different for each allocator.
///
/// @param[in] n Number of objects of type T to allocate
///
/// @return The allocated memory
///
/// @exception std::length_error Integer overflow
/// @exception std::bad_alloc Memory allocation failure
///
T * allocate(const size_t n) const
{
// AlignedAllocator prints a diagnostic message to demonstrate
// what it's doing. Real allocators won't do this.
//std::cout << "Allocating " << n << (n == 1 ? " object" : " objects")
// << " of size " << sizeof(T) << " aligned on " << A << std::endl;
// The return value of allocate(0) is unspecified.
// AlignedAllocator returns NULL in order to avoid depending
// on malloc(0)'s implementation-defined behavior
// (the implementation can define malloc(0) to return NULL,
// in which case the bad_alloc check below would fire).
// All allocators can return NULL in this case.
if (n == 0) return NULL;
// All allocators should contain an integer overflow check.
// The Standardization Committee recommends that std::length_error
// be thrown in the case of integer overflow.
if (n > max_size()) throw std::length_error("AlignedAllocator<T>::allocate() - Integer overflow.");
// AlignedAllocator wraps aligned malloc().
void * const pv = alignedMalloc(n * sizeof(T), A);
// Allocators should throw std::bad_alloc in the case of memory allocation failure.
if(pv == NULL) throw std::bad_alloc();
return static_cast<T *>(pv);
}
void deallocate(T * const p, const size_t /*n*/) const
{
// AlignedAllocator prints a diagnostic message to demonstrate
// what it's doing. Real allocators won't do this.
//std::cout << "Deallocating " << n << (n == 1 ? " object" : " objects")
// << " of size " << sizeof(T) << "." << std::endl;
// AlignedAllocator wraps aligned free().
alignedFree(p);
}
/// The following will be the same for all allocators that ignore hints.
///
template <typename U> T * allocate(const size_t n, const U * /* const hint */) const
{
return allocate(n);
}
// The following will be the same for virtually all allocators.
typedef T *pointer;
typedef const T *const_pointer;
typedef T &reference;
typedef const T &const_reference;
typedef T value_type;
typedef size_t size_type;
typedef ptrdiff_t difference_type;
T *address(T &r) const { return &r; }
const T *address(const T &s) const { return &s; }
size_t max_size() const {
// The following has been carefully written to be independent of
// the definition of size_t and to avoid signed/unsigned warnings.
return (static_cast<size_t>(0) - static_cast<size_t>(1)) / sizeof(T);
}
/// Internal definition for AlignedAllocator.
/// The following must be the same for all allocators.
///
template <typename U> struct rebind { typedef AlignedAllocator<U, A> other; };
bool operator!=(const AlignedAllocator &other) const {
return !(*this == other);
}
void construct(T *const p, const T &t) const {
void *const pv = static_cast<void *>(p);
new (pv) T(t);
}
void destroy(T *const p) const; // Defined below.
/// Returns true if and only if storage allocated from *this
/// can be deallocated from other, and vice versa.
/// Always returns true for stateless allocators.
///
/// @param[in] other The other allocator to be compared
///
/// @return Always true, this is a stateless allocator.
///
bool operator==(const AlignedAllocator &other) const { return true; }
/// Default constructor. It should be empty for stateless allocators.
///
AlignedAllocator() {}
/// Default copy constructor. It should be empty for stateless allocators.
///
AlignedAllocator(const AlignedAllocator &) {}
/// Default rebinding constructor. It should be empty for stateless
/// allocators.
///
template <typename U> AlignedAllocator(const AlignedAllocator<U, A> &) {}
/// Default destructor. It should be empty for stateless allocators.
///
~AlignedAllocator() {}
/// The main allocator routine.
/// The following will be different for each allocator.
///
/// @param[in] n Number of objects of type T to allocate
///
/// @return The allocated memory
///
/// @exception std::length_error Integer overflow
/// @exception std::bad_alloc Memory allocation failure
///
T *allocate(const size_t n) const {
// AlignedAllocator prints a diagnostic message to demonstrate
// what it's doing. Real allocators won't do this.
// std::cout << "Allocating " << n << (n == 1 ? " object" : " objects")
// << " of size " << sizeof(T) << " aligned on " << A << std::endl;
// The return value of allocate(0) is unspecified.
// AlignedAllocator returns NULL in order to avoid depending
// on malloc(0)'s implementation-defined behavior
// (the implementation can define malloc(0) to return NULL,
// in which case the bad_alloc check below would fire).
// All allocators can return NULL in this case.
if (n == 0)
return NULL;
// All allocators should contain an integer overflow check.
// The Standardization Committee recommends that std::length_error
// be thrown in the case of integer overflow.
if (n > max_size())
throw std::length_error(
"AlignedAllocator<T>::allocate() - Integer overflow.");
// AlignedAllocator wraps aligned malloc().
void *const pv = alignedMalloc(n * sizeof(T), A);
// Allocators should throw std::bad_alloc in the case of memory allocation
// failure.
if (pv == NULL)
throw std::bad_alloc();
return static_cast<T *>(pv);
}
void deallocate(T *const p, const size_t /*n*/) const {
// AlignedAllocator prints a diagnostic message to demonstrate
// what it's doing. Real allocators won't do this.
// std::cout << "Deallocating " << n << (n == 1 ? " object" : " objects")
// << " of size " << sizeof(T) << "." << std::endl;
// AlignedAllocator wraps aligned free().
alignedFree(p);
}
/// The following will be the same for all allocators that ignore hints.
///
template <typename U>
T *allocate(const size_t n, const U * /* const hint */) const {
return allocate(n);
}
private:
/// Allocators are not required to be assignable, so
/// all allocators should have a private unimplemented
/// assignment operator. Note that this will trigger the
/// off-by-default (enabled under /Wall) warning C4626
/// "assignment operator could not be generated because a
/// base class assignment operator is inaccessible" within
/// the STL headers, but that warning is useless.
//AlignedAllocator& operator=(const AlignedAllocator& a) {return this;}
//AlignedAllocator& operator=(const AlignedAllocator& a) {return const_cast<AlignedAllocator&>(a);}
AlignedAllocator& operator=(const AlignedAllocator&);
/// Allocators are not required to be assignable, so
/// all allocators should have a private unimplemented
/// assignment operator. Note that this will trigger the
/// off-by-default (enabled under /Wall) warning C4626
/// "assignment operator could not be generated because a
/// base class assignment operator is inaccessible" within
/// the STL headers, but that warning is useless.
// AlignedAllocator& operator=(const AlignedAllocator& a) {return this;}
// AlignedAllocator& operator=(const AlignedAllocator& a) {return
// const_cast<AlignedAllocator&>(a);}
AlignedAllocator &operator=(const AlignedAllocator &);
};
// A compiler bug causes it to believe that p->~T() doesn't reference p.
#ifdef _MSC_VER
#pragma warning(push)
#pragma warning(disable: 4100) // unreferenced formal parameter
#pragma warning(push)
#pragma warning(disable : 4100) // unreferenced formal parameter
#endif
/// The definition of destroy() must be the same for all allocators.
template <typename T, size_t A> inline void AlignedAllocator<T, A>::destroy(T * const p) const
{
p->~T();
template <typename T, size_t A>
inline void AlignedAllocator<T, A>::destroy(T *const p) const {
p->~T();
}
#ifdef _MSC_VER
#pragma warning(pop)
#pragma warning(pop)
#endif
#endif
AlignedMalloc.h
View file @ 654e04d9
......@@ -10,12 +10,13 @@
///
/// @return Pointer to the allocated area
///
extern void* alignedMalloc(size_t aSize, size_t aAlignment);
extern void *alignedMalloc(size_t aSize, size_t aAlignment);
/// Free the memory allocated with alignedMalloc.
///
/// @param[in] aPtr The pointer to the memory (allocated with alignedMalloc) to be freed
/// @param[in] aPtr The pointer to the memory (allocated with alignedMalloc) to
/// be freed
///
extern void alignedFree(void* aPtr);
extern void alignedFree(void *aPtr);
#endif
BayesTest.cpp
View file @ 654e04d9
This diff is collapsed.
BayesTest.h
View file @ 654e04d9
......@@ -8,154 +8,181 @@
#include "ProbabilityMatrixSet.h"
#include "VerbosityLevels.h"
/// The minimum value for class 2 sites probability to be a positive selection site.
/// The minimum value for class 2 sites probability to be a positive selection
/// site.
///
const static double MIN_PROB = 0.50;
const static double ONE_STAR_PROB = 0.95;
const static double MIN_PROB = 0.50;
const static double ONE_STAR_PROB = 0.95;
const static double TWO_STARS_PROB = 0.99;
///@cond Private
/// Helper class to compute BEB_N1D^BEB_DIMS at compile time (that is Y^N)
///
template<unsigned int Y, unsigned int N>
class Pow
{
template <unsigned int Y, unsigned int N> class Pow {
public:
static const int value = Y * Pow<Y, N-1>::value;
static const int value = Y * Pow<Y, N - 1>::value;
};
/// Specialization of the above class to end the recursion
///
template<unsigned int Y>
class Pow<Y, 1>
{
template <unsigned int Y> class Pow<Y, 1> {
public:
static const int value = Y;
static const int value = Y;
};
///@endcond
/// Tests to find the sites under positive selection.
///
/// @author Mario Valle - Swiss National Supercomputing Centre (CSCS)
/// @date 2010-12-22 (initial version)
/// @version 1.1
///
class BayesTest
{
class BayesTest {
public:
/// Constructor.
///
/// @param[in] aForest The forest
/// @param[in] aVerbose The verbosity level
/// @param[in] aNoReduction If true the dependencies computed are for no reduced forests
///
explicit BayesTest(Forest& aForest, unsigned int aVerbose=0, bool aNoReduction=true)
: mForest(aForest), mNumSites(mForest.getNumSites()), mSiteClassProb(BEB_DIMS*mNumSites),
mVerbose(aVerbose), mPriors(mNumSites*BEB_NUM_CAT), mDependencies(mForest, aVerbose), mBEBset(mForest.getNumBranches())
{
// Create the dependency list for forest likelihood computation
mDependencies.computeDependencies(1, aNoReduction);
if(mVerbose >= VERBOSE_ONLY_RESULTS) mDependencies.print("TEST FOR BEB (before optimization)");
mDependencies.optimizeDependencies();
if(mVerbose >= VERBOSE_ONLY_RESULTS) mDependencies.print("TEST FOR BEB");
}
/// Destructor.
///
~BayesTest() {}
/// Bayes Empirical Bayes (BEB) test.
///
/// @param[in] aVars The variables optimized at the end of the H1 run.
/// @param[in] aFgBranch The foreground branch under test.
/// @param[in] aScales The two scales ([0] bg; [1] fg) to rescale the branch lengths. They are computed in H1.
///
void computeBEB(const std::vector<double>& aVars, size_t aFgBranch, const std::vector<double>& aScales);
/// Print the sites under positive selection.
///
/// @param[in] aFgBranch Identifier of the branch marked as foreground branch
///
void printPositiveSelSites(size_t aFgBranch) const;
/// Extract the sites under positive selection and the corresponding probabilities.
///
/// @param[out] aPositiveSelSites Vector of sites under positive selection
/// @param[out] aPositiveSelSitesProb Corresponding probabilities
///
void extractPositiveSelSites(std::vector<unsigned int>& aPositiveSelSites, std::vector<double>& aPositiveSelSitesProb) const;
/// Constructor.
///
/// @param[in] aForest The forest
/// @param[in] aVerbose The verbosity level
/// @param[in] aNoReduction If true the dependencies computed are for no
/// reduced forests
///
explicit BayesTest(Forest &aForest, unsigned int aVerbose = 0,
bool aNoReduction = true)
: mForest(aForest), mNumSites(mForest.getNumSites()),
mSiteClassProb(BEB_DIMS * mNumSites), mVerbose(aVerbose),
mPriors(mNumSites * BEB_NUM_CAT), mDependencies(mForest, aVerbose),
mBEBset(mForest.getNumBranches()) {
// Create the dependency list for forest likelihood computation
mDependencies.computeDependencies(1, aNoReduction);
if (mVerbose >= VERBOSE_ONLY_RESULTS)
mDependencies.print("TEST FOR BEB (before optimization)");
mDependencies.optimizeDependencies();
if (mVerbose >= VERBOSE_ONLY_RESULTS)
mDependencies.print("TEST FOR BEB");
}
/// Destructor.
///
~BayesTest() {}
/// Bayes Empirical Bayes (BEB) test.
///
/// @param[in] aVars The variables optimized at the end of the H1 run.
/// @param[in] aFgBranch The foreground branch under test.
/// @param[in] aScales The two scales ([0] bg; [1] fg) to rescale the branch
/// lengths. They are computed in H1.
///
void computeBEB(const std::vector<double> &aVars, size_t aFgBranch,
const std::vector<double> &aScales);
/// Print the sites under positive selection.
///
/// @param[in] aFgBranch Identifier of the branch marked as foreground branch
///
void printPositiveSelSites(size_t aFgBranch) const;
/// Extract the sites under positive selection and the corresponding
/// probabilities.
///
/// @param[out] aPositiveSelSites Vector of sites under positive selection
/// @param[out] aPositiveSelSitesProb Corresponding probabilities
///
void
extractPositiveSelSites(std::vector<unsigned int> &aPositiveSelSites,
std::vector<double> &aPositiveSelSitesProb) const;
private:
/// This sets up the grid (mPara[][]) according to the priors.
/// It calculates the probability of data at each site given w: f(f_h|w).
/// This is calculated using the branch model (NSsites = 0 model = 2), with
/// BayesEB=2 used to force the use of the correct scale factors in GetPMatBranch().
///
///@verbatim
/// Order of site classes for iw or f(x_h|w):
/// back fore num.sets
/// Branchsite A (121 sets)
/// site class 0: w0 w0 10
/// site class 1: w1=1 w1=1 1
/// site class 2a: w0 w2 100
/// site class 2b: w1=1 w2 10
///@endverbatim
///
/// @param[in] aVars The variables optimized at the end of the H1 run.
/// @param[in] aSiteMultiplicity The site multiplicity vector.
/// @param[in] aFgBranch The foreground branch under test.
/// @param[in] aScales The two scales ([0] bg; [1] fg) to rescale the branch lengths. They are computed in H1.
///
/// @return The computed scale.
///
double getGridParams(const std::vector<double>& aVars, const std::vector<double>& aSiteMultiplicity, size_t aFgBranch, const std::vector<double>& aScales);
/// This gives the indices (ix, iy) and the coordinates (aProbX, aProbY, 1-aProbX-aProbY) for
/// the aTriangleIdx-th triangle, with aTriangleIdx from 0, 1, ..., BEB_N1D*BEB_N1D-1.
///
/// The ternary graph (0-1 on each axis) is partitioned into BEB_N1D*BEB_N1D equal-sized triangles.
/// In the first row (ix=0), there is one triangle (iy=0);
/// In the second row (ix=1), there are 3 triangles (iy=0,1,2);
/// In the i-th row (ix=i), there are 2*i+1 triangles (iy=0,1,...,2*i).
///
/// aProbX rises when ix goes up, but aProbY decreases when iy increases. (aProbX, aProbY) is the
/// centroid in the ij-th small triangle. aProbX and aProbY each takes on 2*BEB_N1D-1 possible values.
///
/// @param[out] aProbX The p0 value on the X axis of the triangular grid.
/// @param[out] aProbY The p1 value on the Y axis of the triangular grid.
/// @param[in] aTriangleIdx The index inside the triangular grid.
///
void getIndexTernary(double* aProbX, double* aProbY, unsigned int aTriangleIdx);
/// This sets up the grid (mPara[][]) according to the priors.
/// It calculates the probability of data at each site given w: f(f_h|w).
/// This is calculated using the branch model (NSsites = 0 model = 2), with
/// BayesEB=2 used to force the use of the correct scale factors in
/// GetPMatBranch().
///
///@verbatim
/// Order of site classes for iw or f(x_h|w):
/// back fore num.sets
/// Branchsite A (121 sets)