Nothing
## This file is part of mvSLOUCH
## This software comes AS IS in the hope that it will be useful WITHOUT ANY WARRANTY,
## NOT even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
## Please understand that there may still be bugs and errors. Use it at your own risk.
## We take no responsibility for any errors or omissions in this package or for any misfortune
## that may befall you or others as a result of its use. Please send comments and report
## bugs to Krzysztof Bartoszek at krzbar@protonmail.ch .
## Function from the slouch package. Author : Jason Pienaar
'fitch.mvsl'<-function (phyltree, niche, deltran = FALSE, acctran = FALSE,
root = NULL)
{
## ----- Added Krzysztof Bartoszek ----------------------------------------
return_tree<-TRUE
if (inherits(phyltree,"phylo")){
if (!is.element("path_from_root",names(phyltree))){phyltree<-phyltree_paths(phyltree);return_tree<-TRUE}
tree.data<-.ape2slouch.mvsl(phyltree)
niche<-c(rep(NA,phyltree$Nnode),as.character(niche))
}else{
stop(sQuote("phyltree"), " must be of class ", sQuote("phylo"))
}
## ------------------------------------------------------------------------
niche <- as.factor(niche)
anc.char <- as.character(tree.data$ancestor)
anc.num <- as.numeric(anc.char)
anc.num[1] <- 0
node.char <- as.character(tree.data$nodes)
node.num <- as.numeric(node.char)
niche.num <- as.numeric(niche)
niches <- levels(niche)
num.code <- 1:length(levels(niche))
translator <- cbind(niches, num.code)
obs.char <- unique(niche.num[!is.na(niche.num)])
N <- length(node.num)
tree.matrix <- matrix(data = c(anc.num, node.num, niche.num),
ncol = 3, )
colnames(tree.matrix) <- c("Ancestors", "Nodes", "States")
Nint <- max(tree.matrix[, 1])
Nnodes <- max(tree.matrix[, 2])
downpass.node.states <- list()
cost <- 0
for (i in 1:Nnodes) {
downpass.node.states[[i]] <- tree.matrix[i, 3]
}
check.node.order <- tree.matrix[, 2] - tree.matrix[, 1]
if (any(check.node.order < 0))
traversal <- 2:Nint
else traversal <- Nint:1
for (i in traversal) {
children <- which(tree.matrix[i, 2] == tree.matrix[,
1])
if (length(children) == 2) {
child1.state <- downpass.node.states[[children[1]]]
child2.state <- downpass.node.states[[children[2]]]
downpass.node.states[[i]] <- intersect(child1.state,
child2.state)
}
if (length(children) > 2) {
child1.state <- downpass.node.states[[children[1]]]
child2.state <- downpass.node.states[[children[2]]]
child3.state <- downpass.node.states[[children[3]]]
tmp <- intersect(child1.state, child2.state)
downpass.node.states[[i]] <- intersect(tmp, child3.state)
}
if (length(downpass.node.states[[i]]) == 0) {
if (length(children) == 2) {
downpass.node.states[[i]] <- union(child1.state,
child2.state)
cost = cost + 1
}
if (length(children) > 2) {
tmp <- union(child1.state, child2.state)
downpass.node.states[[i]] <- union(tmp, child3.state)
cost = cost + 1
}
}
}
if (any(check.node.order < 0)) {
children <- which(tree.matrix[1, 2] == tree.matrix[,
1])
if (length(children) == 2) {
child1.state <- downpass.node.states[[children[1]]]
child2.state <- downpass.node.states[[children[2]]]
downpass.node.states[[1]] <- intersect(child1.state,
child2.state)
}
if (length(children) > 2) {
child1.state <- downpass.node.states[[children[1]]]
child2.state <- downpass.node.states[[children[2]]]
child3.state <- downpass.node.states[[children[3]]]
downpass.node.states[[1]] <- intersect(c(child1.state,
child2.state), child3.state)
}
if (length(downpass.node.states[[1]]) == 0) {
if (length(children) == 2) {
downpass.node.states[[1]] <- union(child1.state,
child2.state)
cost = cost + 1
}
if (length(children) > 2) {
tmp <- union(child1.state, child2.state) ## changed from tmp <- union(child1.state, child.state) by Krzysztof Bartoszek
downpass.node.states[[1]] <- union(tmp, child3.state)
cost = cost + 1
}
}
}
if (any(check.node.order < 0))
traversal <- 2:Nint
else traversal <- Nint:2
pre.traversal <- rev(traversal)
finalpass.node.states <- downpass.node.states
if (!is.null(root))
finalpass.node.states[[1]] <- root
if (length(finalpass.node.states[[1]]) >= 2) {
message("There is an ambiguity at the root node, as given below")
print(translator[as.numeric(finalpass.node.states[[1]])])
message("Set this to one of the ambiguous states using root = state (state needs to be in quotation marks) in the function call before attempting the deltran or acctran reconstructions")
}
for (i in pre.traversal) {
parent <- tree.matrix[i, 2]
ancestor <- tree.matrix[i, 1]
children <- which(tree.matrix[, 1] == tree.matrix[i,
2])
finalpass.node.states[[parent]] <- intersect(downpass.node.states[[parent]],
finalpass.node.states[[ancestor]])
if (setequal(finalpass.node.states[[parent]], finalpass.node.states[[ancestor]]) ==
FALSE) {
if (length(children) == 2) {
child1.state <- downpass.node.states[[children[1]]]
child2.state <- downpass.node.states[[children[2]]]
if (length(intersect(child1.state, child2.state)) !=
0) {
finalpass.node.states[[parent]] <- union(downpass.node.states[[parent]],
(intersect(finalpass.node.states[[ancestor]],
union(child1.state, child2.state))))
}
if (length(intersect(child1.state, child2.state)) ==
0) {
finalpass.node.states[[i]] <- union(downpass.node.states[[parent]],
finalpass.node.states[[ancestor]])
}
if (deltran == TRUE) {
if (length(finalpass.node.states[[i]]) >= 2)
finalpass.node.states[[i]] = finalpass.node.states[[ancestor]]
}
if (acctran == TRUE) {
if (length(finalpass.node.states[[i]]) >= 2)
finalpass.node.states[[i]] = setdiff(finalpass.node.states[[parent]],
finalpass.node.states[[ancestor]])
}
}
if (length(children) > 2) {
child1.state <- downpass.node.states[[children[1]]]
child2.state <- downpass.node.states[[children[2]]]
child3.state <- downpass.node.states[[children[3]]]
if (length(intersect(child3.state, intersect(child1.state,
child2.state))) != 0) {
finalpass.node.states[[parent]] <- union(downpass.node.states[[parent]],
(intersect(finalpass.node.states[[ancestor]],
union(union(child3.state, child1.state),
child2.state))))
}
if (length(intersect(child3.state, intersect(child1.state,
child2.state))) == 0) {
finalpass.node.states[[i]] <- union(downpass.node.states[[parent]],
finalpass.node.states[[ancestor]])
}
}
}
}
dat <- list(treematrix = tree.matrix, Final.states = finalpass.node.states,
downpass.cost = cost, niche.code = translator)
N <- length(dat$Final.states)
count <- 0
ambig <- NA
for (i in 1:N) {
if (length(dat$Final.states[[i]]) >= 2) {
count = count + 1
ambig <- c(ambig, i)
}
}
if (count != 0) {
message("The number of ambiguous nodes are:")
print(count)
message(" and they are at nodes:")
print(ambig[-1])
message("Try deltran OR acctran = TRUE in the function call in order to implement a delayed or accelerated character transformation")
}
f.states <- .make.states.mvsl(dat)
## ----- Added Krzysztof Bartoszek ----------------------------------------
## reorder nodes to be consistent with phylo object edge numbering
f_states_phylo<-rep(NA,nrow(phyltree$edge))
for (i in 1:nrow(phyltree$edge)){## for each edge
i_node_ending<-phyltree$edge[i,2]
c_node_ending<-as.character(i_node_ending)
if (is.element(i_node_ending,phyltree$tip_species_index)){
i_tip_index<-which(phyltree$tip_species_index==i_node_ending)
c_node_ending<-phyltree$tip.label[i_tip_index]
}
f_states_phylo[i]<-as.character(f.states[which(tree.data$species==c_node_ending)])
}
root_regime<-f.states[which(tree.data$species==as.character(phyltree$root_index))]
if (return_tree){
lres<-list(branch_regimes=f_states_phylo,root_regime=root_regime,phyltree=phyltree)##,tree.data=tree.data,f.states=f.states)
}else{
lres<-list(branch_regimes=f_states_phylo,root_regime=root_regime)##,tree.data=tree.data,f.states=f.states)
}
## ------------------------------------------------------------------------
return(lres)
}
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.