Nothing
mytree.symmetric.taxa <-
function (m, waitsp, waitext, complete=TRUE, tiplabel, shiftsp, shiftext, sampling, gsa)
{
# create waiting time functions #
if (is.function(waitsp)){
rnumbsp <- waitsp#parse(text="waitsp()")
} else if (is.character(waitsp)){
rnumbsp <- express.distribution(waitsp)
}
if (is.function(waitext)){
rnumbext <- waitext#parse(text="waitext()")
firstextpar <- "funk"
} else if (is.character(waitext)){
rnumbext <- express.distribution(waitext)
firstextpar <- get.first.par.distribution(waitext)
}
if (is.function(shiftsp$strength)){
rnumbshiftsp <- shiftsp$strength#parse(text="shiftsp$strength()")
} else if (is.character(shiftsp$strength)){
rnumbshiftsp <- express.distribution(shiftsp$strength)
}
if (is.function(shiftext$strength)){
rnumbshiftext <- shiftext$strength#parse(text="shiftext$strength()")
} else if (is.character(shiftext$strength)){
rnumbshiftext <- express.distribution(shiftext$strength)
}
# end create waiting time functions #
labellivingsp=tiplabel[1]
labelextinctsp=tiplabel[2]
shiftsplabel=tiplabel[3]
shiftextlabel=tiplabel[4]
shiftspprob=shiftsp$prob
shiftextprob=shiftext$prob
#testing for shift speciation
testshiftsp <- function (spt){
if (shiftspprob!=0){
if (runif(1,0,1)<shiftspprob) { #cheking if shift happens based on a uniform distribution U[0,1]
shiftspstrength <- rnumbshiftsp()#eval(rnumbshiftsp) #draw shift strength
#print(paste("Hey, a speciation shift happened at node", nextsp - i, "with strenght", shiftspstrength ))
}else{
shiftspstrength <- shiftspm[shiftspm[,"node"]==species,"strength"] #inherit ancestor shift
#print(paste("No speciation shift!",nextsp - i, "is inheriting the shift of",shiftspstrength, "from", species ))
}
}else{
shiftspstrength <- 1
}
spt <- spt*shiftspstrength #update speciation time
shiftspmr <- c((nextsp - i),shiftspstrength ) #link shift with node to be added to shiftspm...
attributes(spt) <- NULL
attributes(shiftspmr) <- NULL
return(list(spt=spt,shiftspmr=shiftspmr))#return updated shift and updated shiftspmatrix
}
#testing for shifts on extinction process
testshiftext <- function (extt){
if (shiftextprob!=0){
if (runif(1,0,1)<shiftextprob) { #cheking if shift happens based on a uniform distribution U[0,1]
shiftextstrength <- rnumbshiftext()#eval(rnumbshiftext) #draw shift strength
#extt <- extt*shiftextstrength #update extinction time
#print(paste("Hey, an extinction shift happened at node", nextsp - i, "with strenght", shiftextstrength ))
}else{
shiftextstrength <- shiftextm[shiftextm[,"node"]==species,"strength"] #inherit ancestor shift
#print(paste("No extinction shift!",nextsp - i, "is inheriting the shift of",shiftextstrength, "from", species ))
#extt <- extt*shiftextstrength #update extt by ancestor shift
}
}else{
shiftextstrength <- 1
}
extt <- extt*shiftextstrength #update extinction time
shiftextmr <- c((nextsp - i),shiftextstrength ) #link shift with node to be added to shiftextm...
attributes(extt) <- NULL
attributes(shiftextmr) <- NULL
return(list(extt=extt,shiftextmr=shiftextmr))#return updated shift and updated shiftextmatrix
}
##SM-E
# tracing back the time for one species until origin
trajectory <- function (trace){
#trace should indicate the edge number that is to be followed until the origin
trajectory <- NULL
while ( length( which(edge[,2] == trace)) ){
atual <- which(edge[,2] == trace)
trajectory <- c(edge.length[atual], trajectory)
trace <- edge[atual,1]
}
return(trajectory)
}
bingo <- FALSE
while (bingo == FALSE) # ######## while bingo
{
stop <- FALSE
stopsearch <- FALSE
#SM-Modified
mytree <-list(edge=NULL, tip.label=NULL, edge.length=NULL, Nnode=NULL, root.edge=NULL, age=NULL,
shiftsp=NULL,shiftext=NULL, shifted.sp.living=NULL,shifted.sp.extinct=NULL,
shifted.ext.living=NULL, shifted.ext.extinct=NULL)
class(mytree) <- "phylo"
edge <- matrix(c(-1,-2), ncol=2)
##SM-S creating shiftsp and shiftext matrixes...
shiftspm <- matrix(c(-1,1,-2,1), byrow=TRUE, ncol=2, dimnames=list(NULL,c("node","strength")))
shiftextm <- matrix(c(-1,1,-2,1), byrow=TRUE, ncol=2, dimnames=list(NULL,c("node","strength")))
##SM-E
leaves <- NULL
realleaves <- NULL
extinct <- NULL
realextinct <- NULL
age <- NULL
timeline <- NULL
extinct <- NULL
tip.label <- NULL
pnodges <- NULL
timeline <- NULL
##SM - S
shiftedspliving <- NULL
shiftedspextinct <- NULL
shiftedextliving <- NULL
shiftedextextinct <- NULL
##SM - E
# initial if in case the (-1,-2) edge get extinct or bigger than age
spt <- rnumbsp()#eval(rnumbsp)
#to remove NaN warnings messages
{
if (firstextpar == 0)
{
extt <- suppressWarnings(rnumbext())
}
else
{
extt <- rnumbext()
}
}
# if to see if the user simulates with extinction = ZERO and avoid error generated by expression-distribution when rate equals zero
{
if (is.nan(extt))
{
extt <- spt +1 # we add so that the sp will always occurs first, and an extinction will never occur
}
}
{
if (spt <= extt)
{
status <- "sp" #occurred an speciation
edge.length <- spt
leaves <- -2
}
else
{
status <- "ext" #occurred an extinction
edge.length <- extt
extinct <- -2
stop <- TRUE
}
}
pnodges <- c(leaves, extinct)
timeline <- edge.length[1]
# for the while (increase of the tree)
while (stop == FALSE)
{
stopsearch <- FALSE
#SM-modified #ATTENTION CHANGE ALL OVER!!!!
positionnextsp <- which(timeline == min(timeline[pnodges%in%leaves]))[1] #ATTENTION, we take the first element since it can happen that we get two equal numbers
nextsp <- min(edge[,2])
species <- pnodges[positionnextsp]
ptime <- timeline[positionnextsp]
pnodges <- pnodges[-positionnextsp]
timeline <- timeline[-positionnextsp]
i <- 1
#print(paste("Now we will start from species ", species))
for (i in 1:2)
{
edge <- rbind( edge, c(species, (nextsp - i)))
spt <- rnumbsp()
##SM -S now we examine for shifts and update spt and shiftspm
testshiftspout <- testshiftsp(spt)
spt <- testshiftspout$spt
shiftspm <- rbind(shiftspm, testshiftspout$shiftspmr)
##SM -E
#to remove NaN warnings messages
{
if (firstextpar == 0)
{
extt <- suppressWarnings(rnumbext())
}
else
{
extt <- rnumbext()
}
}
# if to see if the user simulates with extinction = ZERO and avoid error generated by exp. distribution when rate equals zero
{
if (is.nan(extt))
{
extt <- spt +1000 # we add 1000 so that the sp will always occurs first, and we will never have an extinction
}
}
##SM -S now we examine for shifts and update extt and shiftextm
testshiftextout <- testshiftext(extt)
extt <- testshiftextout$extt
shiftextm <- rbind(shiftextm, testshiftextout$shiftextmr)
##SM -E
{
if (spt <= extt)
{
status <- "sp" #occurred an speciation
edge.length <- c(edge.length, spt)
leaves <- c(leaves, (nextsp - i))
ntime <- spt
}
else
{
status <- "ext" #occurred an extinction
edge.length <- c(edge.length, extt)
extinct <- c(extinct, (nextsp - i))
ntime <- extt
}
}
pnodges <- c(pnodges, (nextsp - i))
timeline <- c(timeline, (ptime + ntime))
}
#SM modified
leaves <- leaves[-which(leaves==species)]
#now we test for the smallest pnodge and see if there is the enough leaving creatures
while (stopsearch == FALSE)
{
{
if (length(pnodges) == 0)
{
stop <- TRUE
stopsearch <- TRUE
}
else
{
{
if (length(pnodges) >= m)
{
#we achieved m species at least
bingo<- TRUE
stop <- TRUE
stopsearch <- TRUE
age <- min(timeline)
realleaves <- pnodges
}
else
{
positionsp <- which(timeline == min(timeline))[1]
#testing fit the smallest tested one is an extinct
{
if (pnodges[positionsp] %in% extinct)
{
#if it is on the extinct, we will remove it from the pnodges and timeline and add it to the realextinct
realextinct <- c(realextinct, pnodges[positionsp])
pnodges <- pnodges[-positionsp]
timeline <- timeline[-positionsp]
}
else
{
stopsearch <- TRUE
}
}
}
}
}
}
}
}
}# ######## while bingo finish
#after this point, all threes that are leaving this last while have at least m leavingspecies....
extinct <- realextinct
#cutting tree at desired age
step <- 1
for (step in 1: length(pnodges))
{
traject <- trajectory(pnodges[step])
traject <- traject[-length(traject)]
edge.length[which(edge[,2]==pnodges[step])] <- (age-sum(traject))
}
# final if... in case of (stop == TRUE) , we write the tree "mytree"
{
if (stop == TRUE)
{
#### replacing to the ape format
prealleaves <- realleaves
{
if (length(realleaves) > 0)
{
realleaves <- c(1:length(realleaves))
i <- 1
for (i in 1:length(realleaves))
{
edge[ which(edge[,2] == prealleaves[i]), 2 ] <- realleaves[i]
## SM- S changing nodes names
shiftspm[ which(shiftspm[,"node"] == prealleaves[i]), "node" ] <- realleaves[i]
shiftextm[ which(shiftextm[,"node"] == prealleaves[i]), "node" ] <- realleaves[i]
## SM- E
}
tip.label <- paste(labellivingsp, realleaves, sep = "")
## SM - S creating vector with shifted species and updating shifted names
#for shifts on speciation
shiftedspliving <- rep(0,length(realleaves)) #preparing vector
data <- shiftspm[shiftspm[,"node"]%in%realleaves,]
shiftedspliving[data[sort.list(data[,"node"]),"strength" ]!=1] <- 1
tip.label[shiftedspliving==1] <- paste(tip.label[shiftedspliving==1], shiftsplabel, sep=" ")
shiftedspliving <- cbind(realleaves,shiftedspliving)
colnames(shiftedspliving) <- c("LivingSpecies", "shift")
#for shifts on extinction
shiftedextliving <- rep(0,length(realleaves))
data <- shiftextm[shiftextm[,"node"]%in%realleaves,]
shiftedextliving[data[sort.list(data[,"node"]),"strength" ]!=1] <- 1
tip.label[shiftedextliving==1] <- paste(tip.label[shiftedextliving==1], shiftextlabel, sep=" ")
shiftedextliving <- cbind(realleaves,shiftedextliving)
colnames(shiftedextliving) <- c("LivingSpecies", "shift")
## SM - E
}
}
pextinct <- extinct
{
if (length(extinct) > 0)
{
extinct <- c((length(realleaves)+1):(length(realleaves)+length(extinct)))
i <- 1
for (i in 1:length(extinct))
{
edge[ which(edge[,2] == pextinct[i]), 2 ] <- extinct[i]
## SM- S changing nodes names
shiftspm[ which(shiftspm[,"node"] == pextinct[i]), "node" ] <- extinct[i]
shiftextm[ which(shiftextm[,"node"] == pextinct[i]), "node" ] <- extinct[i]
## SM- E
}
tip.label.tail <- paste(labelextinctsp, extinct, sep = "")
## SM - S creating vector with shifted species and updating shifted names
#for shifts on speciation
shiftedspextinct <- rep(0,length(extinct))
data <- shiftspm[shiftspm[,"node"]%in%extinct,]
#error handling in case only one species got extinct and data is transformed into a vector and not a matrix...
if(!is(data, "matrix")){
data <- matrix(data, ncol=2)
}
shiftedspextinct[data[sort.list(data[,1]),2 ]!=1] <- 1
tip.label.tail[shiftedspextinct==1] <- paste(tip.label.tail[shiftedspextinct==1], shiftsplabel, sep=" ")
shiftedspextinct <- cbind(extinct,shiftedspextinct)
colnames(shiftedspextinct) <- c("ExtinctSpecies", "shift")
#for shifts on extinction
shiftedextextinct <- rep(0,length(extinct))
data <- shiftextm[shiftextm[,"node"]%in%extinct,]
if(!is(data, "matrix")){
data <- matrix(data, ncol=2)
}
shiftedextextinct[data[sort.list(data[,1]),2 ]!=1] <- 1
tip.label.tail[shiftedextextinct==1] <- paste(tip.label.tail[shiftedextextinct==1], shiftextlabel, sep=" ")
shiftedextextinct <- cbind(extinct,shiftedextextinct)
colnames(shiftedextextinct) <- c("ExtinctSpecies", "shift")
## SM - E
tip.label <- c(tip.label, tip.label.tail )
}
}
#regarding the edges that lead to an extinct or leaving final species, but are not the final edges
potheredges <- levels(as.factor(edge[edge <0]))
otheredges <- rev(seq((max(realleaves, extinct)+1), length.out=length(potheredges)))
#substituting...
i <- 1
for (i in 1:length(potheredges))
{
edge[ edge == potheredges[i] ] <- otheredges[i]
## SM- S changing nodes names
shiftspm[shiftspm[,"node"]==potheredges[i],"node"] <- otheredges[i]
shiftextm[shiftextm[,"node"]==potheredges[i],"node"] <- otheredges[i]
## SM- E
}
mytree$edge <- edge
mytree$tip.label <- tip.label
mytree$edge.length <- edge.length
mytree$Nnode <- length(realleaves) + length(extinct)
mytree$root.edge <- edge.length[1]
root.edge <- edge.length[1]
mytree$age <- age
## SM - S
mytree$shiftsp <- shiftspm
mytree$shiftext <- shiftextm
mytree$shifted.sp.living <- shiftedspliving
mytree$shifted.sp.extinct <- shiftedspextinct
mytree$shifted.ext.living <- shiftedextliving
mytree$shifted.ext.extinct <- shiftedextextinct
## SM - E
}
}
#final handling before plotting with ape
{
if (length(realleaves) == 0)
{
# in case no specie is surviving until final simulation time
mytree <- 0
}
else
{
{
if ( length(realleaves) == 1 & complete == FALSE)
{
#in case only one specie is surviving, even if other speciations events occured in the history
mytree <- 1
}
else
{
{
if (length(realleaves)==1 & length(extinct)==0 & complete==TRUE)
{
#in case only one species is surviving and was the only one that existed
mytree <- 1
}
else
{
#in case non of the above condition is fulfilled, there will be a tree with no initial branch, tree starts at the MRCA
#this is done to be aple to plot and be compatible wth `ape` package
mytree <- collapse.singles(mytree)
#cheking status of 'complete' and take or dont take extincted species out of final tree
{
if (complete == FALSE)
{
mytree<- drop.fossil(mytree)
mytree$root.edge <- root.edge
}
}
}
}
}
}
}
}
if (sampling$frac!=1 & class(mytree)=="phylo" & !gsa){ # do sampling.....
mytree <- sample.mytree(mytree.=mytree, realleaves.=realleaves, extinct.=extinct, sampling.=sampling, complete.=complete)
}
return(mytree)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.