[beast-mcmc] r6745 committed - Made SubtreeJump a 'tunable' operator - destination is picked as the i...

[beast...@googlecode.com]

Revision: 6745
Author: rambaut
Date: Thu Apr 23 15:57:29 2015 UTC
Log: Made SubtreeJump a 'tunable' operator - destination is picked as
the inverse of relatedness in the tree to the power of the tuning
parameter. Large tuning parameters tend to pick related branches, a tuning
parameter of zero picks uniformally.
https://code.google.com/p/beast-mcmc/source/detail?r=6745

Modified:
/trunk/build_tracer.xml
/trunk/src/dr/evomodel/operators/SubtreeJumpOperator.java

=======================================
--- /trunk/build_tracer.xml Tue Jun 24 12:30:44 2014 UTC
+++ /trunk/build_tracer.xml Thu Apr 23 15:57:29 2015 UTC
@@ -129,8 +129,8 @@

</target>

- <property name="version" value="1.6.1pre" />
- <property name="version_number" value="1.6.1" />
+ <property name="version" value="1.6.0" />
+ <property name="version_number" value="1.6.0" />
<property name="release_dir" value="release_tracer" />
<property name="name" value="Tracer" />

@@ -216,11 +216,12 @@
name="${name} v${version}"
mainclass="dr.app.tracer.application.TracerApp"
icon="${Mac_dir}/icons/Tracer.icns"
- jvmversion="1.5+"
+ jvmversion="1.6+"
+ stubfile="${Mac_dir}/universalJavaApplicationStub"
vmoptions="-Xmx1024M"
arguments=""
version="${version}"
- infostring="${name} v${version}, Copyright
2003-2013, Andrew Rambaut, Marc Suchard &amp; Alexei J. Drummond"
+ infostring="${name} v${version}, Copyright
2003-2015, Andrew Rambaut, Marc Suchard &amp; Alexei J. Drummond"
bundleid="tracer" >
<javaproperty name="apple.laf.useScreenMenuBar" value="true"/>
<jarfileset dir="${dist}">
=======================================
--- /trunk/src/dr/evomodel/operators/SubtreeJumpOperator.java Mon Mar 2
17:52:54 2015 UTC
+++ /trunk/src/dr/evomodel/operators/SubtreeJumpOperator.java Thu Apr 23
15:57:29 2015 UTC
@@ -30,10 +30,7 @@
import dr.evomodel.tree.TreeModel;
import dr.evomodelxml.operators.SubtreeJumpOperatorParser;
import dr.evomodelxml.operators.SubtreeSlideOperatorParser;
-import dr.inference.operators.CoercableMCMCOperator;
-import dr.inference.operators.CoercionMode;
-import dr.inference.operators.OperatorFailedException;
-import dr.inference.operators.OperatorUtils;
+import dr.inference.operators.*;
import dr.math.MathUtils;

import java.util.ArrayList;
@@ -45,14 +42,24 @@
* @author Andrew Rambaut
* @version $Id$
*/
-public class SubtreeJumpOperator extends AbstractTreeOperator /*
implements CoercableMCMCOperator */ { // not coercable at the moment.
+public class SubtreeJumpOperator extends AbstractTreeOperator implements
CoercableMCMCOperator {

private TreeModel tree = null;
-// private CoercionMode mode = CoercionMode.DEFAULT;
+ private double size = 1;
+ private CoercionMode mode = CoercionMode.DEFAULT;

- public SubtreeJumpOperator(TreeModel tree, double weight) {
+ /**
+ * Constructor
+ * @param tree
+ * @param weight
+ * @param size a non-negative value used as a power coeficient for the
relatedness weights
+ * @param mode
+ */
+ public SubtreeJumpOperator(TreeModel tree, double weight, double size,
CoercionMode mode) {
this.tree = tree;
setWeight(weight);
+ this.size = size;
+ this.mode = mode;
}

/**
@@ -85,16 +92,18 @@
// get a list of all edges that intersect this height
final List<NodeRef> destinations = getIntersectingEdges(tree,
height);

- // remove the target node and its sibling (shouldn't be there
because their parent's height is exactly equal to the target height).
- destinations.remove(i);
- destinations.remove(CiP);
-
if (destinations.size() < 1) {
throw new OperatorFailedException("no destinations to jump
to");
}
+
+ double[] weights = getDestinationWeights(tree, iP, height,
destinations);
+
+ // remove the target node and its sibling (shouldn't be there
because their parent's height is exactly equal to the target height).
+ destinations.remove(i);
+ destinations.remove(CiP);

// pick uniformly from this list
- final NodeRef j =
destinations.get(MathUtils.nextInt(destinations.size()));
+ final NodeRef j = destinations.get(pickDestination(weights, size));
final NodeRef jP = tree.getParent(j);

tree.beginTreeEdit();
@@ -119,6 +128,20 @@

return logq;
}
+
+ private int pickDestination(double[] weights, double size) {
+ double sum = 0.0;
+ for (int i = 0; i < weights.length; i++) {
+ weights[i] = 1.0 / Math.pow(weights[i], size);
+ sum += weights[i];
+ i++;
+ }
+ for (int i = 0; i < weights.length; i++) {
+ weights[i] /= sum;
+ }
+
+ return MathUtils.randomChoicePDF(weights);
+ }

private List<NodeRef> getIntersectingEdges(Tree tree, double height) {

@@ -131,10 +154,49 @@
intersectingEdges.add(node);
}
}
+ return intersectingEdges;
+ }

- return intersectingEdges;
+ private double[] getDestinationWeights(Tree tree, NodeRef node0,
double height, List<NodeRef> intersectingEdges) {
+ double[] weights = new double[intersectingEdges.size()];
+ double sum = 0.0;
+ for (int i = 0; i < weights.length; i++) {
+ NodeRef node1 = intersectingEdges.get(i);
+ weights[i] =
tree.getNodeHeight(Tree.Utils.getCommonAncestor(tree, node0, node1)) -
height;
+ sum += weights[i];
+ i++;
+ }
+ for (int i = 0; i < weights.length; i++) {
+ weights[i] /= sum;
+ }
+
+ return weights;
+ }
+
+
+ public double getSize() {
+ return size;
+ }
+
+ public void setSize(double size) {
+ this.size = size;
+ }
+
+ public double getCoercableParameter() {
+ return Math.log(getSize());
+ }
+
+ public void setCoercableParameter(double value) {
+ setSize(Math.exp(value));
+ }
+
+ public double getRawParameter() {
+ return getSize();
}

+ public CoercionMode getMode() {
+ return mode;
+ }

public double getTargetAcceptanceProbability() {
return 0.234;
@@ -142,19 +204,18 @@


public String getPerformanceSuggestion() {
- return "";
+ double prob = MCMCOperator.Utils.getAcceptanceProbability(this);
+ double targetProb = getTargetAcceptanceProbability();
+
+ double ws = OperatorUtils.optimizeWindowSize(getSize(),
Double.MAX_VALUE, prob, targetProb);

-// double prob = Utils.getAcceptanceProbability(this);
-// double targetProb = getTargetAcceptanceProbability();
-//
-// double ws = OperatorUtils.optimizeWindowSize(getSize(),
Double.MAX_VALUE, prob, targetProb);
-//
-// if (prob < getMinimumGoodAcceptanceLevel()) {
-// return "Try decreasing size to about " + ws;
-// } else if (prob > getMaximumGoodAcceptanceLevel()) {
-// return "Try increasing size to about " + ws;
-// } else return "";
+ if (prob < getMinimumGoodAcceptanceLevel()) {
+ return "Try decreasing size to about " + ws;
+ } else if (prob > getMaximumGoodAcceptanceLevel()) {
+ return "Try increasing size to about " + ws;
+ } else return "";
}
+

public String getOperatorName() {
return SubtreeJumpOperatorParser.SUBTREE_JUMP + "(" + tree.getId()
+ ")";