Skip to content

Commit dbe5363

Browse files
rbouckaertrbouckaert
committed
debug BetaApproximationLikelihood
1 parent 8ab4d5d commit dbe5363

File tree

2 files changed

+107
-20
lines changed

2 files changed

+107
-20
lines changed

‎examples/testBetaApprox.xml

Lines changed: 93 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,93 @@
1+
<!-- Simulation test 1: species = 4, lineages/species = 2, samples = 1000
2+
u=v=1.0 alpha=2, beta=0.2 lambda=10
3+
-->
4+
5+
<beast version='2.0' namespace='snap:beast.util:beast.core.util:beast.evolution'>
6+
7+
<map name='snapprior'>snap.likelihood.SnAPPrior</map>
8+
<map name='snaptreelikelihood'>snap.likelihood.BetaApproximationLikelihood</map>
9+
10+
<!-- n = 1000 -->
11+
<data spec='snap.Data' id='alignment' dataType='integerdata' statecount='3'>
12+
<sequence taxon='A' totalcount='3'>
13+
2,2,2,2,2,0,0,0,2,2,2,2,2,0,0,0,0,2,0,2,2,0,2,0,0,0,2,2,0,2,2,2,0,0,0,0,0,2,0,2,0,0,2,2,2,0,0,2,0,2,0,2,0,0,0,0,2,0,2,0,2,0,0,0,0,2,0,2,2,0,0,2,2,0,2,2,2,2,2,2,0,0,2,0,0,2,0,2,2,0,2,0,2,2,0,2,0,0,2,0,2,2,0,0,0,2,2,2,0,2,2,2,0,2,2,0,0,0,0,2,0,2,0,2,2,2,2,2,0,0,0,0,2,2,2,2,2,2,0,0,2,2,0,2,0,2,0,0,2,0,0,2,2,0,0,0,0,0,2,2,2,2,2,0,0,0,2,0,2,0,0,0,0,0,0,0,2,2,2,2,0,0,0,2,2,2,0,0,0,0,2,2,2,0,0,2,0,2,0,0,2,2,0,0,2,0,0,0,2,2,0,0,2,0,0,2,2,0,0,0,2,0,0,0,0,0,2,0,0,2,0,0,2,2,2,0,2,0,0,2,0,2,2,2,2,0,2,2,2,2,2,2,0,2,0,2,2,2,2,0,0,2,2,2,0,0,2,2,2,0,2,0,0,0,0,2,0,2,0,2,2,2,0,0,0,0,2,0,2,2,2,0,0,2,2,2,2,2,0,2,2,0,0,0,0,0,0,0,2,2,2,2,2,0,0,0,2,0,2,0,0,2,0,0,2,0,2,2,2,0,0,2,0,0,0,0,2,2,2,2,2,2,0,0,2,2,0,2,0,2,2,0,0,0,2,2,0,0,2,0,0,2,2,0,0,2,2,2,2,2,0,0,2,2,2,2,0,0,2,0,2,0,2,2,0,0,2,2,2,0,2,0,2,2,2,2,0,0,2,2,0,2,2,0,0,0,2,2,2,0,0,2,0,0,0,0,2,0,0,0,2,0,2,0,2,0,0,0,2,2,2,0,2,0,2,0,0,2,0,2,0,0,2,2,2,0,0,2,2,2,2,2,2,0,0,2,0,0,2,0,0,0,2,2,2,0,0,0,2,0,2,2,2,0,0,0,2,0,0,0,2,0,0,0,0,2,0,0,0,2,2,2,0,2,0,2,2,0,2,0,0,0,2,2,0,2,2,2,0,2,2,2,0,2,2,0,2,0,2,0,0,2,2,0,0,0,0,2,2,2,2,0,2,0,0,2,0,0,0,2,0,0,2,2,0,2,2,0,2,0,1,0,2,0,0,0,2,2,2,0,2,0,2,2,2,0,2,0,2,0,2,2,0,2,2,2,0,2,2,0,0,0,0,2,2,0,0,0,0,0,0,2,2,0,0,2,0,0,1,2,2,0,2,0,0,0,0,2,0,2,2,0,0,2,0,0,0,0,2,2,2,0,2,0,0,0,0,0,0,2,2,0,0,0,2,0,2,0,2,0,2,0,2,0,0,2,0,0,2,0,2,0,0,2,0,2,2,0,2,0,0,2,0,0,0,2,2,0,0,0,0,2,2,0,0,2,0,2,2,0,0,0,2,0,0,2,2,0,2,0,0,2,0,0,0,0,2,2,2,0,0,0,0,0,0,0,0,0,2,2,2,2,2,2,0,2,2,0,2,0,0,0,0,0,2,2,0,0,2,0,0,2,2,0,0,0,0,2,0,2,2,0,0,2,0,2,2,2,2,2,2,2,0,0,0,0,2,0,0,2,2,2,2,2,2,0,2,2,2,0,2,0,0,0,0,2,2,0,0,2,2,0,0,2,0,0,0,0,0,2,2,0,0,0,0,0,0,0,2,0,0,0,0,0,0,2,2,0,2,0,0,0,2,0,2,0,0,2,2,0,0,0,0,2,2,2,0,2,2,0,2,0,2,0,2,0,0,2,0,2,2,0,0,0,0,0,2,0,0,0,2,2,0,2,2,0,0,2,2,0,2,2,2,2,0,2,0,2,0,2,0,0,2,0,2,0,2,0,2,2,0,0,0,0,0,2,0,2,2,0,2,0,2,2,2,2,0,0,2,0,2,2,0,0,0,0,2,2,0,0,0,0,2,0,2,2,2,2,2,2,2,2,0,0,0,0,0,2,2,0,2,2,0,2,2,0,0,0,2,0,0,2,0,2,2,0,0,2,2,2,0,2,2,0,2,2,0,0,2,0,2,2,2,0,0,2,2,0,2,2,2,2,2,1,2,2,0,2,2,0,0,2,0,2,0,2,2,2,2,2,0,0,0,0,0,2,2,0,2,2,
14+
</sequence>
15+
<sequence taxon='B' totalcount='3'>
16+
2,2,2,2,2,0,0,0,2,2,2,2,2,0,0,0,0,2,2,2,2,0,2,0,0,0,2,2,0,2,2,2,0,0,0,0,0,2,0,2,0,0,2,2,2,0,0,2,0,2,0,2,0,0,0,0,2,0,2,0,2,0,0,0,2,2,0,2,2,0,0,2,2,0,2,2,2,2,2,2,0,0,2,0,0,2,0,2,2,0,2,0,2,0,0,2,0,0,2,0,2,2,0,0,0,2,2,2,0,2,2,2,0,2,2,0,0,0,0,2,0,2,0,2,2,2,2,2,2,0,0,0,2,2,2,2,2,2,0,0,2,2,0,2,0,2,0,0,2,0,0,2,2,0,0,0,0,0,2,2,2,2,2,0,0,0,2,0,2,0,0,0,0,0,0,0,2,2,2,2,0,0,0,2,2,2,0,0,0,0,2,2,2,0,0,2,0,2,0,0,2,2,0,0,2,0,0,0,2,2,0,0,2,0,0,2,2,0,0,0,0,0,0,0,0,0,2,0,0,2,0,0,2,2,2,0,2,0,0,2,0,2,2,2,2,0,2,2,2,2,2,2,0,2,0,2,2,2,2,0,0,2,2,2,0,0,2,2,2,0,2,0,0,0,0,2,0,2,0,2,2,2,0,0,0,0,2,0,2,2,0,0,0,2,2,2,2,0,0,2,2,0,0,0,0,0,0,0,2,2,2,2,0,0,0,0,2,0,2,0,0,2,0,0,2,2,2,2,2,0,0,2,0,0,0,0,2,2,2,2,2,2,0,0,2,2,0,2,0,2,2,0,0,2,2,2,0,0,2,0,0,2,2,0,0,2,2,2,2,2,0,0,2,2,2,2,0,0,2,0,2,2,2,2,0,0,2,2,2,0,2,0,2,2,0,2,0,0,2,2,0,2,2,0,0,0,2,2,2,0,0,2,0,0,0,0,2,0,0,0,2,0,2,0,0,0,0,0,2,2,2,0,2,0,2,0,0,2,0,2,0,0,2,2,2,0,0,2,2,2,2,2,2,2,0,2,0,0,2,2,0,0,2,2,2,0,0,0,2,0,2,2,2,0,0,0,2,0,0,0,2,0,0,0,0,2,2,0,0,2,2,2,2,2,2,2,2,0,2,0,0,0,2,2,0,2,2,2,0,2,2,0,0,2,2,0,2,0,2,0,0,2,2,0,0,0,0,2,2,2,2,0,2,0,0,2,0,0,0,2,0,0,2,0,0,2,2,0,0,0,0,0,2,0,0,0,2,2,2,0,2,0,2,2,2,0,2,0,2,0,2,2,0,2,2,2,0,2,0,0,0,0,2,2,2,0,0,0,0,0,0,2,2,0,0,2,0,0,0,2,2,0,2,0,0,0,0,2,0,2,2,0,0,2,0,0,0,0,2,2,2,0,2,0,0,0,0,0,0,2,2,0,0,0,2,0,2,2,2,0,2,0,2,0,0,2,0,0,0,0,2,0,0,2,0,2,2,0,0,0,0,2,0,0,0,2,2,0,0,0,0,2,2,0,0,2,0,2,2,0,0,0,2,0,0,2,2,0,2,0,0,2,0,0,0,0,2,2,2,0,0,0,0,0,0,0,0,0,2,2,2,2,2,2,0,2,2,0,2,0,0,0,0,0,2,2,0,0,2,0,0,2,2,0,0,0,0,2,0,2,0,0,0,0,0,0,2,2,2,2,2,2,0,0,0,0,0,0,0,2,2,2,2,2,2,0,2,2,2,2,2,0,0,0,0,2,2,0,2,2,2,0,0,2,0,0,0,0,0,2,2,0,0,0,0,0,0,0,2,0,0,0,0,0,0,2,2,0,2,0,0,0,2,0,2,0,0,2,2,1,0,0,0,2,2,2,0,2,2,0,0,0,2,0,2,0,0,2,0,0,2,0,0,0,0,0,2,0,0,0,2,2,0,2,2,0,0,0,0,0,2,2,2,2,0,2,0,2,2,2,0,0,2,0,2,0,2,0,2,2,0,0,0,0,0,2,0,2,2,0,2,0,2,2,2,2,0,0,2,0,2,2,0,0,0,0,2,2,0,0,0,0,2,0,2,2,2,2,2,2,2,2,0,0,0,0,0,2,2,0,2,2,0,2,2,0,0,0,2,0,0,2,0,2,2,0,0,2,2,2,0,2,2,0,2,2,0,0,2,0,2,2,2,0,0,2,2,0,2,2,2,2,2,0,0,2,0,2,2,0,1,2,0,2,0,2,2,0,2,2,0,1,0,0,0,2,2,0,2,2,
17+
</sequence>
18+
<sequence taxon='C' totalcount='3'>
19+
1,2,2,2,0,0,0,2,2,2,2,2,2,0,0,0,0,2,2,2,2,0,2,0,0,0,2,2,0,2,2,0,0,0,0,0,0,2,0,2,0,0,2,2,0,0,0,2,0,2,0,2,0,0,0,0,2,0,2,0,2,0,0,0,0,2,0,2,0,0,0,2,2,0,2,2,2,2,2,2,0,0,2,2,0,2,0,2,2,0,2,0,2,2,0,2,0,0,2,0,2,2,0,0,0,2,2,2,0,0,2,2,0,2,2,0,0,0,0,2,0,2,0,2,2,2,2,2,0,0,0,0,2,2,2,2,0,2,0,0,2,2,0,0,0,2,2,0,2,0,0,2,2,0,0,0,0,0,2,0,2,2,2,0,0,0,2,0,2,0,0,0,0,2,0,0,2,2,2,2,0,0,0,2,2,2,0,0,2,0,2,2,2,0,0,2,0,2,0,0,0,2,0,0,0,0,0,0,0,2,0,0,2,2,0,2,2,0,0,0,0,2,0,0,0,0,2,0,0,2,0,0,2,2,2,0,0,2,0,2,0,2,2,2,2,0,2,2,2,2,0,2,0,0,0,2,2,2,2,0,0,2,0,2,0,0,2,2,2,0,2,2,0,0,2,2,0,2,0,2,2,2,0,0,0,0,2,0,2,2,2,0,0,2,2,2,2,2,0,0,2,0,2,0,0,0,0,0,2,2,2,2,0,0,0,0,2,0,2,0,0,2,0,2,2,2,2,2,2,0,0,2,0,0,0,0,2,2,2,2,2,2,0,0,2,2,0,2,0,2,2,0,0,0,2,2,0,0,2,0,2,2,2,0,0,2,2,2,2,2,0,0,2,2,2,2,2,0,2,0,2,2,2,2,0,0,2,2,0,0,2,0,2,0,2,2,0,0,2,2,0,2,2,0,0,0,2,2,2,0,0,2,0,0,0,0,2,0,0,0,2,0,2,0,0,0,0,0,2,2,2,0,2,0,2,0,0,2,0,2,0,0,2,2,2,0,0,2,0,2,2,2,0,2,0,0,0,0,2,0,0,0,2,2,2,0,0,0,2,0,2,2,2,2,0,0,2,0,0,0,2,0,0,0,0,2,2,0,0,0,2,2,2,2,2,2,2,0,2,0,0,0,2,2,0,2,2,2,0,2,2,0,0,0,2,0,2,0,2,0,0,2,2,0,0,0,0,2,2,2,2,0,2,0,0,2,0,0,0,2,0,0,2,2,0,2,2,2,2,0,0,0,2,0,0,0,2,2,2,0,2,2,2,2,2,0,2,0,2,0,2,2,0,2,2,2,0,2,2,0,0,0,2,2,0,0,0,0,0,0,0,2,2,0,0,2,0,0,0,2,0,0,2,0,0,0,0,2,0,2,2,0,0,2,0,0,0,0,2,2,2,2,2,0,0,0,0,0,0,2,2,0,0,0,2,0,2,2,2,0,2,0,2,0,0,2,0,0,2,0,2,0,0,2,0,2,2,0,2,0,0,2,0,0,0,2,2,0,0,0,0,2,2,0,2,2,0,2,2,0,0,0,2,0,0,2,2,0,2,0,0,2,0,0,0,0,2,0,2,0,0,0,0,0,0,0,0,0,2,2,2,0,2,2,0,2,2,2,2,0,0,0,0,0,2,2,0,0,2,0,0,2,2,0,0,0,0,2,2,2,0,0,0,0,0,0,2,2,2,2,2,2,0,0,0,0,0,0,0,2,2,2,2,2,2,0,2,2,2,2,2,0,0,2,2,0,2,0,0,2,2,0,0,2,0,0,0,0,0,2,2,0,0,0,0,0,2,0,2,0,0,0,0,0,0,2,0,0,2,0,2,0,2,0,2,0,0,2,2,0,0,0,2,2,2,2,0,2,2,0,0,0,0,0,2,0,0,2,0,0,2,0,0,0,2,0,2,0,0,0,2,2,0,2,2,0,0,2,0,0,2,2,2,2,0,2,0,2,0,2,0,0,2,0,2,0,2,0,2,2,0,0,0,0,0,2,0,2,2,0,2,0,2,2,2,2,0,0,0,0,2,2,0,0,0,2,2,2,0,0,0,0,2,0,2,2,2,2,2,2,2,0,0,0,2,0,0,2,2,0,2,2,0,2,2,0,0,0,2,0,0,2,0,2,2,0,0,2,2,0,0,0,2,0,2,2,2,0,2,0,0,2,2,0,0,2,2,0,2,2,2,0,2,0,0,2,0,2,2,0,0,2,0,2,0,2,2,0,2,2,0,0,0,0,0,2,2,0,2,2,
20+
</sequence>
21+
<sequence taxon='D' totalcount='3'>
22+
2,2,2,2,0,0,0,0,2,2,2,2,2,0,0,0,0,2,2,2,2,0,2,0,0,0,2,2,0,2,2,0,0,0,0,0,0,2,0,2,0,0,2,2,0,0,0,2,0,2,0,2,0,0,0,0,2,0,2,0,2,0,0,0,0,2,0,2,0,0,0,2,2,0,2,2,2,2,2,2,0,0,2,2,0,2,0,2,2,0,2,0,2,2,0,2,0,0,2,0,2,2,0,0,0,2,2,2,0,0,2,2,0,2,2,0,0,0,1,2,0,2,0,2,2,2,2,2,0,0,0,0,2,2,2,2,0,2,0,0,2,2,0,2,0,2,0,0,2,0,0,2,2,0,0,0,0,0,2,0,2,2,2,0,0,0,2,0,2,0,0,0,0,2,0,0,2,2,2,2,1,0,0,2,2,2,0,0,2,0,2,2,2,0,0,2,0,2,0,0,1,2,0,0,0,0,0,0,0,2,0,0,2,2,0,2,2,0,0,0,0,2,0,0,0,0,2,0,0,2,0,0,2,2,2,0,0,2,0,2,0,2,2,2,2,0,2,2,2,2,0,2,0,0,0,2,2,2,2,0,0,2,0,2,0,0,2,2,2,0,2,2,0,0,2,2,0,2,0,2,2,2,0,0,0,0,2,0,2,2,2,0,0,2,2,2,2,2,0,0,2,0,2,0,0,0,0,0,2,2,2,2,0,0,0,0,2,0,2,0,0,2,0,2,2,2,2,2,2,0,0,2,0,0,0,0,2,2,2,2,2,2,0,0,2,2,0,2,0,2,2,0,0,0,2,2,0,0,2,0,2,2,2,0,0,2,2,2,2,2,0,0,2,2,2,2,0,0,2,0,2,2,2,2,0,0,2,2,0,0,2,0,2,0,2,2,0,0,2,2,0,2,2,0,0,0,2,2,2,0,0,2,0,0,0,0,2,0,0,0,2,0,2,0,0,0,0,0,2,2,2,0,2,0,2,0,0,2,2,2,0,0,2,2,2,0,0,2,0,2,2,2,0,2,0,0,0,0,2,0,0,0,2,2,2,0,0,0,2,0,2,2,2,2,0,0,2,0,0,0,2,0,0,0,0,2,2,0,0,2,2,2,2,2,2,2,2,0,2,0,0,0,2,2,0,2,0,2,0,2,2,0,0,0,2,0,2,0,2,0,0,2,2,0,0,0,0,2,2,2,2,0,2,0,1,2,0,0,0,2,0,0,2,2,0,2,2,2,2,0,0,0,2,0,0,0,2,2,2,0,2,2,2,2,2,0,2,0,2,0,2,2,0,2,2,2,0,2,2,0,0,0,2,2,0,0,0,0,0,0,0,2,2,0,0,2,0,0,0,2,0,0,2,0,0,0,0,2,0,2,2,0,0,2,0,0,0,0,2,1,2,2,2,0,0,0,0,0,0,2,2,0,0,0,2,0,2,2,2,0,2,0,2,0,0,2,0,0,2,0,2,0,0,2,0,2,2,0,2,0,0,2,0,0,0,2,2,0,0,0,0,2,2,0,2,2,0,2,2,0,0,0,2,0,0,2,2,0,2,0,0,2,0,0,0,0,2,0,2,0,0,0,0,0,0,0,0,0,2,2,2,2,2,2,0,2,2,2,2,0,0,0,0,0,2,2,0,0,2,0,0,2,2,0,0,0,0,2,2,2,0,0,0,0,0,0,2,2,2,2,2,2,0,0,0,0,0,0,0,2,2,2,2,2,2,0,2,2,2,2,2,0,0,2,2,0,2,0,0,2,2,0,0,2,0,0,0,0,0,2,2,0,0,1,0,0,2,0,2,0,0,0,0,0,0,2,0,0,2,0,2,0,2,0,2,0,0,2,2,0,0,0,2,2,2,2,0,2,2,0,0,0,2,0,2,0,0,1,0,0,2,0,0,0,2,0,2,0,2,1,2,2,0,2,2,0,0,2,0,0,2,2,2,2,0,2,0,2,0,2,0,0,2,0,2,0,2,0,2,2,0,0,0,0,0,2,0,2,2,0,2,0,2,2,2,2,0,0,0,0,2,2,0,0,0,2,2,2,0,0,0,0,2,0,2,2,2,2,2,2,2,0,0,0,2,0,0,2,2,0,2,2,0,2,2,1,0,0,2,0,0,2,0,2,2,0,0,2,2,0,0,2,2,0,2,2,2,0,1,0,0,2,2,0,0,2,2,0,2,2,2,0,2,0,0,2,0,2,2,0,0,2,0,2,0,2,2,0,2,2,0,0,0,0,0,2,2,0,2,2,
23+
</sequence>
24+
25+
26+
</data>
27+
28+
<run id="mcmc" spec='beast.core.MCMC' chainLength="1000000" preBurnin="0">
29+
<state storeEvery='1000'>
30+
<parameter name='stateNode' id='coalescenceRate' dimension='7' value='10'/>
31+
32+
<!-- the following values may or may not be appropriate. Read the manual for further directions -->
33+
<parameter name='stateNode' id='v' value='1.0' lower='0.0'/>
34+
<parameter name='stateNode' id='u' value='1.0' lower='0.0'/>
35+
<parameter name='stateNode' id='alpha' value='2' lower='0.0'/>
36+
<parameter name='stateNode' id='beta' value='200' lower='0.0'/>
37+
<parameter name='stateNode' id='lambda' value='10' lower='0.0'/>
38+
<stateNode clusterType="upgma" estimate="true" id="tree" nodetype="snap.NodeData" spec="beast.util.ClusterTree" taxa='@alignment'/>
39+
</state>
40+
41+
<distribution id="posterior" spec='CompoundDistribution'>
42+
<snapprior name='distribution' id='prior' alpha='@alpha' beta='@beta' lambda='@lambda' tree='@tree' coalescenceRate='@coalescenceRate'/>
43+
<snaptreelikelihood name='distribution' id='treeLikelihood' data='@alignment' tree='@tree' coalescenceRate='@coalescenceRate'/>
44+
</distribution>
45+
46+
<!--
47+
<stateDistribution idref='prior'/>
48+
-->
49+
50+
<operator spec='operators.NodeSwapper' weight='0.5' tree='@tree'/>
51+
<operator spec='operators.NodeBudger' weight='5' size='0.5' tree='@tree'/>
52+
<operator spec='operators.ScaleOperator' scaleFactor='0.25' weight='1' tree='@tree'/>
53+
<operator spec='operators.GammaMover' scale='0.5' weight='0.5' coalescenceRate='@coalescenceRate'/>
54+
<operator spec='operators.RateMixer' scaleFactors='0.25' weight='1' tree='@tree' coalescenceRate='@coalescenceRate'/>
55+
56+
<!-- estimate u and v -->
57+
<operator spec='operators.MutationMover' window='0.1' weight='0.25' u='@u' v='@v'/>
58+
59+
60+
<logger logEvery="1000">
61+
<log idref="u"/>
62+
<log idref="v"/>
63+
<log idref="prior"/>
64+
<log idref="treeLikelihood"/>
65+
<log idref="posterior"/>
66+
<log spec='snap.ThetaLogger' coalescenceRate='@coalescenceRate'/>
67+
<log spec='beast.evolution.tree.TreeHeightLogger' tree='@tree'/>
68+
</logger>
69+
<logger logEvery="1000" fileName="testBetaApproximation.log">
70+
<model idref='posterior'/>
71+
<log idref="u"/>
72+
<log idref="v"/>
73+
<log idref="prior"/>
74+
<log idref="treeLikelihood"/>
75+
<log idref="posterior"/>
76+
<log idref='coalescenceRate'/>
77+
<log spec='snap.ThetaLogger' coalescenceRate='@coalescenceRate'/>
78+
<log spec='beast.evolution.tree.TreeHeightLogger' tree='@tree'/>
79+
<log spec='TreeLengthLogger' tree='@tree'/>
80+
</logger>
81+
<logger logEvery="1000" fileName="testBetaApproximation.trees" mode='tree'>
82+
<log spec='beast.evolution.tree.TreeWithMetaDataLogger' tree="@tree" metadata='@coalescenceRate'/>
83+
</logger>
84+
85+
<!--
86+
<logger logEvery="1000" fileName="test2.$(seed).trees" mode='tree'>
87+
<log spec='snap.CoalescentUnitTreeLogger' tree="@tree" coalescenceRate='@coalescenceRate'/>
88+
</logger>
89+
-->
90+
</run>
91+
92+
</beast>
93+

‎src/snap/likelihood/BetaApproximationLikelihood.java

Lines changed: 14 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -1,32 +1,31 @@
11
package snap.likelihood;
22

33

4+
45
import snap.Data;
56
import beast.core.Citation;
67
import beast.core.Description;
7-
import beast.evolution.alignment.Alignment;
8-
import beast.evolution.branchratemodel.BranchRateModel;
9-
import beast.evolution.branchratemodel.StrictClockModel;
8+
import beast.core.Input;
9+
import beast.core.Input.Validate;
10+
import beast.core.parameter.RealParameter;
1011
import beast.evolution.likelihood.GenericTreeLikelihood;
11-
import beast.evolution.sitemodel.SiteModel;
1212
import beast.evolution.tree.Node;
1313
import beast.evolution.tree.Tree;
1414
import beast.math.Binomial;
1515
import beast.math.GammaFunction;
1616

17-
import org.apache.commons.math.distribution.BetaDistribution;
1817

1918
@Description("Implementation of Siren et al. Beta SNP approximation using Hiscott et al. integrator")
2019
@Citation("")
2120
public class BetaApproximationLikelihood extends GenericTreeLikelihood {
21+
public Input<RealParameter> coalescenceRateInput = new Input<RealParameter>("coalescenceRate", "population size parameter with one value for each node in the tree");
2222

23+
public BetaApproximationLikelihood() {
24+
siteModelInput.setRule(Validate.OPTIONAL);
25+
}
2326

2427
Tree tree;
25-
<<<<<<< HEAD
26-
Double rootAlpha; //Parameter for root distribution.
27-
=======
2828
Double rootTheta; //Parameter for root distribution.
29-
>>>>>>> 998524460c43138ca5124e4ba57b30d99803bf95
3029
Data data;
3130

3231
int nPoints; //Size of mesh used in integration
@@ -55,21 +54,20 @@ public class BetaApproximationLikelihood extends GenericTreeLikelihood {
5554
*/
5655
protected double[] m_branchLengths;
5756
protected double[] storedBranchLengths;
57+
58+
RealParameter coalescenceRate;
5859

5960

6061
@Override
6162
public void initAndValidate() throws Exception {
6263
tree = (Tree) treeInput.get();
63-
data = dataInput.get();
64+
data = (Data) dataInput.get();
6465
patternLogLikelihoods = new double[data.getPatternCount()];
66+
coalescenceRate = coalescenceRateInput.get();
6567

66-
nPoints = 20; //TODO: read this from the xml
68+
nPoints = data.getTaxonCount();//20; //TODO: read this from the xml
6769
int nNodes = tree.getNodeCount();
6870
partialIntegral = new double[nPoints][nNodes];
69-
<<<<<<< HEAD
70-
=======
71-
rootTheta = 0.5; //We should get this from the PARAMETERS.
72-
>>>>>>> 998524460c43138ca5124e4ba57b30d99803bf95
7371

7472
hasDirt = Tree.IS_FILTHY;
7573
}
@@ -82,7 +80,7 @@ public void initAndValidate() throws Exception {
8280
*/
8381
@Override
8482
public double calculateLogP() {
85-
83+
rootTheta = coalescenceRate.getValue(tree.getRoot().getNr()); //We should get this from the PARAMETERS.
8684

8785
double logL=0.0;
8886

@@ -239,11 +237,7 @@ private double transProb(double x,double y, double t) {
239237
}
240238

241239
private double rootProb(double x) {
242-
<<<<<<< HEAD
243-
return betaDensity(x,rootAlpha,rootAlpha);
244-
=======
245240
return betaDensity(x,rootTheta,rootTheta);
246-
>>>>>>> 998524460c43138ca5124e4ba57b30d99803bf95
247241
}
248242

249243

0 commit comments

Comments
 (0)