package dr.evomodel.operators;
import dr.evomodel.treelikelihood.AbstractLikelihoodCore;
import dr.evomodel.treelikelihood.TreeLikelihood;
import dr.inference.operators.MCMCOperator;
import dr.inference.operators.OperatorFailedException;
import dr.inference.operators.SimpleMCMCOperator;
import dr.math.MathUtils;
import dr.xml.*;
/**
* @author Marc A. Suchard
*/
public class TipStateSwapOperator extends SimpleMCMCOperator {
public static final String TIP_STATE_OPERATOR = "tipStateSwapOperator";
public TipStateSwapOperator(TreeLikelihood treeLikelihood, double weight) {
this.treeLikelihood = treeLikelihood;
setWeight(weight);
int patternCount = treeLikelihood.getPatternCount();
states1 = new int[patternCount];
states2 = new int[patternCount];
}
private final int[] states1;
private final int[] states2;
private int index1;
private int index2;
public double doOperation() throws OperatorFailedException {
int tipCount = treeLikelihood.getTreeModel().getExternalNodeCount();
// Choose two tips to swap
index1 = MathUtils.nextInt(tipCount);
index2 = index1;
while (index2 == index1)
index2 = MathUtils.nextInt(tipCount);
swap(index1, index2);
treeLikelihood.makeDirty();
return 0;
}
private void swap(int i, int j) {
AbstractLikelihoodCore likelihoodCore = (AbstractLikelihoodCore) treeLikelihood.getLikelihoodCore();
likelihoodCore.getNodeStates(i, states1);
likelihoodCore.getNodeStates(j, states2);
likelihoodCore.setNodeStates(j, states1);
likelihoodCore.setNodeStates(i, states2);
}
public void reject() {
super.reject();
// There is currently no restore functions for tip states, so manually adjust state
swap(index1, index2);
}
public String getPerformanceSuggestion() {
if (MCMCOperator.Utils.getAcceptanceProbability(this) < getMinimumAcceptanceLevel()) {
return "";
} else if (MCMCOperator.Utils.getAcceptanceProbability(this) > getMaximumAcceptanceLevel()) {
return "";
} else {
return "";
}
}
public String getOperatorName() {
return TIP_STATE_OPERATOR;
}
public static XMLObjectParser PARSER = new AbstractXMLObjectParser() {
public String getParserName() {
return TIP_STATE_OPERATOR;
}
public Object parseXMLObject(XMLObject xo) throws XMLParseException {
TreeLikelihood treeLikelihood =
(TreeLikelihood) xo.getChild(TreeLikelihood.class);
final double weight = xo.getDoubleAttribute("weight");
return new TipStateSwapOperator(treeLikelihood, weight);
}
//************************************************************************
// AbstractXMLObjectParser implementation
//************************************************************************
public String getParserDescription() {
return "This element represents an operator to swap tip states between two random tips.";
}
public Class getReturnType() {
return ExchangeOperator.class;
}
public XMLSyntaxRule[] getSyntaxRules() {
return rules;
}
private final XMLSyntaxRule[] rules = {
AttributeRule.newDoubleRule("weight"),
new ElementRule(TreeLikelihood.class),
};
};
private final TreeLikelihood treeLikelihood;
}