From d48c5663e00817db0488f33d0c9579713ee82cd6 Mon Sep 17 00:00:00 2001 From: Tracy Heath Date: Wed, 16 Apr 2014 17:27:58 -0700 Subject: [PATCH] Major changes and bug fixes --- src/MbMatrix.h | 1 + src/MbRandom.h | 2 +- src/Mcmc.cpp | 51 ++-- src/Model.cpp | 74 ++++-- src/Model.h | 9 +- src/Parameter_speciaton.cpp | 238 +++++++++++++----- src/Parameter_speciaton.h | 14 +- src/Parameter_tree.cpp | 490 ++++++++++++++++++++++++------------ src/Parameter_tree.h | 27 +- src/Parameter_treescale.cpp | 139 +++++++--- src/Parameter_treescale.h | 6 + src/dppdiv.cpp | 30 ++- 12 files changed, 774 insertions(+), 307 deletions(-) diff --git a/src/MbMatrix.h b/src/MbMatrix.h index 9b8dfc2..5bb200d 100644 --- a/src/MbMatrix.h +++ b/src/MbMatrix.h @@ -40,6 +40,7 @@ #include #include #include +#include /*! diff --git a/src/MbRandom.h b/src/MbRandom.h index 3815364..e7eac47 100644 --- a/src/MbRandom.h +++ b/src/MbRandom.h @@ -386,7 +386,7 @@ inline double MbRandom::normalPdf(double mu, double sigma, double x) { */ inline double MbRandom::lnNormalPdf(double mu, double sigma, double x) { - return -0.5 * std::log(2.0 * PI * sigma) - 0.5 * (x - mu) * (x - mu) / (sigma * sigma); + return -0.5 * std::log(2.0 * PI * sigma * sigma) - 0.5 * (x - mu) * (x - mu) / (sigma * sigma); } diff --git a/src/Mcmc.cpp b/src/Mcmc.cpp index 27b8b90..a5c1ed2 100644 --- a/src/Mcmc.cpp +++ b/src/Mcmc.cpp @@ -67,19 +67,16 @@ Mcmc::Mcmc(MbRandom *rp, Model *mp, int nc, int pf, int sf, string ofp, bool wdf void Mcmc::runChain(void) { string pFile = fileNamePref + ".p"; // parameter file name - //string treFile = fileNamePref + ".t"; // tree file name string figTFile = fileNamePref + ".ant.tre"; // write to a file with the nodes colored by their rate classes string dFile = fileNamePref + ".info.out"; string ndFile = fileNamePref + ".nodes.out"; // info about nodes string mtxFile = fileNamePref + ".rates.out"; // info about nodes ofstream pOut(pFile.c_str(), ios::out); - //ofstream tOut(treFile.c_str(), ios::out); ofstream fTOut(figTFile.c_str(), ios::out); ofstream nOut(ndFile.c_str(), ios::out); ofstream mxOut; ofstream dOut; - //writeCalibrationTree(); if(writeInfoFile) dOut.open(dFile.c_str(), ios::out); if(printratef) @@ -134,7 +131,11 @@ void Mcmc::runChain(void) { } else{ modelPtr->updateRejected(); - Tree *t = modelPtr->getActiveTree(); + + // TAH : without this, the lnls were not right for some moves following rejected ones + Tree *t = modelPtr->getActiveTree(); + Treescale *ts = modelPtr->getActiveTreeScale(); + t->setTreeScale(ts->getScaleValue()); t->flipAllCls(); t->flipAllTis(); t->upDateAllCls(); @@ -148,13 +149,37 @@ void Mcmc::runChain(void) { t->setNodeRateValues(); } - if ( n % sampleFrequency == 0 || n == 1) + // sample chain + if ( n % sampleFrequency == 0 || n == 1){ sampleChain(n, pOut, fTOut, nOut, oldLnLikelihood); - + //sampleRtsFChain(n, mxOut); + } + + //Logger & logger = Logger::getInstance(); + //std::ostream &o = logger.debugStream(); + //parm->print(o); + //parm->print(std::cerr); +# if 0 + if( n % 20 == 0 || n == 1) + modelPtr->getActiveTree()->verifyTreeDebug(n, parm->getName()); +# endif + + if(testLnL){ + Tree *t = modelPtr->getActiveTree(); + t->flipAllCls(); + t->flipAllTis(); + t->upDateAllCls(); + t->upDateAllTis(); + modelPtr->upDateRateMatrix(); + modelPtr->setTiProb(); + double chklnl = modelPtr->lnLikelihood(); + cout << " " << n << " " << parm->getName() << " --> OL" << prevlnl << " --- NL" << newLnLikelihood << " (" << chklnl << ", " << oldLnLikelihood << ")\n\n"; + } } - cout << " Markov chain completed." << endl; + Tree *t = modelPtr->getActiveTree(); + int timeEnd = time(NULL); + cout << " Markov chain completed in " << (static_cast(timeEnd - timeSt)) << " seconds" << endl; pOut.close(); - //tOut.close(); fTOut.close(); dOut.close(); nOut.close(); @@ -192,7 +217,6 @@ void Mcmc::sampleChain(int gen, ofstream ¶Out, ofstream &figTOut, if(gen == 1){ paraOut << "Gen\tlnLikelihood\tf(A)\tf(C)\tf(G)\tf(T)"; paraOut << "\tr(AC)\tr(AG)\tr(AT)\tr(CG)\tr(CT)\tr(GT)\tshape\tave rate\tnum rate groups\tconc param\n"; - //treeOut << "#NEXUS\nbegin trees;\n"; figTOut << "#NEXUS\nbegin trees;\n"; nodeOut << "Gen\tlnL"; nodeOut << "\tNetDiv(b-d)\tRelativeDeath(d/b)"; @@ -219,10 +243,12 @@ void Mcmc::sampleChain(int gen, ofstream ¶Out, ofstream &figTOut, } if(treePr == 7){ nodeOut << t->getCalBDSSNodeInfoIndicatorNames(); + nodeOut << "\tnum.tips"; } nodeOut << "\n"; } + // then print stuff paraOut << gen << "\t" << lnl; for(int i=0; igetNumStates(); i++) paraOut << "\t" << f->getFreq(i); @@ -233,10 +259,7 @@ void Mcmc::sampleChain(int gen, ofstream ¶Out, ofstream &figTOut, paraOut << "\t" << nr->getNumRateGroups(); paraOut << "\t" << nr->getConcenParam(); paraOut << "\n"; - - //treeOut << " tree t" << gen << " = "; - //treeOut << t->getTreeDescription() << "\n"; - + figTOut << " tree t" << gen << " = "; figTOut << t->getFigTreeDescription() << "\n"; @@ -278,11 +301,11 @@ void Mcmc::sampleChain(int gen, ofstream ¶Out, ofstream &figTOut, } if(treePr == 7){ nodeOut << t->getCalBDSSNodeInfoIndicatorList(); + nodeOut << "\t" << t->getSumIndicatorV(); } nodeOut << "\n"; if(gen == numCycles){ - //treeOut << "end;\n"; figTOut << "end;\n"; figTOut << "\nbegin figtree;\n"; figTOut << " set appearance.branchColorAttribute=\"rate_cat\";\n"; diff --git a/src/Model.cpp b/src/Model.cpp index 608d4b2..d198606 100644 --- a/src/Model.cpp +++ b/src/Model.cpp @@ -47,17 +47,17 @@ #include #include #include -#include - +#include using namespace std; Model::Model(MbRandom *rp, Alignment *ap, string ts, double pm, double ra, double rb, double hal, double hbe, bool ubl, bool alnm, int offmv, bool rndNo, - string clfn, int nodpr, double bdr, double bda, double fxclkrt, bool roofix, + string clfn, int nodpr, double bdr, double bda, double bds, double fxclkrt, bool roofix, bool sfb, bool ehpc, bool dphpc, int dphpng, bool gamhp, int rmod, bool fxmod, - bool ihp, string tipdfn) { + bool ihp, string tipdfn, bool fxtr) { + // remember pointers to important objects... ranPtr = rp; alignmentPtr = ap; priorMeanN = pm; @@ -75,6 +75,8 @@ Model::Model(MbRandom *rp, Alignment *ap, string ts, double pm, double ra, doubl turnedOffMove = offmv; fixedClockRate = fxclkrt; fixSomeModParams = fxmod; + fixTestRun = fxtr; + estAbsRts = false; double initRootH = 1.0; runIndCalHP = ihp; rHtY = 0.0; @@ -101,6 +103,7 @@ Model::Model(MbRandom *rp, Alignment *ap, string ts, double pm, double ra, doubl cout << "\nStarting with seeds: { " << startS1 << " , " << startS2 << " } \n\n"; + // ...and initialize some important variables numGammaCats = 4; numPatterns = alignmentPtr->getNumChar(); @@ -129,7 +132,7 @@ Model::Model(MbRandom *rp, Alignment *ap, string ts, double pm, double ra, doubl parms[i].push_back( nr ); // restaurant containing node rates parms[i].push_back( conp ); // hyper prior on DPP concentration parameter parms[i].push_back( new Treescale(ranPtr, this, initRootH, rHtY, rHtO, tsPrDist, rtCalib, ehpc) ); // the tree scale prior - parms[i].push_back( new Speciation(ranPtr, this, bdr, bda, initRootH) ); // hyper prior on diversification for cBDP speciation + parms[i].push_back( new Speciation(ranPtr, this, bdr, bda, bds, initRootH) ); // hyper prior on diversification for cBDP speciation parms[i].push_back( excal ); // hyper prior exponential node calibration parameters } numParms = (int)parms[0].size(); @@ -141,16 +144,25 @@ Model::Model(MbRandom *rp, Alignment *ap, string ts, double pm, double ra, doubl parms[0][i]->print(std::cout); + // initialize the probabilities for updating different parameters setUpdateProbabilities(true); if(ehpc) excal->getAllExpHPCalibratedNodes(); + + // allocate and initialize conditional likelihoods initializeConditionalLikelihoods(); + + // allocate the transition probability matrices initializeTransitionProbabilityMatrices(); + + // instantiate the transition probability calculator tiCalculator = new MbTransitionMatrix( getActiveExchangeability()->getRate(), getActiveBasefreq()->getFreq(), true ); setTiProb(); myCurLnL = lnLikelihood(); cout << "lnL = " << myCurLnL << endl; + + } Model::~Model(void) { @@ -267,6 +279,7 @@ ExpCalib* Model::getActiveExpCalib(void) { void Model::initializeConditionalLikelihoods(void) { + // allocate conditional likelihoods int nNodes = 2*alignmentPtr->getNumTaxa()-1; int nChar = alignmentPtr->getNumChar(); int sizeOneNode = nChar * numGammaCats * 4; @@ -281,6 +294,7 @@ void Model::initializeConditionalLikelihoods(void) { clPtr[i][j] = &cls[ i * sizeOneSpace + j * sizeOneNode ]; } + // initialize the tip conditional likelihoods for (int i=0; igetNumTaxa(); i++) { double *cl0 = clPtr[0][i]; @@ -386,6 +400,7 @@ void Model::setTiProb(void) { p->setIsTiDirty(false); } } + // TAH root rate debug. This stuff below is stupid anyway # if ASSIGN_ROOT Node *roo = t->getRoot(); t->setRootRateValue(r->getRateForNodeIndexed(roo->getIdx())); @@ -396,7 +411,10 @@ void Model::setTiProb(Node *p, Shape *s, NodeRate *r) { int activeTi = p->getActiveTi(); int idx = p->getIdx(); - double branchProportion = p->getAnc()->getNodeDepth() - p->getNodeDepth(); + Treescale *ts = getActiveTreeScale(); + double sv = 1.0; //ts->getScaleValue(); // FIXPARM + if(estAbsRts) sv = ts->getScaleValue(); + double branchProportion = (p->getAnc()->getNodeDepth() - p->getNodeDepth()) * sv; double rP = r->getRateForNodeIndexed(p->getIdx()); # if 0 double rA = r->getRateForNodeIndexed(p->getAnc()->getIdx()); @@ -406,6 +424,8 @@ void Model::setTiProb(Node *p, Shape *s, NodeRate *r) { # endif if(rP == 0.0){ + // then this means that this is calculating the lnl for nodes that are not assigned to + // RateGroups in currentRateGroups cerr << "ERROR: Problem rP = 0" << endl; exit(1); } @@ -428,6 +448,9 @@ void Model::setTiProb(Node *p, Shape *s, NodeRate *r) { rt = s->getRate(3); tis[activeTi][idx][3] = tiCalculator->tiProbs( v*rt, tis[activeTi][idx][3] ); # endif + // set node info for printing + //p->setBranchTime(branchProportion); + //p->setRtGrpVal(rP); } void Model::setNodeRateGrpIndxs(void) { @@ -437,8 +460,9 @@ void Model::setNodeRateGrpIndxs(void) { # if ASSIGN_ROOT int rtID = t->getNumNodes() + 2; # else - int rtID = t->getRoot()->getIdx(); + int rtID = t->getRoot()->getIdx(); // TAH root rate debug # endif + //int numRCats = r->getNumRateGroups(); for(int n=0; ngetNumNodes(); n++){ if(n != rtID){ Node *p = t->getNodeByIndex(n); @@ -524,6 +548,12 @@ double Model::readCalibFile(void) { T10 T9 0.9 0.9 */ + /* + TAH + This also initializes the root height parameter + I'm not sure how best to do this + + */ cout << "\nCalibrations:" << endl; bool rootIs = false; @@ -559,6 +589,7 @@ double Model::readCalibFile(void) { fixRootHeight = false; } } + //else if(rooCal->getPriorDistributionType() == 2){ else if(rooCal->getPriorDistributionType() > 1){ fixRootHeight = false; yb = rooCal->getYngTime(); @@ -573,6 +604,7 @@ double Model::readCalibFile(void) { double tmpv; if((*v)->getPriorDistributionType() == 1) tmpv = (*v)->getOldTime(); + //else if((*v)->getPriorDistributionType() == 2) else if((*v)->getPriorDistributionType() > 1) tmpv = (*v)->getYngTime() * 1.1; if(tmpv > yb) @@ -652,10 +684,10 @@ void Model::setUpdateProbabilities(bool initial) { bfp = 0.3; srp = 0.3; shp = 0.3; - ntp = 0.4; + ntp = 0.6; dpp = 0.5; cpa = 0.3; - tsp = 0.3; + tsp = 0.5; spp = 0.4; ehp = 0.0; fcp = 0.0; @@ -686,23 +718,21 @@ void Model::setUpdateProbabilities(bool initial) { else if(turnedOffMove == 5){ dpp = 0.0; cpa = 0.0; - setNodeRateGrpIndxs(); - *parms[0][3] = *parms[1][3]; + setNodeRateGrpIndxs(); // set the rates on the tree + *parms[0][3] = *parms[1][3]; //make sure both trees are the same } - else if(turnedOffMove == 6 || cpfix == true) + else if(turnedOffMove == 6 || cpfix == true) // if this move is turned off then the cp is set using the prior mean number of groups cpa = 0.0; else if(turnedOffMove == 8) spp = 0.0; if(fixRootHeight) tsp = 0.0; - if(treeTimePrior == 1) + if(treeTimePrior == 1) // treeTimePrior = 4 is now BDSS || treeTimePrior == 4) spp = 0.0; if(zeroNodeTimeMove == 1){ ntp = 0.0; cout << "All internal node times are fixed" << endl; } - if(fixedClockRate > 0.0) - cpa = 0.0; if(exponCalibHyperParm){ if(initial) ehp = 0.3; @@ -714,13 +744,25 @@ void Model::setUpdateProbabilities(bool initial) { shp = 0.0; } - if(treeTimePrior > 3) + if(treeTimePrior > 3) // might want to change this spp = 0.5; if(treeTimePrior > 5){ ntp = 0.8; fcp = 0.4; } + if(fixedClockRate > 0.0){ // FIXPARM + cpa = 0.0; + dpp = 0.0; + //spp = 0.0; + } + + // TAH: FIXING 6 Feb + if(fixTestRun){ + spp = 0.0; + //tsp = 0.0; + //ntp = 0.0; + } updateProb.clear(); updateProb.push_back(bfp); // 1 basefreq diff --git a/src/Model.h b/src/Model.h index 2bcfe0f..7981f8e 100644 --- a/src/Model.h +++ b/src/Model.h @@ -63,9 +63,9 @@ class Model { Model(MbRandom *rp, Alignment *ap, std::string ts, double pm, double ra, double rb, double hal, double hbe, bool ubl, bool alnm, int offmv, bool rndNo, std::string clfn, int nodpr, - double bdr, double bda, double fxclkrt, bool roofix, + double bdr, double bda, double bds, double fxclkrt, bool roofix, bool sfb, bool ehpc, bool dphpc, int dphpng, bool gamhp, int rmod, - bool fxmod, bool ihp, std::string tipdfn); + bool fxmod, bool ihp, std::string tipdfn, bool fxtr); ~Model(void); double lnLikelihood(void); double getPriorMeanV(void) { return priorMeanN; } @@ -100,6 +100,9 @@ class Model { bool getExponCalibHyperParm() { return exponCalibHyperParm; } bool getExponDPMCalibHyperParm() { return exponDPMCalibHyperParm; } void setUpdateProbabilities(bool initial); + void setEstAbsRates(bool b) { estAbsRts = b; } + void setFixTestRun(bool b) { fixTestRun = b; } + bool getFixTestRun(void) { return fixTestRun; } private: void initializeConditionalLikelihoods(void); @@ -151,6 +154,8 @@ class Model { bool fixSomeModParams; bool cpfix; bool runIndCalHP; + bool estAbsRts; + bool fixTestRun; }; #endif diff --git a/src/Parameter_speciaton.cpp b/src/Parameter_speciaton.cpp index 4369593..9c5e975 100644 --- a/src/Parameter_speciaton.cpp +++ b/src/Parameter_speciaton.cpp @@ -40,18 +40,25 @@ using namespace std; -Speciation::Speciation(MbRandom *rp, Model *mp, double bdr, double bda, double initRH) : Parameter(rp, mp) { +Speciation::Speciation(MbRandom *rp, Model *mp, double bdr, double bda, double bds, double initRH) : Parameter(rp, mp) { + // Based on BEAST implementation of the birth-death model described in Gernhard (2008) + // Note: the beast implementation doesn't allow for lambda=mu maxdivV = 30000.0; name = "SP"; - relativeDeath = bda; - netDiversificaton = bdr; + relativeDeath = ranPtr->uniformRv(); + netDiversificaton = ranPtr->uniformRv(); probSpeciationS = ranPtr->uniformRv(); - extantSampleRate = ranPtr->betaRv(1.0,1.0); - setAllBDFossParams(); + extantSampleRate = 1.0; treeTimePrior = modelPtr->getTreeTimePriorNum(); + if(mp->getFixTestRun()){ + relativeDeath = bda; + netDiversificaton = bdr; + probSpeciationS = bds; + } + setAllBDFossParams(); if(relativeDeath < 0.0 && netDiversificaton < 0.0){ relativeDeath = ranPtr->uniformRv(); @@ -122,7 +129,8 @@ double Speciation::update(double &oldLnL) { Tree *t = modelPtr->getActiveTree(); double oldtreeprob = getLnTreeProb(t); - double lnProposalRatio = 0.0; + double lnProposalRatio = 0.0; + if(treeTimePrior == 2){ lnProposalRatio += updateNetDivRate(); @@ -135,22 +143,15 @@ double Speciation::update(double &oldLnL) { lnProposalRatio += updateNetDivRate(); } else if(treeTimePrior > 3){ - double *probMove = new double[4]; - probMove[0] = 0.2; - probMove[1] = 0.3; - probMove[2] = 0.3; - probMove[3] = 0.2; - int mvType = ranPtr->categoricalRv(probMove,4); - if(mvType == 0) - lnProposalRatio += updateRelDeathRt(); - else if(mvType == 1) - lnProposalRatio += updateNetDivRate(); - else if(mvType == 2) - lnProposalRatio += updateBDSSFossilProbS(); - else if(mvType == 3) - lnProposalRatio += updateBDSSSampleProbRho(); - delete [] probMove; - setAllBDFossParams(); + + updateRelDeathRt(t); + updateNetDivRate(t); + updateBDSSFossilProbS(t); + modelPtr->setLnLGood(true); + modelPtr->setMyCurrLnl(oldLnL); + return 0.0; + + } double newtreeprob = getLnTreeProb(t); double lnPriorRatio = (newtreeprob - oldtreeprob); @@ -165,26 +166,67 @@ double Speciation::updateRelDeathRt(void) { double rdwindow = 0.2; double oldRD = relativeDeath; - double u; - double newRD; - u = ranPtr->uniformRv(-0.5,0.5) * (rdwindow); - newRD = oldRD + u; + double newRD = getNewValSWindoMv(oldRD, 0.0, 0.99999, rdwindow); + relativeDeath = newRD; + return 0.0; + +} + +double Speciation::updateRelDeathRt(Tree *t) { + + double oldtreeprob = getLnTreeProb(t); + double lnPropR = 0.0; + double rdwindow = 0.2; + double oldRD = relativeDeath; + double newRD = getNewValSWindoMv(oldRD, 0.0, 0.99999, rdwindow); + relativeDeath = newRD; + double newtreeprob = getLnTreeProb(t); + double lnPriorRatio = (newtreeprob - oldtreeprob); + double lnR = lnPriorRatio + lnPropR; + double r = modelPtr->safeExponentiation(lnR); + if(ranPtr->uniformRv() < r){ + setAllBDFossParams(); + } + else{ + relativeDeath = oldRD; + setAllBDFossParams(); + } + return 0.0; + +} + +double Speciation::getNewValScaleMv(double &nv, double ov, double vmin, double vmax, double tv){ + + double rv = ranPtr->uniformRv(); + double c = tv * (rv - 0.5); + double newcv = ov * exp(c); + bool validV = false; + do{ + if(newcv < vmin) + newcv = vmin * vmin / newcv; + else if(newcv > vmax) + newcv = vmax * vmax / newcv; + else + validV = true; + } while(!validV); + nv = newcv; + return c; +} + +double Speciation::getNewValSWindoMv(double ov, double vmin, double vmax, double tv){ + double nv = 0.0; + double u = ranPtr->uniformRv(-0.5,0.5) * (tv); + nv = ov + u; bool validV = false; do{ - if(newRD < 0.0) - newRD = 0.0 - newRD; - else if(newRD > 0.9999) - newRD = (2 * 0.9999) - newRD; + if(nv < vmin) + nv = 2.0 * vmin - nv; + else if(nv > vmax) + nv = (2.0 * vmax) - nv; else validV = true; }while(!validV); - relativeDeath = newRD; - - double betA = 1.0; - double betB = 1.0; - double nv = ranPtr->lnBetaPdf(betA, betB, newRD); - double dv = ranPtr->lnBetaPdf(betA, betB, oldRD); - return nv - dv; //0.0; + return nv; } double Speciation::updateNetDivRate(void) { @@ -193,26 +235,36 @@ double Speciation::updateNetDivRate(void) { double oldND = netDiversificaton; double newND; double tuning = log(2.0); - double rv = ranPtr->uniformRv(); - double c = tuning * (rv - 0.5); - double exr = 10.0; - newND = oldND * exp(c); double minV = 0.0001; - double maxV = maxdivV; - bool validV = false; - do{ - if(newND < minV) - newND = minV * minV / newND; - else if(newND > maxV) - newND = maxV * maxV / newND; - else - validV = true; - } while(!validV); + double c = getNewValScaleMv(newND, oldND, minV, maxdivV, tuning); netDiversificaton = newND; lpr = c; - double num = ranPtr->lnExponentialPdf(exr, newND); - double dem = ranPtr->lnExponentialPdf(exr, oldND); - return lpr + (num - dem); + return lpr; +} + +double Speciation::updateNetDivRate(Tree *t) { + + double oldtreeprob = getLnTreeProb(t); + double lpr = 0.0; + double oldND = netDiversificaton; + double newND; + double tuning = log(2.0); + double minV = 0.0001; + double c = getNewValScaleMv(newND, oldND, minV, maxdivV, tuning); + netDiversificaton = newND; + lpr = c; + double newtreeprob = getLnTreeProb(t); + double lnPriorRatio = (newtreeprob - oldtreeprob); + double lnR = lnPriorRatio + lpr; + double r = modelPtr->safeExponentiation(lnR); + if(ranPtr->uniformRv() < r){ + setAllBDFossParams(); + } + else{ + netDiversificaton = oldND; + setAllBDFossParams(); + } + return 0.0; } double Speciation::updateBDSSSampleProbRho(void) { @@ -233,11 +285,7 @@ double Speciation::updateBDSSSampleProbRho(void) { validV = true; }while(!validV); extantSampleRate = newSP; - double betA = 1.0; - double betB = 1.0; - double nv = ranPtr->lnBetaPdf(betA, betB, newSP); - double dv = ranPtr->lnBetaPdf(betA, betB, oldSP); - return nv - dv; //0.0; + return 0.0; } double Speciation::updateBDSSFossilProbS(void) { @@ -258,15 +306,74 @@ double Speciation::updateBDSSFossilProbS(void) { validV = true; }while(!validV); probSpeciationS = newS; - double betA = 1.0; - double betB = 1.0; - double nv = ranPtr->lnBetaPdf(betA, betB, newS); - double dv = ranPtr->lnBetaPdf(betA, betB, oldS); - return nv - dv; //0.0; + return 0.0; +} + +double Speciation::updateBDSSFossilProbS(Tree *t) { + + + double oldtreeprob = getLnTreeProb(t); + double lnPropR = 0.0; + double swindow = 0.2; + double oldS = probSpeciationS; + double newS = getNewValSWindoMv(oldS, 0.0, 0.99999, swindow); + probSpeciationS = newS; + double newtreeprob = getLnTreeProb(t); + double lnPriorRatio = (newtreeprob - oldtreeprob); + double lnR = lnPriorRatio + lnPropR; + double r = modelPtr->safeExponentiation(lnR); + if(ranPtr->uniformRv() < r){ + setAllBDFossParams(); + } + else{ + probSpeciationS = oldS; + setAllBDFossParams(); + } + return 0.0; + +} + +double Speciation::updatePsiRate(Tree *t) { + + double lpr = 0.0; + double oldtreeprob = getLnTreeProb(t); + double oldPsi = fossilRate; + double newPsi; + double tuning = log(2.0); + double rv = ranPtr->uniformRv(); + double c = tuning * (rv - 0.5); + newPsi = oldPsi * exp(c); + double minV = 0.0001; + double maxV = 100.00; + bool validV = false; + do{ + if(newPsi < minV) + newPsi = minV * minV / newPsi; + else if(newPsi > maxV) + newPsi = newPsi * maxV / newPsi; + else + validV = true; + } while(!validV); + fossilRate = newPsi; + probSpeciationS = newPsi / (deathRate + newPsi) ; + lpr = c; + double newtreeprob = getLnTreeProb(t); + double lnPriorRatio = (newtreeprob - oldtreeprob); + double lnR = lnPriorRatio + lpr; + double r = modelPtr->safeExponentiation(lnR); + if(ranPtr->uniformRv() < r){ + setAllBDFossParams(); + } + else{ + fossilRate = oldPsi; + setAllBDFossParams(); + } + return 0.0; } + double Speciation::lnPrior(void) { return 0.0; @@ -325,5 +432,4 @@ void Speciation::setAllBDFossParams(){ fossilRate = (probSpeciationS / (1-probSpeciationS)) * ((relativeDeath * netDiversificaton) / (1 - relativeDeath)); birthRate = netDiversificaton / (1.0 - relativeDeath); deathRate = (relativeDeath * netDiversificaton) / (1 - relativeDeath); - -} \ No newline at end of file +} diff --git a/src/Parameter_speciaton.h b/src/Parameter_speciaton.h index 588d9c6..f09baaa 100644 --- a/src/Parameter_speciaton.h +++ b/src/Parameter_speciaton.h @@ -45,7 +45,7 @@ class Tree; class Speciation : public Parameter { public: - Speciation(MbRandom *rp, Model *mp, double bdr, double bda, double initRH); + Speciation(MbRandom *rp, Model *mp, double bdr, double bda, double bds, double initRH); ~Speciation(void); Speciation &operator=(const Speciation &c); void clone(const Speciation &c); @@ -76,8 +76,11 @@ class Speciation : public Parameter { double deathRate; // mu double fossilRate; // psi double extantSampleRate; // rho = 1 - double fossilStratSampleProb; // + double fossilStratSampleProb; // omega = ? This is for later double probSpeciationS; // \psi / (\mu + \psi) // Need hyperprior + + double getNewValScaleMv(double &nv, double ov, double vmin, double vmax, double tv); + double getNewValSWindoMv(double ov, double vmin, double vmax, double tv); @@ -85,7 +88,12 @@ class Speciation : public Parameter { double updateNetDivRate(); double updateBDSSFossilProbS(); double updateBDSSSampleProbRho(); - + + double updateRelDeathRt(Tree *t); + double updateNetDivRate(Tree *t); + double updateBDSSFossilProbS(Tree *t); + double updatePsiRate(Tree *t); + }; diff --git a/src/Parameter_tree.cpp b/src/Parameter_tree.cpp index 1d2a476..59d124e 100644 --- a/src/Parameter_tree.cpp +++ b/src/Parameter_tree.cpp @@ -90,7 +90,7 @@ Tree::Tree(MbRandom *rp, Model *mp, Alignment *ap, string ts, bool ubl, bool all useInputBLs = ubl; moveAllNodes = allnm; calibNds = clb; - numCalibNds = calibNds.size(); + numCalibNds = (int)calibNds.size(); numAncFossilsk = 0; datedTips = tdt; treeScale = rth; @@ -99,9 +99,14 @@ Tree::Tree(MbRandom *rp, Model *mp, Alignment *ap, string ts, bool ubl, bool all randShufNdMv = rndNods; softBounds = sb; isCalibTree = false; - isTipCals = false; // indicates if there are non-contemp tips + isTipCals = false; expHyperPrCal = exhpc; buildTreeFromNewickDescription(ts); + buildTreeFromNewickDescription(ts); + + + nodeProposal = 2; // proposal type 1=window, 2=scale, 3=slide + if(numCalibNds > 0 && treeTimePrior < 6){ isCalibTree = true; @@ -174,11 +179,10 @@ void Tree::buildTreeFromNewickDescription(string ts) { p->setLft(q); else if (p->getRht() == NULL) p->setRht(q); - else - { + else{ cerr << "ERROR: Problem adding interior node to tree" << endl; exit(1); - } + } } p = q; readingBrlen = false; @@ -193,30 +197,26 @@ void Tree::buildTreeFromNewickDescription(string ts) { readingBrlen = false; } else if ( c == ',' ){ - if (p->getAnc() == NULL) - { + if (p->getAnc() == NULL){ cerr << "ERROR: Problem moving down the tree" << endl; exit(1); - } + } else p = p->getAnc(); readingBrlen = false; } else if ( c == ':' ){ - // we are now reading a branch length... readingBrlen = true; } else if ( c == ';' ){ // we are finished with the tree...check } else{ - // read the taxon name or branch length into a string string s = ""; while ( isValidChar(ts[i]) ) s += ts[i++]; i--; if (readingBrlen == false){ - // add a tip to the tree if ( alignmentPtr->isTaxonPresent(s) == false ){ cerr << "ERROR: Cannot find taxon in alignment" << endl; exit(1); @@ -243,7 +243,6 @@ void Tree::buildTreeFromNewickDescription(string ts) { p->setIsLeaf(true); } else{ - // include the branch length information double x; istringstream buf(s); buf >> x; @@ -252,11 +251,9 @@ void Tree::buildTreeFromNewickDescription(string ts) { } } - // get the down pass sequence getDownPassSequence(); - root->setRtGrpVal(0.5); // TAH root rate debug + root->setRtGrpVal(0.5); - // initialize the node depths } Tree& Tree::operator=(const Tree &t) { @@ -503,7 +500,6 @@ void Tree::initializeNodeDepths(void) { Node *p = potentialSplits[(int)(ranPtr->uniformRv()*potentialSplits.size())]; p->setNodeDepth(nodeTimes[nextTime++]); - // remove the node from the list of nodes that can split for (vector::iterator n=potentialSplits.begin(); n != potentialSplits.end(); n++){ if ( (*n) == p ){ potentialSplits.erase( n ); @@ -511,7 +507,6 @@ void Tree::initializeNodeDepths(void) { } } - // add in the nodes that descend from it if (p->getLft()->getIsLeaf() == false) potentialSplits.push_back(p->getLft()); if (p->getRht()->getIsLeaf() == false) @@ -581,6 +576,19 @@ void Tree::initializeCalibratedNodeDepths(void) { recursiveNodeDepthInitialization(root, nnc, 1.0); } + +void Tree::checkNodeInit(){ + + for(int i=0; igetFossilMRCANodeID(); + Node *p = &nodes[nID]; + double nage = p->getNodeDepth() * treeScale; + double fage = f->getFossilAge(); + cout << nID << " -- " << nage - fage << " --- " << nage << " --- " << fage << endl; + } +} + double Tree::getTemporaryNodeMaxBound(Node *p){ double tempmax = treeScale; @@ -746,10 +754,11 @@ void Tree::print(std::ostream & o) const { } double Tree::update(double &oldLnL) { - + + double lppr = 0.0; double probCalMove = 0.5; - if(treeTimePrior == 6){ + if(treeTimePrior == 6){ recountFossilAttachNums(); if(ranPtr->uniformRv() > probCalMove){ if(moveAllNodes){ @@ -772,42 +781,29 @@ double Tree::update(double &oldLnL) { } recountFossilAttachNums(); } - if(treeTimePrior == 7){ - recountFossilAttachNums(); - double *probMove = new double[4]; - probMove[0] = 0.4; - probMove[1] = 0.4; - probMove[2] = 0.2; - int mvType = ranPtr->categoricalRv(probMove,3); - if(mvType == 0) - lppr = updateAllNodes(oldLnL); - else if(mvType == 1){ - updateFossilBDSSAttachmentTimePhi(); - modelPtr->setLnLGood(true); - modelPtr->setMyCurrLnl(oldLnL); - Tree *t = modelPtr->getActiveTree(); - t->upDateAllCls(); - t->upDateAllTis(); - modelPtr->setTiProb(); - return 0.0; - } - else if(mvType == 2){ + else if(treeTimePrior == 7){ + + Treescale *ts = modelPtr->getActiveTreeScale(); + setTreeScale(ts->getScaleValue()); + + updateAllTGSNodes(oldLnL); + + updateFossilBDSSAttachmentTimePhi(); + modelPtr->setLnLGood(true); + modelPtr->setMyCurrLnl(oldLnL); + modelPtr->setTiProb(); + + for(int i=0; i<5; i++){ updateRJMoveAddDelEdge(); modelPtr->setLnLGood(true); modelPtr->setMyCurrLnl(oldLnL); - Tree *t = modelPtr->getActiveTree(); - t->upDateAllCls(); - t->upDateAllTis(); modelPtr->setTiProb(); - return 0.0; } - recountFossilAttachNums(); + return 0.0; } else{ if(moveAllNodes){ - if(randShufNdMv) - lppr = updateAllNodesRnd(oldLnL); - else lppr = updateAllNodes(oldLnL); + updateAllNodes(oldLnL); } else lppr = updateOneNode(); setAllNodeBranchTimes(); @@ -824,46 +820,68 @@ double Tree::updateFossilBDSSAttachmentTimePhi() { double mu = s->getBDSSExtinctionRateMu(); double fossRate = s->getBDSSFossilSampRatePsi(); double sppSampRate = s->getBDSSSppSampRateRho(); - recountFossilAttachNums(); - for(int i=0; i rndFossIDs; + for(int i=0; i::iterator it=rndFossIDs.begin(); it!=rndFossIDs.end(); it++){ + Fossil *f = fossSpecimens[(*it)]; if(f->getFossilIndicatorVar()){ + + Node *p = &nodes[f->getFossilMRCANodeID()]; - double nodeDepth = p->getNodeDepth(); - double fossDepth = f->getFossilAge() / treeScale; - double oldPhi = f->getFossilSppTime(); - double oldSumLogGammas = getProposedBranchAttachFossils(p, oldPhi); - double newPhi = fossDepth + ranPtr->uniformRv()*(nodeDepth-fossDepth); - f->setFossilSppTime(newPhi); - double newSumLogGammas = getProposedBranchAttachFossils(p, newPhi); + double nodeDepth = p->getNodeDepth() * treeScale; + double fossDepth = f->getFossilAge(); + + double oldPhi = f->getFossilSppTime() * treeScale; + double oldSumLogGammas = getSumLogAllAttachNums(); + + double rv = ranPtr->uniformRv(); + double newPhi, c; + if(nodeProposal == 1){ + double delta = 5.0; + c = doAWindoMove(newPhi, oldPhi, delta, fossDepth, nodeDepth, rv); + } + else if(nodeProposal == 2){ + double tv = tuningVal; + c = doAScaleMove(newPhi, oldPhi, tv, fossDepth, nodeDepth, rv); + } + else if(nodeProposal == 3){ + newPhi = fossDepth + rv*(nodeDepth-fossDepth); + c = 0.0; + } + + + f->setFossilSppTime(newPhi/treeScale); + + double newSumLogGammas = getSumLogAllAttachNums(); + double lnPriorRat = 0.0; double v1 = newSumLogGammas; - double v2 = bdssQFxn(lambda, mu, fossRate, sppSampRate, newPhi*treeScale); - double v3 = bdssQFxn(lambda, mu, fossRate, sppSampRate, oldPhi*treeScale); + double v2 = bdssQFxn(lambda, mu, fossRate, sppSampRate, newPhi); + double v3 = bdssQFxn(lambda, mu, fossRate, sppSampRate, oldPhi); double v4 = oldSumLogGammas; lnPriorRat = (v1 - v2) - (v4 - v3); - double r = modelPtr->safeExponentiation(lnPriorRat); + + + double r = modelPtr->safeExponentiation(lnPriorRat + c); if(ranPtr->uniformRv() < r){ - f->setFossilSppTime(newPhi); - recountFossilAttachNums(); + f->setFossilSppTime(newPhi/treeScale); setNodeOldestAttchBranchTime(p); } else{ - f->setFossilSppTime(oldPhi); - recountFossilAttachNums(); + f->setFossilSppTime(oldPhi/treeScale); setNodeOldestAttchBranchTime(p); } } } - return 0.0; } double Tree::updateRJMoveAddDelEdge() { - recountFossilAttachNums(); double gA = 0.5; if(numAncFossilsk == numCalibNds) gA = 1.0; @@ -895,24 +913,26 @@ void Tree::doAddEdgeMove(){ int mvFoss = pickRandAncestorFossil(); Fossil *f = fossSpecimens[mvFoss]; Node *p = &nodes[f->getFossilMRCANodeID()]; - double oldSumLogGammas = getProposedBranchAttachFossils(p, 0.0); - + double oldSumLogGammas = getSumLogAllAttachNums(); + double alA = 1.0; if(kN == 0) alA = 2.0; else if(k == m) alA = 0.5; - + double cf = p->getNodeDepth() * treeScale; double yf = f->getFossilAge(); double nu = ranPtr->uniformRv() * (cf - yf); - lnHastings = log(alA) + log(k) - log(m - k + 1.0); + lnHastings = log(alA) + (log(k) - log(m - k + 1.0)); lnJacobian = log(cf - yf); double newPhi = yf + nu; double scnewPhi = newPhi / treeScale; f->setFossilSppTime(scnewPhi); - - double newSumLogGammas = getProposedBranchAttachFossils(p, 0.0); + f->setFossilIndicatorVar(1); + + double newSumLogGammas = getSumLogAllAttachNums(); + double v1 = log(bdssP0Fxn(lambda, mu, fossRate, sppSampRate, yf)); double v2 = bdssQFxn(lambda, mu, fossRate, sppSampRate, yf); double v3 = bdssQFxn(lambda, mu, fossRate, sppSampRate, newPhi); @@ -927,20 +947,19 @@ void Tree::doAddEdgeMove(){ if(ranPtr->uniformRv() < r){ f->setFossilIndicatorVar(1); f->setFossilSppTime(scnewPhi); - recountFossilAttachNums(); setNodeOldestAttchBranchTime(p); numAncFossilsk = kN; } else{ f->setFossilIndicatorVar(0); - f->setFossilFossBrGamma(0); f->setFossilSppTime(yf/treeScale); - recountFossilAttachNums(); + numAncFossilsk = k; setNodeOldestAttchBranchTime(p); } } void Tree::doDeleteEdgeMove(){ + int k = numAncFossilsk; int m = numCalibNds; @@ -956,7 +975,7 @@ void Tree::doDeleteEdgeMove(){ int mvFoss = pickRandTipFossil(); Fossil *f = fossSpecimens[mvFoss]; Node *p = &nodes[f->getFossilMRCANodeID()]; - double oldSumLogGammas = getProposedBranchAttachFossils(p, 0.0); + double oldSumLogGammas = getSumLogAllAttachNums(); double alD = 1.0; if(k == 0) @@ -972,9 +991,10 @@ void Tree::doDeleteEdgeMove(){ double scnewPhi = yf / treeScale; f->setFossilIndicatorVar(0); f->setFossilSppTime(scnewPhi); - double newSumLogGammas = getProposedBranchAttachFossils(p, 0.0); + + double newSumLogGammas = getSumLogAllAttachNums(); - lnHastings = log(alD) + log(m-k) - log(k+1); + lnHastings = log(alD) + (log(m-k) - log(k+1)); lnJacobian = -(log(cf - yf)); double v1 = log(bdssP0Fxn(lambda, mu, fossRate, sppSampRate, yf)); @@ -989,16 +1009,13 @@ void Tree::doDeleteEdgeMove(){ if(ranPtr->uniformRv() < r){ f->setFossilIndicatorVar(0); - f->setFossilFossBrGamma(0); f->setFossilSppTime(scnewPhi); - recountFossilAttachNums(); setNodeOldestAttchBranchTime(p); numAncFossilsk = kN; } else{ f->setFossilSppTime(oldScPhi); f->setFossilIndicatorVar(1); - recountFossilAttachNums(); setNodeOldestAttchBranchTime(p); } } @@ -1052,23 +1069,36 @@ double Tree::updateAllNodes(double &oldLnL) { upDateAllTis(); double oldLike = oldLnL; Node *p = NULL; - for(int i = 0; i rndNodeIDs; + for(int i=0; i::iterator it=rndNodeIDs.begin(); it!=rndNodeIDs.end(); it++){ + p = downPassSequence[(*it)]; if(p != root && !p->getIsLeaf()){ - double currDepth = p->getNodeDepth(); + double currDepth = p->getNodeDepth() * treeScale; - double largestTime = getNodeUpperBoundTime(p); - double smallestTime = getNodeLowerBoundTime(p); - + double largestTime = getNodeUpperBoundTime(p) * treeScale; + double smallestTime = getNodeLowerBoundTime(p) * treeScale; if (largestTime > smallestTime){ - double newNodeDepth = smallestTime + ranPtr->uniformRv()*(largestTime-smallestTime); - double lnPrRatio = 0.0; - if(treeTimePrior > 5) - lnPrRatio = lnPriorRatioTGS(newNodeDepth, currDepth, p); - else lnPrRatio = lnPriorRatio(newNodeDepth, currDepth); + double rv = ranPtr->uniformRv(); + double newNodeDepth, c; + if(nodeProposal == 1){ + double delta = 5.0; + c = doAWindoMove(newNodeDepth, currDepth, delta, smallestTime, largestTime, rv); + } + else if(nodeProposal == 2){ + double tv = tuningVal; + c = doAScaleMove(newNodeDepth, currDepth, tv, smallestTime, largestTime, rv); + } + else if(nodeProposal == 3){ + newNodeDepth = smallestTime + rv*(largestTime-smallestTime); + c = 0.0; + } - p->setNodeDepth(newNodeDepth); + double lnPrRatio = lnPriorRatio(newNodeDepth/treeScale, currDepth/treeScale); + p->setNodeDepth(newNodeDepth/treeScale); flipToRootClsTis(p); updateToRootClsTis(p); @@ -1076,38 +1106,129 @@ double Tree::updateAllNodes(double &oldLnL) { double newLnl = modelPtr->lnLikelihood(); double lnLRatio = newLnl - oldLike; - if(p->getIsCalibratedDepth()){ + if(p->getIsCalibratedDepth()){ if(softBounds && p->getNodeCalibPrDist() == 1){ double ycal = p->getNodeYngTime(); double ocal = p->getNodeOldTime(); - lnPrRatio += lnCalibPriorRatio(newNodeDepth*treeScale, currDepth*treeScale, ycal, ocal); + lnPrRatio += lnCalibPriorRatio(newNodeDepth, currDepth, ycal, ocal); } if(p->getNodeCalibPrDist() == 2){ double offst = p->getNodeYngTime(); double nodeExCR = p->getNodeExpCalRate(); - lnPrRatio += lnExpCalibPriorRatio(newNodeDepth*treeScale, currDepth*treeScale, offst, nodeExCR); + lnPrRatio += lnExpCalibPriorRatio(newNodeDepth, currDepth, offst, nodeExCR); } } - double lnR = lnPrRatio + lnLRatio + 0.0; + double lnR = lnPrRatio + lnLRatio + c; double r = modelPtr->safeExponentiation(lnR); if(ranPtr->uniformRv() < r){ oldLike = newLnl; } else{ - p->setNodeDepth(currDepth); + p->setNodeDepth(currDepth/treeScale); flipToRootClsTis(p); updateToRootClsTis(p); } } } - if(treeTimePrior > 5) + } + oldLnL = oldLike; + return 0.0; +} + +double Tree::updateAllTGSNodes(double &oldLnL) { + + upDateAllCls(); + upDateAllTis(); + double oldLike = oldLnL; + Node *p = NULL; + vector rndNodeIDs; + for(int i=0; i::iterator it=rndNodeIDs.begin(); it!=rndNodeIDs.end(); it++){ + p = downPassSequence[(*it)]; + if(p != root && !p->getIsLeaf()){ + double currDepth = p->getNodeDepth() * treeScale; + + double largestTime = getNodeUpperBoundTime(p) * treeScale; + double smallestTime = getNodeLowerBoundTime(p) * treeScale; + + if (largestTime > smallestTime){ + double rv = ranPtr->uniformRv(); + double newNodeDepth, c; + if(nodeProposal == 1){ + double delta = 5.0; + c = doAWindoMove(newNodeDepth, currDepth, delta, smallestTime, largestTime, rv); + } + else if(nodeProposal == 2){ + double tv = tuningVal; + c = doAScaleMove(newNodeDepth, currDepth, tv, smallestTime, largestTime, rv); + } + else if(nodeProposal == 3){ + newNodeDepth = smallestTime + rv*(largestTime-smallestTime); + c = 0.0; + } + + double lnPrRatio = lnPriorRatioTGS(newNodeDepth, currDepth, p); + p->setNodeDepth(newNodeDepth/treeScale); + + flipToRootClsTis(p); + updateToRootClsTis(p); + modelPtr->setTiProb(); + double newLnl = modelPtr->lnLikelihood(); + double lnLRatio = newLnl - oldLike; + + double lnR = lnPrRatio + lnLRatio + c; + double r = modelPtr->safeExponentiation(lnR); + + if(ranPtr->uniformRv() < r){ + oldLike = newLnl; + } + else{ + p->setNodeDepth(currDepth/treeScale); + flipToRootClsTis(p); + updateToRootClsTis(p); + } + } recountFossilAttachNums(); + } } oldLnL = oldLike; return 0.0; } +double Tree::doAScaleMove(double &nv, double cv, double tv, double lb, double hb, double rv){ + + double c = tv * (rv - 0.5); + double newcv = cv * exp(c); + bool validV = false; + do{ + if(newcv < lb) + newcv = lb * lb / newcv; + else if(newcv > hb) + newcv = hb * hb / newcv; + else + validV = true; + } while(!validV); + nv = newcv; + return c; +} + +double Tree::doAWindoMove(double &nv, double cv, double tv, double lb, double hb, double rv){ + + double newCV = cv + (tv * (rv - 0.5)); + do{ + if(newCV < lb) + newCV = 2.0 * lb - newCV; + else if(newCV > hb) + newCV = 2.0 * hb - newCV; + }while(newCV < lb || newCV > hb); + nv = newCV; + return 0.0; +} + + double Tree::updateAllNodesRnd(double &oldLnL) { upDateAllCls(); @@ -1157,8 +1278,9 @@ double Tree::updateAllNodesRnd(double &oldLnL) { double lnR = lnPrRatio + lnLRatio + 0.0; double r = modelPtr->safeExponentiation(lnR); - if(ranPtr->uniformRv() < r) + if(ranPtr->uniformRv() < r){ oldLike = newLnl; + } else{ p->setNodeDepth(currDepth); flipToRootClsTis(p); @@ -1247,12 +1369,14 @@ double Tree::lnPriorRatio(double snh, double soh) { double Tree::lnPriorRatioTGS(double snh, double soh, Node *p) { - double nh = snh*treeScale; - double oh = soh*treeScale; - double oldSumLogGammas = getProposedBranchAttachFossils(p, 0.0); + double nh = snh; + double oh = soh; - p->setNodeDepth(snh); - double newSumLogGammas = getProposedBranchAttachFossils(p, 0.0); + p->setNodeDepth(soh/treeScale); + double oldSumLogGammas = getSumLogAllAttachNums(); + + p->setNodeDepth(snh/treeScale); + double newSumLogGammas = getSumLogAllAttachNums(); Speciation *s = modelPtr->getActiveSpeciation(); s->setAllBDFossParams(); @@ -1287,6 +1411,7 @@ double Tree::bdssQFxn(double b, double d, double psi, double rho, double t){ double c2Val = bdssC2Fxn(b,d,psi,rho); double vX = c1Val * t + 2.0 * log(exp(-c1Val * t) * (1.0 - c2Val) + (1.0 + c2Val)); + return vX; } @@ -1296,11 +1421,17 @@ double Tree::bdssP0Fxn(double b, double d, double psi, double rho, double t){ double c2Val = bdssC2Fxn(b,d,psi,rho); double eCfrac = (exp(-c1Val * t) * (1.0 - c2Val) - (1.0 + c2Val)) / (exp(-c1Val * t) * (1.0 - c2Val) + (1.0 + c2Val)); - double v = ((b + d + psi) + (c1Val * eCfrac)) / (2.0 * b); + double v = 1.0 + ((-(b - d - psi)) + (c1Val * eCfrac)) / (2.0 * b); return v; } +double Tree::bdssP0HatFxn(double b, double d, double rho, double t){ + + double v = 1.0 - ((rho * (b - d)) / ((b*rho) + (((b * (1-rho)) - d) * exp(-(b-d)*t)))); + return v; +} + double Tree::lnCalibPriorRatio(double nh, double oh, double lb, double ub) { @@ -1332,8 +1463,12 @@ double Tree::lnCalibPriorRatio(double nh, double oh, double lb, double ub) { double Tree::lnExpCalibPriorRatio(double nh, double oh, double offSt, double expRate) { + /* + TAH: Prior ratio for offset exponential + */ double numr = 0.0; double dnom = 0.0; + numr = ranPtr->lnExponentialPdf(expRate, nh - offSt); dnom = ranPtr->lnExponentialPdf(expRate, oh - offSt); return numr - dnom; @@ -1397,6 +1532,7 @@ void Tree::updateToRootClsTis(Node *p){ } } +// overloaded function when passing in just the node ID, then this is called by the DPP rate move. void Tree::flipToRootClsTis(int ndID){ Node *q = &nodes[ndID]; @@ -1721,7 +1857,7 @@ string Tree::getCalBDSSNodeInfoParamList(void){ } for(int i=0; igetFossilIndicatorVar()) + if(1) ss << "\t" << f->getFossilSppTime() * treeScale; else ss << "\t" << f->getFossilAge(); } @@ -1951,14 +2087,14 @@ double Tree::getTreeCBDNodePriorProb() { nprb = getTreeCBDNodePriorProb(diff, 0.0); return nprb; } - else if(treeTimePrior == 3){ // Gernhard (2008) conditioned birth death + else if(treeTimePrior == 3){ Speciation *s = modelPtr->getActiveSpeciation(); double diff = s->getNetDiversification(); - double rel = s->getRelativeDeath(); + double rel = s->getRelativeDeath(); nprb = getTreeCBDNodePriorProb(diff, rel); return nprb; } - else if(treeTimePrior == 4 || treeTimePrior == 5){ // Stadler (2010) code from BEAST + else if(treeTimePrior == 4 || treeTimePrior == 5){ Speciation *s = modelPtr->getActiveSpeciation(); Treescale *ts = modelPtr->getActiveTreeScale(); double div = s->getNetDiversification(); @@ -1970,7 +2106,6 @@ double Tree::getTreeCBDNodePriorProb() { return nprb; } else if(treeTimePrior == 6){ - // fossilized birth-death Speciation *s = modelPtr->getActiveSpeciation(); s->setAllBDFossParams(); double lambda = s->getBDSSSpeciationRateLambda(); @@ -1981,7 +2116,6 @@ double Tree::getTreeCBDNodePriorProb() { return nprb; } else if(treeTimePrior == 7){ - // fossilized birth-death Speciation *s = modelPtr->getActiveSpeciation(); s->setAllBDFossParams(); double lambda = s->getBDSSSpeciationRateLambda(); @@ -2001,7 +2135,7 @@ double Tree::getTreeCBDNodePriorProb(double netDiv, double relDeath) { if(treeTimePrior == 1) return 0.0; else if(treeTimePrior == 2){ - double diff = netDiv; + double diff = netDiv; for(int i=0; igetIsLeaf() == false){ @@ -2015,9 +2149,9 @@ double Tree::getTreeCBDNodePriorProb(double netDiv, double relDeath) { } return nprb; } - else if(treeTimePrior == 3){ // Gernhard (2008) conditioned birth death + else if(treeTimePrior == 3){ double diff = netDiv; - double rel = relDeath; + double rel = relDeath; for(int i=0; igetIsLeaf() == false){ @@ -2032,7 +2166,7 @@ double Tree::getTreeCBDNodePriorProb(double netDiv, double relDeath) { } return nprb; } - else if(treeTimePrior == 4 || treeTimePrior == 5){ // Stadler (2010) + else if(treeTimePrior == 4 || treeTimePrior == 5){ Speciation *s = modelPtr->getActiveSpeciation(); Treescale *ts = modelPtr->getActiveTreeScale(); double fossRate = s->getBDSSFossilSampRatePsi(); @@ -2044,7 +2178,6 @@ double Tree::getTreeCBDNodePriorProb(double netDiv, double relDeath) { return nprb; } else if(treeTimePrior == 6){ - // fossilized birth-death Speciation *s = modelPtr->getActiveSpeciation(); s->setAllBDFossParams(); double lambda = s->getBDSSSpeciationRateLambda(); @@ -2055,7 +2188,6 @@ double Tree::getTreeCBDNodePriorProb(double netDiv, double relDeath) { return nprb; } else if(treeTimePrior == 7){ - // fossilized birth-death Speciation *s = modelPtr->getActiveSpeciation(); s->setAllBDFossParams(); double lambda = s->getBDSSSpeciationRateLambda(); @@ -2084,16 +2216,16 @@ double Tree::getTreeBDSSTreeNodePriorProb(double netDiv, double relDeath, double if(p->getIsCalibratedDepth()){ double nh = p->getNodeDepth() * treeScale; nprb += log(fossRate) + bdssQFxn(lambda,mu,fossRate,sppSampRate,nh); - nprb += log(bdssP0Fxn(lambda,mu,fossRate,sppSampRate,nh)); + nprb += log(bdssP0Fxn(lambda,mu,fossRate,sppSampRate,nh)); } } - else if(!p->getIsLeaf()){ + else if(!p->getIsLeaf()){ double nh = p->getNodeDepth() * treeScale; nprb += log(lambda) - bdssQFxn(lambda,mu,fossRate,sppSampRate,nh); } } } - else if(treeTimePrior == 5){ + else if(treeTimePrior == 5){ nprb = ((numTaxa - 2.0) * log(lambda)) - (2.0 * log(1.0 - bdssP0Fxn(lambda, mu, 0.0, sppSampRate, treeScale))); nprb += (numExtantTips * log(4.0 * sppSampRate)) - bdssQFxn(lambda,mu,fossRate,sppSampRate,treeScale); for(int i=0; igetIsLeaf()){ + else if(!p->getIsLeaf()){ double nh = p->getNodeDepth() * treeScale; nprb += -bdssQFxn(lambda,mu,fossRate,sppSampRate,nh); } @@ -2146,30 +2278,31 @@ double Tree::getTreeAncCalBDSSTreeNodePriorProb(double lambda, double mu, double double nprb = 0.0; recountFossilAttachNums(); - - nprb = ( ((numTaxa - 2 + numCalibNds - numAncFossilsk) * log(lambda)) + (numCalibNds * log(fossRate)) ); - nprb -= ( 2.0 * log(1.0 - bdssP0Fxn(lambda, mu, 0.0, sppSampRate, treeScale)) ); - nprb += ( numTaxa * log(4.0 * sppSampRate) ) - bdssQFxn(lambda, mu, fossRate, sppSampRate, treeScale); - - + + nprb = -(2.0 * (log(lambda) + log(1.0 - bdssP0HatFxn(lambda, mu, sppSampRate, treeScale)))); + nprb += (log(4.0 * lambda * sppSampRate)) - bdssQFxn(lambda, mu, fossRate, sppSampRate, treeScale); + for(int i=0; igetIsLeaf()){ double myAge = p->getNodeDepth() * treeScale; - nprb += -bdssQFxn(lambda, mu, fossRate, sppSampRate, myAge); + nprb += log(2.0 * 4.0 * lambda * sppSampRate) - bdssQFxn(lambda, mu, fossRate, sppSampRate, myAge); } } for(int i=0; igetFossilFossBrGamma() ); if(f->getFossilIndicatorVar()){ double fossAge = f->getFossilAge(); double fossPhi = f->getFossilSppTime() * treeScale; - nprb += log(2.0 * f->getFossilFossBrGamma()) + bdssQFxn(lambda, mu, fossRate, sppSampRate, fossAge); - nprb += log( bdssP0Fxn(lambda, mu, fossRate, sppSampRate, fossAge) ); - nprb -= bdssQFxn(lambda, mu, fossRate, sppSampRate, fossPhi); + double fossPr = log(2.0 * lambda) + bdssQFxn(lambda, mu, fossRate, sppSampRate, fossAge); + fossPr += log( bdssP0Fxn(lambda, mu, fossRate, sppSampRate, fossAge) ); + fossPr -= bdssQFxn(lambda, mu, fossRate, sppSampRate, fossPhi); + nprb += fossPr; } } + return nprb; } @@ -2178,24 +2311,24 @@ double Tree::getTreeSpeciationProbability() { if(treeTimePrior == 1) return 0.0; - else if(treeTimePrior == 2){ // Yule + else if(treeTimePrior == 2){ Speciation *s = modelPtr->getActiveSpeciation(); double diff = s->getNetDiversification(); double c1 = (numTaxa - 1) * log(diff); double nps = getTreeCBDNodePriorProb(diff, 0.0); return c1 + nps; } - else if(treeTimePrior == 3){ // Gernhard (2008) conditioned birth death + else if(treeTimePrior == 3){ Speciation *s = modelPtr->getActiveSpeciation(); - double diff = s->getNetDiversification(); + double diff = s->getNetDiversification(); double rel = s->getRelativeDeath(); - double lnC = 0.0; + double lnC = 0.0; double c1 = ((numTaxa - 1) * log(diff)) + (numTaxa * log(1 - rel)); double nps = getTreeCBDNodePriorProb(diff, rel); return lnC + c1 + nps; } else if(treeTimePrior > 3){ - double lnC = 0.0; + double lnC = 0.0; double nps = getTreeCBDNodePriorProb(); return lnC + nps; } @@ -2473,7 +2606,6 @@ void Tree::initializeTGSCalibVariables(){ // depricated int Tree::countDecLinsTimeIntersect(Node *p, double t, double ancAge){ int n = 0; - double myAge = p->getNodeDepth(); if(myAge < t) return 1; @@ -2503,8 +2635,6 @@ double Tree::getNodeLowerBoundTime(Node *p){ } } - - return t; } @@ -2519,8 +2649,6 @@ double Tree::getNodeUpperBoundTime(Node *p){ if(ocal < t) t = ocal; } - - return t; } @@ -2537,7 +2665,7 @@ void Tree::setUPTGSCalibrationFossils() { int nPFoss = p->getNumCalibratingFossils(); p->setNumFCalibratingFossils(nPFoss + 1); - if(nPFoss > 1){ + if(nPFoss > 0){ if(p->getNodeYngTime() < fAge) p->setNodeYngTime(fAge); } @@ -2547,6 +2675,7 @@ void Tree::setUPTGSCalibrationFossils() { (*v)->setNodeIndex(calbNo); p->insertFossilID(calID); + cout << calbNo << " --- " << p->getNodeYngTime() << endl; calID += 1; } } @@ -2587,20 +2716,79 @@ void Tree::initializeFossilSpecVariables(){ -void Tree::recountFossilAttachNums(){ +int Tree::recountFossilAttachNums(){ + int numChanged = 0; for(int i=0; igetFossilFossBrGamma(); double fiPhi = fi->getFossilSppTime(); Node *p = &nodes[fi->getFossilMRCANodeID()]; int fiGamma = getDecendantFossilAttachBranches(p, fiPhi, i); + if(fiGamma != fgammCur) + numChanged += 1; fi->setFossilFossBrGamma(fiGamma); fi->setFossilMRCANodeAge(p->getNodeDepth() * treeScale); } + return numChanged; +} + +void Tree::recountFromNodeFossilAttchNums(Node *p){ + + vector tempFoss; + recursivCreateTempFossVec(tempFoss, p); + + for(int i=0; igetFossilFossBrGamma(); + double fiPhi = fi->getFossilSppTime(); + Node *p = &nodes[fi->getFossilMRCANodeID()]; + int fiGamma = getDecendantFossilAttachBranches(p, fiPhi, i); + fi->setFossilFossBrGamma(fiGamma); + fi->setFossilMRCANodeAge(p->getNodeDepth() * treeScale); + } +} + +int Tree::recursivCreateTempFossVec(vector &v, Node *p){ + + if(p->getIsLeaf()) + return 0; + else{ + for(int i=0; igetNumCalibratingFossils(); i++){ + int idx = p->getIthFossiID(i); + Fossil *f = fossSpecimens[idx]; + v.push_back(f); + } + recursivCreateTempFossVec(v, p->getLft()); + recursivCreateTempFossVec(v, p->getRht()); + return 0; + } } +double Tree::getSumLogAllAttachNums(){ + + double sumLog = 0.0; + recountFossilAttachNums(); + for(int i=0; igetFossilFossBrGamma()); + } + return sumLog; +} + +int Tree::getSumIndicatorV(){ + if(treeTimePrior < 6) + return 0; + int niv = 0; + for(int i=0; igetFossilIndicatorVar(); + } + return niv; +} -int Tree::getFossilLinAttachNumsForFoss(int fID){ + +int Tree::getFossilLinAttachNumsForFoss(int fID){ // TAH: this is never called Fossil *fi = fossSpecimens[fID]; double fiPhi = fi->getFossilSppTime(); @@ -2613,8 +2801,6 @@ int Tree::getFossilLinAttachNumsForFoss(int fID){ } - - int Tree::getDecendantFossilAttachBranches(Node *p, double t, int fID){ int n = 0; @@ -2651,12 +2837,8 @@ double Tree::getProposedBranchAttachFossils(Node *p, double t){ int idx = p->getIthFossiID(i); int a = 0; Fossil *f = fossSpecimens[idx]; - if(f->getFossilIndicatorVar()){ - a = getDecendantFossilAttachBranches(p,f->getFossilSppTime(),idx); - sumLog += log(a); - } - - + a = getDecendantFossilAttachBranches(p,f->getFossilSppTime(),idx); + sumLog += log(a); } sumLog += getProposedBranchAttachFossils(p->getLft(), t); sumLog += getProposedBranchAttachFossils(p->getRht(), t); @@ -2664,7 +2846,7 @@ double Tree::getProposedBranchAttachFossils(Node *p, double t){ return sumLog; } -void Tree::treeScaleUpdateFossilAttchTimes(double sr){ +void Tree::treeScaleUpdateFossilAttchTimes(double sr, double ort, double nrt){ for(int i=0; i recursiveNodeDepthInitialization(Node *p, int &nCont, double maxD); @@ -299,13 +303,14 @@ class Tree : public Parameter { double bdssQFxn(double b, double d, double psi, double rho, double t); // on log scale double bdssP0Fxn(double b, double d, double psi, double rho, double t); + double bdssP0HatFxn(double b, double d, double rho, double t); void initializeTGSCalibVariables(void); void setUPTGSCalibrationFossils(); void initializeFossilSpecVariables(); void initializeFossilAncestorSpecVariables(); - void recountFossilAttachNums(); + int recountFossilAttachNums(); int getFossilLinAttachNumsForFoss(int fID); int getDecendantFossilAttachBranches(Node *p, double t, int fID); @@ -313,8 +318,17 @@ class Tree : public Parameter { void setNodeOldestAttchBranchTime(Node *p); int pickRandAncestorFossil(void); int pickRandTipFossil(void); - - + double getSumLogAllAttachNums(void); + double doAScaleMove(double &nv, double cv, double tv, double lb, double hb, double rv); + double doAWindoMove(double &nv, double cv, double tv, double lb, double hb, double rv); + int getNumDecFossils(Node *p); + + void checkNodeInit(void); + + void recountFromNodeFossilAttchNums(Node *p); + int recursivCreateTempFossVec(std::vector &v, Node *p); + + Alignment *alignmentPtr; int numTaxa; int numNodes; @@ -335,10 +349,15 @@ class Tree : public Parameter { bool softBounds; bool expHyperPrCal; bool isTipCals; + int nodeProposal; std::vector calibNodes; std::vector fossSpecimens; + double tuningVal; + int numMoves; + bool autotune; + }; diff --git a/src/Parameter_treescale.cpp b/src/Parameter_treescale.cpp index 4c7d788..16aa071 100644 --- a/src/Parameter_treescale.cpp +++ b/src/Parameter_treescale.cpp @@ -55,15 +55,16 @@ Treescale::Treescale(MbRandom *rp, Model *mp, double sv, double yb, double ob, i tOrigExpHPRate = 1.0 / (sv * 0.5); treeOriginTime = ranPtr->exponentialRv(tOrigExpHPRate) + sv; - + + if (treeTimePrior > 5) + tsPriorD = 3; + tuning = sv * 0.1; if(tsPriorD == 1){ exponHPCalib = false; if(yngBound > 0.0 && oldBound > 0.0){ isBounded = true; - // TAH: This might be stupid, but trying to get the window size to be - // at least 15% smaller than the calibration bounds if(yngBound != oldBound){ tuning = ((oldBound + yngBound) * 0.5) * 0.1; double windowSize = 2 * tuning; @@ -75,12 +76,6 @@ Treescale::Treescale(MbRandom *rp, Model *mp, double sv, double yb, double ob, i cout << "Calib window = " << calibSize << " Tuning window = " << windowSize << endl; } } - if(0){ //yngBound == 0.0 && oldBound >= 1000000.0){ - retune = true; - //cout << "retune = true" << endl; - //getchar(); - } - //else isBounded = false; } else if(tsPriorD == 2){ isBounded = true; @@ -91,19 +86,17 @@ Treescale::Treescale(MbRandom *rp, Model *mp, double sv, double yb, double ob, i expoRate = 1.0 / (yngBound * 200.0); exponHPCalib = false; } - oldBound = 1000000.0; + oldBound = 100000.0; tuning = ((yngBound + (sv * 1.2)) * 0.5) * 0.1; double windowSize = 2 * tuning; cout << "Exponential rate = " << expoRate << " E[TS] = " << (1 / expoRate) + yngBound << " offset = " << yngBound << " Tuning window = " << windowSize << endl; } - else if(tsPriorD > 2 && treeTimePrior > 5){ // This is a uniform prior on the age of the root + else if(tsPriorD > 2 && treeTimePrior > 5){ isBounded = true; tsPriorD = 1; - - oldBound = 100000.0; - tuning = ((yngBound + (sv * 2.0)) * 0.5) * 0.1; + tuning = ((yngBound + (sv * 1.2)) * 0.5) * 0.1; double windowSize = 2.0 * tuning; double calibSize = oldBound - yngBound; @@ -150,7 +143,7 @@ double Treescale::update(double &oldLnL) { lppr = updateTreeOrigTime(oldLnL); } else{ - lppr = updateTreeScale(oldLnL); + lppr = updateTreeScalePropSE(oldLnL); } return lppr; } @@ -159,14 +152,9 @@ double Treescale::updateTreeScale(double &oldLnL) { Tree *t = modelPtr->getActiveTree(); Node *rt = t->getRoot(); - - - - - // for birth-death prior + double oldtreeprob = getLnTreeProb(t); - // window if(retune && tuning > scaleVal) tuning = scaleVal * 0.5; @@ -191,7 +179,6 @@ double Treescale::updateTreeScale(double &oldLnL) { double oldRH, newRH; oldRH = scaleVal; - // get new value with reflection double u = ranPtr->uniformRv(-0.5,0.5) * (limO - limY); newRH = oldRH + u; while(newRH < lowBound || newRH > hiBound){ @@ -201,12 +188,10 @@ double Treescale::updateTreeScale(double &oldLnL) { newRH = (2 * hiBound) - newRH; } - // scale all node height proportional to root height double scaleRatio = oldRH / newRH; int numNodes = t->getNumNodes(); for(int i=0; igetNodeByIndex(i); - //if(p->getIsLeaf() == false && p != rt){ if(p != rt){ if(p->getIsLeaf() && p->getIsCalibratedDepth() == false){ p->setNodeDepth(0.0); @@ -216,7 +201,7 @@ double Treescale::updateTreeScale(double &oldLnL) { double newP = oldP * scaleRatio; p->setNodeDepth(newP); p->setFossAttchTime(0.0); - if(false){ //(treeTimePrior == 6 && p->getIsCalibratedDepth()){ + if(false){ double oldPhi = p->getFossAttchTime(); double fA = p->getNodeYngTime(); double newPhi = oldPhi * scaleRatio; @@ -233,7 +218,7 @@ double Treescale::updateTreeScale(double &oldLnL) { } scaleVal = newRH; t->setTreeScale(scaleVal); - t->treeScaleUpdateFossilAttchTimes(scaleRatio); + t->treeScaleUpdateFossilAttchTimes(scaleRatio, oldRH, newRH); t->treeUpdateNodeOldestBoundsAttchTimes(); t->setAllNodeBranchTimes(); @@ -245,16 +230,99 @@ double Treescale::updateTreeScale(double &oldLnL) { expoRate = t->getRootCalibExpRate(); lnPriorRatio += lnExponentialTSPriorRatio(newRH, oldRH); } - // TAH: it's uniform - double lnProposalRatio = 0.0; + double lnProposalRatio = 0.0; - // |J| is only for uniform prior on node times double jacobian = 0.0; if(treeTimePrior < 2) jacobian = (log(oldRH) - log(newRH)) * (t->getNumTaxa() - 2); + + + t->flipAllCls(); + t->flipAllTis(); + t->upDateAllCls(); + t->upDateAllTis(); + modelPtr->setTiProb(); + return lnPriorRatio + lnProposalRatio + jacobian; +} +double Treescale::updateTreeScalePropSE(double &oldLnL) { + + Tree *t = modelPtr->getActiveTree(); + Node *rt = t->getRoot(); + + double oldtreeprob = getLnTreeProb(t); + + double rtLB = t->getNodeLowerBoundTime(rt) * scaleVal; + double lowBound = rtLB; + + double hiBound = 100000.0; + if(treeTimePrior == 4){ + if(hiBound > treeOriginTime) + hiBound = treeOriginTime; + } + + if(isBounded){ + if(lowBound < yngBound) + lowBound = yngBound; + if(hiBound > oldBound) + hiBound = oldBound; + } + + double oldRH, newRH; + oldRH = scaleVal; + + double rv = ranPtr->uniformRv(); + double tv = log(2.0); + double c = tv * (rv - 0.5); + newRH = oldRH * exp(c); + bool validV = false; + do{ + if(newRH < lowBound) + newRH = lowBound * lowBound / newRH; + else if(newRH > hiBound) + newRH = hiBound * hiBound / newRH; + else + validV = true; + } while(!validV); + + double scaleRatio = oldRH / newRH; + int numNodes = t->getNumNodes(); + for(int i=0; igetNodeByIndex(i); + if(p != rt){ + if(p->getIsLeaf() && p->getIsCalibratedDepth() == false){ + p->setNodeDepth(0.0); + } + else{ + double oldP = p->getNodeDepth(); + double newP = oldP * scaleRatio; + p->setNodeDepth(newP); + p->setFossAttchTime(0.0); + } + } + } + scaleVal = newRH; + t->setTreeScale(scaleVal); + t->treeScaleUpdateFossilAttchTimes(scaleRatio, oldRH, newRH); + t->treeUpdateNodeOldestBoundsAttchTimes(); + t->setAllNodeBranchTimes(); + + double lnPriorRatio = 0.0; + double newtreeprob = getLnTreeProb(t); + lnPriorRatio += (newtreeprob - oldtreeprob); + if(tsPriorD == 2){ + if(exponHPCalib) + expoRate = t->getRootCalibExpRate(); + lnPriorRatio += lnExponentialTSPriorRatio(newRH, oldRH); + } + double lnProposalRatio = c; + + double jacobian = 0.0; + if(treeTimePrior < 2) + jacobian = (log(oldRH) - log(newRH)) * (t->getNumTaxa() - 2); + t->flipAllCls(); t->flipAllTis(); t->upDateAllCls(); @@ -265,11 +333,11 @@ double Treescale::updateTreeScale(double &oldLnL) { return lnPriorRatio + lnProposalRatio + jacobian; } + double Treescale::updateTreeOrigTime(double &oldLnL) { Tree *t = modelPtr->getActiveTree(); - // for birth-death prior double oldtreeprob = getLnTreeProb(t); @@ -278,12 +346,9 @@ double Treescale::updateTreeOrigTime(double &oldLnL) { double lowBound = scaleVal; double hiBound = limO; - - double oldTO, newTO; oldTO = treeOriginTime; - // get new value with reflection double u = ranPtr->uniformRv(-0.5,0.5) * (limO - limY); newTO = oldTO + u; while(newTO < lowBound || newTO > hiBound){ @@ -292,16 +357,12 @@ double Treescale::updateTreeOrigTime(double &oldLnL) { if(newTO > hiBound) newTO = (2 * hiBound) - newTO; } - - // scale all node height proportional to root height - + treeOriginTime = newTO; double lnPriorRatio = 0.0; double newtreeprob = getLnTreeProb(t); - //cout << lnPriorRatio << endl; lnPriorRatio += (newtreeprob - oldtreeprob); - lnPriorRatio += lnExponentialTreeOrigPriorRatio(newTO, oldTO); double lnProposalRatio = 0.0; @@ -363,7 +424,6 @@ double Treescale::getLnTreeProb(Tree *t) { return nps; } else if(treeTimePrior == 6){ - // calibrated birth death Speciation *s = modelPtr->getActiveSpeciation(); s->setAllBDFossParams(); double nps = t->getTreeCalBDSSTreeNodePriorProb(s->getBDSSSpeciationRateLambda(), @@ -373,7 +433,6 @@ double Treescale::getLnTreeProb(Tree *t) { return nps; } else if(treeTimePrior == 7){ - // calibrated birth death Speciation *s = modelPtr->getActiveSpeciation(); s->setAllBDFossParams(); double nps = t->getTreeAncCalBDSSTreeNodePriorProb(s->getBDSSSpeciationRateLambda(), diff --git a/src/Parameter_treescale.h b/src/Parameter_treescale.h index e6c0416..ee5f644 100644 --- a/src/Parameter_treescale.h +++ b/src/Parameter_treescale.h @@ -67,9 +67,15 @@ class Treescale : public Parameter { bool exponHPCalib; double treeOriginTime; double tOrigExpHPRate; + double numAccepted; + double numTried; double updateTreeScale(double &oldLnL); + double updateTreeScalePropSE(double &oldLnL); + double updateTreeOrigTime(double &oldLnL); + + }; #endif \ No newline at end of file diff --git a/src/dppdiv.cpp b/src/dppdiv.cpp index 068279d..1d4d933 100644 --- a/src/dppdiv.cpp +++ b/src/dppdiv.cpp @@ -125,8 +125,9 @@ int main (int argc, char * const argv[]) { double rateSc = 4.0; // scale param for gamma dist on rates double hyperSh = -1.0; // scale hyperparam for gamma dist on concentration param double hyperSc = -1.0; // scale hyperparam for gamma dist on concentration param - double netDiv = 1.0; // initial diversificaton rate (lambda - mu) - double relDeath = 0.5; // initial relative death rate (mu / lambda) + double netDiv = -1.0; // initial diversificaton rate (lambda - mu) + double relDeath = -1.0; // initial relative death rate (mu / lambda) + double ssbdPrS = -1.0; double fixclokrt = -1.0; // fix the clock rate to this int offmove = 0; // used to turn off one particular move int printFreq = 100; @@ -153,6 +154,8 @@ int main (int argc, char * const argv[]) { int modelType = 1; // 1 = DPP, 2 = strict clock, 3 = uncorrelated-gamma bool fixClockToR = false; bool indHP = false; + bool doAbsRts = false; + bool fixTest = false; if(argc > 1){ for (int i = 1; i < argc; i++){ @@ -210,14 +213,18 @@ int main (int argc, char * const argv[]) { tipDateFN = argv[i+1]; else if(!strcmp(curArg, "-npr")) treeNodePrior = atoi(argv[i+1]); - else if(!strcmp(curArg, "-tgs")) + else if(!strcmp(curArg, "-tgs")) // set treeNodePrior to do calibrated birth-death treeNodePrior = 6; - else if(!strcmp(curArg, "-tga")) + else if(!strcmp(curArg, "-tga")){ // set treeNodePrior to do calibrated birth-death with ancestor fossils treeNodePrior = 7; + doAbsRts = true; + } else if(!strcmp(curArg, "-bdr")) // (lambda - mu) netDiv = atof(argv[i+1]); else if(!strcmp(curArg, "-bda")) // (mu / lambda) relDeath = atof(argv[i+1]); + else if(!strcmp(curArg, "-bds")) // (mu / lambda) + ssbdPrS = atof(argv[i+1]); else if(!strcmp(curArg, "-fix")){ // fix clock fixclokrt = atof(argv[i+1]); fixClockToR = true; @@ -249,6 +256,13 @@ int main (int argc, char * const argv[]) { fixModelPs = true; else if(!strcmp(curArg, "-ihp")) indHP = true; + else if(!strcmp(curArg, "-abs")){ + doAbsRts = true; + rateSh = 1.0; + } + else if(!strcmp(curArg, "-fxtr")){ + fixTest = true; + } else if(!strcmp(curArg, "-h")){ printHelp(false); return 0; @@ -297,10 +311,13 @@ int main (int argc, char * const argv[]) { MbRandom myRandom; myRandom.setSeed(s1, s2); + // set up the model Model myModel(&myRandom, &myAlignment, treeStr, priorMean, rateSh, rateSc, hyperSh, hyperSc, userBLs, moveAllN, offmove, rndNdMv, calibFN, - treeNodePrior, netDiv, relDeath, fixclokrt, rootfix, softbnd, calibHyP, - dpmExpHyp, dpmEHPPrM, gammaExpHP, modelType, fixModelPs, indHP, tipDateFN); + treeNodePrior, netDiv, relDeath, ssbdPrS, fixclokrt, rootfix, softbnd, calibHyP, + dpmExpHyp, dpmEHPPrM, gammaExpHP, modelType, fixModelPs, indHP, tipDateFN, fixTest); + if(doAbsRts) + myModel.setEstAbsRates(true); if(runPrior) myModel.setRunUnderPrior(true); if(justTree){ @@ -312,4 +329,3 @@ int main (int argc, char * const argv[]) { return 0; } -