use of dr.inference.mcmc.MCMCCriterion in project beast-mcmc by beast-dev.
the class MCMCMC method swapChainTemperatures.
private int swapChainTemperatures() {
if (DEBUG) {
System.out.print("Current scores: ");
for (int i = 0; i < chains.length; i++) {
System.out.print("\t");
if (i == coldChain) {
System.out.print("[");
}
System.out.print(chains[i].getCurrentScore());
if (i == coldChain) {
System.out.print("]");
}
}
System.out.println();
}
int newColdChain = coldChain;
int index1 = MathUtils.nextInt(chains.length);
int index2 = MathUtils.nextInt(chains.length);
while (index1 == index2) {
index2 = MathUtils.nextInt(chains.length);
}
double score1 = chains[index1].getCurrentScore();
MCMCCriterion acceptor1 = ((MCMCCriterion) chains[index1].getAcceptor());
double temperature1 = acceptor1.getTemperature();
double score2 = chains[index2].getCurrentScore();
MCMCCriterion acceptor2 = ((MCMCCriterion) chains[index2].getAcceptor());
double temperature2 = acceptor2.getTemperature();
double logRatio = ((score2 - score1) * temperature1) + ((score1 - score2) * temperature2);
boolean swap = (Math.log(MathUtils.nextDouble()) < logRatio);
if (swap) {
if (DEBUG) {
System.out.println("Swapping chain " + index1 + " and chain " + index2);
}
acceptor1.setTemperature(temperature2);
acceptor2.setTemperature(temperature1);
OperatorSchedule schedule1 = schedules[index1];
OperatorSchedule schedule2 = schedules[index2];
for (int i = 0; i < schedule1.getOperatorCount(); i++) {
MCMCOperator operator1 = schedule1.getOperator(i);
MCMCOperator operator2 = schedule2.getOperator(i);
long tmp = operator1.getAcceptCount();
operator1.setAcceptCount(operator2.getAcceptCount());
operator2.setAcceptCount(tmp);
tmp = operator1.getRejectCount();
operator1.setRejectCount(operator2.getRejectCount());
operator2.setRejectCount(tmp);
double tmp2 = operator1.getSumDeviation();
operator1.setSumDeviation(operator2.getSumDeviation());
operator2.setSumDeviation(tmp2);
if (operator1 instanceof CoercableMCMCOperator) {
tmp2 = ((CoercableMCMCOperator) operator1).getCoercableParameter();
((CoercableMCMCOperator) operator1).setCoercableParameter(((CoercableMCMCOperator) operator2).getCoercableParameter());
((CoercableMCMCOperator) operator2).setCoercableParameter(tmp2);
}
}
if (index1 == coldChain) {
newColdChain = index2;
} else if (index2 == coldChain) {
newColdChain = index1;
}
}
return newColdChain;
}
Aggregations