Skip to content
GitLab
Commit 2589ab54 authored by Iakov Davydov's avatar Iakov Davydov
Browse files

spaces to tabs

parent 5e1e24c9
Expand all Hide whitespace changes
Inline Side-by-side
AlignedAllocator.cpp
View file @ 2589ab54
// 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
#ifdef __MTA__
......@@ -28,12 +28,12 @@ void* alignedMalloc(size_t aSize, size_t aAlignment)
#ifdef _MSC_VER
#if 0
--aAlignment;
uintptr_t r = reinterpret_cast<uintptr_t>(malloc(aSize + aAlignment + sizeof(uintptr_t)));
if(!r) return NULL;
uintptr_t t = r + sizeof(uintptr_t);
uintptr_t o = (t + aAlignment) & ~static_cast<uintptr_t>(aAlignment);
reinterpret_cast<uintptr_t*>(o)[-1] = r;
return reinterpret_cast<void*>(o);
uintptr_t r = reinterpret_cast<uintptr_t>(malloc(aSize + aAlignment + sizeof(uintptr_t)));
if(!r) return NULL;
uintptr_t t = r + sizeof(uintptr_t);
uintptr_t o = (t + aAlignment) & ~static_cast<uintptr_t>(aAlignment);
reinterpret_cast<uintptr_t*>(o)[-1] = r;
return reinterpret_cast<void*>(o);
#endif
return _aligned_malloc(aSize, aAlignment);
#else
......@@ -45,14 +45,14 @@ void* alignedMalloc(size_t aSize, size_t aAlignment)
void alignedFree(void* aPtr)
{
if(!aPtr) return;
if(!aPtr) return;
#ifdef _MSC_VER
#if 0
free(reinterpret_cast<void*>(reinterpret_cast<uintptr_t*>(aPtr)[-1]));
free(reinterpret_cast<void*>(reinterpret_cast<uintptr_t*>(aPtr)[-1]));
#endif
_aligned_free(aPtr);
#else
free(aPtr);
free(aPtr);
#endif
}
......@@ -60,40 +60,40 @@ void alignedFree(void* aPtr)
#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()
{
using namespace std;
using namespace std;
cout << "Constructing l:" << endl;
cout << "Constructing l:" << endl;
vector<double, AlignedAllocator<double, 8> > l;
vector<double, AlignedAllocator<double, 8> > l;
l.reserve(10);
cout << endl << "l.push_back(1729):" << endl;
cout << endl << "l.push_back(1729):" << endl;
l.push_back(1729.);
l.push_back(1729.);
cout << endl << "l.push_back(2161):" << endl;
cout << endl << "l.push_back(2161):" << endl;
l.push_back(2161.);
l.push_back(2161.);
cout << endl;
cout << endl;
double* p = &l[0];
int x = reinterpret_cast<int>(p);
cout << "Aligned on 16: " << x%16 << endl;
cout << "Aligned on 8: " << x%8 << endl;
cout << "Aligned on 4: " << x%4 << endl;
cout << endl;
cout << "Aligned on 8: " << x%8 << endl;
cout << "Aligned on 4: " << x%4 << endl;
cout << endl;
for (vector<double, AlignedAllocator<double, 8> >::const_iterator i = l.begin(); i != l.end(); ++i) {
cout << "Element: " << *i << endl;
}
for (vector<double, AlignedAllocator<double, 8> >::const_iterator i = l.begin(); i != l.end(); ++i) {
cout << "Element: " << *i << endl;
}
cout << endl << "Destroying l:" << endl;
cout << endl << "Destroying l:" << endl;
}
#endif
AlignedAllocator.h
View file @ 2589ab54
......@@ -3,8 +3,8 @@
#define ALIGNED_ALLOCATOR_H
// The following headers are required for all allocators.
#include <cstddef> // Required for size_t and ptrdiff_t and NULL
//#include <new> // Required for placement new and std::bad_alloc
#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
#include "AlignedMalloc.h"
......@@ -13,89 +13,89 @@
/// 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;
///
/// @author Mario Valle - Swiss National Supercomputing Centre (CSCS)
/// @date 2010-12-22 (initial version)
/// @version 1.1
/// @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
{
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;
// 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
T * address(T& r) const
{
return &r;
}
return &r;
}
const T * address(const T& s) const
const T * address(const T& s) const
{
return &s;
}
return &s;
}
size_t max_size() const
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);
}
// 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.
/// The following must be the same for all allocators.
///
template <typename U> struct rebind
template <typename U> struct rebind
{
typedef AlignedAllocator<U, A> other;
};
typedef AlignedAllocator<U, A> other;
};
bool operator!=(const AlignedAllocator& other) const
bool operator!=(const AlignedAllocator& other) const
{
return !(*this == other);
}
return !(*this == other);
}
void construct(T * const p, const T& t) const
void construct(T * const p, const T& t) const
{
void * const pv = static_cast<void *>(p);
new (pv) T(t);
}
void * const pv = static_cast<void *>(p);
new (pv) T(t);
}
void destroy(T * const p) const; // Defined below.
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.
/// 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
bool operator==(const AlignedAllocator& other) const
{
return true;
}
return true;
}
/// Default constructor. It should be empty for stateless allocators.
/// Default constructor. It should be empty for stateless allocators.
///
AlignedAllocator() { }
AlignedAllocator() { }
/// Default copy constructor. It should be empty for stateless allocators.
/// Default copy constructor. It should be empty for stateless allocators.
///
AlignedAllocator(const AlignedAllocator&) { }
AlignedAllocator(const AlignedAllocator&) { }
/// Default rebinding constructor. It should be empty for stateless allocators.
/// Default rebinding constructor. It should be empty for stateless allocators.
///
template <typename U> AlignedAllocator(const AlignedAllocator<U, A>&) { }
template <typename U> AlignedAllocator(const AlignedAllocator<U, A>&) { }
/// Default destructor. It should be empty for stateless allocators.
/// Default destructor. It should be empty for stateless allocators.
///
~AlignedAllocator() { }
~AlignedAllocator() { }
/// The main allocator routine.
/// The following will be different for each allocator.
......@@ -107,83 +107,83 @@ public:
/// @exception std::length_error Integer overflow
/// @exception std::bad_alloc Memory allocation failure
///
T * allocate(const size_t n) const
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 << "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 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);
}
// AlignedAllocator wraps aligned free().
alignedFree(p);
}
/// The following will be the same for all allocators that ignore hints.
/// 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
template <typename U> T * allocate(const size_t n, const U * /* const hint */) const
{
return allocate(n);
}
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();
p->~T();
}
#ifdef _MSC_VER
#pragma warning(pop)
#pragma warning(pop)
#endif
#endif
......
BayesTest.cpp
View file @ 2589ab54
......@@ -43,11 +43,11 @@ double BayesTest::getGridParams(const std::vector<double>& aVars, const std::vec
// Calculating f(x_h|w)
// Order of site classes for iw or f(x_h|w):
// back fore #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
// back fore #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
//
for(unsigned int iw=0; iw < BEB_NUM_CAT; ++iw)
{
......@@ -62,15 +62,15 @@ double BayesTest::getGridParams(const std::vector<double>& aVars, const std::vec
omega_fg_is_one = omega_bg_is_one = true;
}
else if(iw < (BEB_N1D+1+BEB_N1D*BEB_N1D)) // class 2a: w0 w2
{
omega_bg = prior_params[2][(iw-BEB_N1D-1) / BEB_N1D];
omega_fg = prior_params[3][(iw-BEB_N1D-1) % BEB_N1D];
{
omega_bg = prior_params[2][(iw-BEB_N1D-1) / BEB_N1D];
omega_fg = prior_params[3][(iw-BEB_N1D-1) % BEB_N1D];
omega_fg_is_one = omega_bg_is_one = false;
}
else // class 2b: w1 w2
{
omega_bg = 1.;
omega_fg = prior_params[3][iw-BEB_N1D-1-BEB_N1D*BEB_N1D];
{
omega_bg = 1.;
omega_fg = prior_params[3][iw-BEB_N1D-1-BEB_N1D*BEB_N1D];
omega_fg_is_one = false; omega_bg_is_one = true;
}
......@@ -111,7 +111,7 @@ double BayesTest::getGridParams(const std::vector<double>& aVars, const std::vec
{
double p = likelihoods[site];
mPriors[iw*mNumSites+site] = (p > 0) ? log(p) : -184.2068074395237; // If p < 0 then the value in codeml.c is: log(1e-80);
}
}
}
double scale = 0.;
......@@ -200,7 +200,7 @@ void BayesTest::computeBEB(const std::vector<double>& aVars, size_t aFgBranch, c
}
}
// Calculate marginal prob of data, fX, and postpara[]. scale2 is scale.
// Calculate marginal prob of data, fX, and postpara[]. scale2 is scale.
if(mVerbose >= VERBOSE_ONLY_RESULTS) std::cout << std::endl << "Calculating f(X), the marginal probability of data." << std::endl;
double fX = 1.;
double scale2 = -1e300;
......@@ -213,7 +213,7 @@ void BayesTest::computeBEB(const std::vector<double>& aVars, size_t aFgBranch, c
double postp0p1[BEB_N1D*BEB_N1D];
for(unsigned int k=0; k < BEB_N1D*BEB_N1D; k++) postp0p1[k] = 1.;
for(unsigned int igrid=0; igrid < BEB_NGRID; ++igrid)
for(unsigned int igrid=0; igrid < BEB_NGRID; ++igrid)
{
// Get one point on the grid
unsigned int ip[BEB_DIMS];
......@@ -297,7 +297,7 @@ void BayesTest::computeBEB(const std::vector<double>& aVars, size_t aFgBranch, c
fX = log(fX)+scale2;
if(mVerbose >= VERBOSE_ONLY_RESULTS) std::cout << "log(fX) = " << (fX+scale1-BEB_DIMS*log(BEB_N1D*1.))
<< " Scales = " << scale1 << " " << scale2 << std::endl;
<< " Scales = " << scale1 << " " << scale2 << std::endl;
// Calculate posterior probabilities for sites. scale1 is scale factor
if(mVerbose >= VERBOSE_ONLY_RESULTS) std::cout << std::endl << "Calculating f(w|X), posterior probs of site classes." << std::endl;
......@@ -308,7 +308,7 @@ void BayesTest::computeBEB(const std::vector<double>& aVars, size_t aFgBranch, c
for(unsigned int j=0; j < BEB_DIMS; ++j) mSiteClassProb[j*mNumSites+site] = 1.;
for(unsigned int igrid=0; igrid < BEB_NGRID; ++igrid)
for(unsigned int igrid=0; igrid < BEB_NGRID; ++igrid)
{
double fh = 0.;
double fhk[BEB_DIMS];
......@@ -387,9 +387,9 @@ void BayesTest::printPositiveSelSites(size_t aFgBranch) const
// Set significance
const char* sig;
if(prob > TWO_STARS_PROB) sig = "**";
if(prob > TWO_STARS_PROB) sig = "**";
else if(prob > ONE_STAR_PROB) sig = "*";
else sig = "";
else sig = "";
std::cout << std::setw(6) << im->first + 1 << ' ' << std::fixed << std::setprecision(6) << prob << sig << std::endl;
}
......
BayesTest.h
View file @ 2589ab54
......@@ -10,7 +10,7 @@
/// The minimum value for class 2 sites probability to be a positive selection site.
///
const static double MIN_PROB = 0.50;
const static double MIN_PROB = 0.50;
const static double ONE_STAR_PROB = 0.95;
const static double TWO_STARS_PROB = 0.99;
......@@ -39,9 +39,9 @@ public:
/// 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
/// @author Mario Valle - Swiss National Supercomputing Centre (CSCS)
/// @date 2010-12-22 (initial version)
/// @version 1.1
///
class BayesTest
{
......@@ -96,18 +96,18 @@ private:
///
///@verbatim
/// Order of site classes for iw or f(x_h|w):
/// back fore num.sets
/// 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
/// 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.
/// @param[in] aScales The two scales ([0] bg; [1] fg) to rescale the branch lengths. They are computed in H1.
///
/// @return The computed scale.
///
......
BranchSiteModel.cpp
View file @ 2589ab54
This diff is collapsed.
BranchSiteModel.h
View file @ 2589ab54
......@@ -16,9 +16,9 @@ static const double THRESHOLD_FOR_LRT = 1.92072941034706202;
/// Common routines for testing the two hypotheses (H0 and H1).
///
/// @author Mario Valle - Swiss National Supercomputing Centre (CSCS)
/// @date 2010-12-23 (initial version)
/// @version 1.1
/// @author Mario Valle - Swiss National Supercomputing Centre (CSCS)
/// @date 2010-12-23 (initial version)
/// @version 1.1
///
class BranchSiteModel
{
......@@ -40,17 +40,17 @@ protected:
/// @param[in] aExtraDebug Extra parameter for testing during development
/// @param[in] aMaxIterations Maximum number of iterations for the maximization
///
BranchSiteModel(Forest& aForest,
size_t aNumBranches,
size_t aNumSites,
BranchSiteModel(Forest& aForest,
size_t aNumBranches,
size_t aNumSites,
unsigned int aSeed,
unsigned int aNumVariables,
bool aOnlyInitialStep,
bool aTrace,
bool aOnlyInitialStep,
bool aTrace,
unsigned int aOptAlgo,
double aDeltaValueForGradient,
double aRelativeError,
bool aNoParallel,
double aDeltaValueForGradient,
double aRelativeError,
bool aNoParallel,
unsigned int aVerbose,
unsigned int aExtraDebug,
unsigned int aMaxIterations,
......@@ -78,8 +78,8 @@ protected:
// mRelativeError(aRelativeError),
mFixedBranchLength(aFixedBranchLength),
mBranches(aNumBranches),
mSeed(aSeed),
mRelativeError(aRelativeError)
mSeed(aSeed),
mRelativeError(aRelativeError)
{
setLimits(aNumBranches, static_cast<size_t>(aNumVariables), aFixedBranchLength);
}
......@@ -345,8 +345,8 @@ protected:
unsigned int mMaxIterations; ///< Maximum number of iterations for the maximization
TreeAndSetsDependencies mDependencies; ///< The dependency list between trees to use in this run
bool mNoParallel; ///< True if no preparation for multithreading should be done
bool mFixedBranchLength; ///< True if branch lengths are kept fixed (not optimised)
std::vector<double> mBranches; ///< Variable with the branch lengths
bool mFixedBranchLength; ///< True if branch lengths are kept fixed (not optimised)
std::vector<double> mBranches; ///< Variable with the branch lengths
private:
unsigned int mSeed; ///< Random number generator seed to be used also by the optimizer
......@@ -357,9 +357,9 @@ private:
/// Null Hypothesis test.
///
/// @author Mario Valle - Swiss National Supercomputing Centre (CSCS)
/// @date 2010-12-23 (initial version)
/// @version 1.1
/// @author Mario Valle - Swiss National Supercomputing Centre (CSCS)
/// @date 2010-12-23 (initial version)
/// @version 1.1
///
///
class BranchSiteModelNullHyp : public BranchSiteModel
......@@ -434,8 +434,8 @@ private:
double combineSiteLikelihoods(void);
private:
TransitionMatrix mQw0; ///< Q matrix for the omega0 case
TransitionMatrix mQ1; ///< Q matrix for the omega1 == 1 case
TransitionMatrix mQw0; ///< Q matrix for the omega0 case
TransitionMatrix mQ1; ///< Q matrix for the omega1 == 1 case
ProbabilityMatrixSetH0 mSet; ///< Set of matrices used for the tree visits
ProbabilityMatrixSetH0 mSetForGradient; ///< Set of matrices used for the tree visits
double mPrevK; ///< Previous k value used to compute matrices
......@@ -448,9 +448,9 @@ private:
/// Alternate Hypothesis test.
///
/// @author Mario Valle - Swiss National Supercomputing Centre (CSCS)
/// @date 2010-12-23 (initial version)
/// @version 1.1
/// @author Mario Valle - Swiss National Supercomputing Centre (CSCS)
/// @date 2010-12-23 (initial version)
/// @version 1.1
///