Nothing
# Squared-change parsimony hidden state prediction for a continuous trait, for tips on a tree.
hsp_squared_change_parsimony = function(tree,
tip_states, # numeric vector of size Ntips. Can include NA.
weighted = TRUE,
check_input = TRUE){
Ntips = length(tree$tip.label);
Nnodes = tree$Nnode;
# basic error checking
if(length(tip_states)!=Ntips) stop(sprintf("ERROR: Length of tip_states (%d) is not the same as the number of tips in the tree (%d)",length(tip_states),Ntips));
if(!is.numeric(tip_states)) stop(sprintf("ERROR: tip_states must be a numeric vector"))
if(check_input){
if((!is.null(names(tip_states))) && any(names(tip_states)!=tree$tip.label)) stop("ERROR: Names in tip_states and tip labels in tree don't match (must be in the same order).")
}
# find known_tips, extract known_subtree and synchronize known_tip_states with known_subtree
known_tips = which(!is.na(tip_states));
if(length(known_tips)==0) stop("ERROR: All tip states are hidden");
extraction = get_subtree_with_tips(tree, only_tips=known_tips, omit_tips=FALSE, collapse_monofurcations=TRUE, force_keep_root=TRUE);
known_subtree = extraction$subtree;
known2all_tips = extraction$new2old_tip
known2all_nodes = extraction$new2old_node
known_tip_states = tip_states[known2all_tips]
# perform ancestral state reconstruction on known_subtree
asr_results = asr_squared_change_parsimony( known_subtree,
tip_states = known_tip_states,
weighted = weighted,
check_input = FALSE);
if(!asr_results$success) return(list(success=FALSE))
# forward-project reconstructions to tips with hidden state
states = rep(0, times=Ntips+Nnodes);
states[known2all_tips] = known_tip_states;
states[known2all_nodes+Ntips] = asr_results$ancestral_states;
states_known = rep(FALSE, times=(Ntips+Nnodes))
states_known[known2all_tips] = TRUE;
states_known[known2all_nodes + Ntips] = TRUE;
states = apply_BM_parsimony_to_missing_clades_CPP( Ntips = Ntips,
Nnodes = Nnodes,
Nedges = nrow(tree$edge),
tree_edge = as.vector(t(tree$edge))-1, # flatten in row-major format and make indices 0-based
states_known = states_known,
states = states);
return(list(success = TRUE,
total_sum_of_squared_changes = asr_results$total_sum_of_squared_changes,
states = states))
}
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.