Tests

sim_phyl

st <- sim_phyl(ct=5)
plot(st$newick,root.edge = T)
axisPhylo(backward = F,root.time = st$t[1])

miss <- st$L[st$L[,3]!=(-1),]
miss_spe <- miss$spec
time <- rbind(data.frame(brtimes = st$L[,2], E=1, spec = st$L[,1]),data.frame(brtimes = miss[,3], E=0, spec = miss[,1]))
time <- time[order(time$brtimes),]
time

##      brtimes E spec
## 1  0.0000000 1   aa
## 2  0.3334744 1   ab
## 3  0.3397169 1   ac
## 4  0.5256167 1   ad
## 5  0.8184213 1   ae
## 28 0.8808321 0   ab
## 6  0.9409274 1   ag
## 7  1.5194814 1   ah
## 29 1.7944359 0   ac
## 8  2.0255450 1   aj
## 9  2.2365991 1   ak
## 10 2.4714988 1   al
## 11 2.4885062 1   am
## 12 2.5404195 1   an
## 13 2.6038842 1   ao
## 14 2.7091034 1   ap
## 15 2.8988630 1   aq
## 16 2.9514716 1   ar
## 32 3.4082254 0   ak
## 17 3.4654121 1   at
## 18 3.5857570 1   au
## 19 3.6410367 1   av
## 20 3.7033597 1   aw
## 21 3.8538806 1   ax
## 33 3.9868879 0   ao
## 22 4.0000750 1   az
## 23 4.2353958 1   ba
## 24 4.4585866 1   bb
## 25 4.5703244 1   bc
## 31 4.6654420 0   ag
## 30 4.7666942 0   ae
## 26 4.8290354 1   bf
## 27 4.8471700 1   bg

time$n <- cumsum(time$E)-cumsum(1-time$E)
time$n <- c(0,time$n[1:length(time$n)-1])
time

##      brtimes E spec  n
## 1  0.0000000 1   aa  0
## 2  0.3334744 1   ab  1
## 3  0.3397169 1   ac  2
## 4  0.5256167 1   ad  3
## 5  0.8184213 1   ae  4
## 28 0.8808321 0   ab  5
## 6  0.9409274 1   ag  4
## 7  1.5194814 1   ah  5
## 29 1.7944359 0   ac  6
## 8  2.0255450 1   aj  5
## 9  2.2365991 1   ak  6
## 10 2.4714988 1   al  7
## 11 2.4885062 1   am  8
## 12 2.5404195 1   an  9
## 13 2.6038842 1   ao 10
## 14 2.7091034 1   ap 11
## 15 2.8988630 1   aq 12
## 16 2.9514716 1   ar 13
## 32 3.4082254 0   ak 14
## 17 3.4654121 1   at 13
## 18 3.5857570 1   au 14
## 19 3.6410367 1   av 15
## 20 3.7033597 1   aw 16
## 21 3.8538806 1   ax 17
## 33 3.9868879 0   ao 18
## 22 4.0000750 1   az 17
## 23 4.2353958 1   ba 18
## 24 4.4585866 1   bb 19
## 25 4.5703244 1   bc 20
## 31 4.6654420 0   ag 21
## 30 4.7666942 0   ae 20
## 26 4.8290354 1   bf 19
## 27 4.8471700 1   bg 20

st$br

##  [1] 0.3334744 0.3397169 0.5256167 0.8184213 0.8808321 0.9409274 1.5194814
##  [8] 1.7944359 2.0255450 2.2365991 2.4714988 2.4885062 2.5404195 2.6038842
## [15] 2.7091034 2.8988630 2.9514716 3.4082254 3.4654121 3.5857570 3.6410367
## [22] 3.7033597 3.8538806 3.9868879 4.0000750 4.2353958 4.4585866 4.5703244
## [29] 4.6654420 4.7666942 4.8290354 4.8471700 5.0000000

st$E

##  [1] 1 1 1 1 0 1 1 0 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 0 1 1 1 1 0 0 1 1

st$E == time$E[2:length(time$E)]

##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [15] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [29] TRUE TRUE TRUE TRUE

st$br[1:(length(st$br)-1)] == time$brtimes[2:length(time$brtimes)]

##  [1]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE
## [12]  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE
## [23]  TRUE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE FALSE

missing = time[which(is.element(time$spec,miss_spe)),]
missing

##      brtimes E spec  n
## 2  0.3334744 1   ab  1
## 3  0.3397169 1   ac  2
## 5  0.8184213 1   ae  4
## 28 0.8808321 0   ab  5
## 6  0.9409274 1   ag  4
## 29 1.7944359 0   ac  6
## 9  2.2365991 1   ak  6
## 13 2.6038842 1   ao 10
## 32 3.4082254 0   ak 14
## 33 3.9868879 0   ao 18
## 31 4.6654420 0   ag 21
## 30 4.7666942 0   ae 20

Phylo2p and update_tree

First, to warm up, we plot a simple tree

obs = '((AB:1,AA:1):1,AC:2):1;'
#library(ape)
obs = read.tree(text=obs)
is.rooted(obs)

## [1] TRUE

plot(obs,show.tip.label = T, root.edge = TRUE)

branching.times(obs)

## 4 5 
## 2 1

Then a longer one

# complete/incomplete trees
comp = '((AB:1,(AA:0.6,AD:0.3):0.4):1,((AE:0.6,AF:0.8):0.2,AC:1.5):0.5):1;'
comp= read.tree(text = comp)
is.rooted(comp)

## [1] TRUE

plot(comp,show.tip.label = T, root.edge = TRUE)

branching.times(comp)

##   7   8   9  10  11 
## 2.0 1.0 0.6 1.5 1.3

Then we drop extinxt species

incomp = drop.fossil(comp)
write.tree(incomp)

## [1] "((AB:1,AA:1):1,AC:2):1;"

is.rooted(incomp)

## [1] TRUE

plot(incomp,show.tip.label = T, root.edge = TRUE)

branching.times(incomp)

## 4 5 
## 2 1

Now we would like to be able to add an extinct species

in_dmea = phylo2p(incomp)
in_dmea$t

## NULL

in_dmea$E

## [1] 1 1

in_dmea$n

## [1] 1 2 3

update_tree(wt=in_dmea$t,t_spe=0.5,t_ext = 1.5 ,E=in_dmea$E,n=in_dmea$n)

## $wt
## [1] 0.5 1.0  NA  NA  NA 0.5
## 
## $E
## [1] 1 0 1 1
## 
## $n
## [1] 1 2 1 2 3

Now , after a little modification on the drop.fossil function, it works properly

in_dmea = phylo2p(obs)
in_dmea$t

## NULL

in_dmea$E

## [1] 1 1

in_dmea$n

## [1] 1 2 3

update_tree(wt=in_dmea$t,t_spe=0.5,t_ext = 1.5 ,E=in_dmea$E,n=in_dmea$n)

## $wt
## [1] 0.5 1.0  NA  NA  NA 0.5
## 
## $E
## [1] 1 0 1 1
## 
## $n
## [1] 1 2 1 2 3

Create_L

st <- sim_phyl()
L <- create_L(t=st$t, E=st$E)
miss <- L[L[,3] != (-1),]
miss_spe <- miss$spec
time <- rbind(data.frame(brtimes = L[,2], E=1, spec = L[,1]),data.frame(brtimes = miss[,3], E=0, spec = miss[,1]))
time <- time[order(time$brtimes),]
time$n <- cumsum(time$E)-cumsum(1-time$E)
all.equal(st$n,time$n)

## [1] TRUE

all.equal(st$n,cumsum(time$E)-cumsum(1-time$E))

## [1] TRUE

all.equal(c(0,st$br),c(time$brtimes,15))

## [1] TRUE

p2phylo

s = sim_phyl(ct=6)
plot(s$newick)

p = phylo2p(s$newick)
p$E == s$E

## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

p$wt == s$t

## [1] FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE  TRUE

#p$wt
#s$t
p$n == s$n

## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

s2 = p2phylo(p)
plot(s2)

voilà!, we found two things:

This algorithm does not label species in a general way yet, but only to models where all species has same probability to speciate.
We found a nice way to ilustrate ‘equivalent’ trees.

s = sim_phyl(ct=2)
plot(s$newick)

s2=p2phylo(phylo2p(s$newick))
plot(s2)

branching.times(s2)

##          8          9         10         11         12         13 
## 1.83283434 0.19032300 0.05646665 1.44325909 0.80553018 0.65776039

s2=p2phylo(phylo2p(s$newick))
plot(s2)

branching.times(s2)

##          8          9         10         11         12         13 
## 1.83283434 1.44325909 0.80553018 0.65776039 0.19032300 0.05646665

s2=p2phylo(phylo2p(s$newick))
plot(s2)

branching.times(s2)

##          8          9         10         11         12         13 
## 1.83283434 0.80553018 0.65776039 0.19032300 1.44325909 0.05646665

s2=p2phylo(phylo2p(s$newick))
plot(s2)

branching.times(s2)

##          8          9         10         11         12         13 
## 1.83283434 1.44325909 0.80553018 0.19032300 0.05646665 0.65776039

s2=p2phylo(phylo2p(s$newick))
plot(s2)

branching.times(s2)

##          8          9         10         11         12         13 
## 1.83283434 1.44325909 0.80553018 0.19032300 0.05646665 0.65776039

s2=p2phylo(phylo2p(s$newick))
plot(s2)

branching.times(s2)

##          8          9         10         11         12         13 
## 1.83283434 1.44325909 0.05646665 0.80553018 0.65776039 0.19032300

s2=p2phylo(phylo2p(s$newick))
plot(s2)

branching.times(s2)

##          8          9         10         11         12         13 
## 1.83283434 1.44325909 0.19032300 0.05646665 0.80553018 0.65776039

s2=p2phylo(phylo2p(s$newick))
plot(s2)

branching.times(s2)

##          8          9         10         11         12         13 
## 1.83283434 1.44325909 0.19032300 0.80553018 0.65776039 0.05646665

plot(s$newick)

branching.times(s$newick)

##          8          9         10         11         12         13 
## 1.83283434 1.44325909 0.80553018 0.19032300 0.65776039 0.05646665