Skip to content

Commit

Permalink
Remove some asserts that are false during unit tests.
Browse files Browse the repository at this point in the history
  • Loading branch information
@bredelings
bredelings committed Mar 29, 2021
1 parent 0c3e670 commit 7ad8cac
Show file tree
Hide file tree
Showing 3 changed files with 20 additions and 24 deletions.
6 changes: 1 addition & 5 deletions src/mcmc.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -150,7 +150,7 @@ void McmcMachinery::initializeMcmcChain(bool useIBD) {


void McmcMachinery::initializeHap() {
assert( currentHap_.size() == 0);
currentHap_.clear();
if ( this->dEploidIO_ -> initialHapWasGiven() ) {
this->currentHap_ = this->dEploidIO_->initialHap;
dout << "given initial hap ?" << endl;
Expand Down Expand Up @@ -196,22 +196,19 @@ double McmcMachinery::rBernoulli(double p) {


void McmcMachinery::initializeExpectedWsaf() {
assert( this->currentExpectedWsaf_.size() == 0);
this->currentExpectedWsaf_ = this->calcExpectedWsaf( this->currentProp_ );
assert( this->currentExpectedWsaf_.size() == this->nLoci_ );
this->cumExpectedWsaf_ = this->currentExpectedWsaf_;
}


void McmcMachinery::initializellk() {
assert( this->currentLLks_.size() == (size_t)0);
this->currentLLks_ = vector <double> (this->nLoci_, 0.0);
assert( this->currentLLks_.size() == this->nLoci_);
}


void McmcMachinery::initializeProp() {
assert( this->currentProp_.size() == (size_t)0 );
this->currentProp_ = ( this->dEploidIO_ -> initialPropWasGiven()) ?
this->dEploidIO_ ->initialProp:
this->titre2prop( this->currentTitre_ );
Expand All @@ -236,7 +233,6 @@ log_double_t McmcMachinery::calcPriorTitre(const vector <double> &tmpTitre) {

void McmcMachinery::initializeTitre() {
/* titre<-rnorm(initial.k, MN_LOG_TITRE, SD_LOG_TITRE); */
assert( currentTitre_.size() == 0);
currentTitre_ = vector <double> (this->kStrain_, 0.0);

if ( this->dEploidIO_->doUpdateProp() ) {
Expand Down
36 changes: 18 additions & 18 deletions src/updateHap.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -129,7 +129,7 @@ void UpdateSingleHap::calcBwdProbs() {
bwdLast[i] = 1.0;
}
(void)normalizeBySum(bwdLast);
assert(bwdProbs_.size() == 0 );
bwdProbs_.clear();
bwdProbs_.push_back(bwdLast);

int j = (this->nLoci_- 1);
Expand Down Expand Up @@ -162,7 +162,7 @@ void UpdateSingleHap::calcFwdBwdProbs() {
this->calcFwdProbs();
this->calcBwdProbs();

assert (this->fwdBwdProbs_.size() == 0);
this->fwdBwdProbs_.clear();
for ( size_t j = 0; j < this->nLoci_; j++ ) {
vector <double> fwdBwdTmp (this->nPanel_, 0.0);
for ( size_t i = 0 ; i < this->nPanel_; i++ ) {
Expand All @@ -176,8 +176,8 @@ void UpdateSingleHap::calcFwdBwdProbs() {

void UpdateSingleHap::calcExpectedWsaf( vector <double> & expectedWsaf, vector <double> &proportion, vector < vector <double> > &haplotypes ) {
//expected.WSAF.0 <- bundle$expected.WSAF - (bundle$prop[ws] * bundle$h[,ws]);
assert ( expectedWsaf0_.size() == 0);
assert ( expectedWsaf1_.size() == 0);
expectedWsaf0_.clear();
expectedWsaf1_.clear();
this->expectedWsaf0_ = vector <double> (expectedWsaf.begin()+this->segmentStartIndex_, expectedWsaf.begin()+(this->segmentStartIndex_+this->nLoci_));
size_t hapIndex = this->segmentStartIndex_;
for ( size_t i = 0; i < expectedWsaf0_.size(); i++ ) {
Expand Down Expand Up @@ -212,7 +212,7 @@ void UpdateSingleHap::buildEmission( double missCopyProb ) {
vector <double> t1u = vecSum(llk0_, log(missCopyProb));
vector <double> t2u = vecSum(llk1_, log(missCopyProb));

assert(emission_.size() == 0 );
emission_.clear();
for ( size_t i = 0; i < this->nLoci_; i++) {
vector <double> tmp ({t1omu[i], t2omu[i], t1u[i], t2u[i]});
double tmaxTmp = max_value(tmp);
Expand All @@ -226,7 +226,7 @@ void UpdateSingleHap::buildEmission( double missCopyProb ) {


void UpdateSingleHap::buildEmissionBasicVersion( double missCopyProb ) {
assert(emission_.size() == 0 );
emission_.clear();
for ( size_t i = 0; i < this->nLoci_; i++) {
vector <double> emissRow ({exp(llk0_[i])*(1.0-missCopyProb) + exp(llk1_[i])*missCopyProb,
exp(llk1_[i])*(1.0-missCopyProb) + exp(llk0_[i])*missCopyProb});
Expand All @@ -238,7 +238,7 @@ void UpdateSingleHap::buildEmissionBasicVersion( double missCopyProb ) {

void UpdateSingleHap::calcFwdProbs() {
size_t hapIndex = this->segmentStartIndex_;
assert ( this->fwdProbs_.size() == 0 );
this->fwdProbs_.clear();
vector <double> fwd1st (this->nPanel_, 0.0);
for ( size_t i = 0 ; i < this->nPanel_; i++) {
fwd1st[i] = this->emission_[0][this->panel_->content_[hapIndex][i]];
Expand Down Expand Up @@ -280,7 +280,7 @@ void UpdateSingleHap::calcHapLLKs( vector <double> &refCount,


void UpdateSingleHap::samplePaths() {
assert ( this->path_.size() == 0 );
this->path_.clear();
// Sample path at the last position
size_t pathTmp = sampleIndexGivenProp ( this->recombRg_, fwdProbs_.back() );
size_t contentIndex = this->segmentStartIndex_ + this->nLoci_ - 1;
Expand Down Expand Up @@ -313,7 +313,7 @@ void UpdateSingleHap::samplePaths() {


void UpdateSingleHap::addMissCopying( double missCopyProb ) {
assert( this->hap_.size() == 0 );
this->hap_.clear();
for ( size_t i = 0; i < this->nLoci_; i++) {
double tmpMax = max_value ( vector <double>({this->llk0_[i], this->llk1_[i]}));
vector <double> emissionTmp ({exp(this->llk0_[i]-tmpMax), exp(this->llk1_[i]-tmpMax)});
Expand All @@ -333,7 +333,7 @@ void UpdateSingleHap::addMissCopying( double missCopyProb ) {


void UpdateSingleHap::sampleHapIndependently( vector <double> &plaf ) {
assert( this->hap_.size() == 0 );
this->hap_.clear();
size_t plafIndex = this->segmentStartIndex_;
for ( size_t i = 0; i < this->nLoci_; i++) {
double tmpMax = max_value ( vector <double> ( {llk0_[i], llk1_[i]} ) ) ;
Expand Down Expand Up @@ -499,7 +499,7 @@ void UpdatePairHap:: buildEmission( double missCopyProb ) {
//exp( tmp.01.1-tmp.max ) + exp( tmp.01.2-tmp.max ) + exp( tmp.01.3-tmp.max ) + exp( tmp.01.4-tmp.max ),
//exp( tmp.11.1-tmp.max ) + exp( tmp.11.2-tmp.max ) + exp( tmp.11.3-tmp.max ) + exp( tmp.11.4-tmp.max ))

assert(this->emission_.size() == 0 );
this->emission_.clear();
for ( size_t i = 0; i < this->nLoci_; i++) {
vector <double> tmp ({tmp_00_1[i], tmp_00_2[i], tmp_00_3[i], tmp_00_4[i],
tmp_01_1[i], tmp_01_2[i], tmp_01_3[i], tmp_01_4[i],
Expand Down Expand Up @@ -542,7 +542,7 @@ vector <double> UpdatePairHap::computeColMarginalDist( vector < vector < double

void UpdatePairHap:: calcFwdProbs( bool forbidCopyFromSame ) {
size_t hapIndex = this->segmentStartIndex_;
assert ( this->fwdProbs_.size() == 0 );
this->fwdProbs_.clear();
vector < vector < double > > fwd1st;
for ( size_t i = 0 ; i < this->nPanel_; i++) { // Row of the matrix
size_t rowObs = (size_t)this->panel_->content_[0][i];
Expand Down Expand Up @@ -598,8 +598,8 @@ vector <size_t> UpdatePairHap::sampleMatrixIndex( vector < vector < double > > &


void UpdatePairHap::samplePaths() {
assert ( this->path1_.size() == 0 );
assert ( this->path2_.size() == 0 );
this->path1_.clear();
this->path2_.clear();

vector <size_t> tmpPath = sampleMatrixIndex(fwdProbs_[nLoci_-1]);
size_t rowI = tmpPath[0];
Expand Down Expand Up @@ -675,8 +675,8 @@ void UpdatePairHap::samplePaths() {


void UpdatePairHap::addMissCopying( double missCopyProb ) {
assert( this->hap1_.size() == 0 );
assert( this->hap2_.size() == 0 );
this->hap1_.clear();
this->hap2_.clear();

for ( size_t i = 0; i < this->nLoci_; i++) {
double tmpMax = max_value ( vector <double>({this->llk00_[i], this->llk01_[i], this->llk10_[i], this->llk11_[i]}));
Expand Down Expand Up @@ -714,8 +714,8 @@ void UpdatePairHap::addMissCopying( double missCopyProb ) {


void UpdatePairHap::sampleHapIndependently(vector <double> &plaf) {
assert( this->hap1_.size() == 0 );
assert( this->hap2_.size() == 0 );
this->hap1_.clear();
this->hap2_.clear();

size_t plafIndex = this->segmentStartIndex_;
for ( size_t i = 0; i < this->nLoci_; i++) {
Expand Down
2 changes: 1 addition & 1 deletion src/utility.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,7 @@ vector <double> calcLLKs(const vector <double> &refCount,
const vector <double> &expectedWsaf,
size_t firstIndex, size_t length,
double fac, double err) {
assert(expectedWsaf.size() == length);
assert(length <= expectedWsaf.size());
vector <double> tmpLLKs(length, 0.0);
size_t index = firstIndex;
for (size_t i = 0; i < length; i++) {
Expand Down

0 comments on commit 7ad8cac

Please sign in to comment.