package dr.evomodel.antigenic.phyloClustering.misc.obsolete; import dr.inference.model.Bounds; import dr.inference.model.Parameter; import dr.inference.operators.AbstractCoercableOperator; import dr.inference.operators.CoercionMode; import dr.inference.operators.MCMCOperator; import dr.inference.operators.OperatorUtils; import dr.math.MathUtils; import java.util.ArrayList; import java.util.List; public class ClusterWalkOperator extends AbstractCoercableOperator { public enum BoundaryCondition { reflecting, absorbing } public ClusterWalkOperator(Parameter parameter, double windowSize, BoundaryCondition bc, double weight, CoercionMode mode) { this(parameter, null, windowSize, bc, weight, mode); } public ClusterWalkOperator(Parameter parameter, Parameter updateIndex, double windowSize, BoundaryCondition bc, double weight, CoercionMode mode) { this(parameter, updateIndex, windowSize, bc, weight, mode, null, null); } public ClusterWalkOperator(Parameter parameter, Parameter updateIndex, double windowSize, BoundaryCondition bc, double weight, CoercionMode mode, Double lowerOperatorBound, Double upperOperatorBound) { super(mode); this.parameter = parameter; this.windowSize = windowSize; this.condition = bc; 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); } } 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() { System.out.println("Walking cluster"); // a random dimension to perturb int index; if (updateMap == null) { index = MathUtils.nextInt(parameter.getDimension()); } else { index = updateMap.get(MathUtils.nextInt(updateMap.size())); } // a random point around old value within windowSize * 2 double draw = (2.0 * MathUtils.nextDouble() - 1.0) * windowSize; double newValue = parameter.getParameterValue(index) + draw; final Bounds<Double> bounds = parameter.getBounds(); final double lower = (lowerOperatorBound == null ? bounds.getLowerLimit(index) : Math.max(bounds.getLowerLimit(index), lowerOperatorBound)); final double upper = (upperOperatorBound == null ? bounds.getUpperLimit(index) : Math.min(bounds.getUpperLimit(index), upperOperatorBound)); if (condition == BoundaryCondition.reflecting) { newValue = reflectValue(newValue, lower, upper); } else if (newValue < lower || newValue > upper) { // throw new OperatorFailedException("proposed value outside boundaries"); return Double.NEGATIVE_INFINITY; } parameter.setParameterValue(index, 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 ClusterWalkOperatorParser.CLUSTER_WALK_OPERATOR + "(" + parameter.getParameterName() + ", " + windowSize + ", " + getWeight() + ")"; } //PRIVATE STUFF private Parameter parameter = null; private double windowSize = 0.01; private List<Integer> updateMap = null; private final BoundaryCondition condition; private final Double lowerOperatorBound; private final Double upperOperatorBound; }