From 65aa7e71865f8db019af3ba131dc3d8f26e551fb Mon Sep 17 00:00:00 2001 From: fil-monti Date: Tue, 15 Oct 2024 13:19:14 -0700 Subject: [PATCH 01/10] Candidate new version that avoids code duplication --- .../GMRFMultilocusSkyrideLikelihood.java | 520 +++++++++--------- 1 file changed, 268 insertions(+), 252 deletions(-) diff --git a/src/dr/evomodel/coalescent/GMRFMultilocusSkyrideLikelihood.java b/src/dr/evomodel/coalescent/GMRFMultilocusSkyrideLikelihood.java index 5a27693636..970f38a354 100644 --- a/src/dr/evomodel/coalescent/GMRFMultilocusSkyrideLikelihood.java +++ b/src/dr/evomodel/coalescent/GMRFMultilocusSkyrideLikelihood.java @@ -66,20 +66,12 @@ public class GMRFMultilocusSkyrideLikelihood extends GMRFSkyrideLikelihood private double[] storedNumCoalEvents; private double[] gridPoints; - // sortedPoints[i][0] is the time of the i-th grid point or sampling or coalescent event - // sortedPoints[i][1] is 0 if the i-th point is a grid point, 1 if it's a sampling point, and 2 if it's a coalescent point - // sortedPoints[i][2] is the number of lineages present in the interval starting at time sortedPoints[i][0] - private final Parameter phiParameter; private final Parameter ploidyFactors; private double[] ploidySums; private double[] storedPloidySums; -// protected SymmTridiagMatrix precMatrix; -// protected SymmTridiagMatrix storedPrecMatrix; - private final SkygridHelper skygridHelper; - // protected List missingCov; private final List covariates; private final List beta; private final List delta; @@ -96,7 +88,6 @@ public class GMRFMultilocusSkyrideLikelihood extends GMRFSkyrideLikelihood private double[] coalescentEventStatisticValues; - // private List treeList; private final List intervalsList; public GMRFMultilocusSkyrideLikelihood(List intervalsList, @@ -104,342 +95,365 @@ public GMRFMultilocusSkyrideLikelihood(List intervalsList, Parameter groupParameter, Parameter precParameter, Parameter lambda, - Parameter beta, + Parameter betaParameter, MatrixParameter dMatrix, boolean timeAwareSmoothing, double cutOff, int numGridPoints, Parameter phi, Parameter ploidyFactorsParameter) { + this(intervalsList, popParameter, groupParameter, precParameter, lambda, betaParameter, + dMatrix, timeAwareSmoothing, null, null, ploidyFactorsParameter, null, null, + null, null, null, null, null, null, cutOff, numGridPoints, phi); + } + + public GMRFMultilocusSkyrideLikelihood(List intervalsList, + Parameter popParameter, + Parameter groupParameter, + Parameter precParameter, + Parameter lambda, + Parameter betaParameter, + MatrixParameter dMatrix, + boolean timeAwareSmoothing, + Parameter specGridPoints, + List covariates, + Parameter ploidyFactorsParameter, + List firstObservedIndexParameter, + List lastObservedIndexParameter, + List covPrecParametersRecent, + List covPrecParametersDistant, + Parameter recentIndices, + Parameter distantIndices, + List betaList, + List deltaList) { + this(intervalsList, popParameter, groupParameter, precParameter, lambda, betaParameter, + dMatrix, timeAwareSmoothing, specGridPoints, covariates, ploidyFactorsParameter, + firstObservedIndexParameter, lastObservedIndexParameter, covPrecParametersRecent, + covPrecParametersDistant, recentIndices, distantIndices, betaList, deltaList, + specGridPoints != null ? specGridPoints.getParameterValue(specGridPoints.getDimension() - 1) : 0, + specGridPoints != null ? specGridPoints.getDimension() : 0, + null); + + } + + private GMRFMultilocusSkyrideLikelihood(List intervalsList, + Parameter popParameter, + Parameter groupParameter, + Parameter precParameter, + Parameter lambda, + Parameter betaParameter, + MatrixParameter dMatrix, + boolean timeAwareSmoothing, + Parameter specGridPoints, + List covariates, + Parameter ploidyFactorsParameter, + List firstObservedIndexParameter, + List lastObservedIndexParameter, + List covPrecParametersRecent, + List covPrecParametersDistant, + Parameter recentIndices, + Parameter distantIndices, + List betaList, + List deltaList, + double cutOff, + int numGridPoints, + Parameter phi) { super(GMRFSkyrideLikelihoodParser.SKYLINE_LIKELIHOOD); - // adding the key word to the the model means the keyword will be logged in the - // header of the logfile. this.addKeyword("skygrid"); - //if (treeList.size() > 1) { - // this.addKeyword("multilocus"); - // } if (intervalsList.size() > 1) { this.addKeyword("multilocus"); } - - this.popSizeParameter = popParameter; - this.groupSizeParameter = groupParameter; - this.precisionParameter = precParameter; - this.lambdaParameter = lambda; - this.betaParameter = beta; - this.dMatrix = dMatrix; - if (dMatrix != null) { - addVariable(dMatrix); - } - this.timeAwareSmoothing = timeAwareSmoothing; - + this.phiParameter = phi; this.cutOff = cutOff; this.numGridPoints = numGridPoints; - this.phiParameter = phi; - this.ploidyFactors = ploidyFactorsParameter; - setupGridPoints(); - - addVariable(popSizeParameter); - addVariable(precisionParameter); - addVariable(lambdaParameter); - if (betaParameter != null) { - addVariable(betaParameter); - skygridHelper = new SkygridCovariateHelper(); + if (specGridPoints != null) { + this.gridPoints = specGridPoints.getParameterValues(); + if (specGridPoints.getDimension() != numGridPoints) { + throw new IllegalArgumentException("Number of grid points not compatible with" + specGridPoints); + } + if (cutOff != gridPoints[numGridPoints - 1]) { + throw new IllegalArgumentException("Cut-off not compatible with grid points"); + } } else { - skygridHelper = new SkygridHelper(); + setupGridPoints(); } - if (phiParameter != null) { - addVariable(phiParameter); - } - addVariable(ploidyFactors); - this.intervalsList = intervalsList; - //this.numTrees = setTree(treeList); + initializeIndices(firstObservedIndexParameter, lastObservedIndexParameter, recentIndices, distantIndices); + initializeBasicParameters(popParameter, groupParameter, precParameter, lambda, + betaParameter, dMatrix, timeAwareSmoothing); - for (IntervalList intervals : intervalsList) { - addModel((Model) intervals); + this.beta = betaList; + this.delta = initializeDelta(deltaList, betaList); + this.intervalsList = intervalsList; + for (IntervalList intervalList : intervalsList) { + addModel((Model) intervalList); } this.numTrees = intervalsList.size(); + this.ploidyFactors = ploidyFactorsParameter; // # of chromosomes per individual + validatePloidyFactors(); //based on intervalsList.size() TODO: What does this check mean? - int correctFieldLength = getCorrectFieldLength(); +// TODO can I put these this.cov* inside setupCovariates() ? (only for second constructor, ok when null for first) + this.covariates = covariates; + this.covPrecParametersRecent = covPrecParametersRecent; + this.covPrecParametersDistant = covPrecParametersDistant; + setupCovariates(); +// TODO end covariates +// TODO I would like to put all these lengths within a common method. But it does not work +// also this look rather redundant. + int correctFieldLength = getCorrectFieldLength(); if (popSizeParameter.getDimension() <= 1) { - // popSize dimension hasn't been set yet, set it here: popSizeParameter.setDimension(correctFieldLength); } - - fieldLength = popSizeParameter.getDimension(); - if (correctFieldLength != fieldLength) { - throw new IllegalArgumentException("Population size parameter should have length " + correctFieldLength); - } + fieldLength = popSizeParameter.getDimension(); // TODO move before if?? + validateFieldLength(correctFieldLength); oldFieldLength = getCorrectOldFieldLength(); - - - if (ploidyFactors.getDimension() != intervalsList.size()) { - throw new IllegalArgumentException("Ploidy factors parameter should have length " + intervalsList.size()); - } +// TODO end lengths // Field length must be set by this point - wrapSetupIntervals(); +// TODO the following method does not work: (so I silenced it and copied the code here) +// initializeSkygridHelper(betaList, betaParameter, deltaList, lastObservedIndexParameter, firstObservedIndexParameter); - coalescentIntervals = new double[oldFieldLength]; - storedCoalescentIntervals = new double[oldFieldLength]; - sufficientStatistics = new double[fieldLength]; - storedSufficientStatistics = new double[fieldLength]; + if (betaList != null || betaParameter != null) { + if (betaList != null) { + for (Parameter betaParam : betaList) { + addVariable(betaParam); + } + } + if (deltaList != null) { + for (Parameter dParam : deltaList) { + addVariable(dParam); + } + } + if (lastObservedIndexParameter != null || firstObservedIndexParameter != null) { + setupGMRFWeightsForMissingCov(); + skygridHelper = new SkygridMissingCovariateHelper(); + } else { + skygridHelper = new SkygridCovariateHelper(); + } + } else { + skygridHelper = new SkygridHelper(); + } - numCoalEvents = new double[fieldLength]; - storedNumCoalEvents = new double[fieldLength]; - ploidySums = new double[fieldLength]; - storedPloidySums = new double[fieldLength]; + wrapSetupIntervals(); // this does nothing! + initializeArrays(); setupGMRFWeights(); - setupSufficientStatistics(); - +// setupSufficientStatistics(); // TODO this is only in first constructor???; check in the Old version if it is called also in the second constructor indirectly addStatistic(new DeltaStatistic()); - initializationReport(); - - // Force all entries in groupSizeParameter = 1 for compatibility with Tracer - if (groupSizeParameter != null) { - for (int i = 0; i < groupSizeParameter.getDimension(); i++) - groupSizeParameter.setParameterValue(i, 1.0); - } + addVariables(); + initializationReport(); // this.coalescentEventStatisticValues = new double[getNumberOfCoalescentEvents()]; - this.covariates = null; - this.beta = null; - this.delta = null; - this.covPrecParametersRecent = null; - this.covPrecParametersDistant = null; + // TODO the following appears only in the first constructor + // // Force all entries in groupSizeParameter = 1 for compatibility with Tracer +// if (groupSizeParameter != null) { +// for (int i = 0; i < groupSizeParameter.getDimension(); i++) +// groupSizeParameter.setParameterValue(i, 1.0); +// } } + private void initializeBasicParameters(Parameter popParameter, + Parameter groupParameter, Parameter precParameter, + Parameter lambda, Parameter betaParameter, + MatrixParameter dMatrix, boolean timeAwareSmoothing) { + this.popSizeParameter = popParameter; + this.groupSizeParameter = groupParameter; + this.precisionParameter = precParameter; + this.lambdaParameter = lambda; + this.betaParameter = betaParameter; ///TODO CHECK THIS BTW THE TWO CONSTRUCTORS + this.dMatrix = dMatrix; + this.timeAwareSmoothing = timeAwareSmoothing; + } - //rewrite this constructor without duplicating so much code - public GMRFMultilocusSkyrideLikelihood(List intervalsList, - Parameter popParameter, - Parameter groupParameter, - Parameter precParameter, - Parameter lambda, - Parameter betaParameter, - MatrixParameter dMatrix, - boolean timeAwareSmoothing, - Parameter specGridPoints, - List covariates, - Parameter ploidyFactorsParameter, - List firstObservedIndexParameter, - List lastObservedIndexParameter, - List covPrecParametersRecent, - List covPrecParametersDistant, - Parameter recentIndices, - Parameter distantIndices, - List betaList, - List deltaList) { +// private void addKeywords() { +// this.addKeyword("skygrid"); +// if (intervalsList.size() > 1) { +// this.addKeyword("multilocus"); +// } +// } - super(GMRFSkyrideLikelihoodParser.SKYLINE_LIKELIHOOD); + private void addVariables() { + if (phiParameter != null) { + addVariable(phiParameter); + } + addVariable(popSizeParameter); + addVariable(precisionParameter); + addVariable(lambdaParameter); + if (betaParameter != null) { + addVariable(betaParameter); + } + addVariable(ploidyFactors); + if (dMatrix != null) { + addVariable(dMatrix); + } + } - // adding the key word to the the model means the keyword will be logged in the - // header of the logfile. - this.addKeyword("skygrid"); - if (intervalsList.size() > 1) { - this.addKeyword("multilocus"); + + private void validateFieldLength(int correctFieldLength) { + if (correctFieldLength != fieldLength) { + throw new IllegalArgumentException("Population size parameter should have length " + correctFieldLength); } + } - this.gridPoints = specGridPoints.getParameterValues(); - this.numGridPoints = gridPoints.length; - this.cutOff = gridPoints[numGridPoints - 1]; + private void validatePloidyFactors() { + if (ploidyFactors.getDimension() != intervalsList.size()) { + throw new IllegalArgumentException("Ploidy factors parameter should have length " + intervalsList.size()); + } + } - if (firstObservedIndexParameter != null) { - this.firstObservedIndex = new int[firstObservedIndexParameter.size()]; - for (int i = 0; i < firstObservedIndexParameter.size(); i++) { - this.firstObservedIndex[i] = (int) firstObservedIndexParameter.get(i).getParameterValue(0); + private List initializeDelta(List deltaList, List betaList) { + if (deltaList != null) { + return deltaList; + } + if (betaList != null) { // only for the second constructor + List newDelta = new ArrayList<>(); + for (int i = 0; i < betaList.size(); i++) { + Parameter deltaParam = new Parameter.Default(1.0); + deltaParam.setParameterValue(0, 1.0); + newDelta.add(deltaParam); } + return newDelta; // Return initialized newDelta here + } else { + return null; // Return null if both deltaList and betaList are null (first constructor) + } + } - if (recentIndices != null) { - // indices specify which covariates require default unobserved covariate data prior - this.recIndices = new int[firstObservedIndexParameter.size()]; - for (int i = 0; i < firstObservedIndexParameter.size(); i++) { - this.recIndices[i] = (int) recentIndices.getParameterValue(i); - } - } else { - // If specific covariates not specified by indices, need default unobserved covariate data prior for all covariates - this.recIndices = new int[firstObservedIndexParameter.size()]; - for (int i = 0; i < firstObservedIndexParameter.size(); i++) { - this.recIndices[i] = i + 1; - } - } - } + private void initializeIndices(List firstObservedIndexParameter, + List lastObservedIndexParameter, + Parameter recentIndices, + Parameter distantIndices) { + if (firstObservedIndexParameter != null) { + initializeFirstObservedIndex(firstObservedIndexParameter, recentIndices); + } if (lastObservedIndexParameter != null) { - this.lastObservedIndex = new int[lastObservedIndexParameter.size()]; - for (int i = 0; i < lastObservedIndexParameter.size(); i++) { - this.lastObservedIndex[i] = (int) lastObservedIndexParameter.get(i).getParameterValue(0); - } - - if (distantIndices != null) { - // indices specify which covariates require default unobserved covariate data prior - this.distIndices = new int[lastObservedIndexParameter.size()]; - for (int i = 0; i < lastObservedIndexParameter.size(); i++) { - this.distIndices[i] = (int) distantIndices.getParameterValue(i); - } - } else { - // If specific covariates not specified by indices, need default unobserved covariate data prior for all covariates - this.distIndices = new int[lastObservedIndexParameter.size()]; - for (int i = 0; i < lastObservedIndexParameter.size(); i++) { - this.distIndices[i] = i + 1; - } - } + initializeLastObservedIndex(lastObservedIndexParameter, distantIndices); + } + } + private void initializeFirstObservedIndex(List firstObservedIndexParameter, Parameter recentIndices) { + this.firstObservedIndex = new int[firstObservedIndexParameter.size()]; + for (int i = 0; i < firstObservedIndexParameter.size(); i++) { + this.firstObservedIndex[i] = (int) firstObservedIndexParameter.get(i).getParameterValue(0); } + initializeRecIndices(recentIndices, firstObservedIndexParameter.size()); + } - /*else{ - for(int i=0; i < beta.getDimension(); i++) { - this.lastObservedIndex[i] = popParameter.getDimension(); + private void initializeRecIndices(Parameter recentIndices, int size) { + this.recIndices = new int[size]; + if (recentIndices != null) { + for (int i = 0; i < size; i++) { + this.recIndices[i] = (int) recentIndices.getParameterValue(i); + } + } else { + for (int i = 0; i < size; i++) { + this.recIndices[i] = i + 1; } - }*/ - - this.betaParameter = betaParameter; - if (betaParameter != null) { - addVariable(betaParameter); } + } - this.popSizeParameter = popParameter; - this.groupSizeParameter = groupParameter; - this.precisionParameter = precParameter; - this.lambdaParameter = lambda; - this.beta = betaList; + private void initializeLastObservedIndex(List lastObservedIndexParameter, Parameter distantIndices) { + this.lastObservedIndex = new int[lastObservedIndexParameter.size()]; + for (int i = 0; i < lastObservedIndexParameter.size(); i++) { + this.lastObservedIndex[i] = (int) lastObservedIndexParameter.get(i).getParameterValue(0); + } + initializeDistIndices(distantIndices, lastObservedIndexParameter.size()); + } - if (deltaList != null) { - this.delta = deltaList; + private void initializeDistIndices(Parameter distantIndices, int size) { + this.distIndices = new int[size]; + if (distantIndices != null) { + for (int i = 0; i < size; i++) { + this.distIndices[i] = (int) distantIndices.getParameterValue(i); + } } else { - this.delta = new ArrayList(); - if (betaList != null) { - for (int i = 0; i < betaList.size(); i++) { - Parameter deltaParam = new Parameter.Default(1.0); - deltaParam.setParameterValue(0, 1.0); - this.delta.add(deltaParam); - } + for (int i = 0; i < size; i++) { + this.distIndices[i] = i + 1; } } + } - this.dMatrix = dMatrix; - if (dMatrix != null) { - addVariable(dMatrix); +// private void initializeSkygridHelper(List betaList, Parameter betaParameter, +// List deltaList, +// List lastObservedIndexParameter, +// List firstObservedIndexParameter) { +// if (betaList != null || betaParameter != null) { +// if (betaList != null) { +// for (Parameter betaParam : betaList) { +// addVariable(betaParam); +// } +// } +// if (deltaList != null) { +// for (Parameter dParam : deltaList) { +// addVariable(dParam); +// } +// } +// if (lastObservedIndexParameter != null || firstObservedIndexParameter != null) { +// setupGMRFWeightsForMissingCov(); +// skygridHelper = new SkygridMissingCovariateHelper(); +// } else { +// skygridHelper = new SkygridCovariateHelper(); +// } +// } else { +// skygridHelper = new SkygridHelper(); +// } +// } + + private void addBetaAndDeltaVariables(List betaList) { + if (betaList != null) { + for (Parameter betaParam : betaList) { + addVariable(betaParam); + } } - this.timeAwareSmoothing = timeAwareSmoothing; - this.ploidyFactors = ploidyFactorsParameter; - this.covariates = covariates; + if (delta != null) { + for (Parameter dParam : delta) { + addVariable(dParam); + } + } + } + + + private void setupCovariates() { if (covariates != null) { for (MatrixParameter cov : covariates) { addVariable(cov); } } - this.covPrecParametersRecent = covPrecParametersRecent; if (covPrecParametersRecent != null) { for (Parameter covPrecRecent : covPrecParametersRecent) { addVariable(covPrecRecent); } } - this.covPrecParametersDistant = covPrecParametersDistant; if (covPrecParametersDistant != null) { for (Parameter covPrecDistant : covPrecParametersDistant) { addVariable(covPrecDistant); } } + } - addVariable(popSizeParameter); - addVariable(precisionParameter); - addVariable(lambdaParameter); - - addVariable(ploidyFactors); - - this.intervalsList = intervalsList; - - for (IntervalList intervalList : intervalsList) { - addModel((Model) intervalList); - } - - //this.numTrees = setTree(treeList); - this.numTrees = intervalsList.size(); - - int correctFieldLength = getCorrectFieldLength(); - - if (popSizeParameter.getDimension() <= 1) { - // popSize dimension hasn't been set yet, set it here: - popSizeParameter.setDimension(correctFieldLength); - } - - fieldLength = popSizeParameter.getDimension(); - if (correctFieldLength != fieldLength) { - throw new IllegalArgumentException("Population size parameter should have length " + correctFieldLength); - } - - oldFieldLength = getCorrectOldFieldLength(); - - if (ploidyFactors.getDimension() != intervalsList.size()) { - throw new IllegalArgumentException("Ploidy factor parameter should have length " + intervalsList.size()); - } - - // Field length must be set by this point - - if (betaList != null || betaParameter != null) { - if (betaList != null) { - for (Parameter betaParam : betaList) { - addVariable(betaParam); - } - } - if (deltaList != null) { - for (Parameter dParam : deltaList) { - addVariable(dParam); - } - } - if (lastObservedIndexParameter != null || firstObservedIndexParameter != null) { - setupGMRFWeightsForMissingCov(); - skygridHelper = new SkygridMissingCovariateHelper(); - } else { - skygridHelper = new SkygridCovariateHelper(); - } - } else { - skygridHelper = new SkygridHelper(); - } - - wrapSetupIntervals(); + private void initializeArrays() { coalescentIntervals = new double[oldFieldLength]; storedCoalescentIntervals = new double[oldFieldLength]; + sufficientStatistics = new double[fieldLength]; storedSufficientStatistics = new double[fieldLength]; - numCoalEvents = new double[fieldLength]; storedNumCoalEvents = new double[fieldLength]; ploidySums = new double[fieldLength]; storedPloidySums = new double[fieldLength]; - - setupGMRFWeights(); - - addStatistic(new DeltaStatistic()); - - initializationReport(); - - phiParameter = null; - - this.coalescentEventStatisticValues = new double[getNumberOfCoalescentEvents()]; } -// protected int setTree(List treeList) { -//// treesSet = this; -// this.treeList = treeList; -// makeTreeIntervalList(treeList, true); -// return treeList.size(); -// } - - +// TODO this should be @Override !! there is also in the superclass protected int getCorrectFieldLength() { return numGridPoints + 1; @@ -481,6 +495,7 @@ protected void handleModelChangedEvent(Model model, Object object, int index) { } } + // TODO @Override ???? public void initializationReport() { System.out.println("Creating a GMRF smoothed skyride model for multiple loci (SkyGrid)"); System.out.println("\tPopulation sizes: " + popSizeParameter.getDimension()); @@ -541,6 +556,7 @@ private int moveToNextTimeIndex(int treeIndex, int lastTimeIndex, double[] times // } // } +// TODO @Override ??? protected void setupSufficientStatistics() { Arrays.fill(numCoalEvents, 0); From 10514176228eb18e3366d8a488a5ce75066c0079 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Thu, 24 Oct 2024 12:26:26 -0700 Subject: [PATCH 02/10] added code to get Inverse intensity. Still have to manage edge cases --- .../PiecewiseConstantPopulation.java | 41 ++++++++++++++++++- 1 file changed, 40 insertions(+), 1 deletion(-) diff --git a/src/dr/evolution/coalescent/PiecewiseConstantPopulation.java b/src/dr/evolution/coalescent/PiecewiseConstantPopulation.java index e00d52051a..30e6d42b4b 100644 --- a/src/dr/evolution/coalescent/PiecewiseConstantPopulation.java +++ b/src/dr/evolution/coalescent/PiecewiseConstantPopulation.java @@ -148,9 +148,48 @@ public double getIntensity(double t) { } public double getInverseIntensity(double x) { - throw new RuntimeException("Not implemented!"); +// TODO write down error cases +// if (cif != null) { +// if (!intensitiesKnown) { +// setIntervals(intervals, thetas); +// } +// +// int epoch = Collections.binarySearch(cif, x); +// +// if (epoch < 0) { +// epoch = -epoch - 1; +// +// if (epoch > 0) { +// return endTime.get(epoch - 1) + getInverseIntensity(epoch, x - cif.get(epoch - 1)); +// } else { +// assert epoch == 0; +// return getInverseIntensity(0, x); +// } +// } else { +// return endTime.get(epoch); +// } +// } else { + + double inverseIntensity = 0.0; + int epoch = 0; + double x1 = x; + + while (x1 > getIntensity(epoch)) { + x1 -= getIntensity(epoch); + epoch += 1; + if (epoch == intervals.length) { + break; + } + } + + inverseIntensity += getEpochDuration(epoch) + x1 * getEpochDemographic(epoch + 1); +// inverseIntensity += getInverseIntensity(epoch, x1); + + return inverseIntensity; +// } } + public double getUpperBound(int i) { return 1e9; } From c7ba28687599e49699124e2ab839c14964bcc066 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Thu, 24 Oct 2024 12:27:08 -0700 Subject: [PATCH 03/10] temporary change to allow different model with only gridpoints and thetas --- .../coalescent/demographicmodel/PiecewisePopulationModel.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dr/evomodel/coalescent/demographicmodel/PiecewisePopulationModel.java b/src/dr/evomodel/coalescent/demographicmodel/PiecewisePopulationModel.java index 7f98185911..d30ee10742 100644 --- a/src/dr/evomodel/coalescent/demographicmodel/PiecewisePopulationModel.java +++ b/src/dr/evomodel/coalescent/demographicmodel/PiecewisePopulationModel.java @@ -73,7 +73,7 @@ public PiecewisePopulationModel(String name, Parameter N0Parameter, double[] epo if (isLinear) { piecewiseFunction = new PiecewiseLinearPopulation(epochLengths, new double[N0Parameter.getDimension()], units); } else { - piecewiseFunction = new PiecewiseConstantPopulation(epochLengths, new double[N0Parameter.getDimension()], units); + piecewiseFunction = new PiecewiseConstantPopulation(epochLengths, N0Parameter.getParameterValues(), units); } } From e1e9898d3d92a91522d8821fff3a7690454b1cc8 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Thu, 24 Oct 2024 12:27:19 -0700 Subject: [PATCH 04/10] temporary change to allow different model with only gridpoints and thetas --- .../PiecewisePopulationModelParser.java | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/dr/evomodelxml/coalescent/demographicmodel/PiecewisePopulationModelParser.java b/src/dr/evomodelxml/coalescent/demographicmodel/PiecewisePopulationModelParser.java index 759a4eea65..f198a7ba08 100644 --- a/src/dr/evomodelxml/coalescent/demographicmodel/PiecewisePopulationModelParser.java +++ b/src/dr/evomodelxml/coalescent/demographicmodel/PiecewisePopulationModelParser.java @@ -69,8 +69,13 @@ public Object parseXMLObject(XMLObject xo) throws XMLParseException { return new PiecewisePopulationModel(PIECEWISE_POPULATION, epochSizes, epochWidths, isLinear, units); } else { Parameter populationSize = (Parameter) xo.getElementFirstChild(POPULATION_SIZE); - Parameter growthRates = (Parameter) xo.getElementFirstChild(GROWTH_RATES); - return new PiecewisePopulationModel(PIECEWISE_POPULATION, populationSize, growthRates, epochWidths, units); + if (xo.hasChildNamed(GROWTH_RATES)) { + Parameter growthRates = (Parameter) xo.getElementFirstChild(GROWTH_RATES); + return new PiecewisePopulationModel(PIECEWISE_POPULATION, populationSize, growthRates, epochWidths, units); + } else { + return new PiecewisePopulationModel(PIECEWISE_POPULATION, populationSize, epochWidths, false, units); + + } } } @@ -99,7 +104,7 @@ public XMLSyntaxRule[] getSyntaxRules() { new ElementRule(POPULATION_SIZE, new XMLSyntaxRule[]{new ElementRule(Parameter.class)}), new ElementRule(GROWTH_RATES, - new XMLSyntaxRule[]{new ElementRule(Parameter.class)}) + new XMLSyntaxRule[]{new ElementRule(Parameter.class)}, true) ) ), new ElementRule(EPOCH_WIDTHS, From 94208f410a681d3d82accc96726478e6b3b03247 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Thu, 24 Oct 2024 12:29:37 -0700 Subject: [PATCH 05/10] comment --- src/dr/evolution/coalescent/CoalescentSimulator.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dr/evolution/coalescent/CoalescentSimulator.java b/src/dr/evolution/coalescent/CoalescentSimulator.java index 95abf92e4e..4ba1e82bb9 100644 --- a/src/dr/evolution/coalescent/CoalescentSimulator.java +++ b/src/dr/evolution/coalescent/CoalescentSimulator.java @@ -168,7 +168,7 @@ public SimpleNode[] simulateCoalescent(SimpleNode[] nodes, DemographicFunction d currentHeight = getMinimumInactiveHeight(); setCurrentHeight(currentHeight); } - +// TODO should it enforce this? if(!enforceMaxHeight){ // nextCoalescentHeight = currentHeight + DemographicFunction.Utils.getMedianInterval(demographic, getActiveNodeCount(), currentHeight); nextCoalescentHeight = currentHeight + DemographicFunction.Utils.getSimulatedInterval(demographic, From a7437d9fc66419d4caaadbdf2c9c935bff4e0ec0 Mon Sep 17 00:00:00 2001 From: xji3 Date: Wed, 30 Oct 2024 11:43:33 -0500 Subject: [PATCH 06/10] function to respect grid points for Coalescent Simulator --- .../coalescent/CoalescentSimulator.java | 35 +++++++++++++++++++ .../PiecewiseConstantPopulation.java | 10 ++++++ 2 files changed, 45 insertions(+) diff --git a/src/dr/evolution/coalescent/CoalescentSimulator.java b/src/dr/evolution/coalescent/CoalescentSimulator.java index 4ba1e82bb9..b49c567b8b 100644 --- a/src/dr/evolution/coalescent/CoalescentSimulator.java +++ b/src/dr/evolution/coalescent/CoalescentSimulator.java @@ -197,6 +197,41 @@ private double getMinimumInactiveHeight() { } else return Double.POSITIVE_INFINITY; } + private DemographicBoundary getDemographicFunction(DemographicFunction demographic) { + if (demographic instanceof PiecewiseConstantPopulation) { + return DemographicBoundary.PIECEWISE_CONSTANT; + } else { + return DemographicBoundary.PIECEWISE_CONSTANT; + } + } + + enum DemographicBoundary{ + NONE { + @Override + double getNextBoundary(int currentActiveNodeCount, ArrayList sortedNodeList, + DemographicFunction demographic) { + return getMinimumInactiveNode(currentActiveNodeCount, sortedNodeList); + } + }, + PIECEWISE_CONSTANT { + @Override + double getNextBoundary(int currentActiveNodeCount, ArrayList sortedNodeList, DemographicFunction demographic) { + PiecewiseConstantPopulation piecewisePopulationModel = (PiecewiseConstantPopulation) demographic; + final double currentActiveNodeHeight = sortedNodeList.get(currentActiveNodeCount - 1).getHeight(); + final double nextEpochBoundary = piecewisePopulationModel.getNextEpochBoundary(currentActiveNodeHeight); + final double minInactiveNode = getMinimumInactiveNode(currentActiveNodeCount, sortedNodeList); + return nextEpochBoundary < minInactiveNode ? nextEpochBoundary : minInactiveNode; + } + }; + abstract double getNextBoundary(int currentActiveNodeCount, ArrayList sortedNodeList, + DemographicFunction demographic); + double getMinimumInactiveNode(int currentActiveNodeCount, ArrayList sortedNodeList) { + if (currentActiveNodeCount < sortedNodeList.size()) { + return (sortedNodeList.get(currentActiveNodeCount)).getHeight(); + } else return Double.POSITIVE_INFINITY; + } + } + /** * Set the current height. */ diff --git a/src/dr/evolution/coalescent/PiecewiseConstantPopulation.java b/src/dr/evolution/coalescent/PiecewiseConstantPopulation.java index 30e6d42b4b..55cd72ec41 100644 --- a/src/dr/evolution/coalescent/PiecewiseConstantPopulation.java +++ b/src/dr/evolution/coalescent/PiecewiseConstantPopulation.java @@ -103,6 +103,16 @@ public double getDemographic(double t) { return getDemographic(epoch, t1); } + public double getNextEpochBoundary(double t) { + int epoch = 0; // TODO: XJ - really bad code duplication, but how do I return two values without making another structure... + double accumulatedTime = getEpochDuration(epoch); + while (accumulatedTime < t) { + accumulatedTime += getEpochDuration(epoch); + epoch += 1; + } + return accumulatedTime; + } + /** * Gets the integral of intensity function from time 0 to time t. */ From 7d53afc0400236e82dd4516da4ace5ad4e595cec Mon Sep 17 00:00:00 2001 From: xji3 Date: Wed, 30 Oct 2024 18:05:05 -0500 Subject: [PATCH 07/10] more edits --- .../coalescent/CoalescentSimulator.java | 39 ++++++++++++------- .../coalescent/CoalescentSimulator.java | 7 +++- 2 files changed, 30 insertions(+), 16 deletions(-) diff --git a/src/dr/evolution/coalescent/CoalescentSimulator.java b/src/dr/evolution/coalescent/CoalescentSimulator.java index b49c567b8b..cc05eb9ead 100644 --- a/src/dr/evolution/coalescent/CoalescentSimulator.java +++ b/src/dr/evolution/coalescent/CoalescentSimulator.java @@ -49,8 +49,15 @@ */ public class CoalescentSimulator { - public CoalescentSimulator() {} + public CoalescentSimulator() { + demographicBoundary = DemographicBoundary.NONE; + } + + private final DemographicBoundary demographicBoundary; + public CoalescentSimulator(DemographicFunction demographicFunction) { + demographicBoundary = getDemographicBoundary(demographicFunction); + } /** * Simulates a coalescent tree, given a taxon list. @@ -137,7 +144,7 @@ public SimpleNode[] simulateCoalescent(SimpleNode[] nodes, DemographicFunction d // get at least two tips while (getActiveNodeCount() < 2) { - currentHeight = getMinimumInactiveHeight(); + currentHeight = demographicBoundary.getMinimumInactiveNode(activeNodeCount, nodeList); setCurrentHeight(currentHeight); } @@ -154,8 +161,9 @@ public SimpleNode[] simulateCoalescent(SimpleNode[] nodes, DemographicFunction d while (nextCoalescentHeight < maxHeight && (getNodeCount() > 1)) { - if (nextCoalescentHeight >= getMinimumInactiveHeight()) { - currentHeight = getMinimumInactiveHeight(); + final double nextBoundary = demographicBoundary.getNextBoundary(activeNodeCount, nodeList, demographic); + if (nextCoalescentHeight >= nextBoundary) { + currentHeight = nextBoundary; setCurrentHeight(currentHeight); } else { currentHeight = nextCoalescentHeight; @@ -165,7 +173,7 @@ public SimpleNode[] simulateCoalescent(SimpleNode[] nodes, DemographicFunction d if (getNodeCount() > 1) { // get at least two tips while (getActiveNodeCount() < 2) { - currentHeight = getMinimumInactiveHeight(); + currentHeight = demographicBoundary.getMinimumInactiveNode(activeNodeCount, nodeList); setCurrentHeight(currentHeight); } // TODO should it enforce this? @@ -191,13 +199,13 @@ public SimpleNode[] simulateCoalescent(SimpleNode[] nodes, DemographicFunction d /** * @return the height of youngest inactive node. */ - private double getMinimumInactiveHeight() { - if (activeNodeCount < nodeList.size()) { - return (nodeList.get(activeNodeCount)).getHeight(); - } else return Double.POSITIVE_INFINITY; - } +// private double getMinimumInactiveHeight() { +// if (activeNodeCount < nodeList.size()) { +// return (nodeList.get(activeNodeCount)).getHeight(); +// } else return Double.POSITIVE_INFINITY; +// } - private DemographicBoundary getDemographicFunction(DemographicFunction demographic) { + private DemographicBoundary getDemographicBoundary(DemographicFunction demographic) { if (demographic instanceof PiecewiseConstantPopulation) { return DemographicBoundary.PIECEWISE_CONSTANT; } else { @@ -205,11 +213,11 @@ private DemographicBoundary getDemographicFunction(DemographicFunction demograph } } - enum DemographicBoundary{ + private enum DemographicBoundary{ NONE { @Override double getNextBoundary(int currentActiveNodeCount, ArrayList sortedNodeList, - DemographicFunction demographic) { + DemographicFunction demographic) { return getMinimumInactiveNode(currentActiveNodeCount, sortedNodeList); } }, @@ -225,6 +233,7 @@ enum DemographicBoundary{ }; abstract double getNextBoundary(int currentActiveNodeCount, ArrayList sortedNodeList, DemographicFunction demographic); + double getMinimumInactiveNode(int currentActiveNodeCount, ArrayList sortedNodeList) { if (currentActiveNodeCount < sortedNodeList.size()) { return (sortedNodeList.get(currentActiveNodeCount)).getHeight(); @@ -236,7 +245,7 @@ abstract double getNextBoundary(int currentActiveNodeCount, ArrayList Date: Wed, 6 Nov 2024 21:35:46 -0600 Subject: [PATCH 08/10] first attempt to rescale Coalescent simulated tree by demographic function --- .../app/beast/development_parsers.properties | 1 + .../coalescent/CoalescentRescaler.java | 81 ++++++++++++++++++ .../PiecewisePopulationModel.java | 34 +++++++- .../RescaleAwareDemographic.java | 34 ++++++++ .../coalescent/CoalescentRescalerParser.java | 83 +++++++++++++++++++ 5 files changed, 232 insertions(+), 1 deletion(-) create mode 100644 src/dr/evomodel/coalescent/CoalescentRescaler.java create mode 100644 src/dr/evomodel/coalescent/demographicmodel/RescaleAwareDemographic.java create mode 100644 src/dr/evomodelxml/coalescent/CoalescentRescalerParser.java diff --git a/src/dr/app/beast/development_parsers.properties b/src/dr/app/beast/development_parsers.properties index f5dc2ee5e3..2186be757d 100644 --- a/src/dr/app/beast/development_parsers.properties +++ b/src/dr/app/beast/development_parsers.properties @@ -284,6 +284,7 @@ dr.evomodelxml.continuous.RestrictedPartialsParser dr.evomodelxml.coalescent.BNPRSamplingLikelihoodParser dr.inferencexml.distribution.RandomWalkGeneratorParser dr.inference.operators.RandomWalkGammaPrecisionGibbsOperator +dr.evomodelxml.coalescent.CoalescentRescalerParser # Clamps dr.evomodelxml.tree.AncestralTraitTreeModelParser diff --git a/src/dr/evomodel/coalescent/CoalescentRescaler.java b/src/dr/evomodel/coalescent/CoalescentRescaler.java new file mode 100644 index 0000000000..ddc413ab98 --- /dev/null +++ b/src/dr/evomodel/coalescent/CoalescentRescaler.java @@ -0,0 +1,81 @@ +/* + * CoalescentRescaler.java + * + * Copyright © 2002-2024 the BEAST Development Team + * http://beast.community/about + * + * This file is part of BEAST. + * See the NOTICE file distributed with this work for additional + * information regarding copyright ownership and licensing. + * + * BEAST is free software; you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as + * published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * BEAST is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with BEAST; if not, write to the + * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, + * Boston, MA 02110-1301 USA + * + */ + +package dr.evomodel.coalescent; + +import dr.evolution.coalescent.ConstantPopulation; +import dr.evolution.tree.Tree; +import dr.evomodel.bigfasttree.BigFastTreeIntervals; +import dr.evomodel.coalescent.demographicmodel.ConstantPopulationModel; +import dr.evomodel.coalescent.demographicmodel.RescaleAwareDemographic; +import dr.evomodel.tree.TreeModel; + +/** + * Simulates a set of coalescent intervals given a demographic model. + * + * @author Xiang Ji + * @author Marc Suchard + * @author Filippo Monti + */ +public class CoalescentRescaler { + + private TreeModel tree; + + private ConstantPopulationModel constantPopulationModel; + + private RescaleAwareDemographic rescaleAwareDemographic; + + private BigFastTreeIntervals bigFastTreeIntervals; + + public CoalescentRescaler(TreeModel tree, ConstantPopulationModel constantPopulationModel, RescaleAwareDemographic rescaleAwareDemographic) { + this.tree = tree; + this.constantPopulationModel = constantPopulationModel; + this.rescaleAwareDemographic = rescaleAwareDemographic; + this.bigFastTreeIntervals = new BigFastTreeIntervals("constant.population.size.tree", tree); + } + + public TreeModel rescaleTree() { + double[] intensityTimeProducts = new double[bigFastTreeIntervals.getIntervalCount()]; + double[] intervals = new double[bigFastTreeIntervals.getIntervalCount()]; + for (int i = 0; i < bigFastTreeIntervals.getIntervalCount(); i++) { + intervals[i] = bigFastTreeIntervals.getInterval(i); + intensityTimeProducts[i] = intervals[i] / ((ConstantPopulation) constantPopulationModel.getDemographicFunction()).getN0(); + } + double[] scaledIntervals = rescaleAwareDemographic.rescaleInterval(intervals, intensityTimeProducts); + + double height = 0; + for (int i = 0; i < scaledIntervals.length; i++) { + height += scaledIntervals[i]; + int[] intervalNodeNumbers = bigFastTreeIntervals.getNodeNumbersForInterval(i); + tree.setNodeHeightQuietly(tree.getNode(intervalNodeNumbers[1]), height); + } + tree.pushTreeChangedEvent(); + + return tree; + } + +} diff --git a/src/dr/evomodel/coalescent/demographicmodel/PiecewisePopulationModel.java b/src/dr/evomodel/coalescent/demographicmodel/PiecewisePopulationModel.java index d30ee10742..dea808db5e 100644 --- a/src/dr/evomodel/coalescent/demographicmodel/PiecewisePopulationModel.java +++ b/src/dr/evomodel/coalescent/demographicmodel/PiecewisePopulationModel.java @@ -39,7 +39,7 @@ * @author Alexei Drummond * @author Andrew Rambaut */ -public class PiecewisePopulationModel extends DemographicModel { +public class PiecewisePopulationModel extends DemographicModel implements RescaleAwareDemographic { // // Public stuff @@ -166,6 +166,38 @@ protected void restoreState() { protected void acceptState() { } // no additional state needs accepting + @Override + public double[] rescaleInterval(double[] intervals, double[] targetIntensityIntervalProducts) { + + if (!(piecewiseFunction instanceof PiecewiseConstantPopulation)) { + throw new RuntimeException("Not yet implemented!"); + } + + PiecewiseConstantPopulation piecewiseConstantPopulation = (PiecewiseConstantPopulation) piecewiseFunction; + + double[] scaledIntervals = new double[intervals.length]; + + int currentEpochIndex = 0; + for (int i = 0; i < intervals.length; i++) { + if (intervals[i] == 0) { + scaledIntervals[i] = 0; + } else { + double accumulatedProduct = 0; + double accumulatedIntervalLength = 0; + while(accumulatedProduct + piecewiseConstantPopulation.getEpochDuration(currentEpochIndex) / piecewiseConstantPopulation.getEpochDemographic(currentEpochIndex) < targetIntensityIntervalProducts[i]) { + accumulatedProduct += piecewiseConstantPopulation.getEpochDuration(currentEpochIndex) / piecewiseConstantPopulation.getEpochDemographic(currentEpochIndex); + accumulatedIntervalLength += piecewiseConstantPopulation.getEpochDuration(currentEpochIndex); + currentEpochIndex++; + } + final double remainingProduct = targetIntensityIntervalProducts[i] - accumulatedProduct; + accumulatedIntervalLength += remainingProduct * piecewiseConstantPopulation.getEpochDemographic(currentEpochIndex); + scaledIntervals[i] = accumulatedIntervalLength; + } + } + + return scaledIntervals; + } + public class GrowthRateStatistic extends Statistic.Abstract { public GrowthRateStatistic() { diff --git a/src/dr/evomodel/coalescent/demographicmodel/RescaleAwareDemographic.java b/src/dr/evomodel/coalescent/demographicmodel/RescaleAwareDemographic.java new file mode 100644 index 0000000000..9a9fb7b9d1 --- /dev/null +++ b/src/dr/evomodel/coalescent/demographicmodel/RescaleAwareDemographic.java @@ -0,0 +1,34 @@ +/* + * RescaleAwareDemographic.java + * + * Copyright © 2002-2024 the BEAST Development Team + * http://beast.community/about + * + * This file is part of BEAST. + * See the NOTICE file distributed with this work for additional + * information regarding copyright ownership and licensing. + * + * BEAST is free software; you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as + * published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * BEAST is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with BEAST; if not, write to the + * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, + * Boston, MA 02110-1301 USA + * + */ + +package dr.evomodel.coalescent.demographicmodel; + +public interface RescaleAwareDemographic { + + double[] rescaleInterval(double[] intervals, double[] targetIntensityIntervalProducts); + +} diff --git a/src/dr/evomodelxml/coalescent/CoalescentRescalerParser.java b/src/dr/evomodelxml/coalescent/CoalescentRescalerParser.java new file mode 100644 index 0000000000..e830b6ce06 --- /dev/null +++ b/src/dr/evomodelxml/coalescent/CoalescentRescalerParser.java @@ -0,0 +1,83 @@ +/* + * CoalescentRescalerParser.java + * + * Copyright © 2002-2024 the BEAST Development Team + * http://beast.community/about + * + * This file is part of BEAST. + * See the NOTICE file distributed with this work for additional + * information regarding copyright ownership and licensing. + * + * BEAST is free software; you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as + * published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * BEAST is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with BEAST; if not, write to the + * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, + * Boston, MA 02110-1301 USA + * + */ + +package dr.evomodelxml.coalescent; + +import dr.evolution.coalescent.ConstantPopulation; +import dr.evolution.tree.Tree; +import dr.evolution.util.TaxonList; +import dr.evomodel.coalescent.CoalescentRescaler; +import dr.evomodel.coalescent.demographicmodel.ConstantPopulationModel; +import dr.evomodel.coalescent.demographicmodel.DemographicModel; +import dr.evomodel.coalescent.demographicmodel.RescaleAwareDemographic; +import dr.evomodel.tree.DefaultTreeModel; +import dr.evomodel.tree.TreeModel; +import dr.xml.*; + +public class CoalescentRescalerParser extends AbstractXMLObjectParser { + + public static final String COALESCENT_RESCALER = "coalescentRescaler"; + + @Override + public Object parseXMLObject(XMLObject xo) throws XMLParseException { + TreeModel tree = (TreeModel) xo.getChild(TreeModel.class); + if (tree == null) { + tree = new DefaultTreeModel((Tree) xo.getChild(Tree.class)); + } + ConstantPopulationModel constantPopulation = (ConstantPopulationModel) xo.getChild(ConstantPopulationModel.class); + RescaleAwareDemographic rescaleAwareDemographic = (RescaleAwareDemographic) xo.getChild(RescaleAwareDemographic.class); + CoalescentRescaler coalescentRescaler = new CoalescentRescaler(tree, constantPopulation, rescaleAwareDemographic); + return coalescentRescaler.rescaleTree(); + } + + @Override + public XMLSyntaxRule[] getSyntaxRules() { + return rules; + } + + private final XMLSyntaxRule[] rules = { + new XORRule(new ElementRule(TreeModel.class), + new ElementRule(Tree.class)), + new ElementRule(ConstantPopulationModel.class, 0, Integer.MAX_VALUE), + new ElementRule(RescaleAwareDemographic.class, 0, Integer.MAX_VALUE), + }; + + @Override + public String getParserDescription() { + return ""; + } + + @Override + public Class getReturnType() { + return TreeModel.class; + } + + @Override + public String getParserName() { + return COALESCENT_RESCALER; + } +} From d05608f2beadaaee4d68f7329f81764d07ca313c Mon Sep 17 00:00:00 2001 From: xji3 Date: Thu, 7 Nov 2024 09:12:13 -0600 Subject: [PATCH 09/10] Revert "more edits" This reverts commit 7d53afc0400236e82dd4516da4ace5ad4e595cec. --- .../coalescent/CoalescentSimulator.java | 39 +++++++------------ .../coalescent/CoalescentSimulator.java | 7 +--- 2 files changed, 16 insertions(+), 30 deletions(-) diff --git a/src/dr/evolution/coalescent/CoalescentSimulator.java b/src/dr/evolution/coalescent/CoalescentSimulator.java index cc05eb9ead..b49c567b8b 100644 --- a/src/dr/evolution/coalescent/CoalescentSimulator.java +++ b/src/dr/evolution/coalescent/CoalescentSimulator.java @@ -49,15 +49,8 @@ */ public class CoalescentSimulator { - public CoalescentSimulator() { - demographicBoundary = DemographicBoundary.NONE; - } - - private final DemographicBoundary demographicBoundary; + public CoalescentSimulator() {} - public CoalescentSimulator(DemographicFunction demographicFunction) { - demographicBoundary = getDemographicBoundary(demographicFunction); - } /** * Simulates a coalescent tree, given a taxon list. @@ -144,7 +137,7 @@ public SimpleNode[] simulateCoalescent(SimpleNode[] nodes, DemographicFunction d // get at least two tips while (getActiveNodeCount() < 2) { - currentHeight = demographicBoundary.getMinimumInactiveNode(activeNodeCount, nodeList); + currentHeight = getMinimumInactiveHeight(); setCurrentHeight(currentHeight); } @@ -161,9 +154,8 @@ public SimpleNode[] simulateCoalescent(SimpleNode[] nodes, DemographicFunction d while (nextCoalescentHeight < maxHeight && (getNodeCount() > 1)) { - final double nextBoundary = demographicBoundary.getNextBoundary(activeNodeCount, nodeList, demographic); - if (nextCoalescentHeight >= nextBoundary) { - currentHeight = nextBoundary; + if (nextCoalescentHeight >= getMinimumInactiveHeight()) { + currentHeight = getMinimumInactiveHeight(); setCurrentHeight(currentHeight); } else { currentHeight = nextCoalescentHeight; @@ -173,7 +165,7 @@ public SimpleNode[] simulateCoalescent(SimpleNode[] nodes, DemographicFunction d if (getNodeCount() > 1) { // get at least two tips while (getActiveNodeCount() < 2) { - currentHeight = demographicBoundary.getMinimumInactiveNode(activeNodeCount, nodeList); + currentHeight = getMinimumInactiveHeight(); setCurrentHeight(currentHeight); } // TODO should it enforce this? @@ -199,13 +191,13 @@ public SimpleNode[] simulateCoalescent(SimpleNode[] nodes, DemographicFunction d /** * @return the height of youngest inactive node. */ -// private double getMinimumInactiveHeight() { -// if (activeNodeCount < nodeList.size()) { -// return (nodeList.get(activeNodeCount)).getHeight(); -// } else return Double.POSITIVE_INFINITY; -// } + private double getMinimumInactiveHeight() { + if (activeNodeCount < nodeList.size()) { + return (nodeList.get(activeNodeCount)).getHeight(); + } else return Double.POSITIVE_INFINITY; + } - private DemographicBoundary getDemographicBoundary(DemographicFunction demographic) { + private DemographicBoundary getDemographicFunction(DemographicFunction demographic) { if (demographic instanceof PiecewiseConstantPopulation) { return DemographicBoundary.PIECEWISE_CONSTANT; } else { @@ -213,11 +205,11 @@ private DemographicBoundary getDemographicBoundary(DemographicFunction demograph } } - private enum DemographicBoundary{ + enum DemographicBoundary{ NONE { @Override double getNextBoundary(int currentActiveNodeCount, ArrayList sortedNodeList, - DemographicFunction demographic) { + DemographicFunction demographic) { return getMinimumInactiveNode(currentActiveNodeCount, sortedNodeList); } }, @@ -233,7 +225,6 @@ private enum DemographicBoundary{ }; abstract double getNextBoundary(int currentActiveNodeCount, ArrayList sortedNodeList, DemographicFunction demographic); - double getMinimumInactiveNode(int currentActiveNodeCount, ArrayList sortedNodeList) { if (currentActiveNodeCount < sortedNodeList.size()) { return (sortedNodeList.get(currentActiveNodeCount)).getHeight(); @@ -245,7 +236,7 @@ abstract double getNextBoundary(int currentActiveNodeCount, ArrayList Date: Thu, 7 Nov 2024 09:12:23 -0600 Subject: [PATCH 10/10] Revert "function to respect grid points for Coalescent Simulator" This reverts commit a7437d9fc66419d4caaadbdf2c9c935bff4e0ec0. --- .../coalescent/CoalescentSimulator.java | 35 ------------------- .../PiecewiseConstantPopulation.java | 10 ------ 2 files changed, 45 deletions(-) diff --git a/src/dr/evolution/coalescent/CoalescentSimulator.java b/src/dr/evolution/coalescent/CoalescentSimulator.java index b49c567b8b..4ba1e82bb9 100644 --- a/src/dr/evolution/coalescent/CoalescentSimulator.java +++ b/src/dr/evolution/coalescent/CoalescentSimulator.java @@ -197,41 +197,6 @@ private double getMinimumInactiveHeight() { } else return Double.POSITIVE_INFINITY; } - private DemographicBoundary getDemographicFunction(DemographicFunction demographic) { - if (demographic instanceof PiecewiseConstantPopulation) { - return DemographicBoundary.PIECEWISE_CONSTANT; - } else { - return DemographicBoundary.PIECEWISE_CONSTANT; - } - } - - enum DemographicBoundary{ - NONE { - @Override - double getNextBoundary(int currentActiveNodeCount, ArrayList sortedNodeList, - DemographicFunction demographic) { - return getMinimumInactiveNode(currentActiveNodeCount, sortedNodeList); - } - }, - PIECEWISE_CONSTANT { - @Override - double getNextBoundary(int currentActiveNodeCount, ArrayList sortedNodeList, DemographicFunction demographic) { - PiecewiseConstantPopulation piecewisePopulationModel = (PiecewiseConstantPopulation) demographic; - final double currentActiveNodeHeight = sortedNodeList.get(currentActiveNodeCount - 1).getHeight(); - final double nextEpochBoundary = piecewisePopulationModel.getNextEpochBoundary(currentActiveNodeHeight); - final double minInactiveNode = getMinimumInactiveNode(currentActiveNodeCount, sortedNodeList); - return nextEpochBoundary < minInactiveNode ? nextEpochBoundary : minInactiveNode; - } - }; - abstract double getNextBoundary(int currentActiveNodeCount, ArrayList sortedNodeList, - DemographicFunction demographic); - double getMinimumInactiveNode(int currentActiveNodeCount, ArrayList sortedNodeList) { - if (currentActiveNodeCount < sortedNodeList.size()) { - return (sortedNodeList.get(currentActiveNodeCount)).getHeight(); - } else return Double.POSITIVE_INFINITY; - } - } - /** * Set the current height. */ diff --git a/src/dr/evolution/coalescent/PiecewiseConstantPopulation.java b/src/dr/evolution/coalescent/PiecewiseConstantPopulation.java index 55cd72ec41..30e6d42b4b 100644 --- a/src/dr/evolution/coalescent/PiecewiseConstantPopulation.java +++ b/src/dr/evolution/coalescent/PiecewiseConstantPopulation.java @@ -103,16 +103,6 @@ public double getDemographic(double t) { return getDemographic(epoch, t1); } - public double getNextEpochBoundary(double t) { - int epoch = 0; // TODO: XJ - really bad code duplication, but how do I return two values without making another structure... - double accumulatedTime = getEpochDuration(epoch); - while (accumulatedTime < t) { - accumulatedTime += getEpochDuration(epoch); - epoch += 1; - } - return accumulatedTime; - } - /** * Gets the integral of intensity function from time 0 to time t. */