first version of template that produces XML that runs

rbouckaert · Jul 29, 2020 · 531a67c · 531a67c
1 parent 1183dee
commit 531a67c
Show file tree

Hide file tree

Showing 5 changed files with 98 additions and 50 deletions.
diff --git a/src/beast/app/beauti/StarBeastAlignmentProvider3.java b/src/beast/app/beauti/StarBeastAlignmentProvider3.java
@@ -13,8 +13,8 @@ public List<BEASTInterface> getAlignments(BeautiDoc doc, File[] files) {
 		final List<BEASTInterface> newAlignments = super.getAlignments(doc, files);
 		final int alignmentCount = newAlignments.size();
 
-		doc.autoSetClockRate = false;
-		doc.beauti.autoSetClockRate.setSelected(false);
+		doc.autoSetClockRate = true;
+		doc.beauti.autoSetClockRate.setSelected(true);
 
 		System.out.println(String.format("N_ALIGNMENTS = %d", doc.alignments.size()));
 		// initialize delta exchange operator in order to increase weight to something more sensible

diff --git a/src/starbeast3/StarBeastStartState.java b/src/starbeast3/StarBeastStartState.java
@@ -32,6 +32,7 @@
 import beast.util.ClusterTree;
 import genekernel.GTKPrior;
 import starbeast3.evolution.branchratemodel.BranchRateModelSB3;
+import starbeast3.evolution.branchratemodel.UCRelaxedClockModelSB3;
 
 /**
 * @author Joseph Heled
@@ -99,6 +100,18 @@ public void initAndValidate() {
         super.initAndValidate();
         hasCalibrations = calibratedYule.get() != null;
         genes = genesInput.get();
+
+        if (speciesTreeRatesInput.get() != null) {
+        	BranchRateModelSB3 s = speciesTreeRatesInput.get();
+        	if (s instanceof UCRelaxedClockModelSB3) {
+        		UCRelaxedClockModelSB3 s2 = (UCRelaxedClockModelSB3) s;
+        		RealParameter p = s2.realRatesInput.get();
+        		if (p != null) {
+        			rates = p;
+        			lowerRate = 0.1; // p .getLower();
+        		}
+        	}
+        }
     }
 
     @Override
@@ -157,7 +170,13 @@ public void initStateNodes() {
             }
         }
 
+        if (rates != null) {
+        	// rates.setLower(lowerRate);
+        }
     }
+
+    RealParameter rates;
+    double lowerRate;
 
     private double[] firstMeetings(final Tree gtree, final Map<String, Integer> tipName2Species, final int speciesCount) {
         final Node[] nodes = gtree.listNodesPostOrder(null, null);

diff --git a/src/starbeast3/operators/ParallelMCMCTreeOperator.java b/src/starbeast3/operators/ParallelMCMCTreeOperator.java
@@ -34,34 +34,11 @@
 public class ParallelMCMCTreeOperator extends Operator implements MultiStepOperator {
     final public Input<Boolean> useBactrianOperatorsInput = new Input<>("bactrian", "flag to indicate that bactrian operators should be used where possible", true);
 
-	@Description("Distribution on a tree conditinionally independent from all other distributions given the state of the rest of parameter space")
-	public class TreeDistribution {
-		Tree tree;
-		Distribution treelikelihood;
-		GeneTreeForSpeciesTreeDistribution geneprior;
-
-		public Tree getTree() {return tree;}
-		public void setTree(Tree tree) {this.tree = tree;}
-		public Distribution getTreelikelihood() {return treelikelihood;}
-		public void setTreelikelihood(Distribution treelikelihood) {this.treelikelihood = treelikelihood;}
-		public GeneTreeForSpeciesTreeDistribution getGeneprior() {return geneprior;}
-		public void setGeneprior(GeneTreeForSpeciesTreeDistribution geneprior) {this.geneprior = geneprior;}
-
-		public TreeDistribution(@Param(name="tree", description="tree for which") Tree tree,
-				@Param(name="treelikelihood", description="treelikelihood part of the distribution") Distribution treelikelihood,
-				@Param(name="geneprior", description="prior on the gene tree") GeneTreeForSpeciesTreeDistribution geneprior) {
-			this.tree = tree;
-			this.treelikelihood = treelikelihood;
-			this.geneprior = geneprior;
-		}
-	}
-
-
 	final public Input<Long> chainLengthInput =
             new Input<>("chainLength", "Length of the MCMC chain: each individual ParallelMCMC performs chainLength/nrOfThreads samples",
                     Input.Validate.REQUIRED);
 
-	final public Input<List<TreeDistribution>> distributionInput = new Input<>("distribution", 
+	final public Input<List<ParallelMCMCTreeOperatorTreeDistribution>> distributionInput = new Input<>("distribution", 
 			"Distribution on a tree conditinionally independent from all other distributions given the state of the rest"
 			+ "of parameter space. ",
 			new ArrayList<>());
@@ -79,7 +56,7 @@ public TreeDistribution(@Param(name="tree", description="tree for which") Tree t
 
 	@Override
 	public void initAndValidate() {
-		List<TreeDistribution> distributions = distributionInput.get();
+		List<ParallelMCMCTreeOperatorTreeDistribution> distributions = distributionInput.get();
 	    otherState = otherStateInput.get();
 
 		int nrOfThreads = maxNrOfThreadsInput.get() > 0 ? 
@@ -105,11 +82,11 @@ public int stepCount() {
 		return mcmcs.size();
 	}
 
-	private ParallelMCMC createParallelMCMC(List<TreeDistribution> distributions, long chainLength) {
+	private ParallelMCMC createParallelMCMC(List<ParallelMCMCTreeOperatorTreeDistribution> distributions, long chainLength) {
 		List<Distribution> distrs = new ArrayList<>();
 		List<StateNode> stateNodes = new ArrayList<>();
 		List<Operator> operators = new ArrayList<>();
-		for (TreeDistribution d : distributions) {
+		for (ParallelMCMCTreeOperatorTreeDistribution d : distributions) {
 			distrs.add(d.geneprior);
 			distrs.add(d.treelikelihood);
 

diff --git a/src/starbeast3/operators/ParallelMCMCTreeOperatorTreeDistribution.java b/src/starbeast3/operators/ParallelMCMCTreeOperatorTreeDistribution.java
@@ -0,0 +1,34 @@
+package starbeast3.operators;
+
+import beast.core.BEASTObject;
+import beast.core.Description;
+import beast.core.Distribution;
+import beast.core.Param;
+import beast.evolution.tree.Tree;
+import starbeast3.GeneTreeForSpeciesTreeDistribution;
+
+@Description("Distribution on a tree conditinionally independent from all other distributions given the state of the rest of parameter space")
+public class ParallelMCMCTreeOperatorTreeDistribution extends BEASTObject {
+	Tree tree;
+	Distribution treelikelihood;
+	GeneTreeForSpeciesTreeDistribution geneprior;
+
+	public Tree getTree() {return tree;}
+	public void setTree(Tree tree) {this.tree = tree;}
+	public Distribution getTreelikelihood() {return treelikelihood;}
+	public void setTreelikelihood(Distribution treelikelihood) {this.treelikelihood = treelikelihood;}
+	public GeneTreeForSpeciesTreeDistribution getGeneprior() {return geneprior;}
+	public void setGeneprior(GeneTreeForSpeciesTreeDistribution geneprior) {this.geneprior = geneprior;}
+
+	public ParallelMCMCTreeOperatorTreeDistribution(@Param(name="tree", description="tree for which") Tree tree,
+			@Param(name="treelikelihood", description="treelikelihood part of the distribution") Distribution treelikelihood,
+			@Param(name="geneprior", description="prior on the gene tree") GeneTreeForSpeciesTreeDistribution geneprior) {
+		this.tree = tree;
+		this.treelikelihood = treelikelihood;
+		this.geneprior = geneprior;
+	}
+
+	@Override
+	public void initAndValidate() {
+	}
+}
diff --git a/templates/StarBEAST3.xml b/templates/StarBEAST3.xml
@@ -103,17 +103,23 @@
 			starbeast3.RandomLocalRates.rates,
 			starbeast3.RandomLocalRates.noCache,
 			starbeast3.ConstantPopulations.speciesTree,
-			starbeast3.LinearWithConstantRoot.speciesTree"
+			starbeast3.LinearWithConstantRoot.speciesTree,			
+			starbeast3.GeneTreeForSpeciesTreeDistribution.useTipDates,
+			starbeast3.GeneTreeForSpeciesTreeDistribution.sampling,
+			starbeast3.GeneTreeForSpeciesTreeDistribution.tree,
+			starbeast3.GeneTreeForSpeciesTreeDistribution.populationModel,
+			starbeast3.GeneTreeForSpeciesTreeDistribution.speciesTree,
+			starbeast3.GeneTreeForSpeciesTreeDistribution.treeIntervals"
 		buttonLabelMap="beast.app.beauti.BeautiInitDlg.&gt;&gt; details=Edit parameters">
 
 		<panel spec="BeautiPanelConfig" panelname="Partitions" tiptext="Data Partitions" path="distribution/distribution[id='likelihood']/distribution/data" hasPartitions="none" icon="2220.pngx" forceExpansion="FALSE" type="beast.evolution.alignment.Alignment"/>
 		<mergepoint id="aux-partitions-panels"/>
 
 		<panel spec="BeautiPanelConfig" panelname="Taxon sets" tiptext="Specify taxon sets that make up species" path="operator[id='Reheight.t:Species']/taxonset" hasPartitions="none" icon="1.pngx" forceExpansion="TRUE"/>
 
-		<panel spec="BeautiPanelConfig" panelname="Tip Dates" tiptext="Allows to specify date that a taxon was sampled" path="operator[id='coordinatedUniform.t:Species']/speciesTree" hasPartitions="none" icon="2.png.x" forceExpansion="TRUE"/>
-		<!-- <mergepoint id="aux-tipdates-panels"/> -->
-
+<!-- Tip dates need implementing
+		<panel spec="BeautiPanelConfig" panelname="Tip Dates" tiptext="Allows to specify date that a taxon was sampled" path="state/stateNode[id='Tree.t:Species']" hasPartitions="none" icon="2.png.x" forceExpansion="TRUE"/>
+-->
 		<panel spec="BeautiPanelConfig" panelname="Gene Ploidy" tiptext="The ploidy for each gene (locus)" path="distribution/distribution[id='speciescoalescent']/distribution" hasPartitions="none" icon="2.pngx" forceExpansion="TRUE"/>
 
 		<panel spec="BeautiPanelConfig" panelname="Site Model" tiptext="Site model and substitution model specifications" path="siteModel" hasPartitions="SiteModel" icon="3.pngx" forceExpansion="TRUE"/>
@@ -169,16 +175,9 @@
 			<distr spec="beast.math.distributions.Uniform" lower="0.0" upper="1.0"/>
 		</prior>
 
-    	<distribution
-	    	id="ParallelMCMCTreeOperator.$(n)" 
-    		spec="starbeast3.operators.ParallelMCMCTreeOperator$TreeDistribution"
-    		tree="@Tree.t:$(n)" 
-    		geneprior="@treePrior.t:$(n)" 
-    		treelikelihood="@treeLikelihood.$(n)"/>
-
 	    <distribution data="@$(n)" id="treeLikelihood.$(n)" spec="TreeLikelihood" tree="@Tree.t:$(n)">
             <siteModel gammaCategoryCount="1" id="SiteModel.s:$(n)" spec="SiteModel">
-				<mutationRate spec="parameter.RealParameter" id="mutationRate.s:$(n)" value="1.0" estimate="true"/>
+				<mutationRate spec="parameter.RealParameter" id="mutationRate.s:$(n)" value="1.0" estimate="false"/>
     	        <shape id="gammaShape.s:$(n)" spec="parameter.RealParameter" value="1.0" estimate="false" lower="0"/>
                 <proportionInvariant id="proportionInvariant.s:$(n)" spec="parameter.RealParameter" value="0.0" lower="0" upper="1" estimate="false"/>
         	    <substModel id="hky.s:$(n)" kappa="@kappa.s:$(n)" spec="HKY">
@@ -188,6 +187,13 @@
         	<branchRateModel id="StrictClock.c:$(n)" spec="beast.evolution.branchratemodel.StrictClockModel" clock.rate="@clockRate.c:$(n)"/>
 	    </distribution>
 
+    	<distribution
+	    	id="ParallelMCMCTreeOperator.$(n)" 
+    		spec="starbeast3.operators.ParallelMCMCTreeOperatorTreeDistribution"
+    		tree="@Tree.t:$(n)" 
+    		geneprior="@treePrior.t:$(n)" 
+    		treelikelihood="@treeLikelihood.$(n)"/>
+
 		<log id="TreeHeight.t:$(n)" spec="beast.evolution.tree.TreeHeightLogger" tree="@Tree.t:$(n)"/>
 
 		<operator id="ParallelMCMCRealParameterOperator" spec="starbeast3.operators.ParallelMCMCRealParameterOperator" 
@@ -213,28 +219,28 @@
 		<connect srcID="mutationRate.s:$(n)" targetID="state" inputName="stateNode" if="inposterior(mutationRate.s:$(n)) and mutationRate.s:$(n)/estimate=true"/>
 		<connect srcID="gammaShape.s:$(n)" targetID="state" inputName="stateNode" if="inposterior(gammaShape.s:$(n)) and gammaShape.s:$(n)/estimate=true"/>
 
-		<connect srcID="treePrior.t:$(n)" targetID="speciescoalescent" inputName="distribution" if="inposterior(treePrior.t:$(n)) and Tree.t:$(n)/estimate=true"/>
+		<connect srcID="treePrior.t:$(n)" targetID="speciescoalescent" inputName="distribution" if="inposterior(Tree.t:$(n)) and Tree.t:$(n)/estimate=true"/>
 		<connect srcID="ClockPrior.c:$(n)" targetID="prior" inputName="distribution" if="inposterior(ClockPrior.c:$(n)) and clockRate.c:$(n)/estimate=true"/>
 		<connect srcID="KappaPrior.s:$(n)" targetID="prior" inputName="distribution" if="inposterior(KappaPrior.s:$(n)) and kappa.s:$(n)/estimate=true"/>
 		<connect srcID="GammaShapePrior.s:$(n)" targetID="prior" inputName="distribution" if="inlikelihood(gammaShape.s:$(n)) and gammaShape.s:$(n)/estimate=true"/>
 		<connect srcID="PropInvariantPrior.s:$(n)" targetID="prior" inputName="distribution" if="inlikelihood(proportionInvariant.s:$(n)) and proportionInvariant.s:$(n)/estimate=true"/>
 
 
-		<connect srcID="Tree.t:$(n)" targetID="SBI" inputName="gene" if="inposterior(Tree.t:$(n)) and Tree.t:$(n)/estimate=true"/>
+		<connect srcID="Tree.t:$(n)" targetID="SBI" inputName="gene" if="inlikelihood(Tree.t:$(n)) and Tree.t:$(n)/estimate=true"/>
 		<connect srcID="proportionInvariantScaler.s:$(n)" targetID="mcmc" inputName="operator" if="inposterior(proportionInvariant.s:$(n)) and proportionInvariant.s:$(n)/estimate=true"/>
 		<connect srcID="mutationRateScaler.s:$(n)" targetID="mcmc" inputName="operator" if="inposterior(mutationRate.s:$(n)) and mutationRate.s:$(n)/estimate=true"/>
 		<connect srcID="gammaShapeScaler.s:$(n)" targetID="mcmc" inputName="operator" if="inposterior(gammaShape.s:$(n)) and gammaShape.s:$(n)/estimate=true"/>
 		<connect srcID="ParallelMCMCTreeOperator.$(n)" targetID="ParallelMCMCTreeOperator" inputName="distribution" if="inposterior(Tree.t:$(n)) and Tree.t:$(n)/estimate=true"/>
 
-		<connect srcID="treePrior.t:$(n)" targetID="Reheight.t:Species" inputName="gene" if="inposterior(treePrior.t:$(n)) and treePrior.t:$(n)/estimate=true"/>
-		<connect srcID="treePrior.t:$(n)" targetID="CoordinatedExponential.t:Species" inputName="gene" if="inposterior(treePrior.t:$(n)) and treePrior.t:$(n)/estimate=true"/>
-		<connect srcID="treePrior.t:$(n)" targetID="CoordinatedUniform.t:Species" inputName="gene" if="inposterior(treePrior.t:$(n)) and treePrior.t:$(n)/estimate=true"/>
+		<connect srcID="treePrior.t:$(n)" targetID="Reheight.t:Species" inputName="gene" if="inposterior(Tree.t:$(n)) and Tree.t:$(n)/estimate=true"/>
+		<connect srcID="treePrior.t:$(n)" targetID="CoordinatedExponential.t:Species" inputName="gene" if="inposterior(Tree.t:$(n)) and Tree.t:$(n)/estimate=true"/>
+		<connect srcID="treePrior.t:$(n)" targetID="CoordinatedUniform.t:Species" inputName="gene" if="inposterior(Tree.t:$(n)) and Tree.t:$(n)/estimate=true"/>
 		<connect srcID="Tree.t:$(n)" targetID="updown.all" inputName="down" if="inposterior(Tree.t:$(n)) and Tree.t:$(n)/estimate=true"/>
 		<connect srcID="clockRate.c:$(n)" targetID="updown.all" inputName="up" if="inposterior(clockRate.c:$(n)) and clockRate.c:$(n)/estimate=true"/>
 		<connect srcID="birthRate.t:Species" targetID="updown.all" inputName="up" if="inposterior(birthRate.t:Species) and birthRate.t:Species/estimate=true"/>
 		<connect srcID="birthRate.t:Species" targetID="AdaptableVarianceMultivariateNormalOperatorLogTransform" inputName="f" if="inposterior(birthRate.t:Species) and birthRate.t:Species/estimate=true"/>
 		<connect srcID="treeLikelihood.$(n)" targetID="tracelog" inputName="log" if="inposterior(treeLikelihood.$(n))"/>
-		<connect srcID="treePrior.t:$(n)" targetID="tracelog" inputName="log" if="inposterior(treePrior.t:$(n)) and treePrior.t:$(n)/estimate=true"/>
+		<connect srcID="treePrior.t:$(n)" targetID="tracelog" inputName="log" if="inposterior(Tree.t:$(n)) and Tree.t:$(n)/estimate=true"/>
 		<connect srcID="TreeHeight.t:$(n)" targetID="tracelog" inputName="log" if="inposterior(TreeHeight.t:$(n)) and TreeHeight.t:$(n)/estimate=true"/>
 
 		<connect srcID="kappa.s:$(n)" targetID="tracelog" inputName="log" if="inposterior(kappa.s:$(n)) and kappa.s:$(n)/estimate=true"/>
@@ -248,7 +254,18 @@
 		<connect srcID="kappa.s:$(n)" targetID="AdaptableVarianceMultivariateNormalOperatorLogTransform" inputName="f" if="inposterior(kappa.s:$(n)) and kappa.s:$(n)/estimate=true"/>
 		<connect srcID="Tree.t:$(n)" targetID="state" inputName="stateNode" if="inposterior(Tree.t:$(n)) and Tree.t:$(n)/estimate=true"/>
 
-		<connect srcID="treePrior.t:$(n)" targetID="PopSizeGibbsSampler" inputName="gene" if="inposterior(treePrior.t:$(n)) and treePrior.t:$(n)/estimate=true"/>
+		<connect srcID="treePrior.t:$(n)" targetID="PopSizeGibbsSampler" inputName="gene" if="inposterior(Tree.t:$(n)) and Tree.t:$(n)/estimate=true"/>
+
+			<!--
+			Since "Species" is treated as a separate partition, we need the following
+			hacks to get rid of undesirable Tree.t:Species connections
+			-->
+			<connect srcID="Tree.t:Species" targetID="coordinatedUniform.t:Species" inputName="geneTree" if="Tree.t:$(n)/estimate=XXX"/>
+			<connect srcID="Tree.t:Species" targetID="coordinatedExponential.t:Species" inputName="geneTree" if="Tree.t:$(n)/estimate=XXX"/>
+			<connect srcID="Tree.t:Species" targetID="SBI" inputName="geneTree" if="Tree.t:$(n)/estimate=XXX"/>
+			<!-- end hacks -->
+
+
 
                         <mergepoint id="aux-partitiontemplate"/>
                 </partitiontemplate>
@@ -276,7 +293,7 @@
 			</taxonset>
         </tree>
         <parameter id="birthRate.t:Species" lower="0.0" name="stateNode" upper="500.0" value="213.2021"/>
-        <parameter id="popMean" name="stateNode" value="1"/>
+        <parameter id="popMean" name="stateNode" value="1" lower="0"/>
     </state>
 
     <distribution id="posterior" spec="util.CompoundDistribution">
@@ -285,8 +302,8 @@
                 <parameter id="popSizeTop" name="topPopSize" value="1"/>
 	 		    <populationModel id="speciesTreePopulationModel" spec="starbeast3.evolution.speciation.ConstantPopulations" populationSizes="@popSize" speciesTree="@Tree.t:Species" />
             </distribution>
-
         </distribution>
+
         <distribution id="prior" spec="util.CompoundDistribution">
             <distribution birthDiffRate="@birthRate.t:Species" id="YuleModel.t:Species" spec="beast.evolution.speciation.YuleModel" tree="@Tree.t:Species"/>
             <prior id="YuleBirthRatePrior.t:Species" name="distribution" x="@birthRate.t:Species">
@@ -297,6 +314,7 @@
                 <OneOnX id="OneOnX.01" name="distr"/>
             </prior>
         </distribution>
+
         <distribution id="likelihood" spec="util.CompoundDistribution" useThreads="true">
 
 		</distribution>