From fe31d4bbe80b1bfee1582ff4d7598df4789fbf46 Mon Sep 17 00:00:00 2001 From: GuyBaele Date: Sat, 21 Dec 2024 11:35:44 +0100 Subject: [PATCH 1/6] GMRF transition kernels are not Gibbs operators --- src/dr/app/beauti/types/OperatorType.java | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/dr/app/beauti/types/OperatorType.java b/src/dr/app/beauti/types/OperatorType.java index 6f0e47e7d2..d89f8729fd 100644 --- a/src/dr/app/beauti/types/OperatorType.java +++ b/src/dr/app/beauti/types/OperatorType.java @@ -28,7 +28,6 @@ package dr.app.beauti.types; import dr.evomodel.operators.BitFlipInSubstitutionModelOperator; -import dr.evomodelxml.operators.TreeNodeSlideParser; import dr.inference.operators.RateBitExchangeOperator; import dr.inferencexml.operators.ScaleOperatorParser; @@ -66,8 +65,8 @@ public enum OperatorType { NARROW_EXCHANGE("narrowExchange"), WIDE_EXCHANGE("wideExchange"), EMPIRICAL_TREE_SWAP("empiricalSwap"), - GMRF_GIBBS_OPERATOR("gmrfGibbsOperator"), - SKY_GRID_GIBBS_OPERATOR("gmrfGibbsOperator"), + GMRF_BLOCKUPDATE_OPERATOR("gmrfBlockUpdateOperator"), + SKY_GRID_BLOCKUPDATE_OPERATOR("gmrfBlockUpdateOperator"), SKY_GRID_HMC_OPERATOR("gmrfHMCOperator"), // PRECISION_GMRF_OPERATOR("precisionGMRFOperator"), WILSON_BALDING("wilsonBalding"), From b510b75a347b0e60c985d06bbe24bd9cc3bd322e Mon Sep 17 00:00:00 2001 From: GuyBaele Date: Sat, 21 Dec 2024 11:36:06 +0100 Subject: [PATCH 2/6] set tuning to 1 for GMRF transition kernels --- src/dr/app/beauti/options/PartitionTreePrior.java | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/dr/app/beauti/options/PartitionTreePrior.java b/src/dr/app/beauti/options/PartitionTreePrior.java index 3da2ee6c1e..3685518ef8 100644 --- a/src/dr/app/beauti/options/PartitionTreePrior.java +++ b/src/dr/app/beauti/options/PartitionTreePrior.java @@ -28,7 +28,6 @@ package dr.app.beauti.options; import dr.app.beauti.types.*; -import dr.evomodel.coalescent.VariableDemographicModel; import dr.evomodel.speciation.CalibrationPoints; import dr.evomodelxml.coalescent.GMRFSkyrideLikelihoodParser; import dr.evomodelxml.speciation.BirthDeathEpidemiologyModelParser; @@ -271,9 +270,9 @@ public void initModelParametersAndOpererators() { // "demographic.indicators", OperatorType.SCALE_WITH_INDICATORS, 0.5, 2 * demoWeights); createOperatorUsing2Parameters("gmrfGibbsOperator", "gmrfGibbsOperator", "Gibbs sampler for GMRF Skyride", "skyride.logPopSize", - "skyride.precision", OperatorType.GMRF_GIBBS_OPERATOR, -1, 2); + "skyride.precision", OperatorType.GMRF_BLOCKUPDATE_OPERATOR, 1, 2); createOperatorUsing2Parameters("gmrfSkyGridGibbsOperator", "skygrid.logPopSize", "Gibbs sampler for Bayesian SkyGrid", "skygrid.logPopSize", - GMRFSkyrideLikelihoodParser.SKYGRID_PRECISION, OperatorType.SKY_GRID_GIBBS_OPERATOR, -1, 2); + GMRFSkyrideLikelihoodParser.SKYGRID_PRECISION, OperatorType.SKY_GRID_BLOCKUPDATE_OPERATOR, 1, 2); createScaleOperator(GMRFSkyrideLikelihoodParser.SKYGRID_PRECISION, "skygrid precision", 0.75, 1.0); createOperatorUsing2Parameters("gmrfSkyGridHMCOperator", "Multiple", "HMC transition kernel for Bayesian SkyGrid", "skygrid.logPopSize", GMRFSkyrideLikelihoodParser.SKYGRID_PRECISION, OperatorType.SKY_GRID_HMC_OPERATOR, -1, 2); @@ -292,7 +291,7 @@ public void initModelParametersAndOpererators() { createOperator(BirthDeathSerialSamplingModelParser.BDSS + "." + BirthDeathSerialSamplingModelParser.RELATIVE_MU, OperatorType.RANDOM_WALK_LOGIT, demoTuning, 1); createScaleOperator(BirthDeathSerialSamplingModelParser.BDSS + "." - + BirthDeathSerialSamplingModelParser.PSI, demoTuning, 1); // todo random worl op ? + + BirthDeathSerialSamplingModelParser.PSI, demoTuning, 1); // todo random walk op ? createScaleOperator(BirthDeathSerialSamplingModelParser.BDSS + "." + BirthDeathSerialSamplingModelParser.ORIGIN, demoTuning, 1); // createScaleOperator(BirthDeathSerialSamplingModelParser.BDSS + "." From fd162f367c7e25a79a8d3a203a1ed0034a9fdb3b Mon Sep 17 00:00:00 2001 From: GuyBaele Date: Sat, 21 Dec 2024 11:36:53 +0100 Subject: [PATCH 3/6] write scaleFactor to XML for GMRF transition kernels --- .../app/beauti/generator/OperatorsGenerator.java | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/dr/app/beauti/generator/OperatorsGenerator.java b/src/dr/app/beauti/generator/OperatorsGenerator.java index 4234751be5..ede15bb1b4 100644 --- a/src/dr/app/beauti/generator/OperatorsGenerator.java +++ b/src/dr/app/beauti/generator/OperatorsGenerator.java @@ -231,11 +231,11 @@ private void writeOperator(Operator operator, XMLWriter writer) { case SCALE_WITH_INDICATORS: writeScaleWithIndicatorsOperator(operator, writer); break; - case GMRF_GIBBS_OPERATOR: - writeGMRFGibbsOperator(operator, prefix, writer); + case GMRF_BLOCKUPDATE_OPERATOR: + writeGMRFBlockUpdateOperator(operator, prefix, writer); break; - case SKY_GRID_GIBBS_OPERATOR: - writeSkyGridGibbsOperator(operator, prefix, writer); + case SKY_GRID_BLOCKUPDATE_OPERATOR: + writeSkyGridBlockUpdateOperator(operator, prefix, writer); break; case SKY_GRID_HMC_OPERATOR: writeSkyGridHMCOperator(operator, prefix, writer); @@ -518,12 +518,11 @@ private void writeSampleNonActiveOperator(Operator operator, XMLWriter writer) { writer.writeCloseTag(SampleNonActiveGibbsOperatorParser.SAMPLE_NONACTIVE_GIBBS_OPERATOR); } - private void writeSkyGridGibbsOperator(Operator operator, String treePriorPrefix, XMLWriter writer) { + private void writeSkyGridBlockUpdateOperator(Operator operator, String treePriorPrefix, XMLWriter writer) { writer.writeOpenTag( GMRFSkyrideBlockUpdateOperatorParser.GRID_BLOCK_UPDATE_OPERATOR, new Attribute[] { -// This is a Gibbs operator so shouldn't have a tuning parameter? -// new Attribute.Default(GMRFSkyrideBlockUpdateOperatorParser.SCALE_FACTOR, operator.getTuning()), + new Attribute.Default(GMRFSkyrideBlockUpdateOperatorParser.SCALE_FACTOR, operator.getTuning()), getWeightAttribute(operator.getWeight()) } ); @@ -688,7 +687,7 @@ private void writeShrinkageClockHMCOperator(Operator operator, String prefix, XM writer.writeCloseTag(HamiltonianMonteCarloOperatorParser.HMC_OPERATOR); } - private void writeGMRFGibbsOperator(Operator operator, String treePriorPrefix, XMLWriter writer) { + private void writeGMRFBlockUpdateOperator(Operator operator, String treePriorPrefix, XMLWriter writer) { writer.writeOpenTag( GMRFSkyrideBlockUpdateOperatorParser.BLOCK_UPDATE_OPERATOR, new Attribute[]{ From 41a001c5381ead17c7151a4cdfd0a9dc1709f878 Mon Sep 17 00:00:00 2001 From: GuyBaele Date: Tue, 24 Dec 2024 12:04:01 +0100 Subject: [PATCH 4/6] location scale gradient generator code was in the wrong place --- src/dr/app/beauti/generator/ClockModelGenerator.java | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/dr/app/beauti/generator/ClockModelGenerator.java b/src/dr/app/beauti/generator/ClockModelGenerator.java index c19446180c..2c2dc44e07 100644 --- a/src/dr/app/beauti/generator/ClockModelGenerator.java +++ b/src/dr/app/beauti/generator/ClockModelGenerator.java @@ -340,6 +340,9 @@ public void writeBranchRatesModel(PartitionClockModel clockModel, XMLWriter writ writer.writeIDref(DefaultTreeModel.TREE_MODEL, treePrefix + DefaultTreeModel.TREE_MODEL); writer.writeCloseTag(CTMCScalePriorParser.MODEL_NAME); + } + + if (generateScaleGradient){ //location gradient writer.writeOpenTag(LocationScaleGradientParser.NAME, new Attribute[]{ new Attribute.Default<>(XMLParser.ID, prefix + LocationGradient.LOCATION_GRADIENT), @@ -352,9 +355,6 @@ public void writeBranchRatesModel(PartitionClockModel clockModel, XMLWriter writ writer.writeCloseTag(LocationScaleGradientParser.LOCATION); writer.writeCloseTag(LocationScaleGradientParser.NAME); - } - - if (generateScaleGradient){ //scale gradient writer.writeOpenTag(LocationScaleGradientParser.NAME, new Attribute[]{ new Attribute.Default<>(XMLParser.ID, prefix + ScaleGradient.SCALE_GRADIENT), From fb0d1ded52c150537898045caee6e3e3ca70a829 Mon Sep 17 00:00:00 2001 From: Plemey Date: Sat, 28 Dec 2024 15:15:15 +0100 Subject: [PATCH 5/6] stat for Simon --- .../ContinuousDiffusionStatistic.java | 40 ++++++++++++++++++- 1 file changed, 39 insertions(+), 1 deletion(-) diff --git a/src/dr/evomodel/continuous/ContinuousDiffusionStatistic.java b/src/dr/evomodel/continuous/ContinuousDiffusionStatistic.java index f72aeac688..4fa6325c2f 100644 --- a/src/dr/evomodel/continuous/ContinuousDiffusionStatistic.java +++ b/src/dr/evomodel/continuous/ContinuousDiffusionStatistic.java @@ -74,6 +74,7 @@ public class ContinuousDiffusionStatistic extends Statistic.Abstract { public static final String SPEARMAN = "spearman"; public static final String CORRELATION_COEFFICIENT = "correlationCoefficient"; public static final String DISTANCE_TIME_CORRELATION = "distanceTimeCorrelation"; + public static final String SQUAREDDISTANCE_TIME4_CORRELATION = "squaredDistanceTimeFourCorrelation"; public static final String R_SQUARED = "Rsquared"; public static final String STATISTIC = "statistic"; public static final String TRAIT = "trait"; @@ -458,6 +459,18 @@ public double getStatisticValue(int dim) { Regression r = new Regression(convertDoubles(times),convertDoubles(distances)); return r.getCorrelationCoefficient(); } + } else if (summaryStat == summaryStatistic.SQUAREDDISTANCE_TIME4_CORRELATION) { + List squareddistances = squareElements(distances); + List timesFour = elementsTimesFour(times); + if (summaryMode == Mode.SPEARMAN) { + return getSpearmanRho(convertDoubles(timesFour),convertDoubles(squareddistances)); + } else if (summaryMode == Mode.R_SQUARED) { + Regression r = new Regression(convertDoubles(timesFour), convertDoubles(squareddistances)); + return r.getRSquared(); + } else { + Regression r = new Regression(convertDoubles(timesFour),convertDoubles(squareddistances)); + return r.getCorrelationCoefficient(); + } } else { return treeLength; } @@ -490,6 +503,22 @@ private double[] toArray(List list) { return returnArray; } + private static List squareElements(List inputList) { + List squaredList = new ArrayList<>(); + for (Double number : inputList) { + squaredList.add(number * number); + } + return squaredList; + } + + private static List elementsTimesFour (List inputList) { + List returnList = new ArrayList<>(); + for (Double number : inputList) { + returnList.add(number * 4); + } + return returnList; + } + private double[] imputeValue(double[] nodeValue, double[] parentValue, double time, double nodeHeight, double parentHeight, double[] precisionArray, double rate, boolean trueNoise) { final double scaledTimeChild = (time - nodeHeight) * rate; @@ -932,7 +961,8 @@ enum summaryStatistic { WAVEFRONT_DISTANCE, WAVEFRONT_DISTANCE_PHYLO, WAVEFRONT_RATE, - DISTANCE_TIME_CORRELATION + DISTANCE_TIME_CORRELATION, + SQUAREDDISTANCE_TIME4_CORRELATION } enum BranchSet { @@ -1023,6 +1053,12 @@ public Object parseXMLObject(XMLObject xo) throws XMLParseException { System.err.println(name+": mode = "+mode+" ignored for "+DISTANCE_TIME_CORRELATION+", reverting to correlation coefficient mode"); statMode = Mode.CORRELATION_COEFFICIENT; } + } else if (statistic.equals(SQUAREDDISTANCE_TIME4_CORRELATION)) { + summaryStat = summaryStatistic.SQUAREDDISTANCE_TIME4_CORRELATION; + if (mode.equals(AVERAGE) || mode.equals(WEIGHTED_AVERAGE) || mode.equals(COEFFICIENT_OF_VARIATION) || mode.equals(MEDIAN)){ + System.err.println(name+": mode = "+mode+" ignored for "+SQUAREDDISTANCE_TIME4_CORRELATION+", reverting to correlation coefficient mode"); + statMode = Mode.CORRELATION_COEFFICIENT; + } } else if (statistic.equals(WAVEFRONT_DISTANCE)) { summaryStat = summaryStatistic.WAVEFRONT_DISTANCE; if (!mode.equals(WEIGHTED_AVERAGE)) { @@ -1065,6 +1101,8 @@ public Object parseXMLObject(XMLObject xo) throws XMLParseException { summaryStat = summaryStatistic.DIFFUSION_COEFFICIENT; } else if (statistic.equals(DISTANCE_TIME_CORRELATION)) { summaryStat = summaryStatistic.DISTANCE_TIME_CORRELATION; + } else if (statistic.equals(SQUAREDDISTANCE_TIME4_CORRELATION)) { + summaryStat = summaryStatistic.SQUAREDDISTANCE_TIME4_CORRELATION; } else { System.err.println(name+": unknown statistic: "+statistic+". Reverting to diffusion rate."); summaryStat = summaryStatistic.DIFFUSION_RATE; From f37ecb6a49f0be86da0fcce5cd4a576ea1809624 Mon Sep 17 00:00:00 2001 From: Plemey Date: Mon, 30 Dec 2024 19:11:44 +0100 Subject: [PATCH 6/6] simplifying stat for Simon --- .../continuous/ContinuousDiffusionStatistic.java | 15 +++------------ 1 file changed, 3 insertions(+), 12 deletions(-) diff --git a/src/dr/evomodel/continuous/ContinuousDiffusionStatistic.java b/src/dr/evomodel/continuous/ContinuousDiffusionStatistic.java index 4fa6325c2f..227b7c8aa5 100644 --- a/src/dr/evomodel/continuous/ContinuousDiffusionStatistic.java +++ b/src/dr/evomodel/continuous/ContinuousDiffusionStatistic.java @@ -461,14 +461,13 @@ public double getStatisticValue(int dim) { } } else if (summaryStat == summaryStatistic.SQUAREDDISTANCE_TIME4_CORRELATION) { List squareddistances = squareElements(distances); - List timesFour = elementsTimesFour(times); if (summaryMode == Mode.SPEARMAN) { - return getSpearmanRho(convertDoubles(timesFour),convertDoubles(squareddistances)); + return getSpearmanRho(convertDoubles(times),convertDoubles(squareddistances)); } else if (summaryMode == Mode.R_SQUARED) { - Regression r = new Regression(convertDoubles(timesFour), convertDoubles(squareddistances)); + Regression r = new Regression(convertDoubles(times), convertDoubles(squareddistances)); return r.getRSquared(); } else { - Regression r = new Regression(convertDoubles(timesFour),convertDoubles(squareddistances)); + Regression r = new Regression(convertDoubles(times),convertDoubles(squareddistances)); return r.getCorrelationCoefficient(); } } else { @@ -511,14 +510,6 @@ private static List squareElements(List inputList) { return squaredList; } - private static List elementsTimesFour (List inputList) { - List returnList = new ArrayList<>(); - for (Double number : inputList) { - returnList.add(number * 4); - } - return returnList; - } - private double[] imputeValue(double[] nodeValue, double[] parentValue, double time, double nodeHeight, double parentHeight, double[] precisionArray, double rate, boolean trueNoise) { final double scaledTimeChild = (time - nodeHeight) * rate;