Skip to content

Commit

Permalink
Merge branch 'hmc-clock' of https://github.com/beast-dev/beast-mcmc i…
Browse files Browse the repository at this point in the history 
…nto hmc-clock
  • Loading branch information
@xji3
xji3 committed Dec 6, 2023
2 parents 88e85b2 + 8126ce9 commit c81ca40
Show file tree
Hide file tree
Showing 26 changed files with 2,521 additions and 95 deletions.
199 changes: 172 additions & 27 deletions ci/TestXML/testGamGpLikelihoodAndGradient.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,7 @@
<state code="C"/>
<!--
<state code="D"/>
-->
-->
</generalDataType>

<attributePatterns id="loc.pattern" attribute="loc">
Expand Down Expand Up @@ -92,7 +92,112 @@
<parameter id="loc.clock.rate" value="1E-4" lower="0.0"/>
</rate>
</strictClockBranchRates>

<gamGpModel id="gamGpModel"> <!-- returns density of realized field-->
<realizedField>
<parameter id="loc.coefficients0" value="1 0 -1 0 0 0"/> <!-- 1 2 3 4 5put other numbers-->
</realizedField>
<gaussianNoise>
<parameter id="loc.noise" value="1"/>
</gaussianNoise>
<kernel>
<parameter id="loc.kernel" value="1 3 4" lower="0.0"/>
</kernel>
<independentVariables>
<designMatrix id="loc.designMatrix0">
<parameter value="1 0 0
0 0 0"/>

</designMatrix>
</independentVariables>

<independentVariables>
<designMatrix id="loc.designMatrix1">
<parameter value="0 1 0
0 0 0"/>
</designMatrix>
</independentVariables>
<independentVariables>
<designMatrix id="loc.designMatrix2">
<parameter value="0 0 1
0 0 0"/>
</designMatrix>
</independentVariables>
<independentVariables>
<designMatrix id="loc.designMatrix3">
<parameter value="0 0 0
1 0 0"/>
</designMatrix>
</independentVariables>
<independentVariables>
<designMatrix id="loc.designMatrix4">
<parameter value="0 0 0
0 1 0"/>
</designMatrix>
</independentVariables>

<independentVariables>
<designMatrix id="loc.designMatrix5">
<parameter value="0 0 0
0 0 1"/>
</designMatrix>
</independentVariables>

<!-- nothing -->
</gamGpModel>


<randomField id="test.field">
<data>
<parameter idref="loc.coefficients0"/>
</data>
<distribution>
<gaussianProcessField id="test.gp" dim="6">
<orderVariance>
<parameter id="orderVariance" value="1.0"/>
</orderVariance>
<gaussianNoise>
<parameter idref="loc.noise"/>
</gaussianNoise>
<basis>
<designMatrix idref="loc.designMatrix0"/>
<kernel type="dotProduct"/>
</basis>
<basis>
<designMatrix idref="loc.designMatrix1"/>
<kernel type="dotProduct"/>
</basis>
<basis>
<designMatrix idref="loc.designMatrix2"/>
<kernel type="dotProduct"/>
</basis>
<basis>
<designMatrix idref="loc.designMatrix3"/>
<kernel type="dotProduct"/>
</basis>
<basis>
<designMatrix idref="loc.designMatrix4"/>
<kernel type="dotProduct"/>
</basis>
<basis>
<designMatrix idref="loc.designMatrix5"/>
<kernel type="dotProduct"/>
</basis>
</gaussianProcessField>
</distribution>
</randomField>

<randomFieldGradient id="test.gradient">
<randomField idref="test.field"/>
<parameter idref="loc.coefficients0"/>
</randomFieldGradient>

<report>
<gamGpModel idref="gamGpModel"/>
<randomField idref="test.field"/>
<randomFieldGradient idref="test.gradient"/>
</report>

<glmSubstitutionModel id="loc.model" normalize="true">
<generalDataType idref="loc.dataType"/>
<rootFrequencies>
Expand All @@ -103,23 +208,10 @@
</frequencies>
</frequencyModel>
</rootFrequencies>
<gamGpModel id="gamGpModel">
<independentVariables>
<parameter id="loc.coefficients0" value="0 0 0 0 0 0"/>
<designMatrix id="loc.designMatrix0">
<parameter value="1 1 1
1 1 1"/>
</designMatrix>
</independentVariables>

<independentVariables>
<parameter id="loc.coefficients2" value="0 0 0 0 0 0"/>
<designMatrix id="loc.designMatrix2">
<parameter value="0 1 0
0 0 0"/>
</designMatrix>
</independentVariables>
</gamGpModel>
<!-- <gamGpModel idref="gamGpModel"/> -->
<logRates>
<parameter idref="loc.coefficients0"/>
</logRates>
</glmSubstitutionModel>

<siteModel id="loc.siteModel">
Expand All @@ -137,24 +229,77 @@
<strictClockBranchRates idref="loc.branchRates"/>
</treeDataLikelihood>

<!-- START BY UNCOMMENTING THE REPORT BELOW -->

<!--

<report>
<gamGpModel idref="gamGpModel"/>
<treeDataLikelihood idref="treeLikelihood"/>
</report>
-->

<!--
<gamGpSubstitutionModelGradient id="gradient1" traitName="loc" effects="fixed">
<approximateLogCtmcRateGradient id="gradient1" traitName="loc">
<treeDataLikelihood idref="treeLikelihood"/>
<glmSubstitutionModel idref="loc.model"/>
</gamGpSubstitutionModelGradient>
</approximateLogCtmcRateGradient>

<report>
<glmSubstitutionModelGradient idref="gradient1"/>
<approximateLogCtmcRateGradient idref="gradient1"/>
</report>
-->

<jointGradient id="joint.gradient">
<approximateLogCtmcRateGradient idref="gradient1"/>
<randomFieldGradient idref="test.gradient"/>
</jointGradient>

<report>
Joint gradient
<jointGradient idref="joint.gradient"/>
</report>

<operators id="operators" optimizationSchedule="log">
<randomWalkOperator windowSize="0.1" weight="1">
<parameter idref="loc.coefficients0"/>
</randomWalkOperator>
<!-- <scaleOperator scaleFactor="0.75" weight="1">-->
<!-- <parameter idref="loc.kernel"/>-->
<!-- </scaleOperator>-->
</operators>
<!-- scale operator random walk positive CHECK! -->

<mcmc id="mcmc" chainLength="100000" autoOptimize="true">
<joint id="joint"> <!-- joint, from now on a list of density -->
<prior id="prior">
<!-- <gamGpModel idref="gamGpModel"/> -->
<randomField idref="test.field"/>
</prior>
<likelihood id="likelihood">
<treeLikelihood idref="treeLikelihood"/>
</likelihood>
</joint>
<operators idref="operators"/>
<log id="screenLog" logEvery="200"> <!-- keep track of the chain,printed to the screen-->
<column label="Joint" dp="4" width="12"> <!-- "label" is the col name; "dp" is the number of decimals and "width" is the btw cols space -->
<joint idref="joint"/>
</column>
<column label="Prior" dp="4" width="12">
<prior idref="prior"/>
</column>
<!-- <column label="Likelihood" dp="4" width="12">-->
<!-- <likelihood idref="likelihood"/>-->
<!-- </column>-->
<column label="Root Height" sf="6" width="12">
<parameter idref="treeModel.rootHeight"/>
</column>
<column label="Field" sf="6" width="12">
<parameter idref="loc.coefficients0"/>
</column>
</log>
<log id="fileLog" logEvery="10" fileName="test.log">
<joint idref="joint"/>
<prior idref="prior"/>
<!-- <likelihood idref="likelihood"/>-->
<parameter idref="loc.coefficients0"/>
</log>
</mcmc>

<traceAnalysis fileName="test.log" stdError="true"/>

</beast>
54 changes: 47 additions & 7 deletions ci/TestXML/testTimeVaryingRates.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,8 @@
<newick id="dummyTree" usingDates="false" usingHeights="true">
((A:1,B:1):1,C:3)
</newick>

<!--
<timeVaryingRates id="dummyRates1">
<newick idref="dummyTree"/>
<rate>
Expand Down Expand Up @@ -90,7 +91,9 @@
((A:[&#38;rate=1.0]1.0,B:[&#38;rate=1.0]1.0):[&#38;rate=1.0]1.0,C:[&#38;rate=1.0]3.0);
</expected>
</assertEqual>
-->

<!--
<alignment id="dummyAlignment" dataType="nucleotide">
<sequence>
<taxon idref="A"/>
Expand Down Expand Up @@ -186,6 +189,7 @@
<report>
<branchRateGradient idref="dummyGradient3"/>
</report>
-->

<!-- The list of taxa to be analysed (can also include dates/ages). -->
<!-- ntax=17 -->
Expand Down Expand Up @@ -372,16 +376,20 @@
<report>
<transformedVectorSumTransform idref="clock.rate"/>
</report>

<!--
<timeVaryingRates id="branchRates">
<treeModel idref="treeModel"/>
<rate>
<!-- <parameter id="clock.rate" value="1E-2 2E-2 3E-3 4E-2 5E-3 6E-4" lower="0.0"/>-->
<!~~ <parameter id="clock.rate" value="1E-2 2E-2 3E-3 4E-2 5E-3 6E-4" lower="0.0"/>~~>
<parameter idref="clock.rate"/>
</rate>
<gridPoints>
<parameter id="grid" value="100 200 300 400 500" />
</gridPoints>
</timeVaryingRates>
-->

<!-- The HKY substitution model (Hasegawa, Kishino & Yano, 1985) -->
<HKYModel id="hky">
Expand All @@ -403,8 +411,7 @@
<HKYModel idref="hky"/>
</substitutionModel>
</siteModel>


<!--
<bayesianBridge id="timeVaryingRatePrior">
<parameter idref="tmp1"/>
<globalScale>
Expand Down Expand Up @@ -438,7 +445,7 @@
<compoundGradient idref="grad.bdss.prior"/>
</report>
<!-- Likelihood for tree given sequence data -->
<!~~ Likelihood for tree given sequence data ~~>
<treeDataLikelihood id="treeLikelihood" useAmbiguities="false" usePreOrder="true">
<patterns idref="patterns"/>
<treeModel idref="treeModel"/>
Expand Down Expand Up @@ -468,10 +475,10 @@
<jointGradient id="jointGradient">
<gradientWrtIncrements1D id="test" incrementTransformType="log">
<branchRateGradient idref="gradient"/> <!-- wrt clock.rate -->
<branchRateGradient idref="gradient"/> <!~~ wrt clock.rate ~~>
<compoundParameter idref="clockRateIncrements"/>
</gradientWrtIncrements1D>
<compoundGradient idref="grad.bdss.prior"/> <!-- wrt log.inverse.scan [tmp0, tmp1] -->
<compoundGradient idref="grad.bdss.prior"/> <!~~ wrt log.inverse.scan [tmp0, tmp1] ~~>
</jointGradient>
<report id="jointGradientReport">
Expand All @@ -489,4 +496,37 @@
<report idref="jointGradientReport"/>
</expected>
</assertEqual>
-->

<!-- functionalForm="piecewiseLogConstant" -->

<transform id="test" type="log" dim="10"/>

<timeVaryingRates id="logBranchRates" >
<treeModel idref="treeModel"/>
<rate>
<parameter id="log.clock.rate" value="-3 -2 -1 0 1 2"/>
</rate>
<gridPoints>
<parameter id="log.grid" value="100 200 300 400 500" />
</gridPoints>
<transform id="test2" type="log"/>
</timeVaryingRates>

<treeDataLikelihood id="logTreeLikelihood" useAmbiguities="false" usePreOrder="true">
<patterns idref="patterns"/>
<treeModel idref="treeModel"/>
<siteModel idref="siteModel"/>
<timeVaryingRates idref="logBranchRates"/>
</treeDataLikelihood>

<branchRateGradient id="logGradient" traitName="Sequence" useHessian="false">
<treeDataLikelihood idref="logTreeLikelihood"/>
</branchRateGradient>

<report>
<branchRateGradient idref="logGradient"/>
</report>


</beast>
4 changes: 4 additions & 0 deletions src/dr/app/beast/development_parsers.properties
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,6 +80,9 @@ dr.inferencexml.distribution.RandomFieldGradientParser
dr.inferencexml.distribution.GaussianMarkovRandomFieldParser
dr.inferencexml.distribution.BayesianBridgeMarkovRandomFieldParser
dr.inferencexml.distribution.BaselineIncrementFieldParser
dr.inferencexml.distribution.GaussianProcessFieldParser
dr.inferencexml.distribution.GaussianProcessKernelParser
dr.inferencexml.distribution.GaussianProcessPredictionParser

# GAUSSIAN PROCESS
dr.evomodelxml.coalescent.GPSkytrackAnalysisParser
Expand Down Expand Up @@ -262,6 +265,7 @@ dr.evomodelxml.treelikelihood.thorneytreelikelihood.UniformSubtreePruneRegraftPa
dr.inferencexml.operators.MaskMoveOperatorParser
dr.evomodelxml.continuous.hmc.GlmSubstitutionModelGradientParser
dr.evomodelxml.continuous.hmc.GamGpSubstitutionModelGradientParser
dr.evomodelxml.continuous.hmc.ApproximateLogCtmcRateGradientParser
dr.evomodelxml.substmodel.BirthDeathSubstitutionModelParser

# Uncertain attributes:
Expand Down
3 changes: 2 additions & 1 deletion src/dr/evomodel/coalescent/GMRFSkygridLikelihood.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -541,7 +541,8 @@ protected void setupSufficientStatistics() {
}
currentGridIndex = minGridIndex;

lastCoalescentTime = currentTime + intervalsList.get(i).getTotalDuration();
lastCoalescentTime = // currentTime +
intervalsList.get(i).getTotalDuration();

// theLastTime = lastCoalescentTime;

Expand Down
4 changes: 4 additions & 0 deletions src/dr/evomodel/substmodel/GlmSubstitutionModel.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,10 @@ public GlmSubstitutionModel(String name, DataType dataType, FrequencyModel rootF

}

public LogAdditiveCtmcRateProvider getRateProvider() {
return glm;
}

public GeneralizedLinearModel getGeneralizedLinearModel() {
if (glm instanceof GeneralizedLinearModel) {
return (GeneralizedLinearModel) glm;
Expand Down
Loading

0 comments on commit c81ca40

Please sign in to comment.