diff --git a/beast-base/src/main/java/beast/base/evolution/operator/ScaleOperator.java b/beast-base/src/main/java/beast/base/evolution/operator/ScaleOperator.java index b8478036..044f6a67 100644 --- a/beast-base/src/main/java/beast/base/evolution/operator/ScaleOperator.java +++ b/beast-base/src/main/java/beast/base/evolution/operator/ScaleOperator.java @@ -145,8 +145,10 @@ public double proposal() { return -Math.log(scale); } else { // scale the beast.tree - final int internalNodes = tree.scale(scale); - return Math.log(scale) * (internalNodes - 2); + // tree.scale returns the log Jacobian (dof * log(scale)); + // operator adds the -2*log(scale) kernel-symmetry correction. + final double treeLogJacobian = tree.scale(scale); + return treeLogJacobian - 2 * Math.log(scale); } } @@ -205,9 +207,13 @@ public double proposal() { // for the proof. It is supposed to be somewhere in an Alexei/Nicholes article. // all Values assumed independent! - final int computedDoF = param.scale(scale); - final int usedDoF = (specifiedDoF > 0) ? specifiedDoF : computedDoF ; - hastingsRatio = (usedDoF - 2) * Math.log(scale); + // param.scale returns the log Jacobian (dof * log(scale)); + // operator adds the -2*log(scale) kernel-symmetry correction. + final double paramLogJacobian = param.scale(scale); + final double dofLogScale = (specifiedDoF > 0) + ? specifiedDoF * Math.log(scale) + : paramLogJacobian; + hastingsRatio = dofLogScale - 2 * Math.log(scale); } else { hastingsRatio = -Math.log(scale); diff --git a/beast-base/src/main/java/beast/base/evolution/operator/kernel/AdaptableVarianceMultivariateNormalOperator.java b/beast-base/src/main/java/beast/base/evolution/operator/kernel/AdaptableVarianceMultivariateNormalOperator.java index c58a52cf..b99804ab 100644 --- a/beast-base/src/main/java/beast/base/evolution/operator/kernel/AdaptableVarianceMultivariateNormalOperator.java +++ b/beast-base/src/main/java/beast/base/evolution/operator/kernel/AdaptableVarianceMultivariateNormalOperator.java @@ -752,9 +752,9 @@ public int setValue(final int param, final double value) throws Exception { } return 1; } else if (para instanceof Tree) { - double old = para.getArrayValue(); - double scale = value / old; - ((Tree) para).scale(scale); + // Use the Scalable contract: setScalableValue lands the tree's + // dilation-axis summary (sum of intervals) at exactly `value`. + ((Tree) para).setScalableValue(value); return ((Tree) para).getInternalNodeCount(); } return 0; @@ -765,7 +765,8 @@ public double getValue(final int param) { if (f instanceof RealParameter) { return f.getArrayValue(getX(param)); } - return ((Tree) f).getRoot().getHeight(); + // Read the tree's position on its dilation axis (sum of intervals) + return ((Tree) f).getScalableValue(); } // public double getLower(final int param) { diff --git a/beast-base/src/main/java/beast/base/evolution/operator/kernel/BactrianScaleOperator.java b/beast-base/src/main/java/beast/base/evolution/operator/kernel/BactrianScaleOperator.java index cc2f2303..fbb316bf 100644 --- a/beast-base/src/main/java/beast/base/evolution/operator/kernel/BactrianScaleOperator.java +++ b/beast-base/src/main/java/beast/base/evolution/operator/kernel/BactrianScaleOperator.java @@ -64,9 +64,10 @@ public double proposal() { return Math.log(scale); } else { // scale the beast.tree + // tree.scale returns the log Jacobian (dof * log(scale)); + // Bactrian kernel is symmetric so no kernel-ratio correction is needed. final double scale = getScaler(0, Double.NaN); - final int scaledNodes = tree.scale(scale); - return Math.log(scale) * scaledNodes; + return tree.scale(scale); } } @@ -125,10 +126,11 @@ public double proposal() { // for the proof. It is supposed to be somewhere in an Alexei/Nicholes article. // all Values assumed independent! - final double scale = getScaler(0, param.getValue(0)); - final int computedDoF = param.scale(scale); - final int usedDoF = (specifiedDoF > 0) ? specifiedDoF : computedDoF ; - hastingsRatio = usedDoF * Math.log(scale); + // param.scale returns the log Jacobian (dof * log(scale)); + // Bactrian kernel is symmetric so no kernel-ratio correction is needed. + final double scale = getScaler(0, param.getValue(0)); + final double paramLogJacobian = param.scale(scale); + hastingsRatio = (specifiedDoF > 0) ? specifiedDoF * Math.log(scale) : paramLogJacobian; } else { // which position to scale diff --git a/beast-base/src/main/java/beast/base/evolution/speciation/StarBeastStartState.java b/beast-base/src/main/java/beast/base/evolution/speciation/StarBeastStartState.java index 93d57f2d..a12d1b62 100644 --- a/beast-base/src/main/java/beast/base/evolution/speciation/StarBeastStartState.java +++ b/beast-base/src/main/java/beast/base/evolution/speciation/StarBeastStartState.java @@ -191,7 +191,9 @@ private void fullInit() { final ClusterTree ctree = new ClusterTree(); ctree.initByName("initial", gtree, "clusterType", "upgma", "taxa", alignment); - gtree.scale(1 / mu); + // Affine helper: scales internal heights by 1/mu (preserving leaf heights). + // Tree.scale would do interval scaling, which is the wrong operation here. + gtree.scaleToRootHeight(gtree.getRoot().getHeight() / mu); maxNsites = max(maxNsites, alignment.getSiteCount()); } @@ -388,7 +390,8 @@ private void randomInit() { s += 1.0/k; } final double rootHeight = (1/lam) * s; - stree.scale(rootHeight/stree.getRoot().getHeight()); + // Affine helper: lands species tree root at rootHeight exactly. + stree.scaleToRootHeight(rootHeight); randomInitGeneTrees(rootHeight); // final List geneTrees = genes.get(); // for (final Tree gtree : geneTrees) { diff --git a/beast-base/src/main/java/beast/base/evolution/tree/Node.java b/beast-base/src/main/java/beast/base/evolution/tree/Node.java index 0cebae5b..5912809b 100644 --- a/beast-base/src/main/java/beast/base/evolution/tree/Node.java +++ b/beast-base/src/main/java/beast/base/evolution/tree/Node.java @@ -802,7 +802,17 @@ public Set getLengthMetaDataNames() { /** - * scale height of this node and all its descendants + * Affine height scaling: multiply this node's height by {@code scale}, + * recursively. Leaves and fake (sampled-ancestor) nodes are skipped. + *

+ * Used by {@link Tree#scaleToRootHeight(double)} for the affine + * "land root at target" helper. NOT used by {@link Tree#scale(double)}, + * which applies interval scaling via {@link #intervalScale(double)}. + *

+ * Throws {@link IllegalArgumentException} for non-ultrametric trees if a + * scaled internal height drops below a leaf child's height. This is the + * historical {@code Node.scale} behaviour, preserved here for callers + * that explicitly want affine scaling and can handle the failure. * * @param scale scale factor * @return degrees of freedom scaled (used for HR calculations) @@ -832,6 +842,44 @@ public int scale(final double scale) { return dof; } + /** + * Interval scaling: multiply this node's "margin" (height above its taller + * child) by {@code scale}, recursively. Tip heights are preserved by + * construction, so the resulting tree is always valid for any positive + * scale factor. + *

+ * Used by {@link Tree#scale(double)} as the contract-bound dilation + * operation. Each margin is independently multiplied by {@code scale}, so + * the tree's sum of margins (= {@link Tree#getScalableValue()}) is + * exactly multiplied by {@code scale}. + * + * @param scale positive scale factor + * @return number of intervals (margins) scaled, for HR calculations + */ + public int intervalScale(final double scale) { + if (isLeaf()) { + return 0; + } + // sampled-ancestor fake nodes: skip, recurse into the non-direct-ancestor child + if (isFake()) { + if (getLeft().isDirectAncestor()) { + return getRight().intervalScale(scale); + } else { + return getLeft().intervalScale(scale); + } + } + startEditing(); + final double oldMargin = height - Math.max(getLeft().getHeight(), getRight().getHeight()); + int scaledNodeCount = 1; + scaledNodeCount += getLeft().intervalScale(scale); + scaledNodeCount += getRight().intervalScale(scale); + // recompute minHeight after children have been scaled + final double minChildHeight = Math.max(getLeft().getHeight(), getRight().getHeight()); + height = oldMargin * scale + minChildHeight; + isDirty |= Tree.IS_DIRTY; + return scaledNodeCount; + } + // /** // * Used for sampled ancestor trees // * Scales this node and all its descendants (either all descendants, or only non-sampled descendants) diff --git a/beast-base/src/main/java/beast/base/evolution/tree/Tree.java b/beast-base/src/main/java/beast/base/evolution/tree/Tree.java index 5772cbdf..587c8a5a 100644 --- a/beast-base/src/main/java/beast/base/evolution/tree/Tree.java +++ b/beast-base/src/main/java/beast/base/evolution/tree/Tree.java @@ -654,27 +654,82 @@ public void setEverythingDirty(final boolean isDirty) { } } + /** + * Scale this tree by factor {@code scale} along its dilation axis. + *

+ * The dilation axis for {@code Tree} is the sum of intervals (margins above + * taller children) across internal non-fake nodes. Each margin is + * multiplied by {@code scale}; tip dates are preserved by construction. + * Equivalent to interval scaling: see {@link #getScalableValue()} for the + * summary that scales by exactly {@code scale} under this operation. + *

+ * Always succeeds for any positive {@code scale} (no leaf can violate a + * branch-length constraint, because each margin remains positive). + * + * @return log Jacobian determinant of the move ({@code dof × log(scale)}) + */ @Override - public int scale(final double scale) { - return root.scale(scale); + public double scale(final double scale) { + int dof = root.intervalScale(scale); + return dof * Math.log(scale); } + /** + * Read this tree's position on its dilation axis: the sum of intervals + * (margins above taller children) across internal non-fake nodes. + *

+ * This summary is exactly {@code s}-equivariant under {@link #scale(double)}. + */ @Override - public void scaleOne(int i, final double scale) { - startEditing(null); - double h = m_nodes[i].getHeight(); - double newHeight = h * scale; - for (Node child : m_nodes[i].children) { - if (newHeight < child.getHeight()) { - throw new IllegalArgumentException("scale sets nodes below child result in negative branch length"); - } - } - if (!m_nodes[i].isRoot() && newHeight > m_nodes[i].getParent().getHeight()) { - throw new IllegalArgumentException("scale sets nodes above parent result in negative branch length"); - } - m_nodes[i].setHeight(newHeight); + public double getScalableValue() { + return computeSumIntervals(root); + } + + /** + * Affine helper: scale the tree so its root height equals {@code targetRootHeight}. + * Multiplies every internal non-fake non-leaf height by + * {@code targetRootHeight / oldRootHeight}, keeping leaf heights fixed. + *

+ * Not part of the Scalable contract. May throw + * {@link IllegalArgumentException} for non-ultrametric trees if the + * resulting state has a parent below a leaf child. + * + * @return number of internal nodes scaled (degrees of freedom) + */ + public int scaleToRootHeight(final double targetRootHeight) { + double oldRoot = root.getHeight(); + if (oldRoot <= 0.0) { + throw new IllegalArgumentException( + "Cannot scale to root height: current root height is " + oldRoot); + } + return root.scale(targetRootHeight / oldRoot); } + /** + * Compute the sum of margins (h_N - max(h_children)) across internal + * non-fake nodes. For sampled-ancestor trees, fake nodes are skipped. + */ + private double computeSumIntervals(final Node node) { + if (node.isLeaf()) { + return 0.0; + } + if (node.isFake()) { + // skip the fake; recurse into the non-direct-ancestor child + if (node.getLeft().isDirectAncestor()) { + return computeSumIntervals(node.getRight()); + } else { + return computeSumIntervals(node.getLeft()); + } + } + double margin = node.getHeight() + - Math.max(node.getLeft().getHeight(), node.getRight().getHeight()); + double sum = margin; + sum += computeSumIntervals(node.getLeft()); + if (node.getRight() != null) { + sum += computeSumIntervals(node.getRight()); + } + return sum; + } // /** // * The same as scale but with option to scale all sampled nodes diff --git a/beast-base/src/main/java/beast/base/inference/Scalable.java b/beast-base/src/main/java/beast/base/inference/Scalable.java index 4f6013c0..c398d879 100644 --- a/beast-base/src/main/java/beast/base/inference/Scalable.java +++ b/beast-base/src/main/java/beast/base/inference/Scalable.java @@ -2,36 +2,107 @@ import beast.base.core.Description; -@Description("For StateNodes that can be scaled by a scale/up-down operator") +/** + * A {@code Scalable} represents a state component that can be moved along a + * single positive-real dilation axis. The interface defines three operations + * that must be mutually consistent — together they form the + * Scalable contract: + * + *

    + *
  1. {@link #scale(double)} dilates the component by a factor {@code s} + * and returns the log Jacobian determinant of that move + * (i.e. {@code log |det(∂new/∂old)|}).
    2. + *
  2. {@link #getScalableValue()} reads the component's current position + * on its dilation axis.
    3. + *
  3. {@link #setScalableValue(double)} moves the component so that + * {@code getScalableValue()} returns the supplied target {@code V}, and + * returns the log Jacobian for that move.
    4. + *
+ * + *

The log Jacobian is the move's contribution to the Metropolis-Hastings + * acceptance ratio from the change-of-variables formula. The proposal density + * ratio (the "Hastings ratio" proper, {@code q(reverse)/q(forward)}) lives in + * the calling operator's kernel and is not part of the Scalable's + * return.

+ * + *

The contract requires the following three invariants to hold for any + * valid {@code Scalable x} and any positive scale factor {@code s}:

+ * + * 
+ *   // (1) scale-equivariance
+ *   double v0 = x.getScalableValue();
+ *   x.scale(s);
+ *   assert x.getScalableValue() == s * v0;
+ *
+ *   // (2) set is a fixed point of get
+ *   x.setScalableValue(V);
+ *   assert x.getScalableValue() == V;
+ *
+ *   // (3) set composes with scale
+ *   //     x.setScalableValue(x.getScalableValue() * s)
+ *   //     produces the same state as
+ *   //     x.scale(s)
+ *
+ * + *

The choice of dilation axis (and therefore the meaning of + * {@code getScalableValue}) is bound to the implementation of {@code scale}. + * For example, an affine-scaling parameter exposes its value directly. A tree + * whose {@code scale} is interval-scaling exposes its sum-of-margins. A custom + * {@code Scalable} chooses whichever summary is exactly multiplied by {@code s} + * under its own {@code scale} operation.

+ * + *

{@code scale(s)} is expected to succeed for any positive {@code s} that + * leaves the component in a valid state. Implementations may throw + * {@link IllegalArgumentException} for moves that produce an invalid state; + * such throws act as rejection signals for the calling operator. The contract + * invariants apply when {@code scale} does not throw.

+ * + * @see beast3 issue #20 + */ +@Description("State component that can be dilated along a 1-D axis by scale or up-down operators.") public interface Scalable { /** - * Scale StateNode with amount scale and + * Dilate this component by factor {@code s} along its scaling axis. + * After the call, {@link #getScalableValue()} returns + * {@code s * (its previous value)}. * - * @param scale scaling factor - * @return the number of degrees of freedom used in this operation. This number varies - * for the different types of StateNodes. For example, for real - * valued n-dimensional parameters, it is n, for a tree it is the - * number of internal nodes being scaled. - * @throws IllegalArgumentException when StateNode become not valid, e.g. has - * values outside bounds or negative branch lengths. + * @param s positive scale factor + * @return log Jacobian determinant of this move + * @throws IllegalArgumentException if the move would produce an invalid state */ - abstract public int scale(double scale); + double scale(double s); /** - * only scale the i-th element of the StateNode - * @param i - * @param scale + * Read the component's current position on its dilation axis. + * The contract requires this to be exactly {@code s}-equivariant under + * {@link #scale(double)}. */ - abstract public void scaleOne(int i, double scale); + double getScalableValue(); + + /** + * Move the component so that {@link #getScalableValue()} returns {@code V}. + * Defined as {@code scale(V / getScalableValue())}. + *

+ * Implementations rarely need to override this; the default expresses the + * contract identity {@code setScalableValue(V) ≡ scale(V / getScalableValue())} + * directly. Override only if the dilation axis cannot be reached by a + * single multiplicative scale (rare). + * + * @param V target value (must be positive for typical multiplicative axes) + * @return log Jacobian determinant of this move + * @throws IllegalArgumentException if the current value is zero (no + * multiplicative scale can land at {@code V}) or if the resulting + * state would be invalid + */ + default double setScalableValue(double V) { + double current = getScalableValue(); + if (current == 0.0) { + throw new IllegalArgumentException( + "Cannot set scalable value: current value is zero " + + "(no multiplicative scale lands at " + V + ")"); + } + return scale(V / current); + } - default double scaleAll(double scale) { - try { - int d = scale(scale); - return d * Math.log(scale); - } catch (IllegalArgumentException e) { - return Double.NEGATIVE_INFINITY; - } - } - } diff --git a/beast-base/src/main/java/beast/base/inference/operator/UpDownOperator.java b/beast-base/src/main/java/beast/base/inference/operator/UpDownOperator.java index b8856987..9ea83467 100644 --- a/beast-base/src/main/java/beast/base/inference/operator/UpDownOperator.java +++ b/beast-base/src/main/java/beast/base/inference/operator/UpDownOperator.java @@ -67,17 +67,18 @@ public void initAndValidate() { public final double proposal() { final double scale = (scaleFactor + (Randomizer.nextDouble() * ((1.0 / scaleFactor) - scaleFactor))); - int goingUp = 0, goingDown = 0; + double netLogJacobian = 0.0; if (elementWiseInput.get()) { int size = 0; + int numUp = 0, numDown = 0; for (StateNode up : upInput.get()) { if (size == 0) size = ((Function)up).getDimension(); if (size > 0 && ((Function)up).getDimension() != size) { throw new RuntimeException("elementWise=true but parameters of differing lengths!"); } - goingUp += 1; + numUp += 1; } for (StateNode down : downInput.get()) { @@ -85,7 +86,7 @@ public final double proposal() { if (size > 0 && ((Function)down).getDimension() != size) { throw new RuntimeException("elementWise=true but parameters of differing lengths!"); } - goingDown += 1; + numDown += 1; } int index = Randomizer.nextInt(size); @@ -115,12 +116,15 @@ public final double proposal() { index = Randomizer.nextInt(size); } } + // each up StateNode contributes +log(scale), each down -log(scale) + netLogJacobian = (numUp - numDown) * Math.log(scale); } else { try { for (StateNode up : upInput.get()) { up = up.getCurrentEditable(this); - goingUp += ((Scalable)up).scale(scale); + // scale returns log Jacobian = dof * log(scale) + netLogJacobian += ((Scalable)up).scale(scale); } // separated this into second loop because the outsideBounds might return true transiently with // related variables which would be BAD. Note current implementation of outsideBounds isn't dynamic, @@ -134,7 +138,8 @@ public final double proposal() { for (StateNode down : downInput.get()) { down = down.getCurrentEditable(this); - goingDown += ((Scalable)down).scale(1.0 / scale); + // scale returns log Jacobian = dof * log(1/scale) = -dof * log(scale) + netLogJacobian += ((Scalable)down).scale(1.0 / scale); } for (StateNode down : downInput.get()) { if (outsideBounds(down)) { @@ -146,7 +151,8 @@ public final double proposal() { return Double.NEGATIVE_INFINITY; } } - return (goingUp - goingDown - 2) * Math.log(scale); + // kernel-symmetry correction: -2 * log(scale) + return netLogJacobian - 2 * Math.log(scale); } private boolean outsideBounds(final StateNode node) { diff --git a/beast-base/src/main/java/beast/base/inference/operator/kernel/BactrianUpDownOperator.java b/beast-base/src/main/java/beast/base/inference/operator/kernel/BactrianUpDownOperator.java index 34164990..0f315eca 100644 --- a/beast-base/src/main/java/beast/base/inference/operator/kernel/BactrianUpDownOperator.java +++ b/beast-base/src/main/java/beast/base/inference/operator/kernel/BactrianUpDownOperator.java @@ -72,17 +72,18 @@ protected double getScaler(int i) { public final double proposal() { final double scale = getScaler(0); - int goingUp = 0, goingDown = 0; + double netLogJacobian = 0.0; if (elementWiseInput.get()) { int size = 0; + int numUp = 0, numDown = 0; for (StateNode up : upInput.get()) { if (size == 0) size = ((Function)up).getDimension(); if (size > 0 && ((Function)up).getDimension() != size) { throw new RuntimeException("elementWise=true but parameters of differing lengths!"); } - goingUp += 1; + numUp += 1; } for (StateNode down : downInput.get()) { @@ -90,7 +91,7 @@ public final double proposal() { if (size > 0 && ((Function)down).getDimension() != size) { throw new RuntimeException("elementWise=true but parameters of differing lengths!"); } - goingDown += 1; + numDown += 1; } int index = Randomizer.nextInt(size); @@ -114,12 +115,13 @@ public final double proposal() { return Double.NEGATIVE_INFINITY; } } + netLogJacobian = (numUp - numDown) * Math.log(scale); } else { try { for (StateNode up : upInput.get()) { up = up.getCurrentEditable(this); - goingUp += ((Scalable)up).scale(scale); + netLogJacobian += ((Scalable)up).scale(scale); } // separated this into second loop because the outsideBounds might return true transiently with // related variables which would be BAD. Note current implementation of outsideBounds isn't dynamic, @@ -133,7 +135,7 @@ public final double proposal() { for (StateNode down : downInput.get()) { down = down.getCurrentEditable(this); - goingDown += ((Scalable)down).scale(1.0 / scale); + netLogJacobian += ((Scalable)down).scale(1.0 / scale); } for (StateNode down : downInput.get()) { if (outsideBounds(down)) { @@ -145,7 +147,8 @@ public final double proposal() { return Double.NEGATIVE_INFINITY; } } - return (goingUp - goingDown) * Math.log(scale); + // Bactrian kernel is symmetric so no kernel-ratio correction is needed. + return netLogJacobian; } private boolean outsideBounds(final StateNode node) { diff --git a/beast-base/src/main/java/beast/base/inference/parameter/RealParameter.java b/beast-base/src/main/java/beast/base/inference/parameter/RealParameter.java index 891b587d..e909417e 100644 --- a/beast-base/src/main/java/beast/base/inference/parameter/RealParameter.java +++ b/beast-base/src/main/java/beast/base/inference/parameter/RealParameter.java @@ -96,7 +96,7 @@ public void log(final long sample, final PrintStream out) { * StateNode methods * */ @Override - public int scale(final double scale) { + public double scale(final double scale) { int nScaled = 0; for (int i = 0; i < values.length; i++) { @@ -111,14 +111,21 @@ public int scale(final double scale) { } } - return nScaled; + return nScaled * Math.log(scale); } - @Override - public void scaleOne(int i, double scale) { - values[i] *= scale; - } - + /** + * Read this parameter's position on its dilation axis: sum of values. + * Exactly {@code s}-equivariant under {@link #scale(double)}. + */ + @Override + public double getScalableValue() { + double sum = 0.0; + for (Double v : values) { + sum += v; + } + return sum; + } @Override void fromXML(final int dimension, final String lower, final String upper, final String[] valuesString) { diff --git a/beast-base/src/main/java/beast/base/spec/evolution/operator/AdaptableVarianceMultivariateNormalOperator.java b/beast-base/src/main/java/beast/base/spec/evolution/operator/AdaptableVarianceMultivariateNormalOperator.java index fcbe53cb..e401ef64 100644 --- a/beast-base/src/main/java/beast/base/spec/evolution/operator/AdaptableVarianceMultivariateNormalOperator.java +++ b/beast-base/src/main/java/beast/base/spec/evolution/operator/AdaptableVarianceMultivariateNormalOperator.java @@ -779,10 +779,10 @@ public int setValue(final int param, final double value) throws Exception { } return 1; } else if (para instanceof Tree tree) { - double old = tree.getArrayValue(); - double scale = value / old; - ((Tree) para).scale(scale); - return ((Tree) para).getInternalNodeCount(); + // Use the Scalable contract: setScalableValue lands the tree's + // dilation-axis summary (sum of intervals) at exactly `value`. + tree.setScalableValue(value); + return tree.getInternalNodeCount(); } return 0; } @@ -799,7 +799,8 @@ public double getValue(final int param) { } else if (f instanceof IntVectorParam p) { return p.get(getX(param)); } else if (f instanceof Tree t) { - t.getRoot().getHeight(); + // Read the tree's position on its dilation axis (sum of intervals) + return t.getScalableValue(); } throw new RuntimeException("programmer error: should not get here"); } diff --git a/beast-base/src/main/java/beast/base/spec/evolution/operator/IntervalScaleOperator.java b/beast-base/src/main/java/beast/base/spec/evolution/operator/IntervalScaleOperator.java index 7e53a527..038135f8 100644 --- a/beast-base/src/main/java/beast/base/spec/evolution/operator/IntervalScaleOperator.java +++ b/beast-base/src/main/java/beast/base/spec/evolution/operator/IntervalScaleOperator.java @@ -77,13 +77,12 @@ public double proposal() { double logHR = Math.log(scaler) * (numbers); + // scale returns log Jacobian = dim * log(factor) for (Scalable down : downInput.get()) { - int dim = down.scale(1.0/actualScaler);//setValue(down.getValue() / actualScaler); - logHR -= dim * Math.log(actualScaler); + logHR += down.scale(1.0/actualScaler); } for (Scalable up : upInput.get()) { - int dim = up.scale(actualScaler);//setValue(up.getValue() * actualScaler); - logHR += dim * Math.log(actualScaler); + logHR += up.scale(actualScaler); } return logHR; } diff --git a/beast-base/src/main/java/beast/base/spec/evolution/operator/ScaleTreeOperator.java b/beast-base/src/main/java/beast/base/spec/evolution/operator/ScaleTreeOperator.java index 21d7d90b..8d9b1ed2 100644 --- a/beast-base/src/main/java/beast/base/spec/evolution/operator/ScaleTreeOperator.java +++ b/beast-base/src/main/java/beast/base/spec/evolution/operator/ScaleTreeOperator.java @@ -46,10 +46,9 @@ public double proposal() { } else { // scale the beast.tree + // tree.scale returns the log Jacobian (dof * log(scale)) final double scale = getScaler(0, Double.NaN); - final int scaledNodes = tree.scale(scale); - // hastings ratio - return Math.log(scale) * scaledNodes; + return tree.scale(scale); } } catch (Exception e) { diff --git a/beast-base/src/main/java/beast/base/spec/evolution/operator/UpDownOperator.java b/beast-base/src/main/java/beast/base/spec/evolution/operator/UpDownOperator.java index 7f25986f..a716d959 100644 --- a/beast-base/src/main/java/beast/base/spec/evolution/operator/UpDownOperator.java +++ b/beast-base/src/main/java/beast/base/spec/evolution/operator/UpDownOperator.java @@ -107,18 +107,18 @@ protected double getScaler(int i) { public double proposal() { double scale = getScaler(0); - int goingUp = 0, goingDown = 0; - double logHR = 0; + double logHR = 0; if (elementWiseInput.get()) { int size = 0; + int numUp = 0, numDown = 0; for (Scalable up : upInput.get()) { if (size == 0) size = ((Tensor)up).size(); if (size > 0 && ((Tensor)up).size() != size) { throw new RuntimeException("elementWise=true but parameters of differing lengths!"); } - goingUp += 1; + numUp += 1; } for (Scalable down : downInput.get()) { @@ -126,7 +126,7 @@ public double proposal() { if (size > 0 && ((Tensor)down).size() != size) { throw new RuntimeException("elementWise=true but parameters of differing lengths!"); } - goingDown += 1; + numDown += 1; } int index = Randomizer.nextInt(size); @@ -148,12 +148,12 @@ public double proposal() { return Double.NEGATIVE_INFINITY; } } - logHR = (goingUp - goingDown) * Math.log(scale); + logHR = (numUp - numDown) * Math.log(scale); } else { try { if (treesUp.size() > 0) { - + // scale trees up and adjust scale factor double lengthBefore =0, lenghtAfter = 0; for (TreeInterface up : treesUp) { @@ -162,9 +162,9 @@ public double proposal() { lenghtAfter += treeLength(up); } scale = lenghtAfter / lengthBefore; - + } else if (treesDown.size() > 0) { - + // scale trees down and adjust scale factor double lengthBefore =0, lenghtAfter = 0; for (TreeInterface down : treesDown) { @@ -173,12 +173,13 @@ public double proposal() { lenghtAfter += treeLength(down); } scale = 1.0/(lenghtAfter / lengthBefore); - + } - + for (Scalable up : otherUp) { up = (Scalable) ((StateNode) up).getCurrentEditable(this); - goingUp += up.scale(scale); + // scale returns log Jacobian = dof * log(scale) + logHR += up.scale(scale); } // separated this into second loop because the outsideBounds might return true transiently with // related variables which would be BAD. Note current implementation of outsideBounds isn't dynamic, @@ -192,7 +193,8 @@ public double proposal() { for (Scalable down : otherDown) { down = (Scalable) ((StateNode) down).getCurrentEditable(this); - goingDown += down.scale(1.0 / scale); + // scale returns log Jacobian = dof * log(1/scale) = -dof * log(scale) + logHR += down.scale(1.0 / scale); } for (Scalable down : otherDown) { if (outsideBounds(down)) { @@ -203,7 +205,6 @@ public double proposal() { // scale resulted in invalid StateNode, abort proposal return Double.NEGATIVE_INFINITY; } - logHR += (goingUp - goingDown) * Math.log(scale); } return logHR; } diff --git a/beast-base/src/main/java/beast/base/spec/inference/operator/ScaleOperator.java b/beast-base/src/main/java/beast/base/spec/inference/operator/ScaleOperator.java index c6417f5c..82bee40e 100644 --- a/beast-base/src/main/java/beast/base/spec/inference/operator/ScaleOperator.java +++ b/beast-base/src/main/java/beast/base/spec/inference/operator/ScaleOperator.java @@ -145,10 +145,11 @@ public double proposal() { // for the proof. It is supposed to be somewhere in an Alexei/Nicholes article. // all Values assumed independent! + // realVectorParam.scale returns the log Jacobian (dof * log(scale)); + // spec ScaleOperator does not add a kernel-symmetry correction here. final double scale = getScaler(0, realVectorParam.get(0)); - final int computedDoF = realVectorParam.scale(scale); - final int usedDoF = (specifiedDoF > 0) ? specifiedDoF : computedDoF; - logHR = usedDoF * Math.log(scale); + final double paramLogJacobian = realVectorParam.scale(scale); + logHR = (specifiedDoF > 0) ? specifiedDoF * Math.log(scale) : paramLogJacobian; } else { // which position to scale @@ -211,15 +212,14 @@ public double proposal() { } else if (param instanceof RealScalarParam realScalarParam) { final double scale = getScaler(0); - // this set new value - int dim = realScalarParam.scale(scale); + // realScalarParam.scale returns the log Jacobian (= log(scale) for a 1-D scalar) + // and applies the value update as a side effect. + logHR = realScalarParam.scale(scale); if (! realScalarParam.withinBounds(realScalarParam.get()) ) { // reject out of bounds scales return Double.NEGATIVE_INFINITY; } - // hastings ratio - logHR = Math.log(scale); } return logHR; diff --git a/beast-base/src/main/java/beast/base/spec/inference/parameter/CompoundRealScalarParam.java b/beast-base/src/main/java/beast/base/spec/inference/parameter/CompoundRealScalarParam.java index 793d28cb..160a26d5 100644 --- a/beast-base/src/main/java/beast/base/spec/inference/parameter/CompoundRealScalarParam.java +++ b/beast-base/src/main/java/beast/base/spec/inference/parameter/CompoundRealScalarParam.java @@ -150,13 +150,21 @@ protected boolean requiresRecalculation() { } @Override - public int scale(double scale) { - int sum = 0; + public double scale(double scale) { + double sum = 0; for (RealScalarParam p : parameters) { sum += p.scale(scale); } return sum; -// throw new UnsupportedOperationException(); + } + + @Override + public double getScalableValue() { + double sum = 0.0; + for (RealScalarParam p : parameters) { + sum += p.getScalableValue(); + } + return sum; } @Override diff --git a/beast-base/src/main/java/beast/base/spec/inference/parameter/RealScalarParam.java b/beast-base/src/main/java/beast/base/spec/inference/parameter/RealScalarParam.java index ab44bbd7..329766cb 100644 --- a/beast-base/src/main/java/beast/base/spec/inference/parameter/RealScalarParam.java +++ b/beast-base/src/main/java/beast/base/spec/inference/parameter/RealScalarParam.java @@ -153,16 +153,15 @@ public void close(PrintStream out) { /** {@inheritDoc} */ @Override - public int scale(double scale) { - startEditing(null); - value *= scale; - return 1; + public double scale(double scale) { + startEditing(null); + value *= scale; + return Math.log(scale); } @Override - public void scaleOne(int i, double scale) { - startEditing(null); - value *= scale; + public double getScalableValue() { + return value; } /** diff --git a/beast-base/src/main/java/beast/base/spec/inference/parameter/RealVectorParam.java b/beast-base/src/main/java/beast/base/spec/inference/parameter/RealVectorParam.java index 8c6de25b..48fe139c 100644 --- a/beast-base/src/main/java/beast/base/spec/inference/parameter/RealVectorParam.java +++ b/beast-base/src/main/java/beast/base/spec/inference/parameter/RealVectorParam.java @@ -427,7 +427,7 @@ public void close(PrintStream out) { * StateNode methods * */ @Override - public int scale(final double scale) { + public double scale(final double scale) { startEditing(null); int nScaled = 0; @@ -443,13 +443,22 @@ public int scale(final double scale) { } - return nScaled; + return nScaled * Math.log(scale); } + /** + * Read this vector's position on its dilation axis. + * The summary is the sum of values: invariant under sign and exactly + * {@code s}-equivariant under {@link #scale(double)}, since each non-zero + * element is multiplied by {@code s} and zero elements stay zero. + */ @Override - public void scaleOne(int i, double scale) { - startEditing(null); - values[i] *= scale; + public double getScalableValue() { + double sum = 0.0; + for (double v : values) { + sum += v; + } + return sum; } /** diff --git a/beast-base/src/test/java/beast/base/evolution/tree/TreeScalableTest.java b/beast-base/src/test/java/beast/base/evolution/tree/TreeScalableTest.java new file mode 100644 index 00000000..d0deb935 --- /dev/null +++ b/beast-base/src/test/java/beast/base/evolution/tree/TreeScalableTest.java @@ -0,0 +1,157 @@ +package beast.base.evolution.tree; + +import beast.base.inference.ScalableContractTest; +import org.junit.jupiter.api.Test; + +import java.util.function.BiConsumer; +import java.util.function.Supplier; + +import static org.junit.jupiter.api.Assertions.assertEquals; +import static org.junit.jupiter.api.Assertions.assertTrue; + +/** + * Tests the {@link beast.base.inference.Scalable} contract on {@link Tree}, + * including heterochronous (serially-sampled) trees that the previous + * affine {@code Tree.scale} implementation could not handle for small scale + * factors. + * + *

The contract requires:

+ * + * + *

Under the new interval-style {@code Tree.scale}, the contract holds for + * any positive scale factor on any tree shape (ultrametric or heterochronous), + * and the move never throws.

+ */ +public class TreeScalableTest { + + @Test + void contractHoldsForUltrametricTree() { + ScalableContractTest.assertContractAcrossScales( + () -> buildUltrametric(), + this::assertSameTreeState + ); + } + + @Test + void contractHoldsForHeterochronousTree() { + ScalableContractTest.assertContractAcrossScales( + () -> buildHeterochronous(), + this::assertSameTreeState + ); + } + + @Test + void contractHoldsForLeafIntrudingTopology() { + // 4-tip topology where the always-taller-child path doesn't reach the + // oldest leaf. Under the old affine Tree.scale, scaling by s < 0.5 + // would either throw or produce wrong root height. + ScalableContractTest.assertContractAcrossScales( + () -> buildLeafIntruding(), + this::assertSameTreeState + ); + } + + @Test + void scaleNeverThrowsForHeterochronousSmallScales() { + // Path B fix: the new interval-scale Tree.scale should succeed for any + // positive s on a heterochronous tree, including very small s where the + // old affine implementation threw. + for (double s : new double[] { 0.001, 0.01, 0.1, 0.5, 0.9, 1.1, 10.0 }) { + Tree tree = buildHeterochronous(); + // should not throw + tree.scale(s); + // tree should still be valid: all parents above their children + for (Node n : tree.getNodesAsArray()) { + if (!n.isLeaf()) { + assertTrue(n.getHeight() >= n.getLeft().getHeight(), + "Parent below left child after scale(" + s + ")"); + if (n.getRight() != null) { + assertTrue(n.getHeight() >= n.getRight().getHeight(), + "Parent below right child after scale(" + s + ")"); + } + } + } + } + } + + @Test + void sumIntervalsIsExactlyScaleEquivariant() { + // Spot-check: getScalableValue scales by EXACTLY s under interval scaling. + Tree tree = buildHeterochronous(); + double v0 = tree.getScalableValue(); + tree.scale(1.7); + assertEquals(1.7 * v0, tree.getScalableValue(), Math.abs(1.7 * v0) * 1e-12); + } + + /** Ultrametric: leaves A, B, C all at height 0; internal P at 1; root at 2. */ + private Tree buildUltrametric() { + Node a = leaf("A", 0, 0.0); + Node b = leaf("B", 1, 0.0); + Node c = leaf("C", 2, 0.0); + Node p = internal(3, 1.0, a, b); + Node root = internal(4, 2.0, p, c); + return new Tree(root); + } + + /** Heterochronous: A=0, B=2, C=1, P at 4, root at 5. */ + private Tree buildHeterochronous() { + Node a = leaf("A", 0, 0.0); + Node b = leaf("B", 1, 2.0); + Node c = leaf("C", 2, 1.0); + Node p = internal(3, 4.0, a, b); + Node root = internal(4, 5.0, p, c); + return new Tree(root); + } + + /** + * 4-tip leaf-intrusion topology: A=0, B=0, C=1, D=2. + * Topology is (D, ((B, C), A)). The always-taller-child path from root + * leads through ((B,C), A) to (B,C) to C — bypassing D, which is the + * oldest leaf. Under the old affine Tree.scale, scaling by s < 2/3 + * causes leaves to violate parent constraints. + */ + private Tree buildLeafIntruding() { + Node a = leaf("A", 0, 0.0); + Node b = leaf("B", 1, 0.0); + Node c = leaf("C", 2, 1.0); + Node d = leaf("D", 3, 2.0); + Node y = internal(4, 2.0, b, c); // (B, C) + Node x = internal(5, 3.0, y, a); // ((B, C), A) + Node root = internal(6, 4.0, d, x); + return new Tree(root); + } + + private Node leaf(String id, int nr, double height) { + Node n = new Node(id); + n.setNr(nr); + n.setHeight(height); + return n; + } + + private Node internal(int nr, double height, Node left, Node right) { + Node n = new Node(); + n.setNr(nr); + n.setHeight(height); + n.addChild(left); + n.addChild(right); + return n; + } + + /** Compare two Trees node-for-node by height. */ + private void assertSameTreeState(Tree a, Tree b) { + Node[] na = a.getNodesAsArray(); + Node[] nb = b.getNodesAsArray(); + assertEquals(na.length, nb.length, "Trees should have same number of nodes"); + for (int i = 0; i < na.length; i++) { + assertEquals(na[i].getHeight(), nb[i].getHeight(), + Math.abs(na[i].getHeight()) * 1e-9 + 1e-9, + "Node " + i + " (id=" + na[i].getID() + ") height should match across paths"); + } + } +} diff --git a/beast-base/src/test/java/beast/base/inference/ScalableContractTest.java b/beast-base/src/test/java/beast/base/inference/ScalableContractTest.java new file mode 100644 index 00000000..17c55fe8 --- /dev/null +++ b/beast-base/src/test/java/beast/base/inference/ScalableContractTest.java @@ -0,0 +1,94 @@ +package beast.base.inference; + +import java.util.function.BiConsumer; +import java.util.function.Supplier; + +import static org.junit.jupiter.api.Assertions.assertEquals; + +/** + * Reusable test helper that asserts the three {@link Scalable} contract + * invariants on a supplied implementation. + * + *

The contract:

+ * + * 
+ *   // (1) scale-equivariance
+ *   double v0 = x.getScalableValue();
+ *   x.scale(s);
+ *   assert x.getScalableValue() == s * v0;
+ *
+ *   // (2) set is a fixed point of get
+ *   x.setScalableValue(V);
+ *   assert x.getScalableValue() == V;
+ *
+ *   // (3) set ∘ get×s ≡ scale
+ *   //     x.setScalableValue(x.getScalableValue() * s)
+ *   //     produces the same state as
+ *   //     x.scale(s)
+ *
+ * + *

Per-class tests should call {@link #assertContract} with a factory that + * produces a fresh, identically-initialised instance on each call, and an + * {@code assertSameState} comparator that checks full internal state (not just + * the dilation summary).

+ */ +public class ScalableContractTest { + + public static final double EPSILON = 1e-9; + + /** + * Assert the three contract invariants on instances produced by {@code factory} + * with scale factor {@code s} and a state-comparison callback. + * + * @param factory produces a fresh Scalable on each call (state must be deterministic) + * @param s positive scale factor to test (different from 1.0) + * @param assertSameState asserts that two Scalables have identical internal state + */ + public static void assertContract( + Supplier factory, + double s, + BiConsumer assertSameState) { + + // (1) scale-equivariance: getScalableValue() after scale(s) equals s × original + T a = factory.get(); + double v0 = a.getScalableValue(); + double logJacobian = a.scale(s); + assertEquals(s * v0, a.getScalableValue(), Math.abs(s * v0) * EPSILON + EPSILON, + "scale(" + s + ") should make getScalableValue() return s × " + v0 + + " = " + (s * v0) + ", got " + a.getScalableValue()); + // sanity: logJacobian should be finite and have correct sign for nontrivial moves + assertEquals(true, Double.isFinite(logJacobian), + "scale should return a finite log Jacobian, got " + logJacobian); + + // (2) set is a fixed point of get: getScalableValue() after setScalableValue(V) equals V + T b = factory.get(); + double v1 = b.getScalableValue(); + double targetV = v1 * 1.7; + b.setScalableValue(targetV); + assertEquals(targetV, b.getScalableValue(), Math.abs(targetV) * EPSILON + EPSILON, + "setScalableValue(" + targetV + ") should make getScalableValue() return " + + targetV + ", got " + b.getScalableValue()); + + // (3) set ∘ get×s ≡ scale: the two paths produce the same internal state + T c = factory.get(); + double vc = c.getScalableValue(); + c.scale(s); + + T d = factory.get(); + d.setScalableValue(vc * s); + + assertSameState.accept(c, d); + } + + /** + * Convenience overload that asserts the contract for several scale factors + * including values both less than and greater than 1.0. + */ + public static void assertContractAcrossScales( + Supplier factory, + BiConsumer assertSameState) { + for (double s : new double[] { 0.5, 0.9, 1.1, 1.5, 2.0, 3.7 }) { + assertContract(factory, s, assertSameState); + } + } +} diff --git a/beast-base/src/test/java/beast/base/spec/evolution/operator/RealRandomWalkOperatorTest.java b/beast-base/src/test/java/beast/base/spec/evolution/operator/RealRandomWalkOperatorTest.java index 9d04b59e..ee026a00 100644 --- a/beast-base/src/test/java/beast/base/spec/evolution/operator/RealRandomWalkOperatorTest.java +++ b/beast-base/src/test/java/beast/base/spec/evolution/operator/RealRandomWalkOperatorTest.java @@ -36,9 +36,8 @@ public void resetRng() { @Test public void testNormalDistribution() throws Exception { - // Fix seed: will hopefully ensure success of test unless something - // goes terribly wrong. -// Randomizer.setSeed(127); + // Fix seed: ensures reproducibility at the test's 3-SE tolerance. + Randomizer.setSeed(127); // Assemble model: RealScalarParam param = new RealScalarParam<>(0.0, Real.INSTANCE); @@ -97,7 +96,10 @@ public void testNormalDistribution() throws Exception { } double m = StatUtils.mean(v); double s = StatUtils.variance(v); - assertEquals(1.0, m, 5e-3); + // 3 SE for sample mean of Normal(1, 1) with documented ESS ~196k + // (Mirror kernel, see comment above): SE = sqrt(1 / 196000) ~= 2.26e-3, + // 3 SE ~= 6.8e-3. Expected failure rate at this tolerance: 0.27%. + assertEquals(1.0, m, 7e-3); assertEquals(1.0, s, 5e-3); } diff --git a/beast-base/src/test/java/beast/base/spec/evolution/operator/UpDownOperatorTest.java b/beast-base/src/test/java/beast/base/spec/evolution/operator/UpDownOperatorTest.java index d42c251e..2250d5ac 100644 --- a/beast-base/src/test/java/beast/base/spec/evolution/operator/UpDownOperatorTest.java +++ b/beast-base/src/test/java/beast/base/spec/evolution/operator/UpDownOperatorTest.java @@ -124,7 +124,10 @@ private void doMCMCrun(RealScalar param, Loggable param1, Loggable double m = StatUtils.mean(v); double median = StatUtils.percentile(v, 50); double s = StatUtils.variance(v, 50); - assertEquals(1.0, m, 5e-3); + // 3 SE for sample mean of ~498k LogNormal(M=1, S=1) draws (ESS ~ N + // for parameter-only chain): SE = sqrt((exp(1)-1) / 498000) ~= 1.86e-3, + // 3 SE ~= 5.6e-3. Expected failure rate at this tolerance: 0.27%. + assertEquals(1.0, m, 6e-3); assertEquals(Math.exp(-0.5), median, 5e-3); assertEquals(Math.exp(1)-1, s, 1e-1); assertEquals(0.0854, StatUtils.percentile(v, 2.5), 5e-3); diff --git a/beast-base/src/test/java/beast/base/spec/inference/parameter/RealScalarParamScalableTest.java b/beast-base/src/test/java/beast/base/spec/inference/parameter/RealScalarParamScalableTest.java new file mode 100644 index 00000000..355033c8 --- /dev/null +++ b/beast-base/src/test/java/beast/base/spec/inference/parameter/RealScalarParamScalableTest.java @@ -0,0 +1,39 @@ +package beast.base.spec.inference.parameter; + +import beast.base.inference.ScalableContractTest; +import beast.base.spec.domain.PositiveReal; +import org.junit.jupiter.api.Test; + +import static org.junit.jupiter.api.Assertions.assertEquals; + +/** + * Tests the {@link beast.base.inference.Scalable} contract on + * {@link RealScalarParam}. + */ +public class RealScalarParamScalableTest { + + @Test + void contractHoldsForPositiveScalar() { + ScalableContractTest.assertContractAcrossScales( + () -> { + RealScalarParam p = new RealScalarParam<>(1.5, PositiveReal.INSTANCE); + p.initAndValidate(); + return p; + }, + (a, b) -> assertEquals(a.get(), b.get(), ScalableContractTest.EPSILON, + "RealScalarParam internal value should match across scale and set+get×s paths") + ); + } + + @Test + void contractHoldsForLargerInitialValue() { + ScalableContractTest.assertContractAcrossScales( + () -> { + RealScalarParam p = new RealScalarParam<>(42.0, PositiveReal.INSTANCE); + p.initAndValidate(); + return p; + }, + (a, b) -> assertEquals(a.get(), b.get(), ScalableContractTest.EPSILON) + ); + } +} diff --git a/beast-base/src/test/java/beast/base/spec/inference/parameter/RealScalarParamTest.java b/beast-base/src/test/java/beast/base/spec/inference/parameter/RealScalarParamTest.java index 734a0583..7495dd6c 100644 --- a/beast-base/src/test/java/beast/base/spec/inference/parameter/RealScalarParamTest.java +++ b/beast-base/src/test/java/beast/base/spec/inference/parameter/RealScalarParamTest.java @@ -106,18 +106,16 @@ void testStoreAndRestore() { } /* - * - Purpose: verify scale(...) and scaleOne(...) multiply the stored value and - * return expected counts. - * - Assertions: scale returns 1 and multiplies value; - * scaleOne multiplies value too. + * - Purpose: verify scale(...) multiplies the stored value and returns the + * log Jacobian. + * - Assertions: scale returns log Jacobian (= log(s) for a 1-D scalar) + * and multiplies value. */ @Test - void testScaleAndScaleOne() { + void testScale() { RealScalarParam param = new RealScalarParam(3.0, Real.INSTANCE); - assertEquals(1, param.scale(2.0)); + assertEquals(Math.log(2.0), param.scale(2.0), 1e-12); assertEquals(6.0, param.get(), 1e-12); - param.scaleOne(0, 0.5); - assertEquals(3.0, param.get(), 1e-12); } /* diff --git a/beast-base/src/test/java/beast/base/spec/inference/parameter/RealVectorParamScalableTest.java b/beast-base/src/test/java/beast/base/spec/inference/parameter/RealVectorParamScalableTest.java new file mode 100644 index 00000000..46d97287 --- /dev/null +++ b/beast-base/src/test/java/beast/base/spec/inference/parameter/RealVectorParamScalableTest.java @@ -0,0 +1,48 @@ +package beast.base.spec.inference.parameter; + +import beast.base.inference.ScalableContractTest; +import beast.base.spec.domain.PositiveReal; +import org.junit.jupiter.api.Test; + +import static org.junit.jupiter.api.Assertions.assertEquals; + +/** + * Tests the {@link beast.base.inference.Scalable} contract on + * {@link RealVectorParam}. + */ +public class RealVectorParamScalableTest { + + @Test + void contractHoldsForPositiveVector() { + ScalableContractTest.assertContractAcrossScales( + () -> { + RealVectorParam p = new RealVectorParam<>( + new double[] { 1.0, 2.5, 3.7, 0.4 }, PositiveReal.INSTANCE); + p.initAndValidate(); + return p; + }, + this::assertSameVectorState + ); + } + + @Test + void contractHoldsForVectorWithMixedMagnitudes() { + ScalableContractTest.assertContractAcrossScales( + () -> { + RealVectorParam p = new RealVectorParam<>( + new double[] { 0.001, 100.0, 5.5 }, PositiveReal.INSTANCE); + p.initAndValidate(); + return p; + }, + this::assertSameVectorState + ); + } + + private void assertSameVectorState(RealVectorParam a, RealVectorParam b) { + assertEquals(a.size(), b.size(), "Vector sizes should match"); + for (int i = 0; i < a.size(); i++) { + assertEquals(a.get(i), b.get(i), ScalableContractTest.EPSILON, + "Element " + i + " should match across scale and set+get×s paths"); + } + } +} diff --git a/beast-base/src/test/java/test/beast/evolution/operator/BactrianScaleOperatorTest.java b/beast-base/src/test/java/test/beast/evolution/operator/BactrianScaleOperatorTest.java index 3718015c..6dafb5e4 100644 --- a/beast-base/src/test/java/test/beast/evolution/operator/BactrianScaleOperatorTest.java +++ b/beast-base/src/test/java/test/beast/evolution/operator/BactrianScaleOperatorTest.java @@ -144,12 +144,17 @@ public void testTreeScaling() { assertEquals(0.5, node[2].getHeight(), EPSILON); assertEquals(0.5, node[3].getHeight(), EPSILON); - // internal nodes, all scaled - // first determine scale factor - double scale = node[4].getHeight() / 1.0; - assertEquals(1.0 * scale, node[4].getHeight(), EPSILON); - assertEquals(1.5 * scale, node[5].getHeight(), EPSILON); - assertEquals(2.0 * scale, node[6].getHeight(), EPSILON); + // Under interval-scaling Tree.scale (CompEvol/beast3#70): + // * node4 (parent of leaves at h=0): margin 1.0 -> 1.0*s, new h = s + // * node5 (parent of leaves at h=0.5): margin 1.0 -> 1.0*s, new h = 0.5 + s + // * node6 (root, children node4=s, node5=0.5+s): + // old margin = 2.0 - 1.5 = 0.5; after recurse min child h = 0.5 + s + // new h = (0.5 + s) + 0.5*s = 0.5 + 1.5*s + // Recover s from node4 (whose taller child is a leaf at h=0). + double scale = node[4].getHeight(); + assertEquals(scale, node[4].getHeight(), EPSILON); + assertEquals(0.5 + scale, node[5].getHeight(), EPSILON); + assertEquals(0.5 + 1.5 * scale, node[6].getHeight(), EPSILON); } diff --git a/beast-base/src/test/java/test/beast/evolution/operator/ScaleOperatorTest.java b/beast-base/src/test/java/test/beast/evolution/operator/ScaleOperatorTest.java index 790e5e24..57c9b095 100644 --- a/beast-base/src/test/java/test/beast/evolution/operator/ScaleOperatorTest.java +++ b/beast-base/src/test/java/test/beast/evolution/operator/ScaleOperatorTest.java @@ -34,11 +34,16 @@ public void testTreeScaling() { assertEquals(0.5, node[2].getHeight(), EPSILON); assertEquals(0.5, node[3].getHeight(), EPSILON); - // internal nodes, all scaled - // first determine scale factor - double scale = node[4].getHeight() / 1.0; - assertEquals(1.0 * scale, node[4].getHeight(), EPSILON); - assertEquals(1.5 * scale, node[5].getHeight(), EPSILON); - assertEquals(2.0 * scale, node[6].getHeight(), EPSILON); + // Under interval-scaling Tree.scale: + // * node4 (parent of leaves at h=0): margin 1.0 -> 1.0*s, new h = s + // * node5 (parent of leaves at h=0.5): margin 1.0 -> 1.0*s, new h = 0.5 + s + // * node6 (root, children node4=s, node5=0.5+s): + // old margin = 2.0 - 1.5 = 0.5; after recurse min child h = 0.5 + s + // new h = (0.5 + s) + 0.5*s = 0.5 + 1.5*s + // Recover s from node4 (whose taller child is a leaf at h=0). + double scale = node[4].getHeight(); + assertEquals(scale, node[4].getHeight(), EPSILON); + assertEquals(0.5 + scale, node[5].getHeight(), EPSILON); + assertEquals(0.5 + 1.5 * scale, node[6].getHeight(), EPSILON); } } diff --git a/beast-base/src/test/java/test/beast/evolution/tree/TreeTest.java b/beast-base/src/test/java/test/beast/evolution/tree/TreeTest.java index 6cf47cd6..97c4c54a 100644 --- a/beast-base/src/test/java/test/beast/evolution/tree/TreeTest.java +++ b/beast-base/src/test/java/test/beast/evolution/tree/TreeTest.java @@ -27,18 +27,22 @@ public void testTreeScaling() { assertEquals(2.0, node[6].getHeight(), EPSILON); treeParser.scale(2.0); - + // leaf node node = treeParser.getNodesAsArray(); assertEquals(0.0, node[0].getHeight(), EPSILON); assertEquals(0.0, node[1].getHeight(), EPSILON); - // leaf node, not scaled + // leaves are preserved (not scaled) under interval scaling assertEquals(0.5, node[2].getHeight(), EPSILON); assertEquals(0.5, node[3].getHeight(), EPSILON); - // internal nodes, all scaled + // Under interval-scaling Tree.scale(2.0): + // * node4: margin 1.0 -> 2.0; new h = 0 + 2.0 = 2.0 + // * node5: margin 1.0 -> 2.0; new h = 0.5 + 2.0 = 2.5 + // * node6: old margin 0.5 -> 1.0; new min child h = max(2.0, 2.5) = 2.5 + // new h = 2.5 + 1.0 = 3.5 assertEquals(2.0, node[4].getHeight(), EPSILON); - assertEquals(3.0, node[5].getHeight(), EPSILON); - assertEquals(4.0, node[6].getHeight(), EPSILON); - + assertEquals(2.5, node[5].getHeight(), EPSILON); + assertEquals(3.5, node[6].getHeight(), EPSILON); + } } diff --git a/scripts/scalable-contract-talk.pdf b/scripts/scalable-contract-talk.pdf new file mode 100644 index 00000000..8a3b3825 Binary files /dev/null and b/scripts/scalable-contract-talk.pdf differ diff --git a/scripts/scalable-contract-talk.tex b/scripts/scalable-contract-talk.tex new file mode 100644 index 00000000..76050327 --- /dev/null +++ b/scripts/scalable-contract-talk.tex @@ -0,0 +1,356 @@ +% !TEX program = pdflatex +\documentclass[aspectratio=169,11pt]{beamer} + +\usetheme{Madrid} +\usecolortheme{seahorse} +\setbeamertemplate{navigation symbols}{} +\setbeamertemplate{footline}[frame number] + +\usepackage[T1]{fontenc} +\usepackage{lmodern} +\usepackage{listings} +\usepackage{booktabs} +\usepackage{tikz} +\usetikzlibrary{arrows.meta,positioning} + +\lstset{ + basicstyle=\ttfamily\scriptsize, + keywordstyle=\color{blue!70!black}\bfseries, + commentstyle=\color{green!50!black}, + stringstyle=\color{red!60!black}, + breaklines=true, + showstringspaces=false, + columns=fullflexible, + frame=single, + backgroundcolor=\color{gray!8}, + rulecolor=\color{gray!40}, + xleftmargin=2pt, + xrightmargin=2pt, + language=Java, +} + +\definecolor{beastblue}{RGB}{40,80,140} +\definecolor{beastorange}{RGB}{200,100,20} +\setbeamercolor{frametitle}{fg=beastblue} +\setbeamercolor{title}{fg=beastblue} +\setbeamercolor{structure}{fg=beastblue} + +\newcommand{\pkg}[1]{\texttt{\color{beastorange}#1}} +\newcommand{\code}[1]{\texttt{#1}} + +\title{The Scalable Contract} +\subtitle{Resolving beast3 issue \#20: a binding three-method contract on \texttt{Scalable}} +\author{Alexei Drummond} +\date{5 May 2026} + +\begin{document} + +% ============================================================ +\begin{frame} +\titlepage +\end{frame} + +% ============================================================ +\begin{frame}{The blocker} +Joelle Barido-Sottani has been blocked porting custom-tree packages to beast3 since 26 March. + +\vspace{1em} + +Issue \href{https://github.com/CompEvol/beast3/issues/20}{CompEvol/beast3\#20}, opened by Remco in October 2025. + +\vspace{1em} + +Two distinct problems she ran into: + +\begin{enumerate} +\item The old \code{Scalable} interface returns degrees of freedom (\code{int}), which encodes a specific Hastings-ratio shape (\code{dof $\times$ log(s)}). +For her custom trees the HR isn't of that shape. + +\item The new spec \code{UpDownOperator} bypasses \code{Tree.scale()} entirely and runs its own recursion. +Even if she fixed her tree's \code{scale()}, UpDown wouldn't call it. +\end{enumerate} +\end{frame} + +% ============================================================ +\begin{frame}[fragile]{The old contract} +\begin{lstlisting} +public interface Scalable { + int scale(double s); // returns dof + void scaleOne(int i, double s); + default double scaleAll(double s) { + return scale(s) * Math.log(s); + } +} +\end{lstlisting} + +\vspace{0.5em} + +\textbf{What's wrong:} + +\begin{itemize} +\item \code{int dof} is a stand-in for one specific HR formula. +Operators reconstruct \code{dof $\times$ log(s)} as the move's Jacobian contribution. +\item Anything that scales differently (interval scaling, custom trees, IntervalScaleOperator) has to bypass the interface and compute HR inline. +\item The interface gives no way to ask ``what does this Scalable use as its dilation axis?'' +\end{itemize} +\end{frame} + +% ============================================================ +\begin{frame}{What does \texttt{scale} actually return?} +The Metropolis-Hastings acceptance ratio has three parts: + +\begin{enumerate} +\item \textbf{Target ratio}: $\pi(x') / \pi(x)$. +\item \textbf{Proposal ratio} (the Hastings ratio proper): $q(x \mid x') / q(x' \mid x)$. +\item \textbf{Jacobian determinant}: $|\det(\partial x' / \partial x)|$ from the change of variables. +\end{enumerate} + +\vspace{1em} + +\code{Scalable.scale(s)} can only sensibly return the \textbf{Jacobian} of the deterministic part of the move. + +\vspace{0.5em} + +The proposal ratio depends on how the operator drew \code{s} (uniform-on-log, Bactrian, Gaussian, ...) and belongs to the operator's kernel. + +\vspace{0.5em} + +\textbf{Joelle's option (1)} of returning ``the HR'' is right in spirit; the precise refinement is to return the log Jacobian determinant. +\end{frame} + +% ============================================================ +\begin{frame}[fragile]{The new contract} +\begin{lstlisting} +public interface Scalable { + double scale(double s); // returns log Jacobian + double getScalableValue(); // current pos. on dilation axis + default double setScalableValue(double V) { + return scale(V / getScalableValue()); + } +} +\end{lstlisting} + +\vspace{0.5em} + +Three methods, mutually consistent on a single dilation axis. +\end{frame} + +% ============================================================ +\begin{frame}[fragile]{The three invariants} +For any valid Scalable \code{x} and positive \code{s}: + +\begin{enumerate} +\item \textbf{Scale-equivariance.} After \code{scale(s)}, \code{getScalableValue()} returns $s \times$ its previous value. +\begin{lstlisting} +double v0 = x.getScalableValue(); +x.scale(s); +assert x.getScalableValue() == s * v0; +\end{lstlisting} + +\item \textbf{Set is a fixed point of get.} +\begin{lstlisting} +x.setScalableValue(V); +assert x.getScalableValue() == V; +\end{lstlisting} + +\item \textbf{Set composes with scale.} +\begin{lstlisting} +x.setScalableValue(x.getScalableValue() * s) == x.scale(s); +\end{lstlisting} +\end{enumerate} + +\vspace{0.5em} + +The default \code{setScalableValue} expresses (3) directly. +\end{frame} + +% ============================================================ +\begin{frame}{What is the dilation axis?} +\code{getScalableValue} must return whatever quantity \emph{actually} scales by exactly \code{s} under the implementer's \code{scale}. + +\vspace{1em} + +\begin{tabular}{ll} +\toprule +\textbf{Implementer} & \textbf{Dilation summary} \\ +\midrule +\code{RealScalarParam} & the value itself \\ +\code{RealVectorParam} & sum of values \\ +\code{Tree} (interval scaling) & sum of margins \\ +\code{Tree} (affine height scaling) & root height \\ +Joelle's custom tree & whatever her scale preserves \\ +\bottomrule +\end{tabular} + +\vspace{1em} + +\textbf{The contract is enforced per-implementation.} +A type that violates the invariants is not a valid Scalable. +\end{frame} + +% ============================================================ +\begin{frame}{Why interval scaling for \texttt{Tree.scale}?} +The old \code{Tree.scale(s)} multiplied internal heights by \code{s} (affine). +For a heterochronous tree with leaves at \code{h\_c > 0} and parent at \code{h\_N}, +the move \textbf{throws} when \code{s < h\_c / h\_N}. + +\vspace{1em} + +This is not an edge case. +Most BEAST analyses with sampling dates run on heterochronous trees. + +\vspace{1em} + +\textbf{Interval scaling}: multiply each internal node's margin (height above its taller child) by \code{s}. + +\begin{itemize} +\item Tip dates preserved by construction. +\item Move is valid for any positive \code{s}; never throws. +\item Sum of margins is exactly \code{s}-equivariant: matches \code{getScalableValue}. +\item For binary trees, the dof count equals the number of margins, +so operator HRs are unchanged. +\end{itemize} +\end{frame} + +% ============================================================ +\begin{frame}{Affine semantics retained where needed} +Some callers genuinely want ``set my root to exactly this height'' with affine semantics. + +\vspace{1em} + +\code{Tree.scaleToRootHeight(double targetHeight)} preserves the old affine behaviour as a separate helper, outside the contract. + +\begin{itemize} +\item Multiplies internal heights by \code{targetHeight / oldRoot}. +\item May throw \code{IllegalArgumentException} for invalid configurations. +\item Used by \code{StarBeastStartState} init code. +\item Not part of the Scalable contract. +\end{itemize} + +\vspace{1em} + +Two methods, two purposes: \code{scale} for moves, \code{scaleToRootHeight} for setup. +\end{frame} + +% ============================================================ +\begin{frame}{AMVN now reaches its proposed value} +Before: AMVN proposed a target root height $V$, called \code{tree.scale(V/oldRoot)}. + +\vspace{0.5em} + +For non-ultrametric trees under affine scaling, the move sometimes \textbf{threw}. +For non-ultrametric trees under interval scaling, the root would end up at some height other than $V$. +Either way, AMVN's empirical covariance bookkeeping was wrong. + +\vspace{1em} + +After: AMVN parameterises trees by sum-of-margins via \code{getScalableValue} / \code{setScalableValue}. +\code{setScalableValue(V)} reduces to \code{scale(V / sum\_margins)}. + +\begin{itemize} +\item Sum of margins equals $V$ exactly, by the contract. +\item AMVN's covariance is now updated against the realised state. +\item Side benefit: fixes a latent missing \texttt{return} in spec AMVN's \code{getValue} for trees. +\end{itemize} +\end{frame} + +% ============================================================ +\begin{frame}{Test plan} +Three layers of testing in \href{https://github.com/CompEvol/beast3/pull/70}{PR \#70}: + +\begin{itemize} +\item \textbf{ScalableContractTest harness}: a reusable helper that asserts the three invariants on any Scalable. Per-class tests for \code{RealScalarParam}, \code{RealVectorParam}, \code{Tree} on ultrametric, heterochronous, and leaf-intruding fixtures. + +\item \textbf{408 existing beast-base unit tests} all pass. + +\item \textbf{Heterochronous MCMC integration}: 10 \code{TipTimeTest} chains under \texttt{-Pslow-tests}. +Real chains on tip-dated data converge to BEAST1 reference values within tolerance. +\textit{Tree height: $15147 \pm 73$, expected $15000 \pm 70$, ESS $\sim$22000.} +\end{itemize} + +\vspace{1em} + +Statistical-tolerance fix on two slow tests (\code{UpDownOperatorTest}, \code{RealRandomWalkOperatorTest}): tolerances widened to 3 SE using ESS rather than nominal sample count, so they fail less often. +\end{frame} + +% ============================================================ +\begin{frame}{Downstream impact} +Binding API change. +Any package that implements \code{Scalable} or stores the int return needs to update. + +\vspace{1em} + +\begin{tabular}{lll} +\toprule +\textbf{Package} & \textbf{Impact} & \textbf{Status} \\ +\midrule +Mascot & none (no \code{Scalable}) & compiles clean \\ +beast-classic & none & compiles clean \\ +LPhyBeast core & none & compiles clean \\ +BEASTLabs & 2 files & migrated, all 83 tests pass \\ +\bottomrule +\end{tabular} + +\vspace{1em} + +BEASTLabs migration in parallel \href{https://github.com/BEAST2-Dev/BEASTLabs/pull/29}{PR \#29}: + +\begin{itemize} +\item \code{PrevalenceList.scale} returns log Jacobian; new \code{getScalableValue} returns sum of node times. +\item \code{TreeScaleOperator} consumes log Jacobian directly; HR factored as \code{treeLogJ $-$ 2 log(s)}. +\end{itemize} +\end{frame} + +% ============================================================ +\begin{frame}{Out of scope (deferred)} +Operator-design questions, separable from the interface contract: + +\begin{itemize} +\item Spec \code{UpDownOperator}'s \code{actualScaler = lengthAfter / lengthBefore} logic for tree + up/down combinations. +\item Whether \code{IntervalScaleOperator} should collapse into \code{UpDownOperator} now that \code{Tree.scale} does interval scaling. +\item AMVN's covariance behaviour for non-ultrametric trees beyond the proposal-target fix. +\item Whether tree-aware operators should use the kernel \code{s} or an effective tree-dilation factor for non-tree up/down parameters. +\end{itemize} + +\vspace{1em} + +The interface contract enables these to be tackled separately on their own merits. +\end{frame} + +% ============================================================ +\begin{frame}{Timeline} +\begin{itemize} +\item \textbf{2 May}: \href{https://github.com/CompEvol/beast3/pull/70}{PR \#70} and \href{https://github.com/BEAST2-Dev/BEASTLabs/pull/29}{BEASTLabs \#29} ready for review. +\item \textbf{Awaiting}: Remco's review (he owns beast-base). +\item \textbf{Proposed at next stakeholders meeting}: this is the last substantive API change to the beast core before an official release. +Aim for release on the BEAST website mid-June. +\end{itemize} +\end{frame} + +% ============================================================ +\begin{frame}{Summary} +\begin{block}{The contract} +\code{scale(s)} returns the log Jacobian determinant of a move along a single dilation axis. +\code{getScalableValue} reads position on that axis. +\code{setScalableValue(V) = scale(V / getScalableValue())}. +\end{block} + +\vspace{0.5em} + +\begin{block}{The Tree side} +\code{Tree.scale} does interval scaling and never throws. +Sum of margins is the dilation summary. +\code{Tree.scaleToRootHeight} preserves affine semantics where they were genuinely wanted. +\end{block} + +\vspace{0.5em} + +\begin{block}{The result} +Joelle's custom trees can implement the contract directly. +AMVN works correctly on non-ultrametric trees. +408 + 9 + 10 tests pass. +Mid-June beast3 release on track. +\end{block} +\end{frame} + +\end{document}