/*
* File SubtreeSlide.java
*
* Copyright (C) 2010 Remco Bouckaert remco@cs.auckland.ac.nz
*
* This file is part of BEAST2.
* See the NOTICE file distributed with this work for additional
* information regarding copyright ownership and licensing.
*
* BEAST is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as
* published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
*
* BEAST is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with BEAST; if not, write to the
* Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
* Boston, MA 02110-1301 USA
*/
/*
* SubtreeSlideOperator.java
*
* Copyright (C) 2002-2006 Alexei Drummond and Andrew Rambaut
*
* This file is part of BEAST.
* See the NOTICE file distributed with this work for additional
* information regarding copyright ownership and licensing.
*
* BEAST is free software; you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as
* published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
*
* BEAST is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with BEAST; if not, write to the
* Free Software Foundation, Inc., 51 Franklin St, Fifth Floor,
* Boston, MA 02110-1301 USA
*/
package beast.evolution.operators;
import java.text.DecimalFormat;
import java.util.ArrayList;
import java.util.List;
import beast.core.Description;
import beast.core.Input;
import beast.evolution.tree.Node;
import beast.evolution.tree.Tree;
import beast.util.Randomizer;
/**
* Implements the subtree slide move.
*/
@Description("Moves the height of an internal node along the branch. " +
"If it moves up, it can exceed the root and become a new root. " +
"If it moves down, it may need to make a choice which branch to " +
"slide down into.")
public class SubtreeSlide extends TreeOperator {
final public Input<Double> sizeInput = new Input<>("size", "size of the slide, default 1.0", 1.0);
final public Input<Boolean> gaussianInput = new Input<>("gaussian", "Gaussian (=true=default) or uniform delta", true);
final public Input<Boolean> optimiseInput = new Input<>("optimise", "flag to indicate that the scale factor is automatically changed in order to achieve a good acceptance rate (default true)", true);
final public Input<Double> limitInput = new Input<>("limit", "limit on step size, default disable, " +
"i.e. -1. (when positive, gets multiplied by tree-height/log2(n-taxa).", -1.0);
// shadows size
double size;
private double limit;
@Override
public void initAndValidate() {
size = sizeInput.get();
limit = limitInput.get();
}
/**
* Do a probabilistic subtree slide move.
*
* @return log of Hastings Ratio, or Double.NEGATIVE_INFINITY if proposal should not be accepted *
*/
@Override
public double proposal() {
final Tree tree = treeInput.get(this);
double logq;
Node i;
final boolean markClades = markCladesInput.get();
// 1. choose a random node avoiding root
final int nodeCount = tree.getNodeCount();
do {
i = tree.getNode(Randomizer.nextInt(nodeCount));
} while (i.isRoot());
final Node p = i.getParent();
final Node CiP = getOtherChild(p, i);
final Node PiP = p.getParent();
// 2. choose a delta to move
final double delta = getDelta();
final double oldHeight = p.getHeight();
final double newHeight = oldHeight + delta;
// 3. if the move is up
if (delta > 0) {
// 3.1 if the topology will change
if (PiP != null && PiP.getHeight() < newHeight) {
// find new parent
Node newParent = PiP;
Node newChild = p;
while (newParent.getHeight() < newHeight) {
newChild = newParent;
if( markClades ) newParent.makeDirty(Tree.IS_FILTHY); // JH
newParent = newParent.getParent();
if (newParent == null) break;
}
// the moved node 'p' would become a child of 'newParent'
//
// 3.1.1 if creating a new root
if (newChild.isRoot()) {
replace(p, CiP, newChild);
replace(PiP, p, CiP);
p.setParent(null);
tree.setRoot(p);
}
// 3.1.2 no new root
else {
replace(p, CiP, newChild);
replace(PiP, p, CiP);
replace(newParent, newChild, p);
}
p.setHeight(newHeight);
// 3.1.3 count the hypothetical sources of this destination.
final int possibleSources = intersectingEdges(newChild, oldHeight, null);
//System.out.println("possible sources = " + possibleSources);
logq = -Math.log(possibleSources);
} else {
// just change the node height
p.setHeight(newHeight);
logq = 0.0;
}
}
// 4 if we are sliding the subtree down.
else {
// 4.0 is it a valid move?
if (i.getHeight() > newHeight) {
return Double.NEGATIVE_INFINITY;
}
// 4.1 will the move change the topology
if (CiP.getHeight() > newHeight) {
final List<Node> newChildren = new ArrayList<>();
final int possibleDestinations = intersectingEdges(CiP, newHeight, newChildren);
// if no valid destinations then return a failure
if (newChildren.size() == 0) {
return Double.NEGATIVE_INFINITY;
}
// pick a random parent/child destination edge uniformly from options
final int childIndex = Randomizer.nextInt(newChildren.size());
final Node newChild = newChildren.get(childIndex);
final Node newParent = newChild.getParent();
// 4.1.1 if p was root
if (p.isRoot()) {
// new root is CiP
replace(p, CiP, newChild);
replace(newParent, newChild, p);
CiP.setParent(null);
tree.setRoot(CiP);
} else {
replace(p, CiP, newChild);
replace(PiP, p, CiP);
replace(newParent, newChild, p);
}
p.setHeight(newHeight);
if( markClades ) {
// make dirty the path from the (down) moved node back up to former parent.
Node n = p;
while( n != CiP ) {
n.makeDirty(Tree.IS_FILTHY); // JH
n = n.getParent();
}
}
logq = Math.log(possibleDestinations);
} else {
p.setHeight(newHeight);
logq = 0.0;
}
}
return logq;
}
private double getDelta() {
if (!gaussianInput.get()) {
return (Randomizer.nextDouble() * size) - (size / 2.0);
} else {
return Randomizer.nextGaussian() * size;
}
}
private int intersectingEdges(Node node, double height, List<Node> directChildren) {
final Node parent = node.getParent();
if (parent.getHeight() < height) return 0;
if (node.getHeight() < height) {
if (directChildren != null) directChildren.add(node);
return 1;
}
if (node.isLeaf()) {
// TODO: verify that this makes sense
return 0;
} else {
final int count = intersectingEdges(node.getLeft(), height, directChildren) +
intersectingEdges(node.getRight(), height, directChildren);
return count;
}
}
/**
* automatic parameter tuning *
*/
@Override
public void optimize(final double logAlpha) {
if (optimiseInput.get()) {
double delta = calcDelta(logAlpha);
delta += Math.log(size);
final double f = Math.exp(delta);
// double f = Math.exp(delta);
if( limit > 0 ) {
final Tree tree = treeInput.get();
final double h = tree.getRoot().getHeight();
final double k = Math.log(tree.getLeafNodeCount()) / Math.log(2);
final double lim = (h / k) * limit;
if( f <= lim ) {
size = f;
}
} else {
size = f;
}
}
}
@Override
public double getCoercableParameterValue() {
return size;
}
@Override
public void setCoercableParameterValue(final double value) {
size = value;
}
@Override
public String getPerformanceSuggestion() {
final double prob = m_nNrAccepted / (m_nNrAccepted + m_nNrRejected + 0.0);
final double targetProb = getTargetAcceptanceProbability();
double ratio = prob / targetProb;
if (ratio > 2.0) ratio = 2.0;
if (ratio < 0.5) ratio = 0.5;
final double newDelta = size * ratio;
final DecimalFormat formatter = new DecimalFormat("#.###");
if (prob < 0.10) {
return "Try decreasing size to about " + formatter.format(newDelta);
} else if (prob > 0.40) {
return "Try increasing size to about " + formatter.format(newDelta);
} else return "";
}
}