diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..28791f1 --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +cmake-build-debug/ +data/ \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..f099d73 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,27 @@ +cmake_minimum_required(VERSION 3.14) +project(selection) + +set(CMAKE_CXX_STANDARD 14) +set(CMAKE_CXX_STANDARD_REQUIRED YES) + +include_directories(.) + +add_executable(selection + main.cpp + MbRandom.cpp + MbRandom.h + mcmc.cpp + mcmc.h + measure.cpp + measure.h + param.cpp + param.h + path.cpp + path.h + popsize.cpp + popsize.h + settings.cpp + settings.h) + +find_package(GSL REQUIRED) +target_link_libraries(selection GSL::gsl GSL::gslcblas) \ No newline at end of file diff --git a/mcmc.cpp b/mcmc.cpp index 2047948..d48f04b 100644 --- a/mcmc.cpp +++ b/mcmc.cpp @@ -92,7 +92,7 @@ void mcmc::no_linked_sites(settings& mySettings) { for (int i = 0; i < sample_time_vec.size()-1; i++) { if (sample_time_vec[i]->get_oldest() < sample_time_vec[i]->get_youngest()) { if (!mySettings.get_infer_age()) { - std::cout << "ERROR: Cannot have uncertain times without inferring allele age. Will be fixed in the future" << std::endl; + std::cerr << "ERROR: Cannot have uncertain times without inferring allele age. Will be fixed in the future" << std::endl; exit(1); } sample_time_vec[i]->set_path(curParamPath); @@ -182,14 +182,14 @@ void mcmc::no_linked_sites(settings& mySettings) { curlnL = compute_lnL_sample_only(curPath); double LLRatio = curlnL-oldlnL; - if (curlnL != curlnL || oldlnL != oldlnL) { - std::cout << "ERROR: likelihood is nan!" << std::endl; - std::cout << "Generation = " << gen << std::endl; - std::cout << "Proposal = " << curProp << std::endl; - std::cout << "curlnL = " << curlnL << ", oldlnL = " << oldlnL << std::endl; + if (isnan(curlnL) || isnan(oldlnL)) { + std::cerr << "ERROR: likelihood is NaN!" << std::endl; + std::cerr << "Generation = " << gen << std::endl; + std::cerr << "Proposal = " << curProp << std::endl; + std::cerr << "curlnL = " << curlnL << ", oldlnL = " << oldlnL << std::endl; curPath->print(); for (int t = 0; t < sample_time_vec.size(); t++) { - std::cout << sample_time_vec[t]->get() << std::endl; + std::cerr << sample_time_vec[t]->get() << std::endl; } exit(1); } @@ -243,10 +243,10 @@ void mcmc::no_linked_sites(settings& mySettings) { if (curIdx != -1) { double curTimePath = curPath->get_time(curIdx); if (curTime != curTimePath && curIdx != -1) { - std::cout << "ERROR: sample time index for time " << i << " is lost!" << std::endl; - std::cout << "curTime = " << curTime << std::endl; - std::cout << "curIdx = " << curIdx << std::endl; - std::cout << "curTimePath = " << curTimePath << std::endl; + std::cerr << "ERROR: sample time index for time " << i << " is lost!" << std::endl; + std::cerr << "curTime = " << curTime << std::endl; + std::cerr << "curIdx = " << curIdx << std::endl; + std::cerr << "curTimePath = " << curTimePath << std::endl; exit(1); } } @@ -289,10 +289,10 @@ double mcmc::compute_lnL(wfSamplePath* p, measure* m, wienerMeasure* wm) { sample_prob += p->sampleProb(i); } - if (gir != gir) { - std::cout << "Likelihood is nan at generation " << gen << ". Proposal " << curProp << std::endl; - p->print_traj(std::cout); - p->print_time(std::cout); + if (isnan(gir)) { + std::cerr << "ERROR: Likelihood is NaN at generation " << gen << ". Proposal " << curProp << std::endl; + p->print_traj(std::cerr); + p->print_time(std::cerr); exit(1); } return gir + sample_prob; diff --git a/measure.cpp b/measure.cpp index f224835..14a54eb 100644 --- a/measure.cpp +++ b/measure.cpp @@ -187,9 +187,6 @@ double cbpMeasure::log_girsanov_wf_r(path* p, double alpha1, double alpha2, pops double int_mderiv = 0; i = 0; for (j = 0; j < dconts.size()-1; j++) { - if (i >= p->get_length()) { //THIS IS A FUCKING HACK - break; - } //get the potentials while I'm at it. //first the "beginning" potential Hm_w0 += H_wf_r(p->get_traj(i), p->get_time(i), alpha1, alpha2, rho); @@ -209,8 +206,8 @@ double cbpMeasure::log_girsanov_wf_r(path* p, double alpha1, double alpha2, pops * (p->get_time(i)-p->get_time(i-1)); //then the "end" potential double tmp = H_wf_r(p->get_traj(i), p->get_time(i), alpha1, alpha2, rho, 1); - if (tmp != tmp) { - std::cout << "computation of H_wf_r is nan" << std::endl; + if (isnan(tmp)) { + std::cout << "computation of H_wf_r is NaN" << std::endl; std::cout << "p->get_length() = " << p->get_length() << " i = " << i << " j = " << j << std::endl; } Hm_wt += tmp; @@ -220,9 +217,6 @@ double cbpMeasure::log_girsanov_wf_r(path* p, double alpha1, double alpha2, pops double int_msquare = 0; i = 1; for (j = 0; j < dconts.size()-1; j++) { - if (i == p->get_length()) { //THIS IS A FUCKING HACK - break; - } //integrate over the interval using the trapezoid rule while (p->get_time(i) < dconts[j+1]) { int_msquare += (a2_wf_r(p->get_traj(i),p->get_time(i),alpha1,alpha2,rho)+a2_wf_r(p->get_traj(i-1),p->get_time(i-1),alpha1,alpha2,rho))/2.0 @@ -249,9 +243,6 @@ double cbpMeasure::log_girsanov_wf_r(path* p, double alpha1, double alpha2, pops double int_mtime = 0; i = 1; for (j = 0; j < dconts.size()-1; j++) { - if (i == p->get_length()) { //THIS IS A FUCKING HACK - break; - } //integrate over the interval using the trapezoid rule while (p->get_time(i) < dconts[j+1]) { int_mtime += (dHdt_wf_r(p->get_traj(i),p->get_time(i),alpha1,alpha2,rho)+dHdt_wf_r(p->get_traj(i-1),p->get_time(i-1),alpha1,alpha2,rho))/2.0 @@ -265,9 +256,10 @@ double cbpMeasure::log_girsanov_wf_r(path* p, double alpha1, double alpha2, pops if (!is_bridge) { double gir = (Hm_wt-Hm_w0-1.0/2.0*int_mderiv-1.0/2.0*int_msquare-int_mtime); - if (gir != gir) { - std::cout << "ERROR: log_girsanov_wf_r is nan" << std::endl; - std::cout << "Hm_wt = " << Hm_wt << " Hm_w0 = " << Hm_w0 << " int_mderiv = " << int_mderiv << " int_msquare = " << int_msquare << " int_mtime = " << int_mtime << std::endl; + if (isnan(gir)) { + std::cerr << "ERROR: log_girsanov_wf_r is NaN" << std::endl; + std::cerr << "Hm_wt = " << Hm_wt << " Hm_w0 = " << Hm_w0 << " int_mderiv = " << int_mderiv << " int_msquare = " << int_msquare << " int_mtime = " << int_mtime << std::endl; + exit(1); } return gir; } else { @@ -331,9 +323,6 @@ double cbpMeasure::log_girsanov_wfwf_r(path* p, double alpha1, double alpha1p, d double int_mderiv = 0; i = 0; for (j = 0; j < dconts.size()-1; j++) { - if (i >= p->get_length()) { //THIS IS A FUCKING HACK - break; - } //get the potentials while I'm at it. //first the "beginning" potential Hm_w0 += H_wfwf_r(p->get_traj(i), p->get_time(i), alpha1, alpha1p, alpha2, alpha2p, rho); @@ -359,9 +348,6 @@ double cbpMeasure::log_girsanov_wfwf_r(path* p, double alpha1, double alpha1p, d double int_msquare = 0; i = 1; for (j = 0; j < dconts.size()-1; j++) { - if (i == p->get_length()) { //THIS IS A FUCKING HACK - break; - } //integrate over the interval using the trapezoid rule while (p->get_time(i) < dconts[j+1]) { int_msquare += (a2_wfwf_r(p->get_traj(i),p->get_time(i),alpha1,alpha1p,alpha2,alpha2p,rho)+a2_wfwf_r(p->get_traj(i-1),p->get_time(i-1),alpha1,alpha1p,alpha2,alpha2p,rho))/2.0 @@ -377,9 +363,6 @@ double cbpMeasure::log_girsanov_wfwf_r(path* p, double alpha1, double alpha1p, d double int_mtime = 0; i = 1; for (j = 0; j < dconts.size()-1; j++) { - if (i == p->get_length()) { //THIS IS A FUCKING HACK - break; - } //integrate over the interval using the trapezoid rule while (p->get_time(i) < dconts[j+1]) { int_mtime += (dHdt_wfwf_r(p->get_traj(i),p->get_time(i),alpha1,alpha1p,alpha2,alpha2p,rho)+dHdt_wfwf_r(p->get_traj(i-1),p->get_time(i-1),alpha1,alpha1p,alpha2,alpha2p,rho))/2.0 @@ -393,9 +376,10 @@ double cbpMeasure::log_girsanov_wfwf_r(path* p, double alpha1, double alpha1p, d double gir = Hm_wt-Hm_w0-1.0/2.0*int_mderiv-1.0/2.0*int_msquare-int_mtime; - if (gir != gir) { - std::cout << "ERROR: log_girsanov_wfwf_r is nan" << std::endl; - std::cout << "Hm_wt = " << Hm_wt << " Hm_w0 = " << Hm_w0 << " int_mderiv = " << int_mderiv << " int_msquare = " << int_msquare << " int_mtime = " << int_mtime << std::endl; + if (isnan(gir)) { + std::cerr << "ERROR: log_girsanov_wfwf_r is NaN" << std::endl; + std::cerr << "Hm_wt = " << Hm_wt << " Hm_w0 = " << Hm_w0 << " int_mderiv = " << int_mderiv << " int_msquare = " << int_msquare << " int_mtime = " << int_mtime << std::endl; + exit(1); } return gir; @@ -521,12 +505,13 @@ path* cbpMeasure::prop_bridge(double x0, double xt, double t0, double t, std::ve b4_traj.push_back(0); for (i = 0; i < 4; i++) { b4_traj[j] += pow(bb_paths[i]->get_traj(j),2); - if (b4_traj[j] != b4_traj[j]) { - std::cout << "Failing to propose a BES4 bridge from " << x0 << " to " << xt << " during time interval (" << t0 << ", " << t << ")" << std::endl; - std::cout << "This likely means that the time vector is getting loopy, possibly due to pileup of points" << std::endl; - std::cout << "The " << i << "th Brownian bridge between " << u[i] << " and " << xt*v[i] << " is faulty:" << std::endl; - bb_paths[i]->print_traj(std::cout); - bb_paths[i]->print_time(std::cout << std::setprecision(20)); + if (isnan(b4_traj[j])) { + std::cerr << "ERROR: Failing to propose a BES4 bridge from " << x0 << " to " << xt << " during time interval (" << t0 << ", " << t << ")" << std::endl; + std::cerr << "This likely means that the time vector is getting loopy, possibly due to pileup of points" << std::endl; + std::cerr << "The " << i << "th Brownian bridge between " << u[i] << " and " << xt*v[i] << " is faulty:" << std::endl; + bb_paths[i]->print_traj(std::cerr); + bb_paths[i]->print_time(std::cerr << std::setprecision(20)); + exit(1); } } b4_traj[j] = sqrt(b4_traj[j]); @@ -584,7 +569,7 @@ path* wfMeasure::prop_bridge(double x0, double xt, double t0, double t, std::vec gir = cbp->log_girsanov_wf(test_path, 0, 0); accept_prob = rescale + gir; if (accept_prob > 0) { - std::cerr << "Envelope is not sufficient" << std::endl; + std::cerr << "ERROR: Envelope is not sufficient" << std::endl; exit(1); } u = random->uniformRv(); diff --git a/param.cpp b/param.cpp index 2e75b82..01a25a7 100644 --- a/param.cpp +++ b/param.cpp @@ -194,18 +194,18 @@ double sample_time::propose() { //NEW: refelcted uniform double propRatio = 0; if (curVal > youngest || curVal < oldest) { - std::cout << "ERROR: sample_time proposal is outside of range" << std::endl; - std::cout << "oldest = " << oldest << ", youngest = " << youngest << std::endl; - std::cout << "Allele age = " << curParamPath->get_path()->get_time(0) << std::endl; - std::cout << "oldVal = " << oldVal << ", curVal = " << curVal << std::endl; - std::cout << "Starting curVal = " << startVal << std::endl; - std::cout << "Final curVal = " << curVal << std::endl; - std::cout << "old_idx = " << old_idx << ", cur_idx = " << cur_idx << std::endl; - std::cout << "path->time(old_idx) = " << curParamPath->get_path()->get_time(old_idx) << std::endl; - std::cout << "path->time(cur_idx) = " << curParamPath->get_path()->get_time(cur_idx) << std::endl; - std::cout << "propRatio = " << propRatio << std::endl; + std::cerr << "ERROR: sample_time proposal is outside of range" << std::endl; + std::cerr << "oldest = " << oldest << ", youngest = " << youngest << std::endl; + std::cerr << "Allele age = " << curParamPath->get_path()->get_time(0) << std::endl; + std::cerr << "oldVal = " << oldVal << ", curVal = " << curVal << std::endl; + std::cerr << "Starting curVal = " << startVal << std::endl; + std::cerr << "Final curVal = " << curVal << std::endl; + std::cerr << "old_idx = " << old_idx << ", cur_idx = " << cur_idx << std::endl; + std::cerr << "path->time(old_idx) = " << curParamPath->get_path()->get_time(old_idx) << std::endl; + std::cerr << "path->time(cur_idx) = " << curParamPath->get_path()->get_time(cur_idx) << std::endl; + std::cerr << "propRatio = " << propRatio << std::endl; //std::cin.ignore(); - //exit(1); + exit(1); } return propRatio; } @@ -231,11 +231,12 @@ double param_age::propose() { //NEW: reflected uniform //curVal = reflectedUniform(oldVal, tuning, -INFINITY, topTime); //double propRatio = 0; - if (propRatio != propRatio) { - std::cout << "ERROR: Proposal ratio is nan! Debugging information:" << std::endl; - std::cout << "oldVal: " << oldVal << " curVal: " << curVal << " tuning " << tuning << std::endl; - std::cout << "log(P(theta | theta')) = " << log(random->truncatedHalfNormalPdf(topTime, 0, curVal, tuning, oldVal)) << std::endl; - std::cout << "log(P(theta' | theta)) = " << log(random->truncatedHalfNormalPdf(topTime, 0, oldVal, tuning, curVal)) << std::endl; + if (isnan(propRatio)) { + std::cerr << "ERROR: Proposal ratio is NaN! Debugging information:" << std::endl; + std::cerr << "oldVal: " << oldVal << " curVal: " << curVal << " tuning " << tuning << std::endl; + std::cerr << "log(P(theta | theta')) = " << log(random->truncatedHalfNormalPdf(topTime, 0, curVal, tuning, oldVal)) << std::endl; + std::cerr << "log(P(theta' | theta)) = " << log(random->truncatedHalfNormalPdf(topTime, 0, oldVal, tuning, curVal)) << std::endl; + exit(1); } propRatio += curParamPath->proposeAlleleAge(curVal, oldVal); return propRatio; @@ -296,7 +297,7 @@ double param_path::proposeAlleleAge(double newAge, double oldAge) { int end_index; ((wfSamplePath*)curPath)->set_update_begin(); if (newAge < oldAge) { - end_index = 2*minUpdate; + end_index = std::min(2*minUpdate, int(curPath->get_length()) - 1); } else { end_index = 0; while (curPath->get_time(end_index) < newAge) { @@ -306,9 +307,9 @@ double param_path::proposeAlleleAge(double newAge, double oldAge) { end_index += 2*minUpdate; } if (end_index > curPath->get_length()) { - std::cout << "ERROR: trying to update allele age path past the end of the path!" << std::endl; - std::cout << "path length = " << curPath->get_length() << ", end_index = " << end_index << std::endl; - std::cout << "newAge = " << newAge << std::endl; + std::cerr << "ERROR: trying to update allele age path past the end of the path!" << std::endl; + std::cerr << "path length = " << curPath->get_length() << ", end_index = " << end_index << std::endl; + std::cerr << "newAge = " << newAge << std::endl; curPath->print_time(); exit(1); } @@ -363,11 +364,11 @@ std::vector param_path::make_time_vector(double newAge, int end_index, p // for (int j = 0; j < timesToInclude.size()-1; j++) { // if (!(timesToInclude[j+1]>(timesToInclude[j]+std::numeric_limits::epsilon()))) { -// std::cout << "ERROR: Times to include isn't strictly increasing!" << std::endl; +// std::cerr << "ERROR: Times to include isn't strictly increasing!" << std::endl; // for (int l = 0; l < timesToInclude.size(); l++) { -// std::cout << timesToInclude[l] << " "; +// std::cerr << timesToInclude[l] << " "; // } -// std::cout << std::endl; +// std::cerr << std::endl; // exit(1); // } // } @@ -398,25 +399,25 @@ std::vector param_path::make_time_vector(double newAge, int end_index, p } } - //check that time vector is strictly increasing +// //check that time vector is strictly increasing // for (int j = 0; j < newTimes.size()-1; j++) { // if (!(newTimes[j+1]>newTimes[j])) { -// std::cout << "ERROR: new time vector of length " << newTimes.size() << " not strictly increasing" << std::endl; -// std::cout << "Machine eps is " << std::numeric_limits::epsilon() << std::endl; -// std::cout << "Times to include are" << std::endl; +// std::cerr << "ERROR: new time vector of length " << newTimes.size() << " not strictly increasing" << std::endl; +// std::cerr << "Machine eps is " << std::numeric_limits::epsilon() << std::endl; +// std::cerr << "Times to include are" << std::endl; // for (int l = 0; l < timesToInclude.size(); l++) { -// std::cout << timesToInclude[l] << " "; +// std::cerr << timesToInclude[l] << " "; // } -// std::cout << std::endl; +// std::cerr << std::endl; // for (int l = 0; l < timesToInclude.size() - 1; l++) { -// std::cout << timesToInclude[l+1] << " - " << timesToInclude[l] << " = " << timesToInclude[l+1] - timesToInclude[l] << " "; +// std::cerr << timesToInclude[l+1] << " - " << timesToInclude[l] << " = " << timesToInclude[l+1] - timesToInclude[l] << " "; // } -// std::cout << std::endl; +// std::cerr << std::endl; // for (int l = 0; l < timesToInclude.size() - 1; l++) { -// std::cout << (timesToInclude[l+1]>timesToInclude[l]) << " "; +// std::cerr << (timesToInclude[l+1]>timesToInclude[l]) << " "; // } -// std::cout << std::endl; -// std::cout << "newTimes[" << j << "+1] = " << newTimes[j+1] << ", newTimes[" << j << "] = " << newTimes[j] << std::endl; +// std::cerr << std::endl; +// std::cerr << "newTimes[" << j << "+1] = " << newTimes[j+1] << ", newTimes[" << j << "] = " << newTimes[j] << std::endl; // exit(1); // } // } @@ -500,8 +501,8 @@ double param_path::proposeAgePath(double x0,double xt,double t0,double t, std::v double new_like = myCBP.log_girsanov_wf_r(newPath, a1->get(), a2->get(), rho,0); - if (new_like != new_like) { - std::cerr << "ERROR: new path likelihood is nan! Debugging information sent to stderr:" << std::endl; + if (isnan(new_like)) { + std::cerr << "ERROR: new path likelihood is NaN! Debugging information sent to stderr:" << std::endl; std::cerr << "New path:" << std::endl; std::cerr << "alpha1 = " << a1->get() << " alpha2 = " << a2->get() << std::endl; newPath->print_traj(std::cerr); @@ -512,8 +513,8 @@ double param_path::proposeAgePath(double x0,double xt,double t0,double t, std::v double old_like = myCBP.log_girsanov_wf_r(oldPath, a1->get(), a2->get(), rho,0); - if (old_like != old_like) { - std::cerr << "ERROR: old path likelihood is nan! Debugging information sent to stderr:" << std::endl; + if (isnan(old_like)) { + std::cerr << "ERROR: old path likelihood is NaN! Debugging information sent to stderr:" << std::endl; std::cerr << "alpha1 = " << a1->get() << " alpha2 = " << a2->get() << std::endl; std::cerr << "Old path:" << std::endl; oldPath->print_traj(std::cerr); @@ -526,8 +527,8 @@ double param_path::proposeAgePath(double x0,double xt,double t0,double t, std::v propRatio += -1.0/2.0*xt*xt*(1.0/tNew-1.0/tOld)+2*log(tOld)-2*log(tNew); - if (propRatio != propRatio) { - std::cerr << "ERROR: proposal ratio is nan! Debugging information sent to stderr:" << std::endl; + if (isnan(propRatio)) { + std::cerr << "ERROR: proposal ratio is NaN! Debugging information sent to stderr:" << std::endl; std::cerr << "New path:" << std::endl; newPath->print_traj(std::cerr); newPath->print_time(std::cerr); diff --git a/path.cpp b/path.cpp index 47384e4..afc9311 100644 --- a/path.cpp +++ b/path.cpp @@ -143,15 +143,16 @@ void path::modify(path* p, int i) { //old_time[j] = time[i+j]; time[i+j] = p->get_time(j); if (time[i+j] < time[i+j-1]) { - std::cout << "ERROR: time vector is not sorted!" << std::endl; - std::cout << time[i+j] << " >= " << time[i+j-1] << std::endl; + std::cerr << "ERROR: time vector is not sorted!" << std::endl; + std::cerr << time[i+j] << " >= " << time[i+j-1] << std::endl; exit(1); } } if (trajectory.size() != time.size()) { - std::cout << "ERROR: Path trajectory and time are not same lenght!" << std::endl; - std::cout << "trajectory.size() = " << trajectory.size() << std::endl; - std::cout << "time.size() = " << time.size() << std::endl; + std::cerr << "ERROR: Path trajectory and time are not same lenght!" << std::endl; + std::cerr << "trajectory.size() = " << trajectory.size() << std::endl; + std::cerr << "time.size() = " << time.size() << std::endl; + exit(1); } } else { old_trajectory.resize(0); @@ -168,9 +169,10 @@ void path::reset() { } } if (trajectory.size() != time.size()) { - std::cout << "ERROR: Path trajectory and time are not same lenght!" << std::endl; - std::cout << "trajectory.size() = " << trajectory.size() << std::endl; - std::cout << "time.size() = " << time.size() << std::endl; + std::cerr << "ERROR: Path trajectory and time are not same length!" << std::endl; + std::cerr << "trajectory.size() = " << trajectory.size() << std::endl; + std::cerr << "time.size() = " << time.size() << std::endl; + exit(1); } } @@ -337,11 +339,12 @@ wfSamplePath::wfSamplePath(std::vector& st, popsize* p, wfMeasure* breakPoints.push_back(sample_time_vec[i]->get()); breakPoints.push_back(sample_time_vec[i]->get_oldest()); } - - double max_time = breakPoints[*std::max_element(breakPoints.begin(), breakPoints.end())]; - - std::vector curBreaks = myPop->getBreakTimes(first_nonzero, max_time); - + + double min_time = *std::min_element(breakPoints.begin(), breakPoints.end()); + double max_time = *std::max_element(breakPoints.begin(), breakPoints.end()); + + std::vector curBreaks = myPop->getBreakTimes(std::min(min_time, first_nonzero), max_time); + for (int i = 0; i < curBreaks.size(); i++) { breakPoints.push_back(curBreaks[i]); } @@ -448,7 +451,8 @@ void wfSamplePath::set_allele_age(double a, path* p, int i) { //if it's part of the new trajectory, find where it is search_it = std::lower_bound(time.begin(),time.end(),sample_time_vec[j]->get()); if (search_it == time.end() && sample_time_vec[j]->get() != time[time.size()-1]) { - std::cout << "ERROR: could not find sample time index " << j << " with value " << sample_time_vec[j]->get() << " in time vector!" << std::endl; + std::cerr << "ERROR: could not find sample time index " << j << " with value " << sample_time_vec[j]->get() << " in time vector!" << std::endl; + exit(1); } new_idx = search_it-time.begin(); sample_time_vec[j]->set_idx(new_idx); @@ -466,11 +470,11 @@ void wfSamplePath::set_allele_age(double a, path* p, int i) { void wfSamplePath::resetIntermediate() { //check some things if (update_begin) { - std::cout << std::endl << "ERROR: Trying to reset an intermediate part of the path, but allele age was updated!" << std::endl; + std::cerr << "ERROR: Trying to reset an intermediate part of the path, but allele age was updated!" << std::endl; exit(1); } if (old_index == -1) { - std::cout << std::endl << "ERROR: Trying to reset an intermdiate part of the path, but old_index = -1!" << std::endl; + std::cerr << "ERROR: Trying to reset an intermdiate part of the path, but old_index = -1!" << std::endl; exit(1); } //replace the trajectory @@ -484,7 +488,7 @@ void wfSamplePath::resetIntermediate() { void wfSamplePath::resetBeginning() { //check some things if (!update_begin) { - std::cout << std::endl << "ERROR: Trying to reset the beginning of the path, but allele age wasn't updated!" << std::endl; + std::cerr << "ERROR: Trying to reset the beginning of the path, but allele age wasn't updated!" << std::endl; exit(1); } @@ -555,7 +559,7 @@ void wfSamplePath::updateFirstNonzero() { void path::replace_time(std::vector new_time) { if (new_time.size() != time.size()) { - std::cout << "ERROR: Trying to replace a time vector with one of a different size!" << std::endl; + std::cerr << "ERROR: Trying to replace a time vector with one of a different size!" << std::endl; exit(1); } time = new_time; diff --git a/path.h b/path.h index f08434f..f0aa000 100644 --- a/path.h +++ b/path.h @@ -37,12 +37,12 @@ class path { //element access std::vector get_traj() {return trajectory;}; - double get_traj(int i) {return trajectory[i];}; + double get_traj(int i) {return trajectory.at(i);}; std::vector get_traj(int i, int j); - void set_traj(double x, int i) {trajectory[i] = x;}; + void set_traj(double x, int i) {trajectory.at(i) = x;}; std::vector get_time() {return time;}; std::vector& get_time_ref() {return time;}; - double get_time(int i) {return time[i];}; + double get_time(int i) {return time.at(i);}; std::vector get_time(int i, int j); double get_length() {return trajectory.size();}; double get_length_time() {return time.size();}; diff --git a/popsize.cpp b/popsize.cpp index f566751..9b21926 100644 --- a/popsize.cpp +++ b/popsize.cpp @@ -49,17 +49,17 @@ popsize::popsize(settings& s) { for (int i = 0; i < sizes.size()-1; i++) { //the times better be in decreasing order! if (times[i] <= times[i+1]) { - std::cout << "ERROR: Times are not in decreasing order!" << std::endl; - std::cout << "Time " << i << " <= Time " << i+1 << std::endl; - std::cout << times[i] << " <= " << times[i+1] << std::endl; + std::cerr << "ERROR: Times are not in decreasing order!" << std::endl; + std::cerr << "Time " << i << " <= Time " << i+1 << std::endl; + std::cerr << times[i] << " <= " << times[i+1] << std::endl; exit(1); } } //the last thing better be infinity, and not change if (times[times.size()-1] !=-INFINITY || rates[rates.size()-1] != 0) { - std::cout << "ERROR: Final time point is not infinity OR final epoch not constant" << std::endl; - std::cout << "Size Rate Time" << std::endl; - std::cout << sizes[sizes.size()-1] << " " << rates[rates.size()-1] << " " << times[times.size()-1] << std::endl; + std::cerr << "ERROR: Final time point is not infinity OR final epoch not constant" << std::endl; + std::cerr << "Size Rate Time" << std::endl; + std::cerr << sizes[sizes.size()-1] << " " << rates[rates.size()-1] << " " << times[times.size()-1] << std::endl; exit(1); } computeT(); diff --git a/popsize.h b/popsize.h index 40d50ad..17f5e44 100644 --- a/popsize.h +++ b/popsize.h @@ -36,9 +36,9 @@ class popsize { std::vector getBreakTimes(double t0, double t); //get an entry of the vectors - double getSizes(int i) {return sizes[i];}; - double getRates(int i) {return rates[i];}; - double getTimes(int i) {return times[i];}; + double getSizes(int i) {return sizes.at(i);}; + double getRates(int i) {return rates.at(i);}; + double getTimes(int i) {return times.at(i);}; private: diff --git a/settings.cpp b/settings.cpp index a5b4b73..5b3f619 100644 --- a/settings.cpp +++ b/settings.cpp @@ -43,12 +43,12 @@ settings::settings(int argc, char* const argv[]) { inputFile = ""; mySeed = time(0); output_tsv = 0; - a1prop = 2.0; - a2prop = 2.0; + a1prop = 5.0; + a2prop = 5.0; ageprop = 20.0; - endprop = 5.0; + endprop = 2.0; timeprop = 0.1; - pathprop = 10.0; + pathprop = 20.0; a1start = 25.0; a2start = 50.0; set_gen = 0; @@ -181,8 +181,8 @@ std::vector settings::parse_bridge_pars() { pars.push_back(atof(cur_par.c_str())); } if (pars.size() < 4) { - std::cout << "ERROR: Not enough bridge parameters" << std::endl; - std::cout << "Only " << pars.size() << " specified; 4 are required" << std::endl; + std::cerr << "ERROR: Not enough bridge parameters" << std::endl; + std::cerr << "Only " << pars.size() << " specified; 4 are required" << std::endl; exit(1); } return pars; @@ -200,7 +200,7 @@ void settings::print() { std::cout << "gamma\t" << pars[2] << std::endl; std::cout << "t\t" << pars[3] << std::endl; } else if (mcmc) { - std::cout << "num_gen\t" << num_gen << std::endl; + std::cerr << "num_gen\t" << num_gen << std::endl; if (linked_sites) { //file destinations } else { @@ -211,8 +211,8 @@ void settings::print() { popsize* settings::parse_popsize_file() { if (popFile == "") { - std::cout << "ERROR: No population size history specified! Use -P option" << std::endl; - exit(1); + std::cerr << "ERROR: No population size history specified! Use -P option" << std::endl; + exit(1); } @@ -248,11 +248,11 @@ std::vector settings::parse_input_file(MbRandom* r) { std::istringstream curLine(curLineString); curLine >> curCount >> curSS >> curLowTime >> curHighTime; if (curCount < 0 || curCount > curSS) { - std::cout << "Allele count is not between 0 and sample size: X = " << curCount << ", SS = " << curSS << std::endl; + std::cerr << "ERROR: Allele count is not between 0 and sample size: X = " << curCount << ", SS = " << curSS << std::endl; exit(1); } if (curLowTime > curHighTime) { - std::cout << "Low end of time range higher than high end: t_low = " << curLowTime << ", t_high = " << curHighTime << std::endl; + std::cerr << "ERROR: Low end of time range higher than high end: t_low = " << curLowTime << ", t_high = " << curHighTime << std::endl; exit(1); } //Convert time units @@ -272,7 +272,7 @@ std::vector settings::parse_input_file(MbRandom* r) { //check that most recent time point is fixed int num_sam = sample_time_vec.size(); if (sample_time_vec[num_sam-1]->get_oldest()get_youngest()) { - std::cout << "ERROR: most recent time point must have no uncertainty" << std::endl; + std::cerr << "ERROR: most recent time point must have no uncertainty" << std::endl; exit(1); }