From 6eeaf1158c6ed5d908a876f3f29e02a00d0139e4 Mon Sep 17 00:00:00 2001 From: Joe Zhu Date: Mon, 28 Dec 2020 16:05:15 +0000 Subject: [PATCH] check style --- .ci/checkedList | 2 + .ci/mycppfileList | 6 +- src/panel.cpp | 192 ++++++++++++++++++++++++++-------------------- src/panel.hpp | 92 +++++++++++----------- src/utility.cpp | 137 +++++++++++++++++---------------- src/utility.hpp | 4 +- 6 files changed, 235 insertions(+), 198 deletions(-) diff --git a/.ci/checkedList b/.ci/checkedList index c3896e0a..69110f6a 100644 --- a/.ci/checkedList +++ b/.ci/checkedList @@ -16,3 +16,5 @@ tests/unittest/test_vcfReader.cpp src/chooseK.hpp src/chooseK.cpp tests/unittest/test_dEploidIO.cpp +src/panel.hpp +src/panel.cpp diff --git a/.ci/mycppfileList b/.ci/mycppfileList index f1a58733..9526f982 100644 --- a/.ci/mycppfileList +++ b/.ci/mycppfileList @@ -1,12 +1,8 @@ src/dEploidIO.cpp src/dEploidIO.hpp src/mcmc.cpp -src/panel.cpp -src/panel.hpp src/updateHap.cpp src/updateHap.hpp -src/utility.cpp -src/utility.hpp src/export/dEploidIOExport.cpp src/export/dEploidIOExportPosteriorProb.cpp src/export/writeMcmcRelated.cpp @@ -18,3 +14,5 @@ tests/unittest/test_txtReader.cpp tests/unittest/test_updatePairHap.cpp tests/unittest/test_updateSingleHap.cpp tests/unittest/test_utilities.cpp +src/utility.hpp +src/utility.cpp diff --git a/src/panel.cpp b/src/panel.cpp index 0aec2801..0cebfdae 100644 --- a/src/panel.cpp +++ b/src/panel.cpp @@ -64,7 +64,7 @@ Panel::Panel(vector < double > pRec, } -Panel::Panel(const Panel ©From){ +Panel::Panel(const Panel ©From) { nLoci_ = copyFrom.nLoci_; chrom_ = vector (copyFrom.chrom_.begin(), copyFrom.chrom_.end()); @@ -99,33 +99,36 @@ Panel::Panel(const Panel ©From){ inbreedingPanelSize_ = copyFrom.inbreedingPanelSize_; this->fileName_ = copyFrom.fileName_; - for (size_t i = 0; i < copyFrom.header_.size(); i++){ + for (size_t i = 0; i < copyFrom.header_.size(); i++) { header_.push_back(copyFrom.header_[i]); } } -void Panel::readFromFile( const char inchar[] ) { - this->readFromFileBase( inchar ); - this->setTruePanelSize( this->nInfoLines_ ); - this->setInbreedingPanelSize( this->truePanelSize() ); +void Panel::readFromFile(const char inchar[]) { + this->readFromFileBase(inchar); + this->setTruePanelSize(this->nInfoLines_); + this->setInbreedingPanelSize(this->truePanelSize()); } -void Panel::checkForExceptions( size_t nLoci, string panelFileName ) { - if ( this->content_.size() != nLoci ) { - throw LociNumberUnequal( panelFileName ); +void Panel::checkForExceptions(size_t nLoci, string panelFileName) { + if (this->content_.size() != nLoci) { + throw LociNumberUnequal(panelFileName); } - if ( this->pRec_.size() != nLoci ) { - throw LociNumberUnequal( panelFileName ); + if (this->pRec_.size() != nLoci) { + throw LociNumberUnequal(panelFileName); } return; } -void Panel::computeRecombProbs( double averageCentimorganDistance, double G, bool useConstRecomb, double constRecombProb, bool forbidCopyFromSame ) { +void Panel::computeRecombProbs(double averageCentimorganDistance, + double G, bool useConstRecomb, double constRecombProb, + bool forbidCopyFromSame) { + pRec_.clear(); pRecEachHap_.clear(); pNoRec_.clear(); @@ -133,81 +136,91 @@ void Panel::computeRecombProbs( double averageCentimorganDistance, double G, boo pRecNoRec_.clear(); pNoRecNoRec_.clear(); - assert(pRec_.size() == 0 ); - assert(pRecEachHap_.size() == 0 ); - assert(pNoRec_.size() == 0 ); - assert(pRecRec_.size() == 0 ); - assert(pRecNoRec_.size() == 0 ); - assert(pNoRecNoRec_.size() == 0 ); + assert(pRec_.size() == 0); + assert(pRecEachHap_.size() == 0); + assert(pNoRec_.size() == 0); + assert(pRecRec_.size() == 0); + assert(pRecNoRec_.size() == 0); + assert(pNoRecNoRec_.size() == 0); double averageMorganDistance = averageCentimorganDistance * 100; double geneticDistance; double rho; - double nPanelDouble = (double)this->truePanelSize(); + double nPanelDouble = static_cast(this->truePanelSize()); double nPanlelMinus1 = nPanelDouble - 1.0; - for ( size_t i = 0; i < this->position_.size(); i++) { - for ( size_t j = 1; j < this->position_[i].size(); j++) { - geneticDistance = (double)(this->position_[i][j] - this->position_[i][j-1])/averageMorganDistance ; - //rho = geneticDistance * 2 * Ne; + for (size_t i = 0; i < this->position_.size(); i++) { + for (size_t j = 1; j < this->position_[i].size(); j++) { + geneticDistance = static_cast( + this->position_[i][j] - this->position_[i][j-1]) / + averageMorganDistance; + // rho = geneticDistance * 2 * Ne; rho = geneticDistance * G; - double pRecTmp = ( useConstRecomb ) ? constRecombProb : 1.0 - exp(-rho); - this->pRec_.push_back( pRecTmp ); + double pRecTmp = (useConstRecomb) ? + constRecombProb : 1.0 - exp(-rho); + this->pRec_.push_back(pRecTmp); double pRecEachHapTmp = pRecTmp / nPanelDouble; - this->pRecEachHap_.push_back( pRecTmp / nPanelDouble ); + this->pRecEachHap_.push_back(pRecTmp / nPanelDouble); double pNoRecTmp = 1.0 - pRecTmp; - this->pNoRec_.push_back( pNoRecTmp ); + this->pNoRec_.push_back(pNoRecTmp); - double secondPRecEachHapTmp = ( forbidCopyFromSame ) ? (pRecTmp / nPanlelMinus1) : pRecEachHapTmp; // allowing copy from the same + double secondPRecEachHapTmp = (forbidCopyFromSame) ? + (pRecTmp / nPanlelMinus1) : + pRecEachHapTmp; // allowing copy from the same - this->pRecRec_.push_back ( pRecEachHapTmp * secondPRecEachHapTmp ); - this->pRecNoRec_.push_back ( secondPRecEachHapTmp * pNoRecTmp ); - this->pNoRecNoRec_.push_back ( pNoRecTmp * pNoRecTmp ); + this->pRecRec_.push_back(pRecEachHapTmp * secondPRecEachHapTmp); + this->pRecNoRec_.push_back(secondPRecEachHapTmp * pNoRecTmp); + this->pNoRecNoRec_.push_back(pNoRecTmp * pNoRecTmp); } this->pRec_.push_back(1.0); - this->pRecEachHap_.push_back( 1.0 / nPanelDouble ); + this->pRecEachHap_.push_back(1.0 / nPanelDouble); this->pNoRec_.push_back(0.0); - this->pRecRec_.push_back ( ( ( forbidCopyFromSame ) ? (1.0 / nPanelDouble / nPanlelMinus1) : (1.0 / nPanelDouble / nPanelDouble) ) ); - this->pRecNoRec_.push_back ( 0.0 ); - this->pNoRecNoRec_.push_back ( 0.0 ); + this->pRecRec_.push_back(((forbidCopyFromSame) ? + (1.0 / nPanelDouble / nPanlelMinus1) : + (1.0 / nPanelDouble / nPanelDouble))); + this->pRecNoRec_.push_back(0.0); + this->pNoRecNoRec_.push_back(0.0); } - assert(pRec_.size() == this->nLoci_ ); - assert(pRecEachHap_.size() == this->nLoci_ ); - assert(pNoRec_.size() == this->nLoci_ ); - assert(pRecRec_.size() == this->nLoci_ ); - assert(pRecNoRec_.size() == this->nLoci_ ); - assert(pNoRecNoRec_.size() == this->nLoci_ ); + assert(pRec_.size() == this->nLoci_); + assert(pRecEachHap_.size() == this->nLoci_); + assert(pNoRec_.size() == this->nLoci_); + assert(pRecRec_.size() == this->nLoci_); + assert(pRecNoRec_.size() == this->nLoci_); + assert(pNoRecNoRec_.size() == this->nLoci_); } void Panel::buildExamplePanel1() { this->chrom_ = vector ({"Pf3D7_01_v3"}); - this->position_.push_back( vector ({93157, 94422, 94459, 94487, 95518, 95632, 95641 }) ); + this->position_.push_back(vector ( + {93157, 94422, 94459, 94487, 95518, 95632, 95641})); this->indexOfChromStarts_ = vector ({0}); this->buildExamplePanelContent(); } void Panel::buildExamplePanel2() { - this->chrom_ = vector ({"Pf3D7_01_v3", "Pf3D7_02_v3", "Pf3D7_03_v3"}); - this->position_.push_back( vector ({93157}) ); - this->position_.push_back( vector ({94422, 94459, 94487, 95518, 95632}) ); - this->position_.push_back( vector ({95641 }) ); + this->chrom_ = vector ( + {"Pf3D7_01_v3", "Pf3D7_02_v3", "Pf3D7_03_v3"}); + this->position_.push_back(vector ({93157})); + this->position_.push_back(vector ( + {94422, 94459, 94487, 95518, 95632})); + this->position_.push_back(vector ({95641})); this->indexOfChromStarts_ = vector ({0, 1, 6}); this->buildExamplePanelContent(); } void Panel::buildExamplePanelContent() { - this->content_.push_back( vector ({0,0,0,1}) ); - this->content_.push_back( vector ({0,0,0,1}) ); - this->content_.push_back( vector ({0,0,0,1}) ); - this->content_.push_back( vector ({0,0,0,1}) ); - this->content_.push_back( vector ({0,1,1,0}) ); - this->content_.push_back( vector ({0,0,1,0}) ); - this->content_.push_back( vector ({0,0,1,0}) ); + this->content_.push_back(vector ({0, 0, 0, 1})); + this->content_.push_back(vector ({0, 0, 0, 1})); + this->content_.push_back(vector ({0, 0, 0, 1})); + this->content_.push_back(vector ({0, 0, 0, 1})); + this->content_.push_back(vector ({0, 1, 1, 0})); + this->content_.push_back(vector ({0, 0, 1, 0})); + this->content_.push_back(vector ({0, 0, 1, 0})); this->nLoci_ = this->content_.size(); this->nInfoLines_ = this->content_.back().size(); this->setTruePanelSize(this->nInfoLines_); @@ -218,41 +231,47 @@ void Panel::buildExamplePanelContent() { void Panel::initializeUpdatePanel(size_t inbreedingPanelSizeSetTo) { this->setInbreedingPanelSize(inbreedingPanelSizeSetTo); - // If allows inbreeding, update reference panel by including strain haplotypes + // If allows inbreeding, + // update reference panel by including strain haplotypes dout << "************* Allow inbreeding ************" << endl; dout << "** Initialize inbreeding reference panel **" << endl; - if ( this->truePanelSize() == this->inbreedingPanelSize()) { + if (this->truePanelSize() == this->inbreedingPanelSize()) { return; } for ( size_t siteI = 0; siteI < this->content_.size(); siteI++ ) { - for ( size_t panelStrainJ = this->truePanelSize() ; panelStrainJ < this->inbreedingPanelSize(); panelStrainJ++ ) { - this->content_[siteI].push_back( 1 ); + for (size_t panelStrainJ = this->truePanelSize(); + panelStrainJ < this->inbreedingPanelSize(); panelStrainJ++) { + this->content_[siteI].push_back(1); } assert(inbreedingPanelSizeSetTo == this->content_[siteI].size()); } } -void Panel::updatePanelWithHaps(size_t inbreedingPanelSizeSetTo, size_t excludedStrain, vector < vector > & haps) { +void Panel::updatePanelWithHaps(size_t inbreedingPanelSizeSetTo, + size_t excludedStrain, const vector < vector > & haps) { this->setInbreedingPanelSize(inbreedingPanelSizeSetTo); - // If allows inbreeding, update reference panel by including strain haplotypes + // If allows inbreeding, + // update reference panel by including strain haplotypes dout << "*************** Allow inbreeding **************" << endl; - dout << "*** Update reference panel without strain " << excludedStrain << " ***" << endl; + dout << "*** Update reference panel without strain " << + excludedStrain << " ***" << endl; - if ( this->truePanelSize() == this->inbreedingPanelSize()) { + if (this->truePanelSize() == this->inbreedingPanelSize()) { return; } - for ( size_t siteI = 0; siteI < this->content_.size(); siteI++ ) { + for (size_t siteI = 0; siteI < this->content_.size(); siteI++) { size_t shiftAfter = this->inbreedingPanelSize(); - for ( size_t panelStrainJ = this->truePanelSize() ; panelStrainJ < this->inbreedingPanelSize(); panelStrainJ++ ) { + for (size_t panelStrainJ = this->truePanelSize(); + panelStrainJ < this->inbreedingPanelSize(); panelStrainJ++) { size_t strainIndex = panelStrainJ - this->truePanelSize(); - if ( strainIndex == excludedStrain ) { + if (strainIndex == excludedStrain) { shiftAfter = panelStrainJ; } @@ -265,14 +284,13 @@ void Panel::updatePanelWithHaps(size_t inbreedingPanelSizeSetTo, size_t excluded } -void Panel::trimVec(vector &vec, const vector &idx) { +void Panel::trimVec(vector vec, const vector &idx) { vector ret; - for (auto const& value : idx){ + for (auto const& value : idx) { ret.push_back(vec[value]); } - //return ret; vec.clear(); - for (auto const& value : ret){ + for (auto const& value : ret) { vec.push_back(value); } } @@ -309,32 +327,36 @@ void Panel::findAndKeepMarkersGivenIndexHalf( this->trimVec(this->pNoRecNoRec_, givenIndex); } -void IBDrecombProbs::computeRecombProbs( double averageCentimorganDistance, double G, bool useConstRecomb, double constRecombProb ) { - assert(pRec_.size() == 0 ); - assert(pNoRec_.size() == 0 ); +void IBDrecombProbs::computeRecombProbs(double averageCentimorganDistance, + double G, bool useConstRecomb, double constRecombProb) { + assert(pRec_.size() == 0); + assert(pNoRec_.size() == 0); double averageMorganDistance = averageCentimorganDistance * 100; double geneticDistance; double rho; - for ( size_t i = 0; i < this->position_.size(); i++) { - for ( size_t j = 1; j < this->position_[i].size(); j++) { - geneticDistance = (double)(this->position_[i][j] - this->position_[i][j-1])/averageMorganDistance ; - //rho = geneticDistance * 2 * Ne; + for (size_t i = 0; i < this->position_.size(); i++) { + for (size_t j = 1; j < this->position_[i].size(); j++) { + geneticDistance = + static_cast(this->position_[i][j] - + this->position_[i][j-1]) / averageMorganDistance; + // rho = geneticDistance * 2 * Ne; rho = geneticDistance * G; - double pRecTmp = ( useConstRecomb ) ? constRecombProb : (1.0 - exp(-rho)); - this->pRec_.push_back( pRecTmp ); + double pRecTmp = (useConstRecomb) ? + constRecombProb : (1.0 - exp(-rho)); + this->pRec_.push_back(pRecTmp); double pNoRecTmp = 1.0 - pRecTmp; - this->pNoRec_.push_back( pNoRecTmp ); + this->pNoRec_.push_back(pNoRecTmp); } this->pRec_.push_back(1.0); this->pNoRec_.push_back(0.0); } - assert(pRec_.size() == this->nLoci_ ); - assert(pNoRec_.size() == this->nLoci_ ); + assert(pRec_.size() == this->nLoci_); + assert(pNoRec_.size() == this->nLoci_); } -//vector> outtrans(out[0].size(), - //vector(out.size())); - //for (size_t i = 0; i < out.size(); ++i) - //for (size_t j = 0; j < out[0].size(); ++j) - //outtrans[j][i] = out[i][j]; +// vector> outtrans(out[0].size(), + // vector(out.size())); + // for (size_t i = 0; i < out.size(); ++i) + // for (size_t j = 0; j < out[0].size(); ++j) + // outtrans[j][i] = out[i][j]; diff --git a/src/panel.hpp b/src/panel.hpp index b9bfeff4..0cece9a2 100644 --- a/src/panel.hpp +++ b/src/panel.hpp @@ -26,24 +26,27 @@ #ifndef PANEL #define PANEL +#include +#include + #include "txtReader.hpp" #include "exceptions.hpp" class Panel: public TxtReader{ -#ifdef UNITTEST - friend class TestPanel; - friend class TestInitialHaplotypes; - friend class TestUpdateHap; - friend class TestUpdatePairHap; - friend class TestUpdateSingleHap; -#endif - friend class McmcMachinery; - friend class UpdateSingleHap; - friend class UpdatePairHap; - friend class UpdateHap; - friend class DEploidIO; - friend class InitialHaplotypes; - private: + #ifdef UNITTEST + friend class TestPanel; + friend class TestInitialHaplotypes; + friend class TestUpdateHap; + friend class TestUpdatePairHap; + friend class TestUpdateSingleHap; + #endif + friend class McmcMachinery; + friend class UpdateSingleHap; + friend class UpdatePairHap; + friend class UpdateHap; + friend class DEploidIO; + friend class InitialHaplotypes; + private: // Members vector < double > pRec_; // Used in update single haplotype @@ -55,10 +58,12 @@ class Panel: public TxtReader{ vector < double > pNoRecNoRec_; // pNoRec * pNoRec; size_t truePanelSize_; - void setTruePanelSize ( const size_t setTo ){ this->truePanelSize_ = setTo; } + void setTruePanelSize(const size_t setTo) { + this->truePanelSize_ = setTo; } size_t inbreedingPanelSize_; - void setInbreedingPanelSize ( const size_t setTo ){ this->inbreedingPanelSize_ = setTo; } + void setInbreedingPanelSize(const size_t setTo) { + this->inbreedingPanelSize_ = setTo; } size_t inbreedingPanelSize() const { return this->inbreedingPanelSize_; } size_t truePanelSize() const { return this->truePanelSize_; } @@ -72,57 +77,58 @@ class Panel: public TxtReader{ vector < vector < double > > content, vector < string > header); Panel(const Panel ©From); - - //Panel(const char inchar[] ); + // Panel(const char inchar[]); // Methods - void readFromFile( const char inchar[] ); - void computeRecombProbs( double averageCentimorganDistance, double Ne, bool useConstRecomb, double constRecombProb, bool forbidCopyFromSame ); - void checkForExceptions( size_t nLoci, string panelFileName ); - void initializeUpdatePanel( size_t inbreedingPanelSizeSetTo); - void updatePanelWithHaps( size_t inbreedingPanelSizeSetTo, size_t excludedStrain, vector < vector > & haps); + void readFromFile(const char inchar[]); + void computeRecombProbs(double averageCentimorganDistance, double Ne, + bool useConstRecomb, double constRecombProb, bool forbidCopyFromSame); + void checkForExceptions(size_t nLoci, string panelFileName); + void initializeUpdatePanel(size_t inbreedingPanelSizeSetTo); + void updatePanelWithHaps(size_t inbreedingPanelSizeSetTo, + size_t excludedStrain, const vector < vector > & haps); void print(); void buildExamplePanelContent(); void buildExamplePanel1(); void buildExamplePanel2(); - void trimVec(vector &vec, const vector &idx); - //void findWhoToBeKeptGivenIndex(const vector & givenIndex); + void trimVec(vector vec, const vector &idx); + // void findWhoToBeKeptGivenIndex(const vector & givenIndex); void findAndKeepMarkersGivenIndex(const vector & givenIndex); void findAndKeepMarkersGivenIndexHalf(const vector & givenIndex); - public: - virtual ~Panel() {}; + public: + virtual ~Panel() {} }; class InitialHaplotypes: public Panel{ -#ifdef UNITTEST - friend class TestInitialHaplotypes; -#endif - friend class DEploidIO; - InitialHaplotypes():Panel(){} + #ifdef UNITTEST + friend class TestInitialHaplotypes; + #endif + friend class DEploidIO; + InitialHaplotypes():Panel() {} ~InitialHaplotypes(){} }; class IBDrecombProbs: public VariantIndex{ - friend class IBDpath; + friend class IBDpath; + #ifdef UNITTEST + friend class TestIBDpath; + #endif -#ifdef UNITTEST - friend class TestIBDpath; -#endif - - private: + private: vector < double > pRec_; - vector < double > pNoRec_; // = 1.0 - pRec; + vector < double > pNoRec_; // = 1.0 - pRec; - void computeRecombProbs( double averageCentimorganDistance, double Ne, bool useConstRecomb, double constRecombProb ); + void computeRecombProbs(double averageCentimorganDistance, double Ne, + bool useConstRecomb, double constRecombProb); - public: - IBDrecombProbs():VariantIndex(){}; - IBDrecombProbs(vector < vector < int> > position, size_t nLoci){ + public: + IBDrecombProbs():VariantIndex() {} + IBDrecombProbs(vector < vector < int> > position, size_t nLoci) { this->position_ = position; this->nLoci_ = nLoci; } diff --git a/src/utility.cpp b/src/utility.cpp index 26148afd..cadfb862 100644 --- a/src/utility.cpp +++ b/src/utility.cpp @@ -23,16 +23,18 @@ * */ +#include // std::distance +#include // find + #include "utility.hpp" -#include // std::distance -#include // find -#include "codeCogs/loggammasum.h" // which includes log_gamma.h +#include "codeCogs/loggammasum.h" // which includes log_gamma.h #include "codeCogs/gamma.h" #include "codeCogs/logbeta.h" - -double normal_pdf(double x, double m, double s) { // See http://stackoverflow.com/questions/10847007/using-the-gaussian-probability-density-function-in-c +// See http://stackoverflow.com/questions/10847007/ +// using-the-gaussian-probability-density-function-in-c +double normal_pdf(double x, double m, double s) { static const double inv_sqrt_2pi = 0.3989422804014327; double a = (x - m) / s; @@ -40,37 +42,37 @@ double normal_pdf(double x, double m, double s) { // See http://stackoverflow.co } -double min_value( vector x){ - assert( x.size() > 0 ); +double min_value(vector x) { + assert(x.size() > 0); auto tmpMaxIt = std::min_element(std::begin(x), std::end(x)); return *tmpMaxIt; } -double max_value( vector x){ - assert( x.size() > 0 ); +double max_value(vector x) { + assert(x.size() > 0); auto tmpMaxIt = std::max_element(std::begin(x), std::end(x)); return *tmpMaxIt; } -vector computeCdf (const vector & dist ){ +vector computeCdf(const vector & dist) { vector cdf; double cumsum = 0; - for (double p : dist){ + for (double p : dist) { cumsum += p; - cdf.push_back( cumsum ); + cdf.push_back(cumsum); } - assert( cdf.size() == dist.size() ); - //assert( cdf.back() == 1.0 ); + assert(cdf.size() == dist.size()); + // assert(cdf.back() == 1.0); return cdf; } -double sumOfMat(const vector > & matrix ){ +double sumOfMat(const vector > & matrix) { double tmp = 0.0; - for (auto const& array: matrix){ - for (auto const& value: array){ + for (auto const& array : matrix) { + for (auto const& value : array) { tmp += value; } } @@ -78,26 +80,27 @@ double sumOfMat(const vector > & matrix ){ } -void normalizeBySum(vector & array ){ +void normalizeBySum(vector & array ) { double sumOfArray = sumOfVec(array); - for( vector::iterator it = array.begin(); it != array.end(); ++it) { + for (vector::iterator it = array.begin(); it != array.end(); ++it) { *it /= sumOfArray; } } -void normalizeByMax(vector & array ){ +void normalizeByMax(vector & array ) { double maxOfArray = max_value(array); - for( vector::iterator it = array.begin(); it != array.end(); ++it) { + for (vector::iterator it = array.begin(); it != array.end(); ++it) { *it /= maxOfArray; } } -void normalizeBySumMat(vector > & matrix ){ +void normalizeBySumMat(vector > & matrix) { double tmpsum = sumOfMat(matrix); - for( size_t i = 0; i < matrix.size(); i++ ){ - for( vector::iterator it = matrix[i].begin(); it != matrix[i].end(); ++it) { + for (size_t i = 0; i < matrix.size(); i++) { + for (vector::iterator it = matrix[i].begin(); + it != matrix[i].end(); ++it) { *it /= tmpsum; } } @@ -108,46 +111,50 @@ vector calcLLKs(const vector &refCount, const vector &altCount, const vector &expectedWsaf, size_t firstIndex, size_t length, - double fac, double err){ - assert ( expectedWsaf.size() == length ); - vector tmpLLKs (length, 0.0); + double fac, double err) { + assert(expectedWsaf.size() == length); + vector tmpLLKs(length, 0.0); size_t index = firstIndex; - for ( size_t i = 0; i < length; i++ ){ - assert (expectedWsaf[i] >= 0); - //assert (expectedWsaf[i] <= 1); - tmpLLKs[i] = calcLLK( refCount[index], - altCount[index], - expectedWsaf[i], - err, - fac); + for (size_t i = 0; i < length; i++) { + assert(expectedWsaf[i] >= 0); + // assert (expectedWsaf[i] <= 1); + tmpLLKs[i] = calcLLK(refCount[index], + altCount[index], + expectedWsaf[i], + err, + fac); index++; } return tmpLLKs; } -double calcLLK( double ref, double alt, double unadjustedWsaf, double err, double fac ) { +double calcLLK(double ref, double alt, double unadjustedWsaf, double err, + double fac) { double adjustedWsaf = unadjustedWsaf+err*(1-2*unadjustedWsaf); - double llk = Maths::Special::Gamma::logBeta(alt+adjustedWsaf*fac, ref+(1-adjustedWsaf)*fac) - Maths::Special::Gamma::logBeta(adjustedWsaf*fac,(1-adjustedWsaf)*fac); + double llk = Maths::Special::Gamma::logBeta( + alt+adjustedWsaf*fac, ref+(1-adjustedWsaf)*fac) - + Maths::Special::Gamma::logBeta(adjustedWsaf*fac, (1-adjustedWsaf)*fac); return llk; } -size_t sampleIndexGivenProp ( RandomGenerator* rg, vector proportion ){ +size_t sampleIndexGivenProp(RandomGenerator* rg, vector proportion) { #ifndef NDEBUG - auto biggest = std::max_element(std::begin(proportion), std::end(proportion)); + auto biggest = std::max_element( + std::begin(proportion), std::end(proportion)); return std::distance(proportion.begin(), biggest); #else vector tmpIndex; - for ( size_t i = 0; i < proportion.size(); i++ ){ + for ( size_t i = 0; i < proportion.size(); i++ ) { tmpIndex.push_back(i); } vector tmpCdf = computeCdf(proportion); double u = rg->sample(); size_t i = 0; - for ( ; i < tmpCdf.size() ; i++){ - if ( u < tmpCdf[i] ){ + for ( ; i < tmpCdf.size() ; i++) { + if ( u < tmpCdf[i] ) { break; } } @@ -156,10 +163,10 @@ size_t sampleIndexGivenProp ( RandomGenerator* rg, vector proportion ){ } -vector reshapeMatToVec(vector < vector > &Mat ){ +vector reshapeMatToVec(const vector < vector > &Mat) { vector tmp; - for (auto const& array: Mat){ - for (auto const& value: array){ + for (auto const& array : Mat) { + for (auto const& value : array) { tmp.push_back(value); } } @@ -167,20 +174,21 @@ vector reshapeMatToVec(vector < vector > &Mat ){ } -double betaPdf(double x, double a, double b){ - assert(x>=0 && x<=1); - assert(a>=0); - assert(b>=0); - double p = Maths::Special::Gamma::gamma(a + b) / (Maths::Special::Gamma::gamma(a) * Maths::Special::Gamma::gamma(b)); +double betaPdf(double x, double a, double b) { + assert(x >= 0 && x <= 1); + assert(a >= 0); + assert(b >= 0); + double p = Maths::Special::Gamma::gamma(a + b) / + (Maths::Special::Gamma::gamma(a) * Maths::Special::Gamma::gamma(b)); double q = pow(1 - x, b - 1) * pow(x, a - 1); return p * q; } -double logBetaPdf(double x, double a, double b){ - assert(x>=0 && x<=1); - assert(a>=0); - assert(b>=0); +double logBetaPdf(double x, double a, double b) { + assert(x >= 0 && x <= 1); + assert(a >= 0); + assert(b >= 0); double ret = Maths::Special::Gamma::log_gamma(a+b) - Maths::Special::Gamma::log_gamma(a) - Maths::Special::Gamma::log_gamma(b) + @@ -189,22 +197,23 @@ double logBetaPdf(double x, double a, double b){ } -double binomialPdf(int s, int n, double p){ - assert(p>=0 && p<=1); - double ret=n_choose_k(n, s); - ret *= pow(p, (double)s); - ret *= pow((1.0-p), (double)(n-s)); +double binomialPdf(int s, int n, double p) { + assert(p >= 0 && p <= 1); + double ret = n_choose_k(n, s); + ret *= pow(p, static_cast(s)); + ret *= pow((1.0-p), static_cast(n-s)); return ret; } -//double betaDistConst( double a , double b){ - //double ret = Maths::Special::Gamma::gamma(a + b) / (Maths::Special::Gamma::gamma(a) * Maths::Special::Gamma::gamma(b)); - //return ret; -//} +// double betaDistConst( double a , double b) { + // double ret = Maths::Special::Gamma::gamma(a + b) / + // (Maths::Special::Gamma::gamma(a) * Maths::Special::Gamma::gamma(b)); + // return ret; +// } -double rBeta(double alpha, double beta, RandomGenerator* rg){ +double rBeta(double alpha, double beta, RandomGenerator* rg) { double mxAt = (alpha-1.0)/(alpha+beta-2.0); double mx = betaPdf(mxAt, alpha, beta); double y, u; diff --git a/src/utility.hpp b/src/utility.hpp index ffb12209..a90ebc8e 100644 --- a/src/utility.hpp +++ b/src/utility.hpp @@ -78,7 +78,7 @@ vector vecProd(const vector &vecA, template -T sumOfVec(vector & array ) { +T sumOfVec(const vector & array ) { T tmp = 0; for (auto const& value : array) { tmp += value; @@ -121,7 +121,7 @@ vector calcLLKs(const vector &refCount, double calcLLK(double ref, double alt, double unadjustedWsaf, double err, double fac); size_t sampleIndexGivenProp(RandomGenerator* rg, vector proportion); -vector reshapeMatToVec(vector < vector > &Mat); +vector reshapeMatToVec(const vector < vector > &Mat); double betaPdf(double x, double a, double b); double logBetaPdf(double x, double a, double b); double binomialPdf(int s, int n, double p);