Skip to content

Updates to Ascertainment branch #6

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
ekirving wants to merge 7 commits into
base: ascertainment
Choose a base branch
from
Open

Updates to Ascertainment branch #6

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
ea57c87
further adjustments to the proposal chances, based on ESS of resultin…
ekirving Jul 9, 2019
7ff5814
added cmake config file for debugging with clion
ekirving Jul 12, 2019
dbf6211
replaced hacky checks for iterator overrun with calls to the vector a…
ekirving Jul 11, 2019
60ac2ab
send all error messages to std::cerr and then exit(1)
ekirving Jul 12, 2019
76f236f
improved handling of NaN checks
ekirving Jul 14, 2019
65293b4
bug fix for the missing epoch boundaries issue
ekirving Jul 12, 2019
9788e47
bug fix for index overflow when making an age update near the boundar…
ekirving Jul 15, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions .gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
cmake-build-debug/
data/
27 changes: 27 additions & 0 deletions CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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)
30 changes: 15 additions & 15 deletions mcmc.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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);
}
Expand Down Expand Up @@ -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);
}
}
Expand Down Expand Up @@ -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;
Expand Down
51 changes: 18 additions & 33 deletions measure.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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;
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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 {
Expand Down Expand Up @@ -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);
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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;
Expand Down Expand Up @@ -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]);
Expand Down Expand Up @@ -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();
Expand Down
81 changes: 41 additions & 40 deletions param.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
Expand All @@ -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;
Expand Down Expand Up @@ -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) {
Expand All @@ -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);
}
Expand Down Expand Up @@ -363,11 +364,11 @@ std::vector<double> 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<double>::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);
// }
// }
Expand Down Expand Up @@ -398,25 +399,25 @@ std::vector<double> 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<double>::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<double>::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);
// }
// }
Expand Down Expand Up @@ -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);
Expand All @@ -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);
Expand All @@ -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);
Expand Down
Loading