[beast-mcmc] r6729 committed - Fixed a bug where TT and PT timings were not the same (to 11 decimal p...

[beast...@googlecode.com]

Revision: 6729
Author: mdha...@gmail.com
Date: Thu Apr 9 12:41:26 2015 UTC
Log: Fixed a bug where TT and PT timings were not the same (to 11
decimal places)
https://code.google.com/p/beast-mcmc/source/detail?r=6729

Modified:

/trunk/src/dr/evomodel/epidemiology/casetocase/CaseToCaseTreeLikelihood.java
/trunk/src/dr/evomodel/epidemiology/casetocase/WithinCaseCoalescent.java

=======================================
---
/trunk/src/dr/evomodel/epidemiology/casetocase/CaseToCaseTreeLikelihood.java
Tue Apr 7 10:03:54 2015 UTC
+++
/trunk/src/dr/evomodel/epidemiology/casetocase/CaseToCaseTreeLikelihood.java
Thu Apr 9 12:41:26 2015 UTC
@@ -784,7 +784,7 @@

private double getInfectionTime(double min, double max, AbstractCase
infected){
final double branchLength = max-min;
- return max -
branchLength*infected.getInfectionBranchPosition().getParameterValue(0);
+ return min +
branchLength*(1-infected.getInfectionBranchPosition().getParameterValue(0));
}


@@ -1375,40 +1375,43 @@
}

for(AbstractCase aCase : outbreak.getCases()){
- NodeRef originalNode = getEarliestNodeInPartition(aCase);
- int infectionNodeNo = originalNode.getNumber();
- if(!treeModel.isRoot(originalNode)){
- NodeRef originalParent = treeModel.getParent(originalNode);
- double nodeTime = getNodeTime(originalNode);
- double infectionTime = getInfectionTime(aCase);
- double heightToBreakBranch = getHeight(originalNode) +
(nodeTime - infectionTime);
- FlexibleNode newNode =
(FlexibleNode)outTree.getNode(infectionNodeNo);
- FlexibleNode oldParent =
(FlexibleNode)outTree.getParent(newNode);
+ if(aCase.wasEverInfected()) {
+ NodeRef originalNode = getEarliestNodeInPartition(aCase);
+
+ int infectionNodeNo = originalNode.getNumber();
+ if (!treeModel.isRoot(originalNode)) {
+ NodeRef originalParent =
treeModel.getParent(originalNode);
+ double nodeTime = getNodeTime(originalNode);
+ double infectionTime = getInfectionTime(aCase);
+ double heightToBreakBranch = getHeight(originalNode) +
(nodeTime - infectionTime);
+ FlexibleNode newNode = (FlexibleNode)
outTree.getNode(infectionNodeNo);
+ FlexibleNode oldParent = (FlexibleNode)
outTree.getParent(newNode);

- outTree.beginTreeEdit();
- outTree.removeChild(oldParent, newNode);
- FlexibleNode infectionNode = new FlexibleNode();
- infectionNode.setHeight(heightToBreakBranch);
- infectionNode.setLength(oldParent.getHeight() -
heightToBreakBranch);
- infectionNode.setAttribute(PARTITIONS_KEY,
getNodePartition(treeModel, originalParent));
- newNode.setLength(nodeTime - infectionTime);
+ outTree.beginTreeEdit();
+ outTree.removeChild(oldParent, newNode);
+ FlexibleNode infectionNode = new FlexibleNode();
+ infectionNode.setHeight(heightToBreakBranch);
+ infectionNode.setLength(oldParent.getHeight() -
heightToBreakBranch);
+ infectionNode.setAttribute(PARTITIONS_KEY,
getNodePartition(treeModel, originalParent));
+ newNode.setLength(nodeTime - infectionTime);

- outTree.addChild(oldParent, infectionNode);
- outTree.addChild(infectionNode, newNode);
- outTree.endTreeEdit();
- } else {
- double nodeTime = getNodeTime(originalNode);
- double infectionTime = getInfectionTime(aCase);
- double heightToInstallRoot = getHeight(originalNode) +
(nodeTime - infectionTime);
- FlexibleNode newNode =
(FlexibleNode)outTree.getNode(infectionNodeNo);
- outTree.beginTreeEdit();
- FlexibleNode infectionNode = new FlexibleNode();
- infectionNode.setHeight(heightToInstallRoot);
- infectionNode.setAttribute(PARTITIONS_KEY, "The_Ether");
- outTree.addChild(infectionNode, newNode);
- newNode.setLength(heightToInstallRoot -
getHeight(originalNode));
- outTree.setRoot(infectionNode);
- outTree.endTreeEdit();
+ outTree.addChild(oldParent, infectionNode);
+ outTree.addChild(infectionNode, newNode);
+ outTree.endTreeEdit();
+ } else {
+ double nodeTime = getNodeTime(originalNode);
+ double infectionTime = getInfectionTime(aCase);
+ double heightToInstallRoot = getHeight(originalNode) +
(nodeTime - infectionTime);
+ FlexibleNode newNode = (FlexibleNode)
outTree.getNode(infectionNodeNo);
+ outTree.beginTreeEdit();
+ FlexibleNode infectionNode = new FlexibleNode();
+ infectionNode.setHeight(heightToInstallRoot);
+
infectionNode.setAttribute(PARTITIONS_KEY, "The_Ether");
+ outTree.addChild(infectionNode, newNode);
+ newNode.setLength(heightToInstallRoot -
getHeight(originalNode));
+ outTree.setRoot(infectionNode);
+ outTree.endTreeEdit();
+ }
}
}

=======================================
---
/trunk/src/dr/evomodel/epidemiology/casetocase/WithinCaseCoalescent.java
Thu Apr 2 15:11:48 2015 UTC
+++
/trunk/src/dr/evomodel/epidemiology/casetocase/WithinCaseCoalescent.java
Thu Apr 9 12:41:26 2015 UTC
@@ -85,7 +85,7 @@

if(DEBUG){

- super.debugOutputTree("bleh.nex", false);
+ super.debugOutputTree("bleh.nex", true);
}

double logL = 0;
@@ -155,18 +155,20 @@
if (recalculateCoalescentFlags[number]) {
Treelet treelet = partitionsAsTrees.get(aCase);

- if (children.size() != 0) {
- SpecifiedZeroCoalescent coalescent = new
SpecifiedZeroCoalescent(treelet, demoModel,
- treelet.getZeroHeight(), mode ==
Mode.TRUNCATE);
- partitionTreeLogLikelihoods[number] =
coalescent.calculateLogLikelihood();
- coalescencesLogLikelihood +=
partitionTreeLogLikelihoods[number];
- if (DEBUG &&
partitionTreeLogLikelihoods[number] == Double.POSITIVE_INFINITY) {
- debugOutputTree("infCoalescent.nex",
false);
- debugTreelet(treelet, aCase
+ "_partition.nex");
- }
- } else {
- partitionTreeLogLikelihoods[number] = 0.0;
+
+
+ if (children.size() != 0) {
+ SpecifiedZeroCoalescent coalescent = new
SpecifiedZeroCoalescent(treelet, demoModel,
+ treelet.getZeroHeight(), mode ==
Mode.TRUNCATE);
+ partitionTreeLogLikelihoods[number] =
coalescent.calculateLogLikelihood();
+ coalescencesLogLikelihood +=
partitionTreeLogLikelihoods[number];
+ if (DEBUG && partitionTreeLogLikelihoods[number]
== Double.POSITIVE_INFINITY) {
+ debugOutputTree("infCoalescent.nex", false);
+ debugTreelet(treelet, aCase
+ "_partition.nex");
}
+ } else {
+ partitionTreeLogLikelihoods[number] = 0.0;
+ }
recalculateCoalescentFlags[number] = false;
} else {
coalescencesLogLikelihood +=
partitionTreeLogLikelihoods[number];
@@ -293,33 +295,46 @@
AbstractCase aCase = outbreak.getCase(i);
if(aCase.wasEverInfected() &&
partitionsAsTrees.get(aCase)==null){

- if(aCase.wasEverInfected()) {
+ NodeRef partitionRoot = getEarliestNodeInPartition(aCase);

- NodeRef partitionRoot =
getEarliestNodeInPartition(aCase);

- double infectionTime =
getInfectionTime(branchMap.get(partitionRoot.getNumber()));
- double rootTime = getNodeTime(partitionRoot);
+ double infectionTime =
getInfectionTime(branchMap.get(partitionRoot.getNumber()));
+ double rootTime = getNodeTime(partitionRoot);

- FlexibleNode newRoot = new FlexibleNode();
+ FlexibleNode newRoot = new FlexibleNode();

- FlexibleTree littleTree = new FlexibleTree(newRoot);
- littleTree.beginTreeEdit();
+ FlexibleTree littleTree = new FlexibleTree(newRoot);
+ littleTree.beginTreeEdit();

- if (!treeModel.isExternal(partitionRoot)) {
- for (int j = 0; j <
treeModel.getChildCount(partitionRoot); j++) {
- copyPartitionToTreelet(littleTree,
treeModel.getChild(partitionRoot, j), newRoot, aCase);
- }
+ if (!treeModel.isExternal(partitionRoot)) {
+ for (int j = 0; j <
treeModel.getChildCount(partitionRoot); j++) {
+ copyPartitionToTreelet(littleTree,
treeModel.getChild(partitionRoot, j), newRoot, aCase);
}
+ }

- littleTree.endTreeEdit();
+ littleTree.endTreeEdit();

- littleTree.resolveTree();
+ littleTree.resolveTree();

+ double sampleTipHeight = 0;

- partitionsAsTrees.put(aCase, new Treelet(littleTree,
- littleTree.getRootHeight() + rootTime -
infectionTime));
+ if(littleTree.getExternalNodeCount()>1) {
+ for (int j = 0; j < littleTree.getExternalNodeCount();
j++) {
+ NodeRef node = littleTree.getExternalNode(j);
+ if
(!littleTree.getNodeTaxon(node).getId().startsWith("Transmission_")) {
+ sampleTipHeight =
littleTree.getNodeHeight(node);
+ break;
+ }

+ }
}
+
+ Treelet treelet = new Treelet(littleTree,
+ sampleTipHeight +
(aCase.examTime-getInfectionTime(aCase)));
+
+ partitionsAsTrees.put(aCase, treelet);
+
+
}
}
}
@@ -446,6 +461,7 @@
}

public double calculateLogLikelihood() {
+
return calculatePartitionTreeLogLikelihood(getIntervals(),
getDemographicFunction(), 0, zeroHeight,
truncate);
}
@@ -458,8 +474,6 @@

double logL = 0.0;

-
-
double startTime = -zeroHeight;
final int n = intervals.getIntervalCount();

@@ -533,6 +547,7 @@
startTime = finishTime;
} else {
if(!(demographicFunction instanceof LinearGrowth)){
+
throw new RuntimeException("Function must have zero
population at t=0 if truncate=false");
}