diff --git a/src/main/java/coupledMCMC/CoupledMCMC.java b/src/main/java/coupledMCMC/CoupledMCMC.java index 659e672..ed182f3 100644 --- a/src/main/java/coupledMCMC/CoupledMCMC.java +++ b/src/main/java/coupledMCMC/CoupledMCMC.java @@ -774,24 +774,24 @@ void swapOperatorTuning(HeatedChain mcmc1, HeatedChain mcmc2) { for (int i = 0; i < operatorList1.size(); i++){ - Operator operator1 = operatorList1.get(i); - Operator operator2 = operatorList2.get(i); - - int m_nNrRejected = operator1.get_m_nNrRejected(); - int m_nNrAccepted = operator1.get_m_nNrAccepted(); - int m_nNrRejectedForCorrection = operator1.get_m_nNrRejectedForCorrection(); - int m_nNrAcceptedForCorrection = operator1.get_m_nNrAcceptedForCorrection(); - - operator1.setAcceptedRejected(operator2.get_m_nNrRejected(), operator2.get_m_nNrAccepted(), operator2.get_m_nNrRejectedForCorrection(), operator2.get_m_nNrAcceptedForCorrection()); - operator2.setAcceptedRejected(m_nNrAccepted, m_nNrRejected, m_nNrAcceptedForCorrection, m_nNrRejectedForCorrection); - - - double coercableParameterValue = operator1.getCoercableParameterValue(); - operator1.setCoercableParameterValue(operator2.getCoercableParameterValue()); - operator2.setCoercableParameterValue(coercableParameterValue); - + swapTuning(operatorList1.get(i), operatorList2.get(i)); } } + + // operators tune themselves from their acceptance rate, so the counts must travel with the tuned value + static void swapTuning(Operator operator1, Operator operator2) { + int accepted = operator1.get_m_nNrAccepted(); + int rejected = operator1.get_m_nNrRejected(); + int acceptedForCorrection = operator1.get_m_nNrAcceptedForCorrection(); + int rejectedForCorrection = operator1.get_m_nNrRejectedForCorrection(); + operator1.setAcceptedRejected(operator2.get_m_nNrAccepted(), operator2.get_m_nNrRejected(), + operator2.get_m_nNrAcceptedForCorrection(), operator2.get_m_nNrRejectedForCorrection()); + operator2.setAcceptedRejected(accepted, rejected, acceptedForCorrection, rejectedForCorrection); + + double coercableParameterValue = operator1.getCoercableParameterValue(); + operator1.setCoercableParameterValue(operator2.getCoercableParameterValue()); + operator2.setCoercableParameterValue(coercableParameterValue); + } /* makes the schedule of when to swap which chains */ @SuppressWarnings("unchecked") diff --git a/src/test/java/coupledMCMC/SwapOperatorTuningTest.java b/src/test/java/coupledMCMC/SwapOperatorTuningTest.java new file mode 100644 index 0000000..cac6b08 --- /dev/null +++ b/src/test/java/coupledMCMC/SwapOperatorTuningTest.java @@ -0,0 +1,49 @@ +package coupledMCMC; + +import beast.base.spec.domain.Real; +import beast.base.spec.inference.operator.RealRandomWalkOperator; +import beast.base.spec.inference.parameter.RealScalarParam; +import org.junit.jupiter.api.Test; + +import java.util.List; + +import static org.junit.jupiter.api.Assertions.assertEquals; + +public class SwapOperatorTuningTest { + + private static RealRandomWalkOperator operator(double windowSize) { + RealScalarParam x = new RealScalarParam<>(0.0, Real.INSTANCE); + RealRandomWalkOperator op = new RealRandomWalkOperator(); + op.initByName("scalar", x, "windowSize", windowSize, "weight", 1.0); + return op; + } + + private static List counts(RealRandomWalkOperator op) { + return List.of(op.get_m_nNrAccepted(), op.get_m_nNrRejected(), + op.get_m_nNrAcceptedForCorrection(), op.get_m_nNrRejectedForCorrection()); + } + + @Test + void swapTuningExchangesAcceptAndRejectCountsWithoutMixingThemUp() { + RealRandomWalkOperator a = operator(1.0); + RealRandomWalkOperator b = operator(2.0); + a.setAcceptedRejected(3, 7, 2, 6); + b.setAcceptedRejected(5, 1, 4, 0); + + CoupledMCMC.swapTuning(a, b); + + assertEquals(List.of(5, 1, 4, 0), counts(a)); + assertEquals(List.of(3, 7, 2, 6), counts(b)); + } + + @Test + void swapTuningExchangesTheTunedParameterToo() { + RealRandomWalkOperator a = operator(1.0); + RealRandomWalkOperator b = operator(2.0); + + CoupledMCMC.swapTuning(a, b); + + assertEquals(2.0, a.getCoercableParameterValue()); + assertEquals(1.0, b.getCoercableParameterValue()); + } +}