/*
* RandomWalkOperator.java
*
* Copyright (c) 2002-2015 Alexei Drummond, Andrew Rambaut and Marc Suchard
*
* 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 dr.inference.operators;
import dr.inference.model.Bounds;
import dr.inference.model.Parameter;
import dr.inferencexml.operators.RandomWalkOperatorParser;
import dr.math.MathUtils;
import dr.util.Transform;
import java.util.ArrayList;
import java.util.List;
/**
* A generic random walk operator for use with a multi-dimensional parameters.
*
* @author Alexei Drummond
* @author Andrew Rambaut
* @version $Id: RandomWalkOperator.java,v 1.16 2005/06/14 10:40:34 rambaut Exp $
*/
public class RandomWalkOperator extends AbstractCoercableOperator {
public enum BoundaryCondition {
rejecting,
reflecting,
absorbing,
logit
}
public RandomWalkOperator(Parameter parameter, double windowSize, BoundaryCondition bc, double weight, CoercionMode mode) {
this(parameter, null, windowSize, bc, weight, mode);
}
public RandomWalkOperator(Parameter parameter, Parameter updateIndex, double windowSize, BoundaryCondition bc,
double weight, CoercionMode mode) {
this(parameter, updateIndex, windowSize, bc, weight, mode, null, null);
}
public RandomWalkOperator(Parameter parameter, Parameter updateIndex, double windowSize, BoundaryCondition boundaryCondition,
double weight, CoercionMode mode, Double lowerOperatorBound, Double upperOperatorBound) {
super(mode);
this.parameter = parameter;
this.windowSize = windowSize;
this.boundaryCondition = boundaryCondition;
setWeight(weight);
if (updateIndex != null) {
updateMap = new ArrayList<Integer>();
for (int i = 0; i < updateIndex.getDimension(); i++) {
if (updateIndex.getParameterValue(i) == 1.0)
updateMap.add(i);
}
updateMapSize=updateMap.size();
}
this.lowerOperatorBound = lowerOperatorBound;
this.upperOperatorBound = upperOperatorBound;
}
/**
* @return the parameter this operator acts on.
*/
public Parameter getParameter() {
return parameter;
}
public final double getWindowSize() {
return windowSize;
}
/**
* change the parameter and return the hastings ratio.
*/
public final double doOperation() {
// a random dimension to perturb
if (parameter.getDimension() <= 0) {
throw new RuntimeException("Illegal Dimension");
}
int dim;
if (updateMap == null) {
dim = MathUtils.nextInt(parameter.getDimension());
} else {
dim = updateMap.get(MathUtils.nextInt(updateMapSize));
}
// a random point around old value within windowSize * 2
double draw = (2.0 * MathUtils.nextDouble() - 1.0) * windowSize;
double oldValue = parameter.getParameterValue(dim);
final Bounds<Double> bounds = parameter.getBounds();
final double lower = (lowerOperatorBound == null ? bounds.getLowerLimit(dim) : Math.max(bounds.getLowerLimit(dim), lowerOperatorBound));
final double upper = (upperOperatorBound == null ? bounds.getUpperLimit(dim) : Math.min(bounds.getUpperLimit(dim), upperOperatorBound));
if (boundaryCondition == BoundaryCondition.logit) {
// scale oldValue to [0,1]
double x1 = (oldValue - lower) / (upper - lower);
// logit transform it, add the draw, inverse transform it
double x2 = Transform.LOGIT.inverse(Transform.LOGIT.transform(x1) + draw);
// parameter takes new value scaled back into interval [lower, upper]
parameter.setParameterValue(dim, (x2 * (upper - lower)) + lower);
// HR is the ratio of Jacobians for the before and after values in interval [0,1]
return Transform.LOGIT.getLogJacobian(x1) - Transform.LOGIT.getLogJacobian(x2);
} else {
double newValue = oldValue + draw;
if (boundaryCondition == BoundaryCondition.reflecting) {
newValue = reflectValue(newValue, lower, upper);
} else if (boundaryCondition == BoundaryCondition.absorbing && (newValue < lower || newValue > upper)) {
return 0.0;
} else if (boundaryCondition == BoundaryCondition.rejecting && (newValue < lower || newValue > upper)) {
return Double.NEGATIVE_INFINITY;
}
parameter.setParameterValue(dim, newValue);
return 0.0;
}
}
public double reflectValue(double value, double lower, double upper) {
double newValue = value;
if (value < lower) {
if (Double.isInfinite(upper)) {
// we are only going to reflect once as the upper bound is at infinity...
newValue = lower + (lower - value);
} else {
double remainder = lower - value;
double widths = Math.floor(remainder / (upper - lower));
remainder -= (upper - lower) * widths;
// even reflections
if (widths % 2 == 0) {
newValue = lower + remainder;
// odd reflections
} else {
newValue = upper - remainder;
}
}
} else if (value > upper) {
if (Double.isInfinite(lower)) {
// we are only going to reflect once as the lower bound is at -infinity...
newValue = upper - (newValue - upper);
} else {
double remainder = value - upper;
double widths = Math.floor(remainder / (upper - lower));
remainder -= (upper - lower) * widths;
// even reflections
if (widths % 2 == 0) {
newValue = upper - remainder;
// odd reflections
} else {
newValue = lower + remainder;
}
}
}
return newValue;
}
public double reflectValueLoop(double value, double lower, double upper) {
double newValue = value;
while (newValue < lower || newValue > upper) {
if (newValue < lower) {
newValue = lower + (lower - newValue);
}
if (newValue > upper) {
newValue = upper - (newValue - upper);
}
}
return newValue;
}
//MCMCOperator INTERFACE
public final String getOperatorName() {
return parameter.getParameterName();
}
public double getCoercableParameter() {
return Math.log(windowSize);
}
public void setCoercableParameter(double value) {
windowSize = Math.exp(value);
}
public double getRawParameter() {
return windowSize;
}
public double getTargetAcceptanceProbability() {
return 0.234;
}
public double getMinimumAcceptanceLevel() {
return 0.1;
}
public double getMaximumAcceptanceLevel() {
return 0.4;
}
public double getMinimumGoodAcceptanceLevel() {
return 0.20;
}
public double getMaximumGoodAcceptanceLevel() {
return 0.30;
}
public final String getPerformanceSuggestion() {
double prob = MCMCOperator.Utils.getAcceptanceProbability(this);
double targetProb = getTargetAcceptanceProbability();
double ws = OperatorUtils.optimizeWindowSize(windowSize, parameter.getParameterValue(0) * 2.0, prob, targetProb);
if (prob < getMinimumGoodAcceptanceLevel()) {
return "Try decreasing windowSize to about " + ws;
} else if (prob > getMaximumGoodAcceptanceLevel()) {
return "Try increasing windowSize to about " + ws;
} else return "";
}
public String toString() {
return RandomWalkOperatorParser.RANDOM_WALK_OPERATOR + "(" + parameter.getParameterName() + ", " + windowSize + ", " + getWeight() + ")";
}
//PRIVATE STUFF
private Parameter parameter = null;
private double windowSize = 0.01;
private List<Integer> updateMap = null;
private int updateMapSize;
private final BoundaryCondition boundaryCondition;
private final Double lowerOperatorBound;
private final Double upperOperatorBound;
}