use of dr.evomodel.epidemiology.casetocase.AbstractCase in project beast-mcmc by beast-dev.
the class TransmissionTreeOperator method doOperation.
public double doOperation() {
TreeModel tree = c2cLikelihood.getTreeModel();
BranchMapModel branchMap = c2cLikelihood.getBranchMap();
AbstractCase[] newBranchMap = branchMap.getArrayCopy();
int[] oldParents = getParentsArray(tree);
double[] oldHeights = getHeightsArray(tree);
double hr = innerOperator.doOperation();
int[] newParents = getParentsArray(tree);
ArrayList<Integer> changedNodes = new ArrayList<Integer>();
for (int i = 0; i < tree.getNodeCount(); i++) {
if (oldParents[i] != newParents[i]) {
changedNodes.add(i);
}
}
if (changedNodes.size() != 0) {
if (innerOperator instanceof ExchangeOperator) {
//this is a node swap operator
AbstractCase[] nodePaintings = new AbstractCase[2];
AbstractCase[] parentPaintings = new AbstractCase[2];
for (int i = 0; i < 2; i++) {
nodePaintings[i] = branchMap.get(changedNodes.get(i));
parentPaintings[i] = branchMap.get(oldParents[changedNodes.get(i)]);
}
if (nodePaintings[0] == parentPaintings[0] || nodePaintings[1] == parentPaintings[1]) {
//If this is not true there is nothing to do - the result is already a valid painting
for (int i = 0; i < 2; i++) {
paintUp(tree, nodePaintings[1 - i], nodePaintings[i], branchMap, newBranchMap, tree.getNode(changedNodes.get(i)).getNumber(), newParents);
}
}
} else if (innerOperator instanceof SubtreeSlideOperator || innerOperator instanceof WilsonBalding) {
//this is a node transplantation operator
int movedNode = -1;
int oldChild = -1;
int newChild = -1;
for (int i = 0; i < tree.getNodeCount(); i++) {
// exception later, though.
if (tree.getNodeHeight(tree.getNode(i)) != oldHeights[i]) {
movedNode = i;
break;
}
}
for (int j : changedNodes) {
if (j != movedNode) {
if (tree.getParent(tree.getNode(j)) == tree.getNode(movedNode)) {
newChild = j;
} else {
oldChild = j;
}
}
}
if (movedNode == -1 || oldChild == -1 || newChild == -1) {
// is this a bug or should the move be rejected (i.e., return -Inf HR)?
throw new RuntimeException("Failed to establish relationship between relocated node and others");
}
NodeRef movedNodeObject = tree.getNode(movedNode);
NodeRef oldChildObject = tree.getNode(oldChild);
NodeRef newChildObject = tree.getNode(newChild);
int otherChild = -1;
NodeRef otherChildObject;
//Find the other child of the moved node (the root of the transplanted subtree)
for (int i = 0; i < tree.getChildCount(movedNodeObject); i++) {
if (tree.getChild(movedNodeObject, i) != newChildObject) {
otherChildObject = tree.getChild(movedNodeObject, i);
otherChild = otherChildObject.getNumber();
}
}
//If the child of the moved node is the earliest node with its painting:
if (branchMap.get(otherChild) != branchMap.get(movedNode)) {
newBranchMap[movedNode] = branchMap.get(newChild);
} else {
//Change all paintings up the tree from the old child that used to match the moved node to match
//the old child
paintUp(tree, branchMap.get(movedNode), branchMap.get(oldChild), branchMap, newBranchMap, oldChild, oldParents);
//This may have resulted in the moved node being recoloured wrong.
newBranchMap[movedNode] = branchMap.get(movedNode);
branchMap.setAll(newBranchMap, true);
//Change all paintings up the tree from the moved node that used to match the new child to match
//the moved node
paintUp(tree, branchMap.get(newChild), branchMap.get(movedNode), branchMap, newBranchMap, movedNode, newParents);
}
} else {
//I don't know what this is
throw new UnsupportedOperationException("Operator class " + innerOperator.getOperatorName() + " not yet " + "supported");
}
}
branchMap.setAll(newBranchMap, true);
c2cLikelihood.makeDirty();
return hr;
}
use of dr.evomodel.epidemiology.casetocase.AbstractCase in project beast-mcmc by beast-dev.
the class TransmissionExchangeOperatorA method exchange.
public double exchange() {
TreeModel tree = c2cLikelihood.getTreeModel();
final int nodeCount = tree.getNodeCount();
final NodeRef root = tree.getRoot();
NodeRef i = root;
while (root == i) {
i = tree.getNode(MathUtils.nextInt(nodeCount));
}
ArrayList<NodeRef> candidates = getPossibleExchanges(tree, i);
int candidateCount = candidates.size();
if (candidateCount == 0) {
// throw new OperatorFailedException("No valid exchanges for this node");
return Double.NEGATIVE_INFINITY;
}
NodeRef j = candidates.get(MathUtils.nextInt(candidates.size()));
int jFirstCandidateCount = getPossibleExchanges(tree, j).size();
double HRDenom = (1 / ((double) candidateCount)) + (1 / ((double) jFirstCandidateCount));
NodeRef iP = tree.getParent(i);
NodeRef jP = tree.getParent(j);
if (resampleInfectionTimes) {
BranchMapModel branchMap = c2cLikelihood.getBranchMap();
AbstractCase iCase = branchMap.get(i.getNumber());
AbstractCase jCase = branchMap.get(j.getNumber());
AbstractCase parentCase = branchMap.get(iP.getNumber());
if (iCase != parentCase) {
iCase.setInfectionBranchPosition(MathUtils.nextDouble());
}
if (jCase != parentCase) {
jCase.setInfectionBranchPosition(MathUtils.nextDouble());
}
}
exchangeNodes(tree, i, j, iP, jP);
ArrayList<NodeRef> reverseCandidatesIfirst = getPossibleExchanges(tree, i);
ArrayList<NodeRef> reverseCandidatesJfirst = getPossibleExchanges(tree, j);
double HRNum = (1 / (double) reverseCandidatesIfirst.size()) + (1 / (double) reverseCandidatesJfirst.size());
return Math.log(HRNum / HRDenom);
}
Aggregations