R/BioGeoBEARS_plots_v1.R

Defines functions plot_BioGeoBEARS_results add_statum_boundaries_to_phylo_plot

#######################################################
# Functions to make nice plots
#######################################################


#######################################################
# add_statum_boundaries_to_phylo_plot
# Add stratum boundaries to a plot of a time-scaled phylogeny
# Correcting for the offset in the plot
#######################################################
#' Add stratum boundaries to a phylogeny plot
#' 
#' Adds vertical lines, representing stratum boundaries, 
#' to a plot of a phylogeny. 
#' 
#' \bold{Note:} This function asssumes the phylogeny is plotted with
#' tips to the right. 
#' 
#' \bold{Background:} This function may be useful, because plotting vertical lines
#' onto plots generated with \code{APE}'s \code{\link[ape]{plot.phylo}} plots at the 
#' time-points desired will not work. This is because, confusingly, APE's plotting
#' of phylogenies puts the x-axis minimum at 0, and the max at
#' the height of the tree above the root, plus a buffer for 
#' the tip labels.*
#' 
#' Therefore, \code{add_statum_boundaries_to_phylo_plot} calculates the height of 
#' the highest tip and subtracts the \code{timeperiods}, to get the plotted x-values, then 
#' uses abline to plot the lines. (This information may be helpful if you want to 
#' plot your arbitrary lines on top of phylogenies plotted in other orientations.
#' 
#' *This may change further if the \code{phylo} object
#' \code{tr} has a root edge (a branch below the root), so the function might not 
#' draw the lines in the correct place if there is a root edge.  You can tell if your 
#' tree has a root edge by running \code{names(tr)} and seeing if a \code{root.edge} 
#' is listed. Alternatively, look at \code{tr$root.edge}. 
#'
#' @param tr An \code{APE} \code{phylo} (tree) object.
#' @param timeperiods A vector of one or more timeperiods. The default is 1 timeperiod
#'                    at 1 million years (or 1 time unit).
#' @param lty Line type, e.g. "dashed".  See \code{\link{par}}
#' @param col Color of the line, e.g. "grey50" (default)
#' @param plotlines If TRUE (default), the lines are plotted on the current plot. 
#'                  If FALSE, the function just returns the x-axis positions
#'                  of the lines on the phylogeny plot. 

#' @return \code{line_positions_on_plot} The x-axis positions of the lines.
#' @export
#' @seealso \code{\link[base]{plot}}, \code{\link[base]{par}}, \code{\link[ape]{plot.phylo}}
#' @note Go (BioGeo)BEARS!
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu} 
#' @references
#' 
#' @examples
#' 
#' # Loading the default tree
#' trfn = np(paste(addslash(extdata_dir), "Psychotria_5.2.newick", sep=""))
#' tr = read.tree(trfn)
#' 
#' # Get the tree coordinates (APE 5.0 or higher), i.e. the x and y of each 
#' node.
#' trcoords = plot_phylo3_nodecoords_APE5(tr, plot=FALSE, root.edge=TRUE)
#' 
#' # Set reasonable x-limits (unlike the defaults on APE5.0)
#' xlims = c(min(trcoords$xx), 1.42*max(trcoords$xx))
#' 
#' # Plot the tree
#' trplot = plot(tr, cex=1, x.lim=xlims); axisPhylo()
#' 
#' # Add the stratum boundaries
#' add_statum_boundaries_to_phylo_plot(tr, timeperiods=1, lty="dashed", col="gray50", plotlines=TRUE)
#' 
add_statum_boundaries_to_phylo_plot <- function(tr, timeperiods=1, lty="dashed", col="gray50", plotlines=TRUE)
	{
	# Example code for the programmer to run easily
	SETUP='
	# Loading the default tree
	trfn = np(paste(addslash(extdata_dir), "Psychotria_5.2.newick", sep=""))
	tr = read.tree(trfn)
	
	# Get the tree coordinates (APE 5.0 or higher), i.e. the x and y of each node.
	trcoords = plot_phylo3_nodecoords_APE5(tr, plot=FALSE, root.edge=TRUE)
	
	# Set reasonable x-limits (unlike the defaults on APE5.0)
	xlims = c(min(trcoords$xx), 1.42*max(trcoords$xx))
	
	# Plot the tree
	trplot = plot(tr, cex=1, x.lim=xlims); axisPhylo()
	
	# Add the stratum boundaries
	add_statum_boundaries_to_phylo_plot(tr, timeperiods=1, lty="dashed", col="gray50", plotlines=TRUE)
	' ## END SETUP
	
	ntips = length(tr$tip.label)
	tr_table = prt(tr, printflag=FALSE)
	tr_height = tr_table$time_bp[ntips+1]
	line_positions_on_plot = tr_height-timeperiods
	if (plotlines == TRUE)
		{
		abline(v=line_positions_on_plot, lty=lty, col=col)
		} # END if (plotlines == TRUE)
	return(line_positions_on_plot)
	} # END add_statum_boundaries_to_phylo_plot <- function(tr, lty="dashed")


#######################################################
# plot_BioGeoBEARS_results
#######################################################
#' Plot the results of a BioGeoBEARS run
#' 
#' This function plots on a tree the highest-probability ancestral states (ranges), 
#' splits if desired (these are the ranges/states just after cladogenesis, and are 
#' plotted on the corners of a tree), and/or pie charts at nodes/splits.  A legend 
#' tying the relationship between colors and states/ranges is also optionally plotted.
#'
#' The legend, if desired, is plotted on a separate plot, as it is very difficult 
#' to predict whether or not there will be appropriate space on any given tree plot. 
#' The utility of a classical legend is also debatable, as 
#' \code{plot_BioGeoBEARS_results} plots the colors and state/range names directly 
#' onto the plot.  Any legend will get unwieldy above perhaps 16
#' states/ranges, which is just 4 areas with no constraints (see \code{\link[cladoRcpp]{numstates_from_numareas}}, or type \code{numstates_from_numareas(numareas=4, maxareas=4, include_null_range=TRUE)}.
#'
#' In any case, the human eye can only easily read a few colors on a plot. The 
#' philosophy adopted in \code{BioGeoBEARS} plots is to use bright, primary colors for the 
#' single-area ranges, and then blend these colors for multi-areas ranges. A range 
#' all areas is colored white. Coupled with plotting the letter codes for the 
#' ranges ("A", "AB", "ABC"), this makes for reasonably readable plots. (Of course,
#' some researchers design their own custom plots in Adobe Illustrator or elsewhere.)
#' 
#' Note that \code{plot_BioGeoBEARS_results} assumes
#' that the ancestral states were calculated under the global optimum model (rather than the local optimum, with the model re-optimized for each 
#' possible state at each possible node, as done in e.g. \code{LAGRANGE}), and that these are marginal probabilities, i.e. this is not a joint reconstruction,
#' instead it gives the probabilities of states at each node.  This will not always be readable as a joint reconstruction (it could depict split scenarios
#' that are not possible, for instance.)
#' 
#' @param results_object The results object from \code{\link{bears_optim_run}} (with ancestral states on).
#' @param analysis_titletxt The main title of the plot. If NULL, \code{results_object$inputs$description} is checked.
#' @param addl_params The function will plot the log-likelihood (LnL) and the ML values of the free parameters. If you want additional parameters plotted, list them here.
#' @param plotwhat To plot the ML discrete states, "text".  To plot a piechart of the relative probability of all the states, "pie".
#' @param label.offset Offset for the tree tip labels. If \code{NULL}, program chooses 0.05 x tree height.
#' @param tipcex \code{cex} value for the tiplabels (scaling factor, i.e. 0.5 is half size)
#' @param statecex \code{cex} value for the states (scaling factor, i.e. 0.5 is half size). Used on piecharts if plotwhat="pie".
#' @param splitcex \code{cex} value for the splits (scaling factor, i.e. 0.5 is half size). Used on piecharts if plotwhat="pie".
#' @param titlecex \code{cex} value for the title (scaling factor, i.e. 0.5 is half size). 
#' @param plotsplits If \code{TRUE}, plot states on the corners -- text or pie charts, depending on \code{plotwhat}.
#' @param plotlegend If \code{TRUE}, make a (separate) plot with a legend giving the colors for each state/range, using \code{\link{colors_legend}}.
#' @param legend_ncol The number of columns in the legend.  If \code{NULL} (default), the function calculates \code{floor(sqrt(length(possible_ranges_list_txt) / 2))} 
#' when the number of states is <=64, and \code{sqrt(ceiling(length(possible_ranges_list_txt)))} when > 64. Note that when you have hundreds of states, there is probably 
#' no good way to have a readable legend, and it is easier to just rely upon printing the 
#' character codes for the ML states in the plots, with the colors, and users can then see and trace the common colors/states by eye.
#' @param legend_cex The cex (character expansion size) for the legend.  Defaults to 1, which means the \code{\link[graphics]{legend}} function determines the 
#' size.  The value 2.5 works well for 15 or 16 states/ranges.
#' @param cornercoords_loc The directory location containing the R script \code{plot_phylo3_nodecoords.R}. This function, modified from the APE function
#' \code{\link[ape]{plot.phylo}}, cannot be included directly in the R package as it contains C code that does not pass CRAN's R CMD check. The default, 
#' cornercoords_loc="manual", will not allow split states to be plot.  The R script \code{plot_phylo3_nodecoords.R} is located in the BioGeoBEARS extension data 
#' directory, \code{extdata/a_scripts}.  You should be able to get the full path with \code{list.files(system.file("extdata/a_scripts", package="BioGeoBEARS"), full.names=TRUE)}.
#' @param include_null_range If \code{TRUE} (default), the null range is included in calculation of colors. (Safest for now.)
#' @param tr Tree to plot on. Default \code{NULL}, which means the tree will be read from the file at \code{results_object$inputs$trfn}.
#' @param tipranges Tip geography data. Default \code{NULL}, which means the tree will be read from the file at \code{results_object$inputs$geogfn}.
#' @param if_ties What to do with ties in probability. Currently, 
#' the options are:
#' (1) "takefirst", which takes the first tied state in the 
#' probabilities column (The full probabilities of all states will be
#' shown in e.g. pie charts, of course); (2) "asterisk", which 
#' returns returns the first tied state, but marks it with an 
#' asterisk ("*").
#' @param pie_tip_statecex  \code{cex} value for pies at the tips (scaling factor, i.e. 0.5 is half size).
#' @param juststats If \code{FALSE} (default), plots are plotted. If \code{TRUE}, 
#' no plots are done, 
#' the function just returns the summary statistics.
#' @param root.edge Should the root edge be plotted, if it exists in the tree?  Passed to
#' plot.phylo().  This can be handy if the root state display is getting cut off.
#' @param colors_list_for_states Default \code{NULL} auto-generates colors with 
#' \code{get_colors_for_states_list_0based}. Otherwise, users can specify colors for 
#' each state (remember that e.g. 4 areas can mean 2^4=16 states).
#' @param skiptree Skip the plotting of the tree -- useful for plotting the state labels
#' e.g. on top of a stochastic map. Default \code{FALSE}.
#' @param show.tip.label Same as for APE's \code{plot.phylo}.
#' @param tipcol The tip text color. Default "black".
#' @param dej_params_row Parameters used to generate an SSE simulation. dej_params_row can be 
#' obtained from \code{get_dej_params_row_from_SSEsim_results_processed}.
#' @param plot_max_age The maximum age of the plot (age below the tips). Default
#' is tree height
#' @param skiplabels If \code{TRUE}, skip the nodelabels command, resulting in plotting just the 
#' tree. (Yes, you could have just used \code{FALSE}).  Default \code{FALSE}.
#' @param plot_stratum_lines If \code{TRUE} (default), plot dotted lines at the stratum 
#' boundaries, *if* it's a time-stratified results object.
#' @param simplify_piecharts If \code{TRUE}, just plot one slice for the maximum probability, 
#' and white for "other". This should help with large plots with many ranges, 
#' which can overwhelm some graphics programs (imagine 1000 slices per piechart,
#' times 1000+nodes). Default \code{FALSE}.
#' @param tipboxes_TF Plot the tip boxes (tip states)?  Default \code{TRUE}.
#' @param tiplabel_adj Justification for tiplabel boxes (same as "adj" parameter 
#' of text()). Default is c(0.5). Left justification: c(0). Right: c(1). 
#' Justification top left: c(0,0), etc.
#' @param no.margin Same as in plot.phylo. Default is FALSE (meaning yes, 
#' there will be margins).
#' @param xlims Same as in plot.phylo x.lim. Default is basically c(0,treeheight), so
#' if tiplabels are getting cut off, try e.g. xlims=c(0,1.5*treeheight).
#' @param ylims Same as in plot.phylo y.lim. Default is basically c(0,numtips), so
#' if the spacing between the title and time axis and the tree are too big,
#' try e.g. ylims=c(0+10,numtips-10). Trial and error will get you there. 
#' See: https://stat.ethz.ch/pipermail/r-sig-phylo/2013-March/002540.html
#' @return NULL
#' @export
#' @seealso \code{\link{get_leftright_nodes_matrix_from_results}}, \code{\link{corner_coords}}, \code{\link[ape]{plot.phylo}}, \code{\link[ape]{plot.phylo}}, \code{\link[ape]{tiplabels}}, \code{\link[graphics]{legend}}, \code{\link[base]{floor}}, \code{\link[base]{ceiling}}, \code{\link[base]{floor}}, \code{\link[cladoRcpp]{numstates_from_numareas}}, \code{\link[base]{system.file}}, \code{\link[base]{list.files}}
#' @note Go (BioGeo)BEARS!
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu} 
#' @references
#' \url{https://code.google.com/p/lagrange/}
#' @examples
#' test=1
#' # See example script at PhyloWiki for example
plot_BioGeoBEARS_results <- function(results_object, analysis_titletxt=NULL, addl_params=list(), plotwhat="text", label.offset=NULL, tipcex=0.8, statecex=0.7, splitcex=0.6, titlecex=0.8, plotsplits=TRUE, plotlegend=FALSE, legend_ncol=NULL, legend_cex=1, cornercoords_loc="auto", tr=NULL, tipranges=NULL, if_ties="takefirst", pie_tip_statecex=0.7, juststats=FALSE, xlab="Millions of years ago", root.edge=TRUE, colors_list_for_states=NULL, skiptree=FALSE, show.tip.label=TRUE, tipcol="black", dej_params_row=NULL, plot_max_age=NULL, skiplabels=FALSE, plot_stratum_lines=TRUE, include_null_range=NULL, plot_null_range=FALSE, simplify_piecharts=FALSE, tipboxes_TF=TRUE, tiplabel_adj=c(0.5), no.margin=FALSE, xlims=NULL, ylims=NULL)
	{
	
	junk='
	# manual_ranges_txt=NULL, 
	# @manual_ranges_txt If you dont want to use the default text for each range, produced
	# by areas_list_to_states_list_new(), specify the list here.

	
	scriptdir = "/Dropbox/_njm/__packages/BioGeoBEARS_setup/inst/extdata/a_scripts/"
	plot_BioGeoBEARS_results(results_object, analysis_titletxt=NULL, addl_params=list(), plotwhat="text", label.offset=NULL, tipcex=0.8, statecex=0.8, splitcex=0.8, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=NULL, tipranges=NULL)
	
	# Defaults
	addl_params=list("j"); plotwhat="text"; label.offset=0.45; tipcex=0.7; statecex=0.7; splitcex=0.6; titlecex=0.8; plotsplits=TRUE; cornercoords_loc=scriptdir; include_null_range=TRUE; tr=tr; tipranges=tipranges; juststats = FALSE; plotlegend=FALSE; 	xlab="Millions of years ago"; if_ties="takefirst"
	
	
	# Setup
results_object = resDEC
analysis_titletxt ="BioGeoBEARS DEC on Mariana M1v4_unconstrained"
addl_params=list("j"); plotwhat="text"; label.offset=0.45; tipcex=0.7; statecex=0.7; splitcex=0.6; titlecex=0.8; plotsplits=TRUE; cornercoords_loc=scriptdir; include_null_range=TRUE; tr=tr; tipranges=tipranges
juststats=FALSE; plotlegend=FALSE; 	xlab="Millions of years ago"; if_ties="takefirst"
show.tip.label=TRUE
tipcol="black"; dej_params_row=NULL; plot_max_age=NULL; skiplabels=FALSE; 
colors_list_for_states=NULL
skiptree=FALSE
include_null_range=NULL
plot_stratum_lines=TRUE
	plot_null_range = FALSE
	' # endjunk
	
	# Default; no longer used
	if (is.null(include_null_range) == TRUE)
		{
		include_null_range = results_object$inputs$include_null_range
		}
	
	# Force this in, if user-specified
	results_object$inputs$include_null_range = include_null_range
	
	#######################################################
	# User modifications to border color (externally, using
	# 'par(fg=NA)' or whatever
	#######################################################
	# border color (for modifying this)
	tmp_fg = par("fg")
	par(fg="black")	# set to default for most things

	
	#######################################################
	# Plot ancestral states - DEC
	#######################################################


	# Setup
	#results_object = resDEC_strat
	BioGeoBEARS_run_object = results_object$inputs

	# Read the tree from file, if needed
	if (is.null(tr))
		{
		#tr = read.tree(BioGeoBEARS_run_object$trfn)
		tr = check_trfn(trfn=BioGeoBEARS_run_object$trfn)
		}
	tr_pruningwise = reorder(tr, "pruningwise")
	
	# Basic tree info
	tips = 1:length(tr_pruningwise$tip.label)
	nodes = (length(tr_pruningwise$tip.label)+1):(length(tr_pruningwise$tip.label)+tr_pruningwise$Nnode)


	
	# Read the tipranges from file, if needed.
	if (is.null(tipranges))
		{
		# Get geographic ranges at tips
		if (BioGeoBEARS_run_object$use_detection_model == FALSE)
			{
			tipranges = getranges_from_LagrangePHYLIP(lgdata_fn=np(BioGeoBEARS_run_object$geogfn))
			}
		if (BioGeoBEARS_run_object$use_detection_model == TRUE)
			{
			if (BioGeoBEARS_run_object$use_detection_model == TRUE)
				{
				tipranges = tipranges_from_detects_fn(detects_fn=BioGeoBEARS_run_object$detects_fn)
				} # END if (inputs$use_detection_model == TRUE)
			} # END if (BioGeoBEARS_run_object$use_detection_model == TRUE)
		} # END if (is.null(tipranges))
	
	# Basic areas info
	areas = getareas_from_tipranges_object(tipranges)
	areas

	numareas = length(areas)
	numareas
	
	if (!is.na(results_object$inputs$max_range_size))
		{
		max_range_size = results_object$inputs$max_range_size
		} else {
		max_range_size = length(areas)
		}
	max_range_size


	if (is.null(results_object$inputs$states_list))
		{
		numstates = numstates_from_numareas(numareas=length(areas), maxareas=max_range_size, include_null_range=results_object$inputs$include_null_range)
		numstates
		states_list_areaLetters = areas_list_to_states_list_new(areas, maxareas=max_range_size, include_null_range=results_object$inputs$include_null_range)
		#states_list
		states_list_0based_index = rcpp_areas_list_to_states_list(areas, maxareas=max_range_size, include_null_range=results_object$inputs$include_null_range)
		#states_list_0based_index
		} else {
		states_list_0based_index = results_object$inputs$states_list
		#states_list = states_list_indexes_to_areastxt(states_list=states_list_0based_index, areanames=areas, counting_base=0, concat=FALSE, sep="")
		}


	# calculate the ML marginal probs of states at the base of each branch
	# above each split (local, non-joint probs, global model)
	# 2014-05-15_NJM: Used to add:
	# results_object$ML_marginal_prob_each_split_at_branch_bottom_BELOW_node = ML_marginal_prob_each_split_at_branch_bottom_BELOW_node / rowSums(ML_marginal_prob_each_split_at_branch_bottom_BELOW_node)
	# ... but this is now totally pointless since this is done automatically
	#results_object = get_MLsplitprobs_from_results(results_object)
	#names(results_object)

	# Extract ML parameter values, and LnL
	# This will work with optim, optimx2012, or optimx2013
	
	# Handy summary outputs
	param_ests = extract_params_from_BioGeoBEARS_results_object(results_object, returnwhat="table", addl_params=addl_params, paramsstr_digits=4)

	
	# If you want to skip the plotting and just extract
	# the parameter values
	if (juststats == TRUE)
		{
		return(param_ests)		
		} else {
		paramstr = extract_params_from_BioGeoBEARS_results_object(results_object, returnwhat="string", addl_params=addl_params, paramsstr_digits=4)
		} # if (juststats == TRUE)

	# Get the parameter names
	param_names = extract_params_from_BioGeoBEARS_results_object(results_object, returnwhat="param_names", addl_params=addl_params, paramsstr_digits=4)

		
	# Major title (short description)
	if (is.null(analysis_titletxt))
		{
		tmptxt = results_object$inputs$description
		if (any(is.null(tmptxt), tmptxt=="", tmptxt=="defaults", tmptxt=="default"))
			{
			analysis_titletxt = ""
			} else {
			analysis_titletxt = results_object$inputs$description
			}
		}
	
	
	if (is.null(dej_params_row))
		{
		# Default text for an inference or stochastic mapping
		analysis_titletxt = paste(analysis_titletxt, "\n", "ancstates: global optim, ", max_range_size, " areas max. ", paramstr, sep="")
		analysis_titletxt
		} else {
		# Text for an SSE simulation
		dej_params_row

		brate_col_TF = names(dej_params_row) == "brate"
		brate_col = (1:length(dej_params_row))[brate_col_TF]
		biogeog_params = dej_params_row[1:(brate_col-1)]
		biogeog_param_names = names(dej_params_row)[1:(brate_col-1)]
		equals_col = "="
		
		tmpcols = cbind(biogeog_param_names, equals_col, unlist(biogeog_params))
		tmpcols
		txtrows = apply(X=tmpcols, MARGIN=1, FUN=paste, sep="", collapse="")
		txtrows
		biogeog_params_txt = paste(txtrows, sep="", collapse="; ")
		biogeog_params_txt
			
		titletxt2 = bquote(paste(.(max_range_size), " areas max., ", .(biogeog_params_txt), "; ", lambda, "=", .(dej_params_row$brate), "; ",  mu, "=", .(dej_params_row$drate), "; ", alpha, "=", .(dej_params_row$rangesize_b_exponent), "; ", omega, "=", .(dej_params_row$rangesize_d_exponent), "", sep=""))
		
		#print(titletxt2)
		
		} # END if (is.null(dej_params_row))




	#######################################################
	# Get the marginal probs of the splits (global ML probs, not local)
	# (These are marginal, rather than joint probs; but not local optima)
	#######################################################
	leftright_nodes_matrix = get_leftright_nodes_matrix_from_results(tr_pruningwise)

	# This gets you the prob. of each state at the left base above each node, and
	# at the right base above each node
	marprobs = results_object$ML_marginal_prob_each_state_at_branch_bottom_below_node
	left_ML_marginals_by_node = marprobs[leftright_nodes_matrix[, 2], ]
	right_ML_marginals_by_node = marprobs[leftright_nodes_matrix[, 1], ]
	right_ML_marginals_by_node

	# If they aren't matrices (because of a 2-species, 1-internal-node tree), fix that
	if (is.null(dim(left_ML_marginals_by_node)))
		{
		left_ML_marginals_by_node = matrix(data=left_ML_marginals_by_node, nrow=1)
		}
	if (is.null(dim(right_ML_marginals_by_node)))
		{
		right_ML_marginals_by_node = matrix(data=right_ML_marginals_by_node, nrow=1)
		}


	#######################################################
	# Extract the outputs ancestral states at nodes, and plot!
	#######################################################
	relprobs_matrix = results_object$ML_marginal_prob_each_state_at_branch_top_AT_node
	
	if (length(nodes) > 1)
		{
		relprobs_matrix_for_internal_states = relprobs_matrix[nodes,]	# subset to just internal nodes
		} else {
		relprobs_matrix_for_internal_states = relprobs_matrix[nodes,]	# subset to just internal nodes
		# Convert back to a matrix
		relprobs_matrix_for_internal_states = matrix(data=relprobs_matrix_for_internal_states, nrow=1, ncol=ncol(relprobs_matrix))
		}
	
	
	relprobs_matrix
	
	if (is.null(states_list_0based_index))
		{
		statenames = areas_list_to_states_list_new(areas, maxareas=max_range_size, include_null_range=results_object$inputs$include_null_range, split_ABC=FALSE)
		ranges_list = as.list(statenames)
		statenames
		} else {
		ranges_list = states_list_0based_to_ranges_txt_list(state_indices_0based=states_list_0based_index, areanames=areas)
		ranges_list
		statenames = unlist(ranges_list)
		statenames
		}


	MLprobs = get_ML_probs(relprobs_matrix)
	MLprobs
	MLstates = get_ML_states_from_relprobs(relprobs_matrix, statenames, returnwhat="states", if_ties=if_ties)


	# Set up colors for each state
	if (is.null(colors_list_for_states))
		{
		# Fix plot_null_range to FALSE (don't want to plot that color)
		colors_matrix = get_colors_for_numareas(length(areas))
		colors_list_for_states = mix_colors_for_states(colors_matrix, states_list_0based_index, plot_null_range=results_object$inputs$include_null_range)
		colors_list_for_states
		} # END if (is.null(colors_list_for_states))

	# Set up colors by possible ranges
	if (is.null(ranges_list))
		{
		possible_ranges_list_txt = areas_list_to_states_list_new(areas,  maxareas=max_range_size, split_ABC=FALSE, include_null_range=results_object$inputs$include_null_range)
		} else {
		possible_ranges_list_txt = ranges_list
		} # if (is.null(ranges_list))
	#possible_ranges_list_txt = areas_list_to_states_list_new(areas,  maxareas=max_range_size, split_ABC=FALSE, include_null_range=include_null_range)

# 	if (plot_null_range == FALSE)
# 		{
# 		possible_ranges_list_txt[possible_ranges_list_txt == "_"] = NULL
# 		possible_ranges_list_txt[possible_ranges_list_txt == ""] = NULL
# 		possible_ranges_list_txt[possible_ranges_list_txt == "NA"] = NULL
# 		possible_ranges_list_txt[is.na(possible_ranges_list_txt)] = NULL
# 		possible_ranges_list_txt[is.null(possible_ranges_list_txt)] = NULL
# 		}

	cols_byNode = rangestxt_to_colors(possible_ranges_list_txt, colors_list_for_states, MLstates)

	# Legend, if desired
	if (plotlegend == TRUE)
		{
		colors_legend(possible_ranges_list_txt, colors_list_for_states, legend_ncol=legend_ncol, legend_cex=legend_cex)
		}
	
	
	# Put in a 0 for root.edge
	if (root.edge == FALSE)
		{
		tr$root.edge = 0
		} # END if (root.edge == FALSE)
	if (root.edge == TRUE)
		{
		if (is.null(tr$root.edge) == TRUE)
			{
			tr$root.edge = 0
			} # END if (is.null(tr$root.edge) == TRUE)
		} # END if (root.edge == TRUE)

	# Default label offset
	if (is.null(label.offset))
		{
		label.offset = 0.007 * (get_max_height_tree(tr) + tr$root.edge)
		}
		
	# Manual changing of xlims for phylogeny plot, with plot_max_age
	if (show.tip.label == TRUE)
		{
		# If plot_max_age *IS NOT* specified
		if (is.null(plot_max_age))
			{
			max_x = 1.25*(get_max_height_tree(tr) + tr$root.edge)
			min_x = 0
			} else {
			# If plot_max_age *IS* specified
			nontree_part_of_x = plot_max_age - (get_max_height_tree(tr) + tr$root.edge)
			max_x = 1.25*(get_max_height_tree(tr) + tr$root.edge)
			min_x = -1 * nontree_part_of_x
			}
		} else {
		# NO tip labels
		if (is.null(plot_max_age))
			{
			# If plot_max_age *IS NOT* specified
			max_x = 1.05*(get_max_height_tree(tr) + tr$root.edge)
			min_x = 0
			} else {
			# If plot_max_age *IS* specified
			nontree_part_of_x = plot_max_age - (get_max_height_tree(tr) + tr$root.edge)
			max_x = 1.05*(get_max_height_tree(tr) + tr$root.edge)
			min_x = -1 * nontree_part_of_x
			}
		} # if (show.tip.label == TRUE)


	###################################################
	# Calculate x-axis ticks (axisPhylo alternative)
	###################################################	
	max_tree_x = 1.0 * (get_max_height_tree(tr) + tr$root.edge)
	
	# Plot min/max
	if (is.null(xlims))
		{
		xlims = c(min_x, max_x)
		} else {
		xlims = xlims
		}
	#print(xlims)

	# Tree min/max
	nodecoords = node_coords(tr, tmplocation=cornercoords_loc, root.edge=root.edge)
	max_tree_x = max(nodecoords$x)
	
	# AxisPhylo() min/max
	if (is.null(plot_max_age))
		{
		xticks_desired_lims = c(0, max_tree_x)
		} else {
		xticks_desired_lims = c(0, plot_max_age)
		}

	xticks_desired = pretty(xticks_desired_lims)
	
	# Translate into plot coordinates
	xaxis_ticks_locs = max_tree_x - xticks_desired
	
	#print(xticks_desired)
	#print(xaxis_ticks_locs)
	###################################################	
		
	# Skip tree plotting if it has already been done:
	if (skiptree != TRUE)
		{
		#######################################################
		# Otherwise, plot the tree!!
		#######################################################
		plot(tr_pruningwise, x.lim=xlims, y.lim=ylims, show.tip.label=FALSE, label.offset=label.offset, cex=tipcex, no.margin=no.margin, root.edge=root.edge)
		
		if (show.tip.label == TRUE)
			{
			tiplabels_to_plot = sapply(X=tr_pruningwise$tip.label, FUN=substr, start=1, stop=30)
			if (skiplabels == FALSE)
				{
				tiplabels(text=tiplabels_to_plot, tip=tips, cex=tipcex, adj=0, bg="white", frame="n", pos=4, offset=label.offset, col=tipcol)	# pos=4 means labels go to the right of coords
				} # END if (skiplabels == FALSE)
			} # END if (show.tip.label == TRUE)
		
		#axisPhylo(cex.axis=titlecex)
		axis(side=1,  at=xaxis_ticks_locs, labels=xticks_desired)
		# Plot the title
		mtext(text=xlab, side=1, line=2, cex=titlecex)
		}
	
	# Add states / piecharts
	if (plotwhat == "text")
		{
		# Use statecex for pie chart size at nodes AND states at tips
		par(fg=tmp_fg)	# so that user can change border externally
		
		if (skiplabels == FALSE)
			{
			nodelabels(text=MLstates[nodes], node=nodes, bg=cols_byNode[nodes], cex=statecex)		
			tiplabels(text=MLstates[tips], tip=tips, bg=cols_byNode[tips], cex=statecex, adj=tiplabel_adj)
			} # END if (skiplabels == FALSE)
		
		par(fg="black")	# set to default for most things
		}
	if (plotwhat == "pie")
		{
		# Use statecex for pie chart size at nodes BUT for states at tips,
		# use "pie_tip_statecex"
		par(fg=tmp_fg)	# so that user can change border externally

		if (skiplabels == FALSE)
			{
			# DOSIMPLIFY PIE CHARTS
			if (simplify_piecharts == TRUE)
				{
				# columns to keep in the final
				colnums_to_keep_in_probs = NULL
				
				# Probs table
				probs = results_object$ML_marginal_prob_each_state_at_branch_top_AT_node
				probs2 = probs
				maxprob = rep(0, nrow(probs))
				other = rep(0, nrow(probs))
				num_to_keep = 1
				cat("\nSince simplify_piecharts==TRUE, reducing prob pie charts to (most probable, other)...\n")
				for (i in 1:nrow(probs))
					{
					cat(i, " ", sep="")
					tmprow = probs[i,]
					positions_highest_prob_to_lowest = rev(order(tmprow))
					# If there are ties, we take the first ones
					positions_to_keep = positions_highest_prob_to_lowest[1:num_to_keep]
					colnums_to_keep_in_probs = c(colnums_to_keep_in_probs, positions_to_keep)
					keepTF = rep(FALSE, length(tmprow))
					keepTF[positions_to_keep] = TRUE
	
					# Sum the others
					otherTF = keepTF == FALSE
					other[i] = sum(tmprow[otherTF])
					tmprow[otherTF] = 0
					probs2[i,] = tmprow
					} # END for (i in 1:nrow(probs))
				cat("\n")
				
				colnums_to_keep_in_probs_in_order = sort(unique(colnums_to_keep_in_probs))
				probs3 = cbind(probs2[,colnums_to_keep_in_probs_in_order], other)
				# Subset to just internal nodes
				probs3 = probs3[nodes,]
				
				newcols = c(colors_list_for_states[colnums_to_keep_in_probs_in_order], "white")

				# DO SIMPLIFY PIE CHARTS
				nodelabels(pie=probs3, node=nodes, piecol=newcols, cex=statecex)
				} else {
				# DON'T SIMPLIFY PIE CHARTS
				nodelabels(pie=relprobs_matrix_for_internal_states, node=nodes, piecol=colors_list_for_states, cex=statecex)
				} # END if (simplify_piecharts == TRUE)

			# Plot the tiplabels, if desired
			if (tipboxes_TF == TRUE)
				{
				tiplabels(text=MLstates[tips], tip=tips, bg=cols_byNode[tips], cex=pie_tip_statecex, adj=tiplabel_adj)
				} # END if (tipboxes_TF = TRUE)
			} # END if (skiplabels == FALSE)

		par(fg="black")	# set to default for most things
		}
	
	if (skiptree != TRUE)
		{
		if (titlecex > 0)
			{
			#par(ps = 12, cex = titlecex, cex.main = titlecex)
			par(cex.main = titlecex)
			title(analysis_titletxt)
			if (!is.null(dej_params_row))
				{
				# Subtitle for SSE simulation plots
				title(titletxt2, line=1)
				#print(titletxt2)
				}
			#par("font.main") = 2
			}
		}

	if (plotsplits == TRUE)
		{
		# Error check; users must specify the location of the function "plot_phylo3_nodecoords"
		if (cornercoords_loc == "manual")
			{
			stoptxt = cat("\nNOTE: To plot splits, this function needs to access the function 'plot_phylo3_nodecoords'.\n",
							"The function is modified from an APE function, and cannot be directly included in the package,\n",
							"due to some C code that does not meet CRAN standards. To solve this, give plot_BioGeoBEARS_results\n",
							"a 'cornercoords_loc' string that gives the directory of plot_phylo3_nodecoords.R.  Typically this\n",
							"can be found via: ", 'tmp=np(system.file("extdata/a_scripts", package="BioGeoBEARS"))\n',
							"then: list.files(tmp); print(tmp)\n", sep="")
			plotsplits = FALSE
			}
		}

	if (plotsplits == TRUE)
		{
		#######################################################
		# Also add the splits to the plot
		#######################################################
		# First, get the corner coordinates
		coords_df = corner_coords(tr, tmplocation=cornercoords_loc, root.edge=root.edge)
		coords_df

		# LEFT SPLITS
		relprobs_matrix = left_ML_marginals_by_node
		
		if (plotwhat == "text")
			{
			MLprobs = get_ML_probs(relprobs_matrix)
			MLprobs
			MLstates = get_ML_states_from_relprobs(relprobs_matrix, statenames, returnwhat="states", if_ties=if_ties)
			MLstates
			length(MLstates)
		
			# Set up colors
			#possible_ranges_list_txt = areas_list_to_states_list_new(areas,  maxareas=max_range_size, split_ABC=FALSE, include_null_range=include_null_range)
			cols_byNode = rangestxt_to_colors(possible_ranges_list_txt, colors_list_for_states, MLstates)
			
			par(fg=tmp_fg)	# so that user can change border externally
			if (skiplabels == FALSE)
				{			
				cornerlabels(text=MLstates, coords=coords_df$leftcorns, bg=cols_byNode, cex=splitcex)
				} # END if (skiplabels == FALSE)

			par(fg="black")	# set to default for most things
			}
		
		if (plotwhat == "pie")
			{
			par(fg=tmp_fg)	# so that user can change border externally
			cornerpies(pievals=relprobs_matrix, coords=coords_df$leftcorns, piecol=colors_list_for_states, cex=splitcex)
			par(fg="black")	# set to default for most things
			}



		# RIGHT SPLITS
		relprobs_matrix = right_ML_marginals_by_node

		if (plotwhat == "text")
			{
			MLprobs = get_ML_probs(relprobs_matrix)
			MLprobs
			MLstates = get_ML_states_from_relprobs(relprobs_matrix, statenames, returnwhat="states", if_ties=if_ties)
			MLstates
			length(MLstates)

			# Set up colors
			#possible_ranges_list_txt = areas_list_to_states_list_new(areas,  maxareas=max_range_size, split_ABC=FALSE, include_null_range=include_null_range)
			cols_byNode = rangestxt_to_colors(possible_ranges_list_txt, colors_list_for_states, MLstates)

			par(fg=tmp_fg)	# so that user can change border externally
			if (skiplabels == FALSE)
				{
				cornerlabels(text=MLstates, coords=coords_df$rightcorns, bg=cols_byNode, cex=splitcex)
				} # END if (skiplabels == FALSE)
			par(fg="black")	# set to default for most things
			}

		if (plotwhat == "pie")
			{
			par(fg=tmp_fg)	# so that user can change border externally
			cornerpies(pievals=relprobs_matrix, coords=coords_df$rightcorns, piecol=colors_list_for_states, cex=splitcex)			
			par(fg="black")	# set to default for most things
			}
		}

	# Plot vertical dashed lines for timeperiods
	# Is it time-stratified? Plot lines if desired
	if ( ((is.null(BioGeoBEARS_run_object$timeperiods) == FALSE)) && (plot_stratum_lines==TRUE) )
		{
		timeperiods = BioGeoBEARS_run_object$timeperiods
		line_positions_on_plot = add_statum_boundaries_to_phylo_plot(tr, timeperiods=timeperiods, lty="dashed", col="gray50", plotlines=TRUE)
		} # END plot vertical dashed lines for timeperiods



	# Handy summary outputs
	param_ests = matrix(data=param_ests, nrow=1)
	param_ests = adf2(param_ests)
	
	param_ests = dfnums_to_numeric(param_ests)
	names(param_ests) = c("LnL", "nparams", param_names)
	
	return(param_ests)
	} # END plot_BioGeoBEARS_results






#######################################################
# Colors and legends
#######################################################


#######################################################
# get_colors_for_numareas
#######################################################
#' Get colors for a certain number of single areas
#' 
#' Get colors for a given number of single areas. I have
#' chosen default colors that are very bright/primary.
#' 
#' @param numareas The number of areas
#' @param use_rainbow If TRUE, force use of \code{rainbow()}
#' @return \code{colors_matrix} The colors for the single areas, 1 column per area
#' @export
#' @seealso \code{\link[stats]{optim}}
#' @note Go (BioGeo)BEARS!
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu} 
#' @references
#' \url{http://phylo.wikidot.com/matzke-2013-international-biogeography-society-poster}
#' @examples
#' testval=1
#' 
get_colors_for_numareas = function(numareas, use_rainbow=FALSE)
	{
	# default
	area_colors = sapply(X=(1:numareas)/numareas, FUN=gray)
	
	if (use_rainbow == TRUE)
		{
		area_colors = rainbow(numareas)
		colors_matrix = col2rgb(area_colors)
		return(colors_matrix)
		}
	
	if (numareas == 2)
		{
		area_colors = c("blue", "green")
		}
	if (numareas == 3)
		{
		area_colors = c("blue", "green", "red")
		}
	if (numareas == 4)
		{
		area_colors = c("blue", "green3", "yellow", "red")
		}
	if (numareas == 5)
		{
		area_colors = c("blue", "cyan", "green3", "yellow", "red")
		}
	if (numareas == 6)
		{
		area_colors = c("blue", "cyan", "green3", "yellow", "red", "violet")
		}
	if (numareas == 7)
		{
		area_colors = c("blue", "cyan", "green3", "yellow", "orange", "red", "violet")
		}
	if (numareas > 7)
		{
		area_colors = col2rgb(rainbow(numareas))
		return(area_colors)
		}
	
	colors_matrix = col2rgb(area_colors)
	return(colors_matrix)	
	}





#######################################################
# mix_colors_for_states
#######################################################
#' Mix colors logically to produce colors for multi-area ranges
#' 
#' Mixex colors logically to produce colors for multi-area ranges. 
#' Because there are so many possible ranges/states in most analyses,
#' and because the human eye can only easily read <10 colors on a
#' map (particularly if you are anomalous trichromatic, like me!)
#' devising a reasonable color setup is nontrivial. The idea here
#' is that, when a range occupies multiple colors, the colors of
#' the individual areas will be mixed. If a range contains all
#' areas, it is colored white. 
#'
#' Sometimes people try to make legends listing all of the dozens
#' of colors used, even for the multi-area "mixed" colors. However,
#' as a graphical matter, I would recommend just giving the colors
#' of the single areas, and any particularly common/important
#' multi-area ranges. If the ancestral ranges are labeled with e.g.
#' "A", "ABC", etc. (as is typical), there is limited 
#' value in a separate legend for the colors anyway.
#' 
#' @param colors_matrix A column with a color for each single area
#' @param states_list_0based_index States list giving areas, 0-based
#' @param plot_null_range If FALSE, null ranges are excluded (however coded). 
#' Default FALSE, and should be false unless you *really* want to plot null range.
#' @return \code{colors_list_for_states} The colors for the ML states
#' @export
#' @seealso \code{\link[stats]{optim}}
#' @note Go (BioGeo)BEARS!
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu} 
#' @references
#' \url{http://phylo.wikidot.com/matzke-2013-international-biogeography-society-poster}
#' @examples
#' testval=1
#' 
mix_colors_for_states <- function(colors_matrix, states_list_0based_index, plot_null_range=FALSE)
	{
	# Eliminate the null range, if present
	if (plot_null_range == FALSE)
		{
		states_list_0based_index[states_list_0based_index == "_"] = NULL
		states_list_0based_index[states_list_0based_index == ""] = NULL
		states_list_0based_index[states_list_0based_index == "NA"] = NULL
		states_list_0based_index[is.na(states_list_0based_index)] = NULL
		states_list_0based_index[is.null(states_list_0based_index)] = NULL
		}
	
	
	colors_list_for_states = rep(NA, length(states_list_0based_index))   #matrix(data=NA, nrow=3, ncol=length(states_list_0based_index))
	
	# There is a column for each single area
	numareas = ncol(colors_matrix)
	
	
	# Combine the colors of the input areas, and, if 2 areas or more, divide by (numareas-0.5)
	# (assures no duplication of colors)
	for (i in 1:length(states_list_0based_index))
		{
		tmpareas_in_state_0based_index = states_list_0based_index[[i]]
		tmp_numareas = length(tmpareas_in_state_0based_index)
		
		# Check for null etc.
		# Avoid warning by checking for length = 1 first...
		if (length(tmpareas_in_state_0based_index) == 1)
			{
			if (tmpareas_in_state_0based_index == "_" || tmpareas_in_state_0based_index == "" || is.na(tmpareas_in_state_0based_index) || is.null(tmpareas_in_state_0based_index))
				{
				# NULL/empty range gets black
				colors_list_for_states[i] = rgb(red=0, green=0, blue=0, maxColorValue=255)
				next()
				} # END (tmpareas_in_state_0based_index == "_" ...
			} # END if (length(tmpareas_in_state_0based_index == "_" ||) == 1)
		
		
		
		# Fill in the colors for this area
		
		# Make it white, if its all areas
		if (tmp_numareas == numareas)
			{
			# input white
			colors_list_for_states[i] = rgb(red=255, green=255, blue=255, maxColorValue=255)
			next()
			}
		
		# Otherwise
		sum_color_nums = rep(0, 3)
		for (j in 1:tmp_numareas)
			{
			tmparea_1based_index = 1 + tmpareas_in_state_0based_index[j]
			
			color_nums = colors_matrix[ ,tmparea_1based_index]
			sum_color_nums = sum_color_nums + color_nums
			}
		
		if (tmp_numareas >= 2)
			{
			divide_by = (tmp_numareas - 0.0)
			} else {
			divide_by = 1
			}
		
		sum_color_nums = round(sum_color_nums / divide_by)
		
		# Save the colors as hex format
		colors_list_for_states[i] = rgb(
			red=sum_color_nums[1], 
			green=sum_color_nums[2], 
			blue=sum_color_nums[3], 
			maxColorValue=255)
		}
	
	return(colors_list_for_states)
	}



# Convert a list of ranges text (KOM, MH, KOMIH, etc.) 
#######################################################
# rangestxt_to_colors
#######################################################
#' Convert a list of ranges texts (KOM, MH, KOMIH, etc.) to colors
#' 
#' Converts a list of ranges texts (KOM, MH, KOMIH, etc.) to colors.
#' 
#' @param possible_ranges_list_txt A list of the allowed ranges/states
#' @param colors_list_for_states The corresponding colors
#' @param MLstates The ML states for the internal nodes
#' @return \code{MLcolors} The colors for the ML states
#' @export
#' @seealso \code{\link[stats]{optim}}
#' @note Go (BioGeo)BEARS!
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu} 
#' @references
#' \url{http://phylo.wikidot.com/matzke-2013-international-biogeography-society-poster}
#' @examples
#' testval=1
#' 
#' possible_ranges_list_txt = c("A", "AB", "ABC")
#' colors_list_for_states = c("red", "orange", "white")
#' MLstates = c("ABC", "AB", "AB", "A", "A", "A", "ABC")
#' rangestxt_to_colors(possible_ranges_list_txt, colors_list_for_states, MLstates)
#' 
rangestxt_to_colors <- function(possible_ranges_list_txt, colors_list_for_states, MLstates)
	{
	setup='
	possible_ranges_list_txt = c("A", "AB", "ABC")
	colors_list_for_states = c("red", "orange", "white")
	MLstates = c("ABC", "AB", "AB", "A", "A", "A", "ABC")
	rangestxt_to_colors(possible_ranges_list_txt, colors_list_for_states, MLstates)
	'
	
	if (length(possible_ranges_list_txt) != length(colors_list_for_states))
		{
		cat("\nlength(possible_ranges_list_txt) = ", length(possible_ranges_list_txt))
		cat("\nlength(colors_list_for_states) = ", length(colors_list_for_states))
		cat("\nlength(MLstates) = ", length(MLstates))
		
		
		stop("Error: possible_ranges_list_txt and colors_list_for_states must be the same length")
		}
	
	ranges_colors = 1:length(MLstates)
	
	# Nums
	#nums = 1:length(possible_ranges_list_txt)
	

	match_indices = get_indices_where_list1_occurs_in_list2(list1=MLstates, list2=possible_ranges_list_txt)
	MLcolors = colors_list_for_states[match_indices]
	
	return(MLcolors)
	}




#######################################################
# colors_legend
#######################################################
#' Plot a colors legend for geographic ranges
#' 
#' Plots a colors legend for geographic ranges.
#' 
#' @param possible_ranges_list_txt A list of the allowed ranges/states
#' @param colors_list_for_states The corresponding colors
#' @param legend_ncol The number of columns in the legend.  If \code{NULL} (default), the function calculates \code{floor(sqrt(length(possible_ranges_list_txt) / 2))}.
#' Note that when you have hundreds of states, there is probably no good way to have a coherent legend, and it is easier to just rely upon printing the 
#' character codes for the ML states in the plots, with the colors, and users can then see and trace the common colors/states by eye.
#' @param legend_cex The cex (character expansion size) for the legend.  Defaults to 1, which means the \code{\link[graphics]{legend}} function determines the 
#' size.  The value 2.5 works well for 15 or 16 states/ranges.
#' @location The "location" input goes to the legend() command. Default is "top".
#' @make_blank_plot_first If TRUE (default), colors_legend() will first plot a blank plot
#' with the x- and y-axes ranging from 1-10. make_blank_plot_first=FALSE may be useful for
#' adding the legend on top of other plots.
#' @return Nothing
#' @export
#' @seealso \code{\link[graphics]{legend}}, \code{\link[base]{floor}}, \code{\link[base]{ceiling}}, \code{\link[base]{floor}}
#' @note Go (BioGeo)BEARS!
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu} 
#' @references
#' \url{http://phylo.wikidot.com/matzke-2013-international-biogeography-society-poster}
#' @examples
#' testval=1
#' 
colors_legend <- function(possible_ranges_list_txt, colors_list_for_states, legend_ncol=NULL, legend_cex=1, location="top", make_blank_plot_first=TRUE, ...)
	{

	if (is.null(legend_ncol))
		{
		tmp_numstates = length(possible_ranges_list_txt)
		if (tmp_numstates <= 64)
			{
			legend_ncol = floor(sqrt(tmp_numstates / 2))
			legend_ncol
			} else {
			# For huge numbers of states, just do a square
			legend_ncol = ceiling(sqrt(tmp_numstates))
			}
		}
	
	# Plot, no borders (bty="n"), no labels (xlab, ylab), no tick marks (xaxt, yaxt)
	if (make_blank_plot_first == TRUE)
		{
		plot(1:10, 1:10, pch=".", col="white", xaxt="n", yaxt="n", xlab="", ylab="", bty="n")
		}
	#lines(1:4,4:1, col="blue") 
	#legend("top", leg=c("a","b"),col=c("black","blue"), fill=TRUE) 
	legend(x=location, legend=possible_ranges_list_txt, fill=colors_list_for_states, ncol=legend_ncol, title="Legend", cex=legend_cex, ...)#, fill=TRUE) 
	
	} # END colors_legend <- function(possible_ranges_list_txt, colors_list_for_states, legend_ncol=NULL, legend_cex=1, ...)










#######################################################
# map_LGpy_MLsplits_to_tree
#######################################################
#' Take the table of ML splits and node number and map on tree (Python version)
#' 
#' Given a table of splits probabilities from \code{\link{LGpy_splits_fn_to_table}}, map the splits on the tree.
#' 
#' See \code{\link{get_lagrange_nodenums}} for connecting these node numbers to APE node numbers.
#' 
#' @param MLsplits_LGpy A data.frame containing the node numbers, splits, and split probabilities.
#' @param tr An ape phylo object
#' @param tiprange_names The geographic ranges at the tips (i.e. the input data)
#' @return \code{MLsplits_LGpy} A data.frame containing the node numbers, ML splits, and split probabilities; reordered for this plot
#' @export
#' @seealso \code{\link{get_lagrange_nodenums}}, \code{\link{LGpy_splits_fn_to_table}}, \code{\link{LGcpp_splits_fn_to_table}}
#' @note Go (BioGeo)BEARS!
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu} 
#' @references
#' \url{http://phylo.wikidot.com/matzke-2013-international-biogeography-society-poster}
#' \url{https://code.google.com/p/lagrange/}
#' @examples
#' test=1
#' 
map_LGpy_MLsplits_to_tree <- function(MLsplits_LGpy, tr, tiprange_names)
	{
	# Order in LAGRANGE order
	MLsplits_LGpy = MLsplits_LGpy[order(MLsplits_LGpy$nodenum_LGpy),]
	MLsplits_LGpy
	
	# Get the names of the columns in the MLsplits table
	tmpnames = names(MLsplits_LGpy)

	# Add R node numbering, if not already done
	if (("Rnodes" %in% tmpnames) == FALSE)
		{
		downpass_node_matrix = get_lagrange_nodenums(tr)
		#downpass_node_matrix[,2] = 1:18
		#downpass_node_matrix = downpass_node_matrix[order(downpass_node_matrix[,1]), ]
	
		MLsplits_LGpy = cbind(MLsplits_LGpy, downpass_node_matrix)
		names(MLsplits_LGpy) = c(tmpnames, "Rnodes", "LGnodes")
		MLsplits_LGpy
		}
	
	MLsplits_LGpy = MLsplits_LGpy[order(MLsplits_LGpy$Rnodes), ]
	MLsplits_LGpy

	# Plot them	
	par(mfrow=c(2,1))
	plot(tr, label.offset=0.15)
	nodelabels(text=MLsplits_LGpy$splits, node=20:37)
	tiplabels(tiprange_names)
	title("LAGRANGE (python) ML splits")
	axisPhylo()
	mtext(text="million years ago", side=1, line=2)
	
	plot(tr, label.offset=0.15)
	pievals = as.matrix(MLsplits_LGpy[,c("relprob","relprob2")])
	nodelabels(pie=pievals, piecol=c("blue", "white"))
	tiplabels(tiprange_names)
	title("LAGRANGE (python) split probs")
	axisPhylo()
	mtext(text="million years ago", side=1, line=2)

	return(MLsplits_LGpy)
	}



#######################################################
# map_LG_MLsplits_to_tree
#######################################################
#' Take the table of ML splits and node number and map on tree (C++ LAGRANGE version)
#' 
#' Given a table of splits probabilities from \code{\link{LGcpp_splits_fn_to_table}}, map the splits on the tree.
#' 
#' See \code{\link{get_lagrange_nodenums}} for connecting these node numbers to APE node numbers.
#' 
#' @param MLsplits_LGcpp A data.frame containing the node numbers, splits, and split probabilities.
#' @param tr An ape phylo object
#' @param tiprange_names The geographic ranges at the tips (i.e. the input data)
#' @param removechar The character to remove, if needed.
#' @param type The type of LAGRANGE input (default C++)
#' @return \code{MLsplits_LGcpp} A data.frame containing the node numbers, ML splits, and split probabilities; reordered for this plot.
#' @export
#' @seealso \code{\link{get_lagrange_nodenums}}, \code{\link{LGpy_splits_fn_to_table}}, \code{\link{LGcpp_splits_fn_to_table}}
#' @note Go (BioGeo)BEARS!
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu} 
#' @references
#' \url{http://phylo.wikidot.com/matzke-2013-international-biogeography-society-poster}
#' \url{https://code.google.com/p/lagrange/}
#' @examples
#' test=1
#' 
map_LG_MLsplits_to_tree <- function(MLsplits_LGcpp, tr, tiprange_names, removechar=NULL, type="C++")
	{
	defaults='
	MLsplits_LGpy, tr=tr, tiprange_names=tiprange_names, type="python"
	'
	# Order in LAGRANGE order
	MLsplits_LGcpp = order_LGnodes(MLsplits_LGcpp, tr, removechar, type)
	
	
	# Plot them	
	par(mfrow=c(2,1))
	plot(tr, label.offset=0.15)
	nodelabels(text=MLsplits_LGcpp$splits, node=20:37)
	tiplabels(tiprange_names)
	title(paste("LAGRANGE (", type, ") ML splits", sep=""))
	axisPhylo()
	mtext(text="million years ago", side=1, line=2)
	
	plot(tr, label.offset=0.15)
	pievals = as.matrix(MLsplits_LGcpp[,c("relprob","relprob2")])
	nodelabels(pie=pievals, piecol=c("blue", "white"))
	tiplabels(tiprange_names)
	title(paste("LAGRANGE (", type, ") split probs", sep=""))
	axisPhylo()
	mtext(text="million years ago", side=1, line=2)

	return(MLsplits_LGcpp)
	}


#######################################################
# get_statesColors_table
#######################################################
#' Make a color table for each area and their combinations
#' 
#' Given a list of areas, make a color table for the various combinations.
#' 
#' @param areanames A list of the area names.
#' @param plot_null_range If FALSE, the null range is excluded from the state space. 
#' Note: The null range is never plotted, regardless of the setting, but the function 
#' needs to know which, for proper counting/coloring.
#' @return \code{statesColors_table} A table giving the colors for each state.
#' @export
#' @seealso \code{\link{get_lagrange_nodenums}}, \code{\link{LGpy_splits_fn_to_table}}, \code{\link{LGcpp_splits_fn_to_table}}
#' @note Go (BioGeo)BEARS!
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu} 
#' @references
#' \url{http://phylo.wikidot.com/matzke-2013-international-biogeography-society-poster}
#' \url{https://code.google.com/p/lagrange/}
#' @examples
#' test=1
#' 
get_statesColors_table <- function(areanames=c("K","O","M","H"), plot_null_range=FALSE)
	{
	# Make the color matrix for the individual areas
	colors_matrix = get_colors_for_numareas(numareas=length(areanames), use_rainbow=FALSE)
	colors_matrix
	
	# Get the states
	# Here, include_null_range is fixed to FALSE, so that this state is not displayed
	states_list_0based_index = rcpp_areas_list_to_states_list(areas=areanames, include_null_range=plot_null_range)
	
	colors_list_for_states = mix_colors_for_states(colors_matrix, states_list_0based_index, plot_null_range=plot_null_range)
	colors_list_for_states

	possible_ranges_list_txt = states_list_indexes_to_areastxt(states_list=states_list_0based_index, areanames, counting_base=0, sep="")
	possible_ranges_list_txt
	
	statesColors_table = adf2(cbind(possible_ranges_list_txt, colors_list_for_states))
	names(statesColors_table) = c("range", "color")
	return(statesColors_table)
	}


#######################################################
# map_LG_MLsplits_to_tree_corners
#######################################################
#' Map splits to the corners on a phylogeny
#' 
#' Map splits to the "corners" on a phylogeny. Normally, the corners are meaningless,
#' but in a biogeographical analysis with discrete cladogenesis events at splits on 
#' the tree, the provide a convenient place to plot the state probabilities just
#' after speciation. These will not be the same as the state probabilities at the 
#' nodes (unlike most phylogenetic ancestral state estimations), because of the 
#' "instantaneous" changes in the geographic ranges that can occur at speciation. (The
#' exception is models that assume ranges do not change at speciation, like the 
#' Markov-<i>k</i> model of Lewis 2001, or the \code{BAYAREALIKE} model.
#' 
#' @param MLsplits A data.frame containing the node numbers, splits, and split probabilities.
#' @param tr An ape phylo object
#' @param tipranges Tipranges object
#' @param removechar The character to remove, if needed.
#' @param type The type of LAGRANGE input (default C++)
#' @param statesColors_table If not default, a table with a color for each area combination.
#' @param bgcol The background color
#' @param areanames The area names, if different from those in the tipranges object
#' @param newplot Default TRUE; should there be a new plot, or should the splits be added to another plot?
#' @param ... Additional arguments to standard functions
#' @return \code{MLsplits} The splits table, ordered appropriately.
#' @export
#' @seealso \code{\link{get_lagrange_nodenums}}, \code{\link{LGpy_splits_fn_to_table}}, \code{\link{LGcpp_splits_fn_to_table}}
#' @note Go (BioGeo)BEARS!
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu} 
#' @references
#' \url{http://phylo.wikidot.com/matzke-2013-international-biogeography-society-poster}
#' \url{https://code.google.com/p/lagrange/}
#' @examples
#' test=1
#' 
map_LG_MLsplits_to_tree_corners <- function(MLsplits, tr, tipranges, removechar=NULL, type="C++", statesColors_table="default", bgcol="green3", areanames="default", newplot=TRUE, root.edge, ...)
	{
	defaults='
	MLsplits=MLsplits_LGpy
	type="python"
	'
	# Get corner coordinates
	corners_list = corner_coords(tr, root.edge=root.edge)
	leftcorns = corners_list$leftcorns
	rightcorns = corners_list$rightcorns
	
	# Plot splits on corners
	# Ensure correct order
	MLsplits = order_LGnodes(MLsplits, tr, removechar=removechar, type=type)
	MLsplits
	
	
	# Get the ranges at the tips (in the right order)
	tipranges = order_tipranges_by_tree_tips(tipranges, tr)
	tiprange_names = tipranges_to_area_strings(tipranges=tipranges)
	
	if (areanames == "default")
		{
		areanames = getareas_from_tipranges_object(tipranges=tipranges)
		}
	
	# Colors 
	if (statesColors_table == "default")
		{
		statesColors_table = get_statesColors_table(areanames)
		tmp_index = get_indices_where_list1_occurs_in_list2(list1=MLsplits$leftBB, list2=statesColors_table$range)
		left_colors = statesColors_table$color[tmp_index]
		tmp_index = get_indices_where_list1_occurs_in_list2(list1=MLsplits$rightBB, list2=statesColors_table$range)
		right_colors = statesColors_table$color[tmp_index]
		
		tmp_index = get_indices_where_list1_occurs_in_list2(list1=tiprange_names, list2=statesColors_table$range)
		tip_colors = statesColors_table$color[tmp_index]
		} else {
		left_colors = bgcol
		right_colors = bgcol
		tip_colors = bgcol
		}
	
	# Plot them
	if (newplot == TRUE)
		{
		plot(tr, label.offset=0.15)
		axisPhylo()
		mtext(text="million years ago", side=1, line=2)
		}
	cornerlabels(text=MLsplits$leftBB, coords=leftcorns, bg=left_colors, ...)
	cornerlabels(text=MLsplits$rightBB, coords=rightcorns, bg=right_colors, ...)
	tiplabels(text=tiprange_names, bg=tip_colors, ...)

	return(MLsplits)
	}



#######################################################
# map_LG_MLstates_to_tree
#######################################################
#' Map states to the nodes on a phylogeny
#' 
#' Map the most-probable ancestral states to the nodes on a phylogeny. Note that
#' the most-probable ancestral state may still have a probability less than
#' 50%, if the inference for that node is highly uncertain. It is always 
#' necessary to look at, and interpret, the pie charts in addition to plot of
#' most-probable states. See "mistakes not to make" on PhyloWiki.
#' 
#' @param MLstates_LGcpp A data.frame containing the node numbers, states, and states probabilities.
#' @param tr An ape phylo object
#' @param tipranges Tipranges object
#' @param removechar The character to remove, if needed.
#' @param type The type of LAGRANGE input (default C++)
#' @param statesColors_table If not default, a table with a color for each area combination.
#' @param bgcol The background color
#' @param areanames The area names, if different from those in the tipranges object
#' @param newplot Default TRUE; should there be a new plot, or should the splits be added to another plot?
#' @param ... Additional arguments to standard functions
#' @return \code{MLstates_LGcpp} The states table, ordered appropriately.
#' @export
#' @seealso \code{\link{get_lagrange_nodenums}}, \code{\link{LGpy_splits_fn_to_table}}, \code{\link{LGcpp_splits_fn_to_table}}
#' @note Go (BioGeo)BEARS!
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu} 
#' @references
#' \url{http://phylo.wikidot.com/matzke-2013-international-biogeography-society-poster}
#' \url{https://code.google.com/p/lagrange/}
#' @examples
#' test=1
#' 
map_LG_MLstates_to_tree <- function(MLstates_LGcpp, tr, tipranges, removechar=NULL, type="C++", statesColors_table="default", bgcol="green3", areanames="default", newplot=TRUE, ...)
	{
	# Order in LAGRANGE order
	MLstates_LGcpp = order_LGnodes(MLstates_LGcpp, tr, removechar, type=type, type2="states")
	
	# Get the ranges at the tips (in the right order)
	tipranges = order_tipranges_by_tree_tips(tipranges, tr)
	tiprange_names = tipranges_to_area_strings(tipranges=tipranges)
	
	if (areanames == "default")
		{
		areanames = getareas_from_tipranges_object(tipranges=tipranges)
		}
	
	# Colors 
	if (statesColors_table == "default")
		{
		statesColors_table = get_statesColors_table(areanames)
		tmp_index = get_indices_where_list1_occurs_in_list2(list1=MLstates_LGcpp$states, list2=statesColors_table$range)
		state_colors = statesColors_table$color[tmp_index]
		
		tmp_index = get_indices_where_list1_occurs_in_list2(list1=tiprange_names, list2=statesColors_table$range)
		tip_colors = statesColors_table$color[tmp_index]
		} else {
		left_colors = bgcol
		right_colors = bgcol
		tip_colors = bgcol
		}

	# Plot them
	if (newplot == TRUE)
		{
		plot(tr, label.offset=0.15)
		axisPhylo()
		mtext(text="million years ago", side=1, line=2)
		}

	nodelabels(text=MLstates_LGcpp$states, node=20:37, bg=state_colors, ...)
	tiplabels(tiprange_names, bg=tip_colors, ...)
	#title(paste("LAGRANGE (", type, ") ML states", sep=""))
	
# 	plot(tr, label.offset=0.15)
# 	pievals = as.matrix(MLstates_LGcpp[,c("relprob","relprob2")])
# 	nodelabels(pie=pievals, piecol=c("blue", "white"))
# 	tiplabels(tiprange_names)
# 	title(paste("LAGRANGE (", type, ") state probs", sep=""))
# 	axisPhylo()
# 	mtext(text="million years ago", side=1, line=2)

	return(MLstates_LGcpp)
	}






#######################################################
# order_LGnodes
#######################################################
#' Order LAGRANGE-numbered nodes so that they can be plotted in R
#' 
#' Node numbers between ancestral states are different for ape's \code{phylo}
#' objects, C++ Lagrange,
#' Python Lagrange, and DIVA. This function reorders the Lagrange
#' splits table, for plotting in R. 
#' 
#' @param MLsplits_LGcpp A data.frame containing the node numbers, splits, and split probabilities.
#' @param tr An ape phylo object
#' @param removechar The character to remove, if needed.
#' @param type The type of LAGRANGE input (default C++). Inputs other than C++ assume
#' the Python lagrange ordering.
#' @param type2 "splits" or "states"
#' @return \code{MLsplits} The splits table, ordered appropriately.
#' @export
#' @seealso \code{\link{get_lagrange_nodenums}}, \code{\link{LGpy_splits_fn_to_table}}, \code{\link{LGcpp_splits_fn_to_table}}
#' @note Go (BioGeo)BEARS!
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu} 
#' @references
#' \url{http://phylo.wikidot.com/matzke-2013-international-biogeography-society-poster}
#' \url{https://code.google.com/p/lagrange/}
#' @examples
#' test=1
#' 
order_LGnodes <- function(MLsplits_LGcpp, tr=NULL, removechar=NULL, type="C++", type2="splits")
	{
	# Order in LAGRANGE order
	if (type == "C++")
		{
		MLsplits_LGcpp = MLsplits_LGcpp[order(MLsplits_LGcpp$nodenum_LGcpp),]
		} else {
		MLsplits_LGcpp = MLsplits_LGcpp[order(MLsplits_LGcpp$nodenum_LGpy),]
		}
	MLsplits_LGcpp
	
	# Get the names of the columns in the MLsplits table
	tmpnames = names(MLsplits_LGcpp)

	# Add R node numbering, if not already done
	if (("Rnodes" %in% tmpnames) == FALSE)
		{
		downpass_node_matrix = get_lagrange_nodenums(tr)
		#downpass_node_matrix[,2] = 1:18
		#downpass_node_matrix = downpass_node_matrix[order(downpass_node_matrix[,1]), ]
	
		MLsplits_LGcpp = cbind(MLsplits_LGcpp, downpass_node_matrix)
		names(MLsplits_LGcpp) = c(tmpnames, "Rnodes", "LGnodes")
		MLsplits_LGcpp
		}
	
	MLsplits_LGcpp = MLsplits_LGcpp[order(MLsplits_LGcpp$Rnodes), ]
	MLsplits_LGcpp
	
	# Change e.g. O_M -> OM, K_O_M_H -> KOMH
	if (!is.null(removechar))	
		{
		if (type2 == "splits")
			{
			MLsplits_LGcpp$splits = gsub(pattern=removechar, replacement="", x=MLsplits_LGcpp$splits)
			MLsplits_LGcpp$leftBB = gsub(pattern=removechar, replacement="", x=MLsplits_LGcpp$leftBB)
			MLsplits_LGcpp$rightBB = gsub(pattern=removechar, replacement="", x=MLsplits_LGcpp$rightBB)
			} else {
			MLsplits_LGcpp$states = gsub(pattern=removechar, replacement="", x=MLsplits_LGcpp$states)
			}
		}
	
	return(MLsplits_LGcpp)
	}



#######################################################
# cornerlabels
#######################################################
#' Make labels for plotting ranges on corners
#' 
#' This function makes labels for plotting ranges on corners.
#' 
#' @param text The text to put at the corners.
#' @param coords The coordinates at which to plot the labels
#' @param bg The background color
#' @param col The text color
#' @param adj Position adjustment; default \code{adj=c(0.5,0.5)}
#' @param ... Additional arguments to standard functions
#' @return nothing
#' @export
#' @seealso \code{\link{cornerpies}}, \code{\link{corner_coords}}, \code{\link{get_lagrange_nodenums}}, 
#' \code{\link{LGpy_splits_fn_to_table}}, \code{\link{LGcpp_splits_fn_to_table}}
#' @note Go (BioGeo)BEARS!
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu} 
#' @references
#' \url{http://phylo.wikidot.com/matzke-2013-international-biogeography-society-poster}
#' \url{https://code.google.com/p/lagrange/}
#' @examples
#' test=1
#' 
cornerlabels <- function(text, coords, bg="green3", col="black", adj=c(0.5,0.5), ...)
	{
	args <- list(...)
    CEX <- if ("cex" %in% names(args)) 
    	{
	    args$cex
	    } else {
	    par("cex")
	    }
	   
	# Draw the rectangles	   
	width <- strwidth(text, units = "inches", cex = CEX)
	height <- strheight(text, units = "inches", cex = CEX)

	width <- xinch(width)
	height <- yinch(height)

	XX = coords$x
	YY = coords$y	
	xl <- XX - width * adj[1] - xinch(0.03)
	xr <- xl + width + xinch(0.03)
	yb <- YY - height * adj[2] - yinch(0.02)
	yt <- yb + height + yinch(0.05)
	rect(xl, yb, xr, yt, col = bg)
	
	# Write the text
	text(XX, YY, text, adj = adj, col = col, ...)
	
	return()
	}



#######################################################
# cornerpies
#######################################################
#' Make pie charts for plotting ranges on corners
#' 
#' This function makes pie charts for plotting ranges on corners.  It makes use of 
#' \code{ape:::floating.pie.asp} (copied to ape_floating_pie_asp in 
#' BioGeoBEARS) to plot the pie charts on the corners.
#' 
#' To get the corner coordinates, use \code{\link{corner_coords}}.  Please note the 
#' special input required in that function to get it to access a corner-coordinates 
#' function in the extensions data (\code{extdata}) directory.
#' 
#' @param pievals The matrix (numnodes x numstates) of probabilities to plot.
#' @param coords The coordinates at which to plot the labels.
#' @param piecol The color for each possible state.
#' @param adj Position adjustment; default \code{adj=c(0.5,0.5)}
#' @param ... Additional arguments to standard functions
#' @return nothing
#' @export
#' @seealso \code{\link{cornerlabels}}, \code{\link{corner_coords}}, \code{\link{get_lagrange_nodenums}}, 
#' \code{\link{LGpy_splits_fn_to_table}}, \code{\link{LGcpp_splits_fn_to_table}}
#' @note Go (BioGeo)BEARS!
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu} 
#' @references
#' \url{http://phylo.wikidot.com/matzke-2013-international-biogeography-society-poster}
#' \url{https://code.google.com/p/lagrange/}
#' @examples
#' test=1
#' 
cornerpies <- function(pievals, coords, piecol, adj=c(0.5,0.5), ...)
	{
	#require(ape)	# for ape:::floating.pie.asp
	
	# Make sure pievals, i.e. the ancestral state probabilities,
	# are a numeric matrix instead of a data.frame
	pievals = unname(as.matrix(pievals))
	
	args <- list(...)
    CEX <- if ("cex" %in% names(args)) 
    	{
	    args$cex
	    } else {
	    par("cex")
	    }
	   
# 	# Draw the rectangles	   
# 	width <- strwidth(text, units = "inches", cex = CEX)
# 	height <- strheight(text, units = "inches", cex = CEX)
# 
# 	width <- xinch(width)
# 	height <- yinch(height)

	XX = coords$x
	YY = coords$y	
# 	xl <- XX - width * adj[1] - xinch(0.03)
# 	xr <- xl + width + xinch(0.03)
# 	yb <- YY - height * adj[2] - yinch(0.02)
# 	yt <- yb + height + yinch(0.05)
# 	rect(xl, yb, xr, yt, col = bg)
# 	
# 	# Write the text
# 	text(XX, YY, text, adj = adj, col = col, ...)
	
	if (is.vector(pie))
		{
		pie <- cbind(pie, 1 - pie)
		}
	xrad <- CEX * diff(par("usr")[1:2])/50
	xrad <- rep(xrad, nrow(pievals))
	XX <- XX + adj[1] - 0.5
	YY <- YY + adj[2] - 0.5
	
	# Loop through the nodes
	for (i in 1:nrow(pievals))
		{
		if (any(is.na(pievals[i, ]))) 
			{
			next()
			}
		ape_floating_pie_asp(XX[i], YY[i], pievals[i, ], radius = xrad[i], col = piecol)
		}
	
	
	return()
	}





#######################################################
# add_corners
#######################################################
#' Iterate up through a plotted tree, getting the coordinates of the corners 
#' 
#' What it says.
#'
#' @param startnode The node to start at (this is a recursive function)
#' @param tr A tree object in \code{\link[ape]{phylo}} format.
#' @param nodecoords The accumulating list of node coordinates
#' @param corners_list The accumulating list of corners
#' @return \code{corners_list} 
#' @export
#' @seealso \code{\link[ape]{phylo}}, \code{\link{get_nodenums}}
#' @note Go (BioGeo)BEARS!
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu} 
#' @references
#' \url{http://phylo.wikidot.com/matzke-2013-international-biogeography-society-poster}
#' @examples
#' blah=1
#' 
add_corners <- function(startnode, tr, nodecoords, corners_list)
	{
	# Get the daughters of this node 
	daughtersTF = tr$edge[,1] == startnode
	
	# Run through them if they exist
	if (sum(daughtersTF) > 0)
		{
		# Get the daughter nodenums
		daughters = tr$edge[daughtersTF,2]

		# Get and store node number
		nums = 1:nrow(corners_list$nodevals)
		filled_TF = !is.na(corners_list$nodevals)
		
		if (sum(filled_TF) > 0)
			{
			row_we_are_at = 1 + max(nums[filled_TF])
			} else {
			row_we_are_at = 1
			}
		
		# Store the location
		corners_list$nodevals[row_we_are_at, 1] = startnode
		
		# Get left corner
		# x coordinate
		corners_list$corner_coords_left[row_we_are_at,1] = nodecoords[startnode,1]
		# y coordinate
		corners_list$corner_coords_left[row_we_are_at,2] = nodecoords[daughters[2],2]

		# Get right corner
		# x coordinate
		corners_list$corner_coords_right[row_we_are_at,1] = nodecoords[startnode,1]
		# y coordinate
		corners_list$corner_coords_right[row_we_are_at,2] = nodecoords[daughters[1],2]
		
		# Then propagate up
		for (d in daughters)
			{
			#startnode = d
			corners_list = add_corners(startnode=d, tr, nodecoords, corners_list)
			}
		
		return(corners_list)
		} else {
		return(corners_list)		
		}
	return(corners_list)
	}

# Get the corner coordinates
#######################################################
# corner_coords
#######################################################
#' Get the corner coordinates
#' 
#' Gets the coordinates of the corners when the tree is plotted.
#'
#' Because this function needs to use a modified version of the APE plot.phylo
#' function, and for complex reasons APE's .C functions cannot be used 
#' elsewhere without causing problems with R CMD check, this function is
#' left up to user specification.  Basically, the user puts in
#' the name of the function, which is available in the extension data
#' (\code{extdata/a_scripts}) directory of the package.  The defaults work on the 
#' developer's machine, other users may have to e.g. change "manual" to \code{tmplocation},
#' where \code{tmplocation} is specified as in the example.
#'
#' @param tr A tree object in \code{\link[ape]{phylo}} format.
#' @param coords_fun The name of the function to use to get node coordinates. Default: 
#' "plot_phylo3_nodecoords". 
#' @param tmplocation The directory location containing the R 
#' script \code{plot_phylo3_nodecoords.R}. This function, modified from the \code{\link[ape]{ape}} function
#' \code{\link[ape]{plot.phylo}}, cannot be included directly in the R package as it contains C code that 
#' does not pass CRAN's R CMD check. Default is "auto". Option "manual" will throw an error 
#' check unless your path structure matches the developer's.  Most users
#' should probably use the \code{\link[base]{system.file}} command in the examples, below. 
#' Using cornercoords_loc="manual", will not allow split states to be plot. The R script \code{plot_phylo3_nodecoords.R} 
#' is located in the BioGeoBEARS extension data 
#' directory, \code{extdata/a_scripts}.  You should be able to get the full path with 
#' \code{list.files(system.file("extdata/a_scripts", package="BioGeoBEARS"), full.names=TRUE)}.
#' @param root.edge Plot the root.edge, if it exists? Default TRUE.
#' @return \code{corners_list} 
#' @export
#' @seealso \code{\link[ape]{phylo}}, \code{\link{get_nodenums}}
#' @note Go (BioGeo)BEARS!
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu} 
#' @references
#' \url{http://phylo.wikidot.com/matzke-2013-international-biogeography-society-poster}
#' @examples
#' 
#' # Set location like this if you don't have plot_phylo3_nodecoords
#' # hardcoded/sourced elsehwhere
#' # tmplocation = np(system.file("extdata/a_scripts", package="BioGeoBEARS"))
#' # 
#' \dontrun{
#' extdata_dir = np(system.file("extdata", package="BioGeoBEARS"))
#' trfn = np(paste(extdata_dir, "/Psychotria_5.2.newick", sep=""))
#' tr = read.tree(trfn)
#' tmplocation = np(system.file("extdata/a_scripts", package="BioGeoBEARS"))
#' corner_coords(tr, coords_fun="plot_phylo3_nodecoords", tmplocation=tmplocation)
#' }
#' 
#'
#' 
corner_coords <- function(tr, coords_fun="plot_phylo3_nodecoords", tmplocation="auto", root.edge=TRUE)
	{
	defaults='
	coords_fun="plot_phylo3_nodecoords"
	tmplocation="manual"
	'
	
	# The plot_phylo3_nodecoords causes problems for CRAN, since APE's .C functions
	# aren't exportable, or something.
	# So, we'll source this script from extdata

	# coords_fun="plot_phylo3_nodecoords"; location="manual"

	# trfn = "/Dropbox/_njm/__packages/BioGeoBEARS_setup/
	# inst/extdata/examples/Psychotria_M0/LGcpp/Psychotria_5.2.newick"
	
	# tr = read.tree(trfn)
	
	# 2019-08-01_NJM fix for some scripts that defaulted to "manual"
	if (tmplocation == "manual")
		{
		scriptdir = "/Dropbox/_njm/__packages/BioGeoBEARS_setup/inst/extdata/a_scripts/"
		} else if (tmplocation == "auto") {
		extdata_dir = np(system.file("extdata", package="BioGeoBEARS"))
		scriptdir = slashslash(paste(extdata_dir, "a_scripts" , sep="/"))
		} else {
		scriptdir = tmplocation
		}


	file_to_source = np(paste(addslash(scriptdir), coords_fun, ".R", sep=""))
	source(file=file_to_source)
	
	# Get the coordinates of the vertices in the tree
	
	# Initialize
	trcoords = matrix(data=c(1,1), nrow=1, ncol=2)
	
	# Set up the command as a string
	# trcoords = plot_phylo3_nodecoords(tr, plot=FALSE)

	# 2017-11-02_edit for new versino of APE
	if (packageVersion("ape") < 5.0)
		{
		# Set up the command as a string
		# trcoords = plot_phylo3_nodecoords(tr, plot=FALSE)
		cmdstr = paste("trcoords = ", coords_fun, "(tr, plot=FALSE, root.edge=root.edge)", sep="")
		eval(parse(text=cmdstr))
	
		# X and Y coords for nodes, 1-37
		nodecoords = cbind(trcoords$xx, trcoords$yy)
		} else {
		# Temporary plot, not to screen (hopefully)
		trcoords = plot_phylo3_nodecoords_APE5(tr, plot=FALSE, root.edge=root.edge)
		
		# X and Y coords for nodes, 1-37
		nodecoords = cbind(trcoords$xx, trcoords$yy)
		} # END if (packageVersion("ape") < 5.0)
	
	# Go through the edge matrix from the root, take the x coord of the node,
	# and the y coord of the descendant
	rootnode = get_nodenum_structural_root(tr)
	
	# Get corner coords
	corners_list = NULL
	corners_list$corner_coords_left = matrix(NA, ncol=2, nrow=tr$Nnode)
	corners_list$corner_coords_right = matrix(NA, ncol=2, nrow=tr$Nnode)
	corners_list$nodevals = matrix(NA, ncol=1, nrow=tr$Nnode)
	
	corners_list = add_corners(startnode=rootnode, tr, nodecoords, corners_list)
	
	
	left = corners_list$corner_coords_left
	right = corners_list$corner_coords_right
	node = corners_list$nodevals
	
	leftcorns = adf(cbind(node, left))
	names(leftcorns) = c("node", "x", "y")
	row.names(leftcorns) = node

	rightcorns = adf(cbind(node, right))
	names(rightcorns) = c("node", "x", "y")
	row.names(rightcorns) = node
	
	corners_list = NULL
	corners_list$leftcorns = leftcorns
	corners_list$rightcorns = rightcorns
	
	return(corners_list)
	}



# Get the node coordinates
#######################################################
# node_coords
#######################################################
#' Get the node coordinates
#' 
#' Gets the coordinates of the nodes when the tree is plotted.
#'
#' Because this function needs to use a modified version of the APE plot.phylo
#' function, and for complex reasons APE's .C functions cannot be used 
#' elsewhere without causing problems with R CMD check, this function is
#' left up to user specification.  Basically, the user puts in
#' the name of the function, which is available in the extension data
#' (\code{extdata/a_scripts}) directory of the package.  The defaults work on the 
#' developer's machine, other users may have to e.g. change "manual" to \code{tmplocation},
#' where \code{tmplocation} is specified as in the example.
#'
#' @param tr A tree object in \code{\link[ape]{phylo}} format.
#' @param coords_fun The name of the function to use to get node coordinates. Default: 
#' "plot_phylo3_nodecoords". 
#' @param tmplocation The directory location containing the R 
#' script \code{plot_phylo3_nodecoords.R}. This function, modified from the \code{\link[ape]{ape}} function
#' \code{\link[ape]{plot.phylo}}, cannot be included directly in the R package as it contains C code that 
#' does not pass CRAN's R CMD check. Default is "auto". Option "manual" will throw an error 
#' check unless your path structure matches the developer's.  Most users
#' should probably use the \code{\link[base]{system.file}} command in the examples, below. 
#' Using cornercoords_loc="manual", will not allow split states to be plot. The R script \code{plot_phylo3_nodecoords.R} 
#' is located in the BioGeoBEARS extension data 
#' directory, \code{extdata/a_scripts}.  You should be able to get the full path with 
#' \code{list.files(system.file("extdata/a_scripts", package="BioGeoBEARS"), full.names=TRUE)}.
#' @param root.edge Plot the root.edge, if it exists? Default TRUE.
#' @return \code{nodecoords}, a data.frame with nodecoords$x and nodecoords$y specifying
#' the plot coordinates of each node.
#' @export
#' @seealso \code{\link[ape]{phylo}}, \code{\link{get_nodenums}}
#' @note Go (BioGeo)BEARS!
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu} 
#' @references
#' \url{http://phylo.wikidot.com/matzke-2013-international-biogeography-society-poster}
#' @examples
#' 
#' # Set location like this if you don't have plot_phylo3_nodecoords
#' # hardcoded/sourced elsehwhere
#' # tmplocation = np(system.file("extdata/a_scripts", package="BioGeoBEARS"))
#' # 
#' \dontrun{
#' extdata_dir = np(system.file("extdata", package="BioGeoBEARS"))
#' trfn = np(paste(extdata_dir, "/Psychotria_5.2.newick", sep=""))
#' tr = read.tree(trfn)
#' tmplocation = np(system.file("extdata/a_scripts", package="BioGeoBEARS"))
#' nodes_xy = node_coords(tr, coords_fun="plot_phylo3_nodecoords", tmplocation=tmplocation)
#' nodes_xy
#' }
#' 
#'
#' 
node_coords <- function(tr, coords_fun="plot_phylo3_nodecoords", tmplocation="auto", root.edge=TRUE)
	{
	defaults='
	coords_fun="plot_phylo3_nodecoords"
	tmplocation="manual"
	'
	
	# The plot_phylo3_nodecoords causes problems for CRAN, since APE's .C functions
	# aren't exportable, or something.
	# So, we'll source this script from extdata

	# coords_fun="plot_phylo3_nodecoords"; location="manual"

	# trfn = "/Dropbox/_njm/__packages/BioGeoBEARS_setup/
	# inst/extdata/examples/Psychotria_M0/LGcpp/Psychotria_5.2.newick"
	
	# tr = read.tree(trfn)

	# 2019-08-01_NJM fix for some scripts that defaulted to "manual"
	if (tmplocation == "manual")
		{
		scriptdir = "/Dropbox/_njm/__packages/BioGeoBEARS_setup/inst/extdata/a_scripts/"
		} else if (tmplocation == "auto") {
		extdata_dir = np(system.file("extdata", package="BioGeoBEARS"))
		scriptdir = slashslash(paste(extdata_dir, "a_scripts" , sep="/"))
		} else {
		scriptdir = tmplocation
		}
	file_to_source = np(paste(addslash(scriptdir), coords_fun, ".R", sep=""))
	source(file=file_to_source)
	
	# Get the coordinates of the vertices in the tree
	
	# Initialize
	trcoords = matrix(data=c(1,1), nrow=1, ncol=2)
	
	# Set up the command as a string
	# trcoords = plot_phylo3_nodecoords(tr, plot=FALSE)

	# 2017-11-02_edit for new versino of APE
	if (packageVersion("ape") < 5.0)
		{
		# Set up the command as a string
		# trcoords = plot_phylo3_nodecoords(tr, plot=FALSE)
		cmdstr = paste("trcoords = ", coords_fun, "(tr, plot=FALSE, root.edge=root.edge)", sep="")
		eval(parse(text=cmdstr))
	
		# X and Y coords for nodes, 1-37
		nodecoords = cbind(trcoords$xx, trcoords$yy)
		} else {
		# Temporary plot, not to screen (hopefully)
		trcoords = plot_phylo3_nodecoords_APE5(tr, plot=FALSE, root.edge=root.edge)
		
		# X and Y coords for nodes, 1-37
		nodecoords = cbind(trcoords$xx, trcoords$yy)
		} # END if (packageVersion("ape") < 5.0)


	nodes_xy = adf(nodecoords)
	names(nodes_xy) = c("x", "y")
	row.names(nodes_xy) = 1:nrow(nodecoords)
	
	return(nodes_xy)
	}






#######################################################
# PLOT PROBABILITIES ON A PER-AREA BASIS
#######################################################

#' Plot per-area occupation probabilities in a barchart at each node
#'
#' Although it is traditional to plot the most-probable ancestral 
#' states at each node, or plot the piechart of all state probabilities
#' at each node, there are other imaginable options. This function plots
#' a barplot at each node, where the probability of each area being
#' occupied has been calculated by summing across all of the state/range
#' probabilities. The bar heights represent the integrated probability of
#' each individual area being occupied. Particularly for complex analyses
#' with high uncertainty, this may be a useful plot.
#'
#' @param tr An ape phylo object.
#' @param res A \code{BioGeoBEARS_results_object} that is the result of a 
#' \code{\link{bears_optim_run}}. (typically \code{res}, 
#' \code{resDEC}, \code{resDECj}, etc.)
#' @param areas The list of area abbreviations, e.g. from 
#'              \code{getareas_from_tipranges_object(tipranges)}
#' @param states_list_0based A list of states, where each
#' list item is a vector of 0-based area numbers.
#' @param titletxt A title for the plot. Default \code{""}.
#' @param cols_each_area The color for each indvidual area. Auto-determined if 
#'                       \code{NULL} (default).
#' @param barwidth_proportion The proportional width of the bars (with 
#'        respect to total plot width, I think). Default 0.02.
#' @param barheight_proportion The proportional width of the bars (with 
#'        respect to total plot height, I think). Default 0.025.
#' @param offset_nodenums Node numbers on which to offset the barplot
#'        (for improved readability).
#' @param offset_xvals A multiplier on barwidth to move the offset_nodenums 
#'        left or right. Each xval corresponds to a value in offset_nodenums.
#' @param offset_yvals A multiplier on barheight to move the offset_nodenums 
#'        up or down. Each yval corresponds to a value in offset_nodenums.
#' @param root.edge An input for \code{node_coords}, which is passed to 
#'        \code{plot_phylo3_nodecoords}. Default \code{TRUE}.
#' @param border The color of the border of the boxes holding areas. Default is
#' the string "default", which means that the borders will be "gray50" if 
#' there is no color specified (that is, cols_each_area=NULL, which means bars 
#' will be "gray70"), and borders will be "black" if color is 
#' specified.  By changing the box borders and the tree color (see 
#' parameters \code{border} or  \code{trcol} in \code{plot_per_area_probs}, 
#' or \code{edge.color} in \code{ape::plot.phylo}).
#' the user can emphasize or de-emphasize one or the other.
#' @param trcol The color of the lines in the phylogeny when plotted.  
#' Default is "black", but "gray60" looks good if you want to 
#' de-emphasize the tree with respect to e.g. node areas.
#' @param plot_rangesizes If \code{FALSE} (default), plot the barplot with each area 
#'        probability. If \code{TRUE}, collapse to the probabilities of each range size and
#'        plot that.
#' @return probs_each_area The probabilities of each area being occupied, at each node.
#' @export
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu}
#' @examples
#' test=1
plot_per_area_probs <- function(tr, res, areas, states_list_0based, titletxt="", cols_each_area=NULL, barwidth_proportion=0.02, barheight_proportion=0.025, offset_nodenums=NULL, offset_xvals=NULL, offset_yvals=NULL, root.edge=TRUE, border="default", trcol="black", plot_rangesizes=FALSE, nodes_to_plot=NULL)
	{
	defaults='
	areas = getareas_from_tipranges_object(tipranges)
	states_list_0based = rcpp_areas_list_to_states_list(areas=areas, maxareas=max_range_size, include_null_range=TRUE)
	states_list_0based

	cols_each_area=NULL
	offset_nodenums=NULL
	offset_xvals=NULL
	offset_yvals=NULL
	barwidth_proportion=0.02
	barheight_proportion=0.025
	
	border="default"
	trcol="black"
	plot_rangesizes=FALSE
	nodes_to_plot=NULL
	' # END defaults
	
	
	# Error check
	if ((length(cols_each_area) == 1) && (cols_each_area == FALSE))
		{
		errortxt = "\n\nERROR IN plot_per_area_probs():\ncols_each_area=FALSE, but it should be either:\nNULL (which gives grey boxes)\nTRUE (which makes colors auto-generated), or \na list of colors the same length as 'areas'.\n\n"
		cat(errortxt)
		
		stop("\n\nStopping on error.\n\n")
		
		} # END if ((length(cols_each_area) == 1) && (cols_each_area == FALSE))
	
	
	# Get the relative probabilities of each state/range
	relprobs_matrix = res$ML_marginal_prob_each_state_at_branch_top_AT_node
	
	# Collapse to the probabilities of each area
	if (plot_rangesizes == FALSE)
		{
		probs_each_area = infprobs_to_probs_of_each_area(relprobs_matrix, states_list=states_list_0based)
		}
	
	# Collapse to the probabilities of each range size
	if (plot_rangesizes == TRUE)
		{
		rangesizes_by_state = sapply(FUN=length, X=states_list_0based)
		rangesizes = sort(unique(rangesizes_by_state))
		rangesizes
		
		probs_each_area = infprobs_to_probs_of_each_rangesize(relprobs_matrix, states_list=states_list_0based)
		}
	dim(probs_each_area)


	# To get offset_tiplabels:
	ntips = length(tr$tip.label)
	numnodes = tr$Nnode
	tipnums = 1:ntips
	nodenums = (ntips+1):(ntips+numnodes)

	extdata_dir = np(system.file("extdata", package="BioGeoBEARS"))
	tmplocation=slashslash(paste(extdata_dir, "a_scripts" , sep="/"))
	nodes_xy = node_coords(tr, root.edge=root.edge, tmplocation=tmplocation)
	nodes_xy

	# Make bar width a proportion of the width of the plot in x
	#barwidth_proportion = 0.02
	barwidth = (max(nodes_xy$x) - min(nodes_xy$x)) * barwidth_proportion
	barwidth

	#barheight_proportion = 0.025
	barheight = (max(nodes_xy$y) - min(nodes_xy$y)) * barheight_proportion
	barheight

	numareas = ncol(probs_each_area)
	areanums = 1:numareas
	middle = median(areanums)
	middle
	offsets_nodes = (areanums - middle)*barwidth
	offsets_tips = (areanums - 0)*barwidth
	offset_tiplabels = max(offsets_tips) + barwidth/1



	# Plot the tree
	plot(tr, label.offset=offset_tiplabels, root.edge=root.edge, edge.color=trcol)
	# plot(tr, label.offset=offset_tiplabels)
	axisPhylo()
	title(titletxt)

	# Add the areas boxes
	add_per_area_probs_to_nodes(tr, probs_each_area, cols_each_area=cols_each_area, barwidth_proportion=barwidth_proportion, barheight_proportion=barheight_proportion, offset_nodenums=offset_nodenums, offset_xvals=offset_xvals, offset_yvals=offset_yvals, border=border, nodes_to_plot=nodes_to_plot)
	
	return(probs_each_area)
	} # END plot_per_area_probs







#######################################################
# add_per_area_probs_to_nodes
#######################################################

#' Plot per-area occupation probabilities to nodes.
#'
#' Used by \code{\link{plot_per_area_probs}} to plot barcharts of
#' per-area probabilities to each node.
#'
#' @param tr An ape phylo object.
#' @param probs_each_area The probabilities of each area being occupied, at each node.
#' @param cols_each_area The color for each indvidual area. Auto-determined if 
#'                       \code{NULL} (default).
#' @param barwidth_proportion The proportional width of the bars (with 
#'        respect to total plot width, I think). Default 0.02.
#' @param barheight_proportion The proportional width of the bars (with 
#'        respect to total plot height, I think). Default 0.025.
#' @param offset_nodenums Node numbers on which to offset the barplot
#'        (for improved readability).
#' @param offset_xvals A multiplier on barwidth to move the offset_nodenums 
#'        left or right. Each xval corresponds to a value in offset_nodenums.
#' @param offset_yvals A multiplier on barheight to move the offset_nodenums 
#'        up or down. Each yval corresponds to a value in offset_nodenums.
#' @param root.edge An input for \code{node_coords}, which is passed to 
#'        \code{plot_phylo3_nodecoords}. Default \code{TRUE}.
#' @param border The color of the border of the boxes holding areas. Default is
#' the string "default", which means that the borders will be "gray50" if 
#' there is no color specified (that is, cols_each_area=NULL, which means bars 
#' will be "gray70"), and borders will be "black" if color is 
#' specified.  By changing the box borders and the tree color (see 
#' parameter \code{plot_per_area_probs} in \code{plot_per_area_probs}, 
#' or \code{edge.color} in \code{ape::plot.phylo}).
#' the user can emphasize or de-emphasize one or the other.
#' @return \code{NULL}
#' @export
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu}
#' @examples
#' test=1
#'
add_per_area_probs_to_nodes <- function(tr, probs_each_area, cols_each_area=NULL, barwidth_proportion=0.02, barheight_proportion=0.025, offset_nodenums=NULL, offset_xvals=NULL, offset_yvals=NULL, root.edge=TRUE, border="default", nodes_to_plot=NULL)
	{
	defaults='
	cols_each_area=NULL
	offset_nodenums=NULL
	offset_xvals=NULL
	offset_yvals=NULL
	barwidth_proportion=0.02
	barheight_proportion=0.025
	' # END defaults


	# Error check
	if ((length(cols_each_area) == 1) && (cols_each_area == FALSE))
		{
		errortxt = "\n\nERROR IN add_per_area_probs_to_nodes():\n\ncols_each_area=FALSE, but it should be either:\nNULL (which gives grey boxes)\nTRUE (which makes colors auto-generated), or \na list of colors the same length as 'areas'.\n\n"
		cat(errortxt)
		
		stop("\n\nStopping on error.\n\n")
		
		} # END if ((length(cols_each_area) == 1) && (cols_each_area == FALSE))
	

	ntips = length(tr$tip.label)
	numnodes = tr$Nnode
	tipnums = 1:ntips
	nodenums = (ntips+1):(ntips+numnodes)

	extdata_dir = np(system.file("extdata", package="BioGeoBEARS"))
	tmplocation=slashslash(paste(extdata_dir, "a_scripts" , sep="/"))
	nodes_xy = node_coords(tr, root.edge=root.edge, tmplocation=tmplocation)
	nodes_xy


	# Make bar width a proportion of the width of the plot in x
	#barwidth_proportion = 0.02
	barwidth = (max(nodes_xy$x) - min(nodes_xy$x)) * barwidth_proportion
	barwidth

	#barheight_proportion = 0.025
	barheight = (max(nodes_xy$y) - min(nodes_xy$y)) * barheight_proportion
	barheight

	numareas = ncol(probs_each_area)
	areanums = 1:numareas
	middle = median(areanums)
	middle
	offsets_nodes = (areanums - middle)*barwidth
	offsets_tips = (areanums - 0)*barwidth
	offset_tiplabels = max(offsets_tips) + barwidth/1
	


	# If *specific* nodes (internal or tips) are desired, edit:
	if (is.null(nodes_to_plot) == FALSE)
		{
		# tip nodes
		TF = nodes_to_plot <= max(tipnums)
		tipnums_to_use = tipnums[TF]
		
		# internal nodes
		TF = nodes_to_plot > max(tipnums)
		nodenums_to_use = nodenums[TF]
		
		all_nodenums_to_use = c(tipnums_to_use, nodenums_to_use)
		#nodes_xy
		}

	
	
	# Draw the empty boxes

	# xcoords for tips
	xlefts_tips = sapply(X=nodes_xy$x[tipnums], FUN="+", (offsets_tips - barwidth/2))
	xlefts_nodes = sapply(X=nodes_xy$x[nodenums], FUN="+", (offsets_nodes - barwidth/2))
	xrights_tips = sapply(X=nodes_xy$x[tipnums], FUN="+", (offsets_tips + barwidth/2))
	xrights_nodes = sapply(X=nodes_xy$x[nodenums], FUN="+", (offsets_nodes + barwidth/2))
	
	xlefts_matrix = t(cbind(xlefts_tips, xlefts_nodes))
	xrights_matrix = t(cbind(xrights_tips, xrights_nodes))

	ybottoms_per_node = sapply(X=nodes_xy$y, FUN="-", (barheight/2))
	ytops_per_node = sapply(X=nodes_xy$y, FUN="+", (barheight/2))
	
	# Manually modify some positions
	if (is.null(offset_nodenums) == FALSE)
		{
		xlefts_matrix[offset_nodenums,] = xlefts_matrix[offset_nodenums,] + (offset_xvals*barwidth)
		xrights_matrix[offset_nodenums,] = xrights_matrix[offset_nodenums,] + (offset_xvals*barwidth)
		ybottoms_per_node[offset_nodenums] = ybottoms_per_node[offset_nodenums] + (offset_yvals*barheight)
		ytops_per_node[offset_nodenums] = ytops_per_node[offset_nodenums] + (offset_yvals*barheight)
		}
	
	# Convert these into plain lists
	xlefts = c(xlefts_matrix)
	xrights = c(xrights_matrix)

	ybottoms = rep(ybottoms_per_node, times=numareas)
	ytops = rep(ytops_per_node, times=numareas)


	# Plot the box outlines
	#nulls = mapply(FUN=rect, xleft=xlefts, ybottom=ybottoms, xright=xrights, ytop=ytops, MoreArgs=list(col="white", border=border))

	# Then draw black boxes inside these
	# You just have to adjust the top of the black bar, based on prob
	yranges_per_node = ytops_per_node - ybottoms_per_node
	yadd_above_ybottom = yranges_per_node * probs_each_area
	ytops_probs_node = yadd_above_ybottom + ybottoms
	ytops_probs = c(ytops_probs_node)
	

	# IF cols_each_area == TRUE, auto-generate colors
	if (length(cols_each_area) == 1 && (cols_each_area == TRUE))
		{
		tmp_colors_matrix = get_colors_for_numareas(numareas, use_rainbow=FALSE)
		#cols_each_area = c("blue", "green", "yellow", "red")
		cols_each_area = mapply(FUN=rgb, red=tmp_colors_matrix[1,], green=tmp_colors_matrix[2,], blue=tmp_colors_matrix[3,], MoreArgs=list(maxColorValue=255))

		}
	

	# If *specific* nodes (internal or tips) are desired, edit:
	if (is.null(nodes_to_plot) == FALSE)
		{
		# Convert these into plain lists
		xlefts = c(xlefts_matrix[all_nodenums_to_use,])
		xrights = c(xrights_matrix[all_nodenums_to_use,])

		ybottoms = rep(ybottoms_per_node[all_nodenums_to_use], times=numareas)
		ytops = rep(ytops_per_node[all_nodenums_to_use], times=numareas)


		ytops_probs = ytops_probs_node[all_nodenums_to_use,]
		}
	
	
	# Default color is darkgray
	if (is.null(cols_each_area))
		{
		# Auto-generate border, also -- gray50 looks black against lighter gray70
		if (border == "default")
			{
			border = "gray50"
			}

		# Plot the box outlines
		nulls = mapply(FUN=rect, xleft=xlefts, ybottom=ybottoms, xright=xrights, ytop=ytops, MoreArgs=list(col="white", border=border))

		# Draw bars
		nulls = mapply(FUN=rect, xleft=xlefts, ybottom=ybottoms, xright=xrights, ytop=ytops_probs, MoreArgs=list(col="gray70", border=border))
		} else {
		
		
		if ( length(cols_each_area) != length(areas))
			{
			errortxt = paste("\n\nERROR in add_per_area_probs_to_nodes():\n\nif cols_each_area is specified, length(cols_each_area) must equal length(areas).\n\n", sep="")
			
			cat(errortxt)
			
			cat("Your 'areas':\n\n", sep="")
			print(areas)
			
			cat("Your 'cols_each_area':\n\n", sep="")
			print(cols_each_area)
			
			stop("Stopping on error.")
			} # END if ( length(cols_each_area) != length(areas))
		
		# Otherwise, expand colors to each box
		cols_each_area_matrix = matrix(data=cols_each_area, nrow=(ntips+numnodes), ncol=length(cols_each_area), byrow=TRUE)

		# If *specific* nodes (internal or tips) are desired, edit:
		if (is.null(nodes_to_plot) == FALSE)
			{
			cols_each_area_matrix = cols_each_area_matrix[all_nodenums_to_use,]
			}
		
		
		# Plot with different internal box colors
		cols_each_area = c(cols_each_area_matrix)

		# Auto-generate border with colored boxes (black looks best), also
		if (border == "default")
			{
			border = "black"
			}
		
		# Plot the box outlines
		nulls = mapply(FUN=rect, xleft=xlefts, ybottom=ybottoms, xright=xrights, ytop=ytops, MoreArgs=list(col="white", border=border))
		
		# Plot the colored bars
		nulls = mapply(FUN=rect, xleft=xlefts, ybottom=ybottoms, xright=xrights, ytop=ytops_probs, col=cols_each_area, MoreArgs=list(border=border))
		} # END if (is.null(cols_each_area))

	return(NULL)
	}





#######################################################
# add_per_area_probs_to_corners
#######################################################
#' Plot per-area occupation probabilities to corners.
#'
#' Used by \code{\link{plot_per_area_probs}} to plot barcharts of
#' per-area probabilities to each left or right corner.
#'
#' @param tr An ape phylo object.
#' @param probs_each_area The probabilities of each area being occupied, at each node.
#' @param left_or_right Are the "left" or "right" corners to be plotted (run twice to
#'                      do both).
#' @param cols_each_area The color for each indvidual area. Auto-determined if 
#'                       \code{NULL} (default).
#' @param barwidth_proportion The proportional width of the bars (with 
#'        respect to total plot width, I think). Default 0.02.
#' @param barheight_proportion The proportional width of the bars (with 
#'        respect to total plot height, I think). Default 0.025.
#' @param offset_nodenums Node numbers on which to offset the barplot
#'        (for improved readability).
#' @param offset_xvals A multiplier on barwidth to move the offset_nodenums 
#'        left or right. Each xval corresponds to a value in offset_nodenums.
#' @param offset_yvals A multiplier on barheight to move the offset_nodenums 
#'        up or down. Each yval corresponds to a value in offset_nodenums.
#' @param root.edge An input for \code{node_coords}, which is passed to 
#'        \code{plot_phylo3_nodecoords}. Default \code{TRUE}.
#' @param border The color of the border of the boxes holding areas. Default is
#' the string "default", which means that the borders will be "gray50" if 
#' there is no color specified (that is, cols_each_area=NULL, which means bars 
#' will be "gray70"), and borders will be "black" if color is 
#' specified.  By changing the box borders and the tree color (see 
#' parameter \code{plot_per_area_probs} in \code{plot_per_area_probs}, 
#' or \code{edge.color} in \code{ape::plot.phylo}).
#' the user can emphasize or de-emphasize one or the other.
#' @return \code{NULL}
#' @export
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu}
#' @examples
#' test=1
#'
add_per_area_probs_to_corners <- function(tr, probs_each_area, left_or_right, cols_each_area=NULL, barwidth_proportion=0.02, barheight_proportion=0.025, offset_nodenums=NULL, offset_xvals=NULL, offset_yvals=NULL, root.edge=TRUE, border="default")
	{
	defaults='
	cols_each_area=NULL
	offset_nodenums=NULL
	offset_xvals=NULL
	offset_yvals=NULL
	barwidth_proportion=0.02
	barheight_proportion=0.025
	' # END defaults
	
	# Error check
	if ((length(cols_each_area) == 1) && (cols_each_area == FALSE))
		{
		errortxt = "\n\nERROR IN add_per_area_probs_to_corners():\n\ncols_each_area=FALSE, but it should be either:\nNULL (which gives grey boxes)\nTRUE (which makes colors auto-generated), or \na list of colors the same length as 'areas'.\n\n"
		cat(errortxt)
		
		stop("\n\nStopping on error.\n\n")
		
		} # END if ((length(cols_each_area) == 1) && (cols_each_area == FALSE))

	
	ntips = length(tr$tip.label)
	numnodes = tr$Nnode
	tipnums = 1:ntips
	#nodenums = (ntips+1):(ntips+numnodes)
	nodenums = (ntips+1):(ntips+numnodes)
	
	# Get the plot coordinates of the corners ABOVE each node
	extdata_dir = np(system.file("extdata", package="BioGeoBEARS"))
	tmplocation = slashslash(paste(extdata_dir, "a_scripts" , sep="/"))
	corners_list = corner_coords(tr, root.edge=root.edge, tmplocation=tmplocation)
	corners_list
	
	# Get the node numbers of the nodes ABOVE each corner above each node
	leftright_nodes_matrix = get_leftright_nodes_matrix_from_results(tr)
	leftright_nodes_matrix
	
	if (left_or_right == "left")
		{
		corners_xy = corners_list$leftcorns
		nodenums_above_LorR_corner = leftright_nodes_matrix$left
		}
		
	if (left_or_right == "right")
		{
		corners_xy = corners_list$rightcorns
		nodenums_above_LorR_corner = leftright_nodes_matrix$right
		}
	
	
	
	
	# LEFT OR RIGHT SPLITS
	# Probs of each area BELOW the node (nodenum = rownum)
	
	# Make bar width a proportion of the width of the plot in x
	#barwidth_proportion = 0.02
	barwidth = (max(corners_xy$x) - min(corners_xy$x)) * barwidth_proportion
	barwidth

	#barheight_proportion = 0.025
	barheight = (max(corners_xy$y) - min(corners_xy$y)) * barheight_proportion
	barheight

	numareas = ncol(probs_each_area)
	areanums = 1:numareas
	middle = median(areanums)
	middle
	offsets_nodes = (areanums - middle)*barwidth
	offsets_tips = (areanums - 0)*barwidth
	offset_tiplabels = max(offsets_tips) + barwidth/1

	# Draw the empty boxes

	# xcoords for tips
	nodenums_above_LorR_corner
	
	# We just need to plot at the corners above internal nodes, not tips
	#xlefts_tips = sapply(X=corners_xy$x[nodenums], FUN="+", (offsets_tips - barwidth/2))
	rownums = nodenums - ntips
	xlefts_nodes = sapply(X=corners_xy$x[rownums], FUN="+", (offsets_nodes - barwidth/2))
	#xrights_tips = sapply(X=corners_xy$x[tipnums], FUN="+", (offsets_tips + barwidth/2))
	xrights_nodes = sapply(X=corners_xy$x[rownums], FUN="+", (offsets_nodes + barwidth/2))
	
	xlefts_matrix = t(xlefts_nodes)
	xrights_matrix = t(xrights_nodes)

	ybottoms_per_node = sapply(X=corners_xy$y, FUN="-", (barheight/2))
	ytops_per_node = sapply(X=corners_xy$y, FUN="+", (barheight/2))
	
	# Manually modify some positions
	if (is.null(offset_nodenums) == FALSE)
		{
		rownums = match(x=offset_nodenums, table=nodenums)
		xlefts_matrix[rownums,] = xlefts_matrix[rownums,] + (offset_xvals*barwidth)
		xrights_matrix[rownums,] = xrights_matrix[rownums,] + (offset_xvals*barwidth)
		ybottoms_per_node[rownums] = ybottoms_per_node[rownums] + (offset_yvals*barheight)
		ytops_per_node[rownums] = ytops_per_node[rownums] + (offset_yvals*barheight)
		}
	
	# Convert these into plain lists
	xlefts = c(xlefts_matrix)
	xrights = c(xrights_matrix)

	ybottoms = rep(ybottoms_per_node, times=numareas)
	ytops = rep(ytops_per_node, times=numareas)


	# Plot the box outlines
	#nulls = mapply(FUN=rect, xleft=xlefts, ybottom=ybottoms, xright=xrights, ytop=ytops, MoreArgs=list(col="white", border=border))

	# Then draw black boxes inside these
	# You just have to adjust the top of the black bar, based on prob
	yranges_per_node = ytops_per_node - ybottoms_per_node
	yadd_above_ybottom = yranges_per_node * probs_each_area[nodenums_above_LorR_corner,]
	ytops_probs_node = yadd_above_ybottom + ybottoms
	ytops_probs = c(ytops_probs_node)
	

	# IF cols_each_area == TRUE, auto-generate colors
	if (length(cols_each_area) == 1 && (cols_each_area == TRUE))
		{
		tmp_colors_matrix = get_colors_for_numareas(numareas, use_rainbow=FALSE)
		#cols_each_area = c("blue", "green", "yellow", "red")
		cols_each_area = mapply(FUN=rgb, red=tmp_colors_matrix[1,], green=tmp_colors_matrix[2,], blue=tmp_colors_matrix[3,], MoreArgs=list(maxColorValue=255))
		}
	
	
	# Default color is darkgray
	if (is.null(cols_each_area))
		{
		# Auto-generate border, also -- gray50 looks black against lighter gray70
		if (border == "default")
			{
			border = "gray50"
			}
		
		# Plot the box outlines
		nulls = mapply(FUN=rect, xleft=xlefts, ybottom=ybottoms, xright=xrights, ytop=ytops, MoreArgs=list(col="white", border=border))
		
		# Draw bars
		nulls = mapply(FUN=rect, xleft=xlefts, ybottom=ybottoms, xright=xrights, ytop=ytops_probs, MoreArgs=list(col="gray70", border=border))
		} else {
		
		
		if ( length(cols_each_area) != length(areas))
			{
			errortxt = paste("\n\nERROR in add_per_area_probs_to_corners():\n\nif cols_each_area is specified, length(cols_each_area) must equal length(areas).\n\n", sep="")
			cat(errortxt)
			
			cat("Your 'areas':\n\n", sep="")
			print(areas)
			
			cat("Your 'cols_each_area':\n\n", sep="")
			print(cols_each_area)
			
			stop("Stopping on error.")
			} # END if ( length(cols_each_area) != length(area))
		
		# Otherwise, expand colors to each box
		cols_each_area_matrix = matrix(data=cols_each_area, nrow=numnodes, ncol=length(cols_each_area), byrow=TRUE)
		
		
		# Plot with different internal box colors
		cols_each_area = c(cols_each_area_matrix)

		# Auto-generate border with colored boxes (black looks best), also
		if (border == "default")
			{
			border = "black"
			}

		# Plot the box outlines
		nulls = mapply(FUN=rect, xleft=xlefts, ybottom=ybottoms, xright=xrights, ytop=ytops, MoreArgs=list(col="white", border=border))
		
		nulls = mapply(FUN=rect, xleft=xlefts, ybottom=ybottoms, xright=xrights, ytop=ytops_probs, col=cols_each_area, MoreArgs=list(border=border))
		} # END if (is.null(cols_each_area))

	return(NULL)
	}










#######################################################
# plot_BioGeoBEARS_model
#######################################################
#' Graphical display of your anagenetic and cladogenetic biogeography models
#' 
#' This function produces a graphical summary of the model stored in a \code{BioGeoBEARS_run_object}.
#' This could be either an input model, or the result of the ML parameter search.
#'
#' Understanding of phylogenetic methoods in historical biogeography methods is hampered by the
#' difficulty of displaying the models the computer is using.  This function is one attempt to 
#' improve the situation, by plotting the relative weights of the various parameters.
#'
#' Experimental at the moment.
#'
#' @param obj The input object, either a \code{BioGeoBEARS_run_object} (if so, set 
#' \code{obj_is_run_or_results="run"} or an output object from \code{\link{bears_optim_run}} 
#' (if so, specify \code{obj_is_run_or_results="results"}.
#' @param obj_is_run_or_results Specify \code{"run"} or \code{"results"}, as described above 
#' for parameter \code{obj}.
#' @param plotwhat Default is \code{"init"}, which means plotting the starting model parameters. 
#' \code{"est"} plots the estimated model parameters.
#' @param titletxt Additional text for the title of the plot
#' @param statenames State names to pass to \code{\link{plot_cladogenesis_size_probabilities}}. 
#' If \code{NULL} (default), these are auto-generated assuming all areas up to the maximum number 
#' are allowed.
#' @return nada
#' @export
#' @seealso \code{\link{plot_cladogenesis_size_probabilities}}, \code{\link{define_BioGeoBEARS_run}}, \code{\link{define_BioGeoBEARS_model_object}}
#' @note Go (BioGeo)BEARS!
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu} 
#' @references
#' \url{http://phylo.wikidot.com/matzke-2013-international-biogeography-society-poster}
#' @examples
#' blah=1
#' 
plot_BioGeoBEARS_model <- function(obj, obj_is_run_or_results=NULL, plotwhat="init", titletxt="", statenames=NULL)
	{
	runjunk='
	obj = BioGeoBEARS_run_object
	obj_is_run_or_results = "run"
	plotwhat="init"
	titletxt=""
	'
	
	if (is.null(obj_is_run_or_results) == TRUE)
		{
		stoptxt = "\n\nWith plot_BioGeoBEARS_model(), you must specify whether obj_is_run_or_results='run' or 'results'.\n"
		cat(stoptxt)
		stop(stoptxt)
		}
		
	if (obj_is_run_or_results == "run")
		{
		BioGeoBEARS_run_object = obj
		BioGeoBEARS_model_object = BioGeoBEARS_run_object$BioGeoBEARS_model_object
		} else {
		if (obj_is_run_or_results == "results")
			{
			# Extract the correct inputs and outputs
			BioGeoBEARS_run_object = obj$inputs
			BioGeoBEARS_run_object$BioGeoBEARS_model_object = obj$outputs
			BioGeoBEARS_model_object = BioGeoBEARS_run_object$BioGeoBEARS_model_object
			}		
		else {
			stoptxt = "\n\nWith plot_BioGeoBEARS_model(), you must specify whether obj_is_run_or_results='run' or 'results'.\n"
			cat(stoptxt)
			stop(stoptxt)
			}
		}
	
	# Get areas/states/tips
	# Get geographic ranges at tips
	tipranges = getranges_from_LagrangePHYLIP(lgdata_fn=np(BioGeoBEARS_run_object$geogfn))
	
	# Get the list of geographic areas
	areas = getareas_from_tipranges_object(tipranges)
	areas_list = seq(0, length(areas)-1, 1)		# 0-base indexes

	# Change the names to tipranges@df:
	# this doesn't make sense if areas_list is 0-based indexes
	names(tipranges@df) = areas_list
	
	if (is.na(BioGeoBEARS_run_object$max_range_size))
		{
		if (is.null(BioGeoBEARS_run_object$states_list))
			{
			# Maximum range size is all areas
			max_range_size = length(areas)
			} else {
			# If not NA
			# Get max rangesize from states list
			max_range_size = max(sapply(X=BioGeoBEARS_run_object$states_list, FUN=length), na.rm=TRUE)
			}
		} else {
		# Maximum range size hard-coded
		max_range_size = BioGeoBEARS_run_object$max_range_size
		}
	max_numareas = max_range_size

	#######################################################
	# Check that no tips have larger ranges than you allowed
	#######################################################
	tipranges_df_tmp = tipranges@df
	tipranges_df_tmp[tipranges_df_tmp=="?"] = 0
	TF = (rowSums(dfnums_to_numeric(tipranges_df_tmp))) > max_range_size
	if (sum(TF, na.rm=TRUE) > 0)
		{
		cat("\n\nERROR: Tips with ranges too big:\n", sep="")
		print(dfnums_to_numeric(tipranges@df)[TF, ])
		cat("\n\nCheck your input geography file!\n", sep="")
		txt = paste("ERROR: Some tips (listed above) have range sizes larger than ", max_range_size, sep="")
		stop(txt)
		}
	
	# Take the list of areas, and get list of possible states
	# (the user can manually input states if they like)
	if (is.null(BioGeoBEARS_run_object$states_list))
		{
		states_list = rcpp_areas_list_to_states_list(areas=areas, maxareas=max_range_size, include_null_range=BioGeoBEARS_run_object$include_null_range)
		states_list
		#BioGeoBEARS_run_object$states_list = states_list
		#inputs$states_list = states_list
		} else {
		states_list = BioGeoBEARS_run_object$states_list
		#BioGeoBEARS_run_object$states_list = states_list
		#inputs$states_list = states_list
		}





	
	
	# Get everything
	#BioGeoBEARS_model_object = BioGeoBEARS_run_object$BioGeoBEARS_model_object
	#BioGeoBEARS_model_object = params_into_BioGeoBEARS_model_object(BioGeoBEARS_model_object=BioGeoBEARS_model_object, params=params)
	
	# Update linked parameters
	BioGeoBEARS_model_object = calc_linked_params_BioGeoBEARS_model_object(BioGeoBEARS_model_object)

	# Set the dispersal and extinction rate
	d = BioGeoBEARS_model_object@params_table["d",plotwhat]
	e = BioGeoBEARS_model_object@params_table["e",plotwhat]
	a = BioGeoBEARS_model_object@params_table["a",plotwhat]


	#######################################################
	#######################################################
	# Do branch-length exponentiation if desired
	#######################################################
	#######################################################
	b = BioGeoBEARS_model_object@params_table["b",plotwhat]
	# Modify the edge.lengths
	#phy$edge.length = phy$edge.length ^ b


	#######################################################
	#######################################################
	# Do distance-dependence and dispersal multipliers matrix
	#######################################################
	#######################################################
	# Equal dispersal in all directions (unconstrained)
	# Equal extinction probability for all areas
	areas = areas_list
	
	# If there is a distance matrix, use the first one 
	# (non-stratified analysis, here)
	if ( (is.null(BioGeoBEARS_run_object$list_of_distances_mats) == FALSE))
		{
		distances_mat = BioGeoBEARS_run_object$list_of_distances_mats[[1]]
		} else {
		# Default is all areas effectively equidistant
		distances_mat = matrix(1, nrow=length(areas), ncol=length(areas))
		}

	# Get the exponent on distance, apply to distances matrix
	x = BioGeoBEARS_model_object@params_table["x",plotwhat]
	dispersal_multipliers_matrix = distances_mat ^ x

	# Environmental distances
	if ( (is.null(BioGeoBEARS_run_object$list_of_envdistances_mats) == FALSE))
		{
		envdistances_mat = BioGeoBEARS_run_object$list_of_envdistances_mats[[1]]
		} else {
		# Default is all areas effectively equidistant
		envdistances_mat = matrix(1, nrow=length(areas), ncol=length(areas))
		}

	# Get the exponent on environmental distance, apply to distances matrix
	n = BioGeoBEARS_model_object@params_table["n","est"]
	dispersal_multipliers_matrix = dispersal_multipliers_matrix * envdistances_mat ^ n


	# Apply manual dispersal multipliers, if any
	# If there is a manual dispersal multipliers matrix, use the first one 
	# (non-stratified analysis, here)
	if ( (is.null(BioGeoBEARS_run_object$list_of_dispersal_multipliers_mats) == FALSE))
		{
		manual_dispersal_multipliers_matrix = BioGeoBEARS_run_object$list_of_dispersal_multipliers_mats[[1]]
		} else {
		# Default is all areas effectively equidistant
		manual_dispersal_multipliers_matrix = matrix(1, nrow=length(areas), ncol=length(areas))
		}
	
	# Get the exponent on manual dispersal multipliers
	w = BioGeoBEARS_model_object@params_table["w","est"]

	# Apply element-wise
	dispersal_multipliers_matrix = dispersal_multipliers_matrix * manual_dispersal_multipliers_matrix ^ w

	#######################################################
	# multiply parameter d by dispersal_multipliers_matrix
	#######################################################
	dmat_times_d = dispersal_multipliers_matrix * matrix(d, nrow=length(areas), ncol=length(areas))
	amat = dispersal_multipliers_matrix * matrix(a, nrow=length(areas), ncol=length(areas))




	#######################################################
	#######################################################
	# Do area-dependence and extinction multipliers list
	#######################################################
	#######################################################
	if ( (is.null(BioGeoBEARS_run_object$list_of_area_of_areas) == FALSE))
		{
		area_of_areas = BioGeoBEARS_run_object$list_of_area_of_areas[[1]]
		} else {
		# Default is all areas effectively equidistant
		area_of_areas = rep(1, length(areas))
		}
		
	# Get the exponent on extinction, apply to extinction modifiers	
	u = BioGeoBEARS_model_object@params_table["u",plotwhat]
	extinction_modifier_list = area_of_areas ^ (1 * u)
	
	# Apply to extinction rate
	elist = extinction_modifier_list * rep(e, length(areas))



	
	#######################################################
	# Start a big plot with layout()
	#######################################################
	
	#require(plotrix)
	#par(mfrow=c(1,4))
	#Pdecsize_given_ancsize = cbind(maxent01y, maxent01s, maxent01v, maxent01j)
	#color2D.matplot(x=Pdecsize_given_ancsize, c(1,0), c(1,0), c(1,0), show.values=TRUE, show.legend=TRUE, xlab="Ancestral range size", ylab="P(range size) in smaller descendant", axes=FALSE, nslices=100)

	# Layout so plot has 4 columns, 2 main rows and header/footer, and 2 side columns
	numcols = 6
	plotgroup_header = matrix(data=1, nrow=1, ncol=numcols)
	plotgroup_row2 = matrix(data=c(2,3,4,5,6,7), nrow=1, ncol=numcols)

	spmodel_header = matrix(data=8, nrow=1, ncol=numcols)
	plotgroup_row3 = matrix(data=c(9,10,11,12,13,14), nrow=1, ncol=numcols)
	#plotgroup_footer = matrix(data=15, nrow=1, ncol=numcols)
	plotgroup_footer = matrix(data=c(15,15,16,16,16,16), nrow=1, ncol=numcols)
	plotgroups = rbind(plotgroup_header, plotgroup_row2, spmodel_header, plotgroup_row3, plotgroup_footer)
	plotgroups
	
	layout(mat=plotgroups, widths=c(0.2,1,1,1,1,0.3), heights=c(0.2,1,0.1,1,1))
	
	
	#######################################################
	# Header (cell #1)
	#######################################################
	# set the inside box ("i") ahead of time!!  -- removes the 4% extension
	# this friggin' solution took an hour to find,
	# solution here: http://tolstoy.newcastle.edu.au/R/help/06/08/32529.html
	par(mar=c(0,0,0,0), xaxs = "i", yaxs = "i") 
	plot(x=c(0,1), y=c(0,1), pch=".", col="white", xaxt="n", yaxt="n", xlab="", ylab="", bty="n")
	tmptxt = paste("BioGeoBEARS model: ", titletxt, sep="")
	text(x=0.5, y=0.5, tmptxt, cex=1.5, cex.main=2)


	#################################
	# Row 2
	#################################
	# Cell #2 (leftmost, thin)
	plot(0,0, pch=".", col="white", xaxt="n", yaxt="n", xlab="", ylab="", bty="n")

	# cell #3:
	# Table of free vs. fixed params
	plot(x=c(0,1),y=c(0,1), pch=".", col="white", xaxt="n", yaxt="n", xlab="", ylab="", bty="n")
	val_txt = c("d", "e", "a", "x", "n", "w", "u", "b")
	Param = val_txt
	Type = BioGeoBEARS_model_object@params_table[val_txt, "type"]
	Init = BioGeoBEARS_model_object@params_table[val_txt, "init"]
	Est = BioGeoBEARS_model_object@params_table[val_txt, "est"]
	dtf = adf2(cbind(Param, Type, Init, Est))
	dtf$Init = round(as.numeric(dtf$Init), 3)
	dtf$Est = round(as.numeric(dtf$Est), 3)
	row.names(dtf) = NULL
	dtf

	
	# Make the table into a plot
	#require(plotrix)
	#background_colors = matrix(data="#FFEFDB", nrow=1+nrow(dtf), ncol=ncol(dtf))
	addtable2plot(x=0.1, y=0.80, table=dtf, xjust=0, yjust=0)#, bty="o", box.col="black")#, bg=background_colors) # bg doesn't work
	title("Anagenetic\nparameters", font.main=2, cex.main=1, line=-2)
	

	# cell #4
	# Barplot of anagenetic parameters
	par(mar=c(5,3,3,1), xaxs = "r", yaxs = "r") 
	#plot(0,0, pch=".", col="white", xaxt="n", yaxt="n", xlab="", ylab="", bty="n")
	# plot(0,0, pch=".", col="white", xaxt="n", yaxt="n", xlab="", ylab="", bty="n")
	val_txt = c("d", "e", "a", "x", "n", "w", "u")
	plotvals =  BioGeoBEARS_model_object@params_table[val_txt, plotwhat]
	
	# xaxt="s" means plot (anything other than "n")
	# Get x-axis bar centers
	x_axis_bar_centers = barplot(height=plotvals, width=1, space=NULL, names.arg="", legend.text=FALSE, horiz=FALSE, col="white", border="black", xpd=TRUE, xaxt="s", yaxt="s", tck=0, plot=FALSE) 

	# Now plot for realsies
	barplot(height=plotvals, width=1, space=NULL, names.arg="", legend.text=FALSE, horiz=FALSE, col="white", border="black", xpd=TRUE, xaxt="s", yaxt="n", tck=0, yaxp=c(0,1,1)) 
	axis(side=1, at=x_axis_bar_centers, labels=val_txt, tick=FALSE, line=-1)
	title("Anagenetic\nparameters", font.main=2, cex.main=1, line=1)
	axis(side=2, at=NULL, labels=NULL, tick=TRUE, padj=1, las=0, cex.axis=1, tcl=-0.25, line=0, hadj=0.5)
	mtext(text="param. est.", side=2, line=2.1, padj=0.5, adj=0.5, las=3, cex=0.8)#, adj=0.65)
	

	
	# Cell #5 (table of cladogensis params)	
	# Cells #6 (cladogenesis params barplot), #7 (rightmost thin spacer)
	# Cells #8 (cladogenesis header) and #9-14 (cladogenesis rangesize model)
	# Plot the cladogeneis model P(size|event,ancsize) (requires 6 slots)
	plot_cladogenesis_size_probabilities(BioGeoBEARS_run_object, plotwhat=plotwhat)

	
	
	}







#######################################################
# plot_cladogenesis_size_probabilities
#######################################################
#' Graphical display of P(daughter rangesize) for your input or inferred speciation model
#' 
#' This function produces a graphical summary of the daughter rangesize aspect of the 
#' cladogenesis model stored in a \code{BioGeoBEARS_run_object}. This could be either an 
#' input model, or the result of the ML parameter search.
#'
#' The \code{LAGRANGE} DEC model assumes that at cladogenesis events, one daughter species has a 
#' range size of 1 area, and the other daughter either inherits the full ancestral range 
#' (sympatric-subset speciation), inherits the remainder of the ancestral range (vicariance),
#' or as the same range (sympatric-range copying, which is the only option when the ancestor
#' range is of size 1 area.  
#'
#' BioGeoBEARS enables numerous additional models. To see how these are similar or different from
#' the LAGRANGE DEC cladogenesis model, this function can be used.  E.g., comparison of
#' \code{LAGRANGE} DEC to a \code{DIVA}-like model is instructive: see examples. DIVA disallows
#' sympatric-subset speciation (probability 0 under this model), but allows classic vicariance 
#' (a species with 4 areas splitting into 2 daughters, each occupying 2 areas).  LAGRANGE DEC
#' gives 0 probability to a \code{4->(2,2)} history, allowing only \code{4->(3,1)} or 
#' \code{4->(1,3)} histories.
#'
#' Several additional plots relating to the cladogenesis model are also produced.  Best used via
#' \code{\link{plot_BioGeoBEARS_model}}.
#'
#' Experimental at the moment.
#'
#' @param BioGeoBEARS_run_object The input run object.
#' @param plotwhat Default is "input", which means plotting the starting model.
#' @param statenames State names to pass to \code{\link{plot_cladogenesis_size_probabilities}}. 
#' If \code{NULL} (default), these are auto-generated assuming all areas up to the maximum number 
#' are allowed.
#' @return \code{Nothing} 
#' @export
#' @seealso \code{\link{plot_BioGeoBEARS_model}}, \code{\link{define_BioGeoBEARS_run}}, \code{\link{define_BioGeoBEARS_model_object}}
#' @note Go (BioGeo)BEARS!
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu} 
#' @references
#' \url{http://phylo.wikidot.com/matzke-2013-international-biogeography-society-poster}
#' @examples
#' blah=1
#' 
plot_cladogenesis_size_probabilities <- function(BioGeoBEARS_run_object, plotwhat="est", statenames=NULL)
	{
	# Get the model
	BioGeoBEARS_model_object = BioGeoBEARS_run_object$BioGeoBEARS_model_object
	
	# Get areas/states/tips
	# Get geographic ranges at tips
	tipranges = getranges_from_LagrangePHYLIP(lgdata_fn=np(BioGeoBEARS_run_object$geogfn))
	
	# Get the list of geographic areas
	areas = getareas_from_tipranges_object(tipranges)
	areas_list = seq(0, length(areas)-1, 1)		# 0-base indexes
	areanames = areas
	areanames
	
	
	if (is.na(BioGeoBEARS_run_object$max_range_size))
		{
		if (is.null(BioGeoBEARS_run_object$states_list))
			{
			# Maximum range size is all areas
			max_range_size = length(areas)
			} else {
			# If not NA
			# Get max rangesize from states list
			max_range_size = max(sapply(X=BioGeoBEARS_run_object$states_list, FUN=length), na.rm=TRUE)
			}
		} else {
		# Maximum range size hard-coded
		max_range_size = BioGeoBEARS_run_object$max_range_size
		}
	max_numareas = max_range_size




	# Take the list of areas, and get list of possible states
	# (the user can manually input states if they like)
	if (is.null(BioGeoBEARS_run_object$states_list))
		{
		states_list = rcpp_areas_list_to_states_list(areas=areas, maxareas=max_range_size, include_null_range=BioGeoBEARS_run_object$include_null_range)
		states_list
		#BioGeoBEARS_run_object$states_list = states_list
		#inputs$states_list = states_list
		} else {
		states_list = BioGeoBEARS_run_object$states_list
		#BioGeoBEARS_run_object$states_list = states_list
		#inputs$states_list = states_list
		}







	# If there is a distance matrix, use the first one 
	# (non-stratified analysis, here)
	if ( (is.null(BioGeoBEARS_run_object$list_of_distances_mats) == FALSE))
		{
		distances_mat = BioGeoBEARS_run_object$list_of_distances_mats[[1]]
		} else {
		# Default is all areas effectively equidistant
		distances_mat = matrix(1, nrow=length(areas), ncol=length(areas))
		}

	# Get the exponent on distance, apply to distances matrix
	x = BioGeoBEARS_model_object@params_table["x",plotwhat]
	dispersal_multipliers_matrix = distances_mat ^ x

	# Environmental distances
	if ( (is.null(BioGeoBEARS_run_object$list_of_envdistances_mats) == FALSE))
		{
		envdistances_mat = BioGeoBEARS_run_object$list_of_envdistances_mats[[1]]
		} else {
		# Default is all areas effectively equidistant
		envdistances_mat = matrix(1, nrow=length(areas), ncol=length(areas))
		}

	# Get the exponent on environmental distance, apply to distances matrix
	n = BioGeoBEARS_model_object@params_table["n","est"]
	dispersal_multipliers_matrix = dispersal_multipliers_matrix * envdistances_mat ^ n


	# Apply manual dispersal multipliers, if any
	# If there is a manual dispersal multipliers matrix, use the first one 
	# (non-stratified analysis, here)
	if ( (is.null(BioGeoBEARS_run_object$list_of_dispersal_multipliers_mats) == FALSE))
		{
		manual_dispersal_multipliers_matrix = BioGeoBEARS_run_object$list_of_dispersal_multipliers_mats[[1]]
		} else {
		# Default is all areas effectively equidistant
		manual_dispersal_multipliers_matrix = matrix(1, nrow=length(areas), ncol=length(areas))
		}
	
	# Get the exponent on manual dispersal multipliers
	w = BioGeoBEARS_model_object@params_table["w","est"]

	# Apply element-wise
	dispersal_multipliers_matrix = dispersal_multipliers_matrix * manual_dispersal_multipliers_matrix ^ w


	
	# Set up the instantaneous rate matrix (Q matrix)
	# someday we'll have to put "a" (anagenic range-switching) in here...
	#Qmat = rcpp_states_list_to_DEmat(areas_list=areas_list, states_list=states_list, dmat=dmat_times_d, elist=elist, amat=amat, include_null_range=BioGeoBEARS_run_object$include_null_range, normalize_TF=TRUE, makeCOO_TF=force_sparse)

	#######################################################
	# Cladogenic model
	#######################################################
	j = BioGeoBEARS_model_object@params_table["j",plotwhat]
	ysv = BioGeoBEARS_model_object@params_table["ysv",plotwhat]
	ys = BioGeoBEARS_model_object@params_table["ys",plotwhat]
	v = BioGeoBEARS_model_object@params_table["v",plotwhat]
	y = BioGeoBEARS_model_object@params_table["y",plotwhat]
	s = BioGeoBEARS_model_object@params_table["s",plotwhat]
	sum_SPweights = y + s + j + v


	maxent_constraint_01 = BioGeoBEARS_model_object@params_table["mx01",plotwhat]
	
	# Text version of speciation matrix	
	maxent_constraint_01v = BioGeoBEARS_model_object@params_table["mx01v",plotwhat]
	#spPmat = symbolic_to_relprob_matrix_sp(spmat, cellsplit="\\+", mergesym="*", ys=ys, j=j, v=v, maxent_constraint_01=maxent_constraint_01, maxent_constraint_01v=maxent_constraint_01v, max_numareas=max_numareas)
		
	# Set the parameter controlling the size distribution of 
	# the smaller descendant species
	maxent01s_param = BioGeoBEARS_model_object@params_table["mx01s",plotwhat]
	maxent01v_param = BioGeoBEARS_model_object@params_table["mx01v",plotwhat]
	maxent01j_param = BioGeoBEARS_model_object@params_table["mx01j",plotwhat]
	maxent01y_param = BioGeoBEARS_model_object@params_table["mx01y",plotwhat]


	# Cladogenesis model inputs
	spPmat_inputs = NULL

	# Note that this gets the dispersal multipliers matrix, which is applied to 
	# e.g. the j events, NOT the dmat_times_d above which is d*dispersal_multipliers_matrix
	dmat = dispersal_multipliers_matrix
	spPmat_inputs$dmat = dmat

	states_indices = states_list
	# shorten the states_indices by 1 (cutting the 
	# null range state from the speciation matrix)
	if (include_null_range == TRUE)
		{
		states_indices[1] = NULL
		} # END if (include_null_range == TRUE)
	spPmat_inputs$l = states_indices
	spPmat_inputs$s = s
	spPmat_inputs$v = v
	spPmat_inputs$j = j
	spPmat_inputs$y = y
	spPmat_inputs$maxent01s_param = maxent01s_param
	spPmat_inputs$maxent01v_param = maxent01v_param
	spPmat_inputs$maxent01j_param = maxent01j_param
	spPmat_inputs$maxent01y_param = maxent01y_param

	
	
	
	# Calculate the ancsize/decsize speciation model...
	maxent01s = relative_probabilities_of_subsets(max_numareas=max_numareas, maxent_constraint_01=maxent01s_param, NA_val=0)
	maxent01v = relative_probabilities_of_vicariants(max_numareas=max_numareas, maxent_constraint_01v=maxent01v_param, NA_val=0)
	maxent01j = relative_probabilities_of_subsets(max_numareas=max_numareas, maxent_constraint_01=maxent01j_param, NA_val=0)
	maxent01y = relative_probabilities_of_subsets(max_numareas=max_numareas, maxent_constraint_01=maxent01y_param, NA_val=0)

	# Matrix of probs for each ancsize
	maxprob_as_function_of_ancsize_and_decsize = mapply(FUN=max, maxent01s, maxent01s, maxent01s, maxent01s, MoreArgs=list(na.rm=TRUE))
	maxprob_as_function_of_ancsize_and_decsize = matrix(data=maxprob_as_function_of_ancsize_and_decsize, nrow=nrow(maxent01s), ncol=ncol(maxent01s))
	maxprob_as_function_of_ancsize_and_decsize[maxprob_as_function_of_ancsize_and_decsize > 0] = 1
	maxprob_as_function_of_ancsize_and_decsize[maxprob_as_function_of_ancsize_and_decsize <= 0] = 0

	# Now, go through, and make a list of the max minsize for each decsize
	max_minsize_as_function_of_ancsize = apply(X=maxprob_as_function_of_ancsize_and_decsize, MARGIN=1, FUN=maxsize)
	
	rangesizes = 1:nrow(maxent01s)
	rangesizes_at_y = rev(rangesizes - 0.5)
	rangesizes_at_x = rangesizes - 0.5

	maxent01y = adf2(round(maxent01y,2))
	#rownames(maxent01y) = paste("ancsize=", rangesizes, sep="")
	rownames(maxent01y) = rangesizes
	colnames(maxent01y) = rangesizes
	maxent01s = adf2(round(maxent01s,2))
	#rownames(maxent01s) = paste("ancsize=", rangesizes, sep="")
	rownames(maxent01s) = rangesizes
	colnames(maxent01s) = rangesizes
	maxent01v = adf2(round(maxent01v,2))
	#rownames(maxent01v) = paste("ancsize=", rangesizes, sep="")
	rownames(maxent01v) = rangesizes
	colnames(maxent01v) = rangesizes
	maxent01j = adf2(round(maxent01j,2))
	#rownames(maxent01j) = paste("ancsize=", rangesizes, sep="")
	rownames(maxent01j) = rangesizes
	colnames(maxent01j) = rangesizes

	# Set to 0, if the parameter is 0
	y_cellcolors = NA
	s_cellcolors = NA
	v_cellcolors = NA
	j_cellcolors = NA
	if (y == 0)
		{
		maxent01y[maxent01y != 0] = 0
		y_cellcolors = matrix(data="white", nrow=nrow(maxent01y), ncol=ncol(maxent01y))
		}
	if (s == 0)
		{
		maxent01s[maxent01s != 0] = 0
		s_cellcolors = matrix(data="white", nrow=nrow(maxent01s), ncol=ncol(maxent01s))
		}
	if (v == 0)
		{
		maxent01v[maxent01v != 0] = 0
		v_cellcolors = matrix(data="white", nrow=nrow(maxent01v), ncol=ncol(maxent01v))
		}
	if (j == 0)
		{
		maxent01j[maxent01j != 0] = 0
		j_cellcolors = matrix(data="white", nrow=nrow(maxent01j), ncol=ncol(maxent01j))
		}



	# Cell #5
	# Return to no margins
	par(mar=c(0,0,0,0), xaxs = "i", yaxs = "i") 
	
	# cells # 5, 6
	#plot(0,0, pch=".", col="white", xaxt="n", yaxt="n", xlab="", ylab="", bty="n")
	# Table of free vs. fixed params
	plot(x=c(0,1),y=c(0,1), pch=".", col="white", xaxt="n", yaxt="n", xlab="", ylab="", bty="n")
	val_txt = c("ysv", "ys", "y", "s", "v", "mx01", "mx01y", "mx01s", "mx01v", "mx01j")
	Param = val_txt
	Type = BioGeoBEARS_model_object@params_table[val_txt, "type"]
	Init = BioGeoBEARS_model_object@params_table[val_txt, "init"]
	Est = BioGeoBEARS_model_object@params_table[val_txt, "est"]
	dtf = adf2(cbind(Param, Type, Init, Est))
	dtf = dfnums_to_numeric(dtf)
	row.names(dtf) = NULL
	dtf$Init = round(dtf$Init, 2)
	dtf$Est = round(dtf$Est, 2)
	dtf
	
	# Make the table into a plot
	#require(plotrix)
	plotrix::addtable2plot(x=0.1, y=0.80, table=dtf, xjust=0, yjust=0, cex=1)
	title("Cladogenetic\nparameters", font.main=2, cex.main=1, line=-2)


	# Row #2, right side: speciation model per-event probabilities
	# Rightmost square cell
	#plot(0,0, pch=".", col="white", xaxt="n", yaxt="n", xlab="", ylab="", bty="n")
	par(mar=c(5,3,3,1), xaxs = "r", yaxs = "r") 
	
	
	# If you just want to plot the default (unscaled) y,s,v,j values:
	plot_unscaled_ysvj = FALSE
	if (plot_unscaled_ysvj)
		{
		barplot_yaxis = c(0, 0.5, 1)
		barplot_yaxis_txt = c(0, 0.5, 1)
		ylims = c(0,1.35)
		
		val_txt = c("y", "s", "v", "j")
		plotvals =  BioGeoBEARS_model_object@params_table[val_txt, plotwhat]
		#plotvals = c(y,s,v,j)
		
		# xaxt="s" means plot (anything other than "n")
	
		# Get x-axis bar centers
		x_axis_bar_centers = barplot(height=plotvals, width=1, space=NULL, names.arg="", legend.text=FALSE, horiz=FALSE, col="white", border="black", ylim=ylims, xpd=TRUE, xaxt="s", yaxt="n", tck=0, yaxp=c(0,1,1), plot=FALSE) 

		# Now plot for realsies
		barplot(height=plotvals, width=1, space=NULL, names.arg="", legend.text=FALSE, horiz=FALSE, col="white", border="black", ylim=ylims, xpd=TRUE, xaxt="s", yaxt="n", tck=0, yaxp=c(0,1,1)) 
		axis(side=1, at=x_axis_bar_centers, labels=c("y", "s", "v", "j"), tick=FALSE, line=-1)
		title("Relative prob. of\neach type of\ncladogenesis event", font.main=2, cex.main=1)
		axis(side=2, at=barplot_yaxis, labels=barplot_yaxis_txt, tick=TRUE, padj=0.5, las=1, cex.axis=1, tcl=-0.25, line=0, hadj=0.5)
		mtext(text="relative prob.", side=2, line=2.1, padj=0.5, adj=0.5, las=3, cex=0.8)#, adj=0.65)
		
		} else {
		# If you ARE rescaling, e.g. so ysvj add up to 1
		barplot_yaxis = c(0, 0.5, 1)
		barplot_yaxis_txt = c(0, 0.5, 1)
		ylims = c(0, 1.15)

		params_table = BioGeoBEARS_model_object@params_table
		tmpvals = get_clado_perEvent_weights(params_table, plotwhat=plotwhat)
		
		plotvals = c(tmpvals$y, tmpvals$s, tmpvals$v, tmpvals$j)
		# xaxt="s" means plot (anything other than "n")
	
		# Get x-axis bar centers
		x_axis_bar_centers = barplot(height=plotvals, width=1, space=NULL, names.arg="", legend.text=FALSE, horiz=FALSE, col="white", border="black", ylim=ylims, xpd=TRUE, xaxt="s", yaxt="n", tck=0, yaxp=c(0,1,1), plot=FALSE) 

		# Now plot for realsies
		barplot(height=plotvals, width=1, space=NULL, names.arg="", legend.text=FALSE, horiz=FALSE, col="white", border="black", ylim=ylims, xpd=TRUE, xaxt="s", yaxt="n", tck=0, yaxp=c(0,1,1)) 
		axis(side=1, at=x_axis_bar_centers, labels=c("y", "s", "v", "j"), tick=FALSE, line=-1)
		title("Relative prob. of\neach type of\ncladogenesis event", font.main=2, cex.main=1)
		axis(side=2, at=barplot_yaxis, labels=barplot_yaxis_txt, tick=TRUE, padj=0.5, las=1, cex.axis=1, tcl=-0.25, line=0, hadj=0.5)
		mtext(text="relative prob.", side=2, line=2.1, padj=0.5, adj=0.5, las=3, cex=0.8)#, adj=0.65)
		
		} # end parplot in right square cell of row #2.

	# Plot text of the numbers, above the bars, if desired
	plot_text = TRUE
	if (plot_text == TRUE)
		{
		# pos=3 means above
		text(x=x_axis_bar_centers, y=plotvals, labels=round(plotvals,2), pos=3, offset=0.2, cex=0.9)
		}
	
	
	# Rightmost thin cell
	par(mar=c(0,0,0,0), xaxs = "i", yaxs = "i") 
	plot(0,0, pch=".", col="white", xaxt="n", yaxt="n", xlab="", ylab="", bty="n")


	# Row 3: title for speciation model
	par(mar=c(0,0,0,0), xaxs = "i", yaxs = "i") 
	plot(x=c(0,1), y=c(0,1), pch=".", col="white", xaxt="n", yaxt="n", xlab="", ylab="", bty="n")
	tmptxt = paste("Prob(smaller daughter range size, given ancestor range size)", sep="")
	#text(x=0.5, y=0.5, tmptxt, cex=1.2)
	title(tmptxt, cex=1.5, line=-1)

	# Row 4 left thin column
	plot(x=c(0,1), y=c(0,1), pch=".", col="white", xaxt="n", yaxt="n", xlab="", ylab="", bty="n")
	#mtext(text="ancsize", side=2, line=-3, hadj=0.5, padj=1, las=3, cex=1)
	mtext(text="ancsize", side=2, line=-3, adj=0.65, padj=1, las=3, cex=0.8)
	
	# Plot maxent
	par(mar=c(5,3,3,1), xaxs = "r", yaxs = "r") 
	
	# for color2D.matplot
	#require(plotrix)	# for color2D.matplot
	
	plotrix::color2D.matplot(x=maxent01y, c(1,0), c(1,0), c(1,0), cellcolors=y_cellcolors, show.values=TRUE, show.legend=FALSE, xlab="", ylab="", axes=FALSE, nslices=100, xaxt="n", yaxt="n")
	axis(side=2, at=rangesizes_at_y, labels=rangesizes, tick=FALSE, padj=0.5, las=1, cex.axis=1.5, line=-0.8)
	axis(side=3, at=rangesizes_at_x, labels=rangesizes, tick=FALSE, hadj=0.5, cex.axis=1.5, line=-0.8)			
	title("y:Sympatric (copying)", line=2, font.main=2, cex.main=1.1)
	mtext("decsize", side=1, cex=0.8, line=0.2, font.main=1, cex.main=0.9)
	
	plotrix::color2D.matplot(x=maxent01s, c(1,0), c(1,0), c(1,0), cellcolors=s_cellcolors, show.values=TRUE, show.legend=FALSE, xlab="", ylab="", axes=FALSE, nslices=100, xaxt="n", yaxt="n")
	axis(side=2, at=rangesizes_at_y, labels=rangesizes, tick=FALSE, padj=0.5, las=1, cex.axis=1.5, line=-0.8)
	axis(side=3, at=rangesizes_at_x, labels=rangesizes, tick=FALSE, hadj=0.5, cex.axis=1.5, line=-0.8)			
	title("s:Sympatric (subset)", line=2, font.main=2, cex.main=1.1)
	mtext("decsize", side=1, cex=0.8, line=0.2, font.main=1, cex.main=0.9)
	
	plotrix::color2D.matplot(x=maxent01v, c(1,0), c(1,0), c(1,0), cellcolors=v_cellcolors, show.values=TRUE, show.legend=FALSE, xlab="", ylab="", axes=FALSE, nslices=100, xaxt="n", yaxt="n")
	axis(side=2, at=rangesizes_at_y, labels=rangesizes, tick=FALSE, padj=0.5, las=1, cex.axis=1.5, line=-0.8)
	axis(side=3, at=rangesizes_at_x, labels=rangesizes, tick=FALSE, hadj=0.5, cex.axis=1.5, line=-0.8)			
	title("v:Vicariance", line=2, font.main=2, cex.main=1.1)
	mtext("decsize", side=1, cex=0.8, line=0.2, font.main=1, cex.main=0.9)
	
	plotrix::color2D.matplot(x=maxent01j, c(1,0), c(1,0), c(1,0), cellcolors=j_cellcolors, show.values=TRUE, show.legend=FALSE, xlab="", ylab="", axes=FALSE, nslices=100, xaxt="n", yaxt="n")
	axis(side=2, at=rangesizes_at_y, labels=rangesizes, tick=FALSE, padj=0.5, las=1, cex.axis=1.5, line=-0.8)
	axis(side=3, at=rangesizes_at_x, labels=rangesizes, tick=FALSE, hadj=0.5, cex.axis=1.5, line=-0.8)			
	title("j:Founder-event (jump)", line=2, font.main=2, cex.main=1.1)
	mtext("decsize", side=1, cex=0.8, line=0.2, font.main=1, cex.main=0.9)
	
	
	par(mar=c(0,0,0,0), xaxs = "i", yaxs = "i") 
	# Row 2 right thin column
	plot(x=c(0,1), y=c(0,1), pch=".", col="white", xaxt="n", yaxt="n", xlab="", ylab="", bty="n")
	testcol = rev(plotrix::color.gradient(c(0,1),c(0,1),c(0,1),nslices=100))
	plotrix::color.legend(xl=0.15, yb=0.3, xr=0.4, yt=0.8, legend=c(0,1), align="rb", rect.col=testcol, gradient="y")	# gradient in x-axis
	


	
	#######################################################
	# Bottom row: depict the cladogenesis process in some fashion.
	# Let's just do the rowSums of the cladogenesis matrix
	#######################################################
	# Rangesize of each state
	tmpstates = states_list
	if ((length(tmpstates[[1]]) == 1) && (is.na(tmpstates[[1]]) == TRUE))
		{
		tmpstates[[1]] = NULL
		}
	rangesizes = sapply(X=tmpstates, FUN=length, simplify=TRUE)
	rangesizes

	# Footer
	#par(mar=c(0,0,0,0), xaxs = "i", yaxs = "i") 
	par(mar=c(5,3,3,1), xaxs = "r", yaxs = "r") 

	tmpca_1 = rep(1, (ncol(tipranges@df)-1))
	tmpcb_1 = rep(1, (ncol(tipranges@df)-1))

	COO_weights_columnar = rcpp_calc_anclikes_sp_COOweights_faster(Rcpp_leftprobs=tmpca_1, Rcpp_rightprobs=tmpcb_1, l=tmpstates, s=s, v=v, j=j, y=y, dmat=dmat, maxent01s=maxent01s, maxent01v=maxent01v, maxent01j=maxent01j, maxent01y=maxent01y, max_minsize_as_function_of_ancsize=max_minsize_as_function_of_ancsize, printmat=FALSE)
	COO_weights_columnar

	# This gives 15 states
	Rsp_rowsums = rcpp_calc_rowsums_for_COOweights_columnar(COO_weights_columnar=COO_weights_columnar)
	Rsp_rowsums

	# Do a count of nonzeros
	COO_weights_columnar_count = COO_weights_columnar
	COO_weights_columnar_count[[4]][COO_weights_columnar_count[[4]] > 0] = 1
	Rsp_rowsums_count = rcpp_calc_rowsums_for_COOweights_columnar(COO_weights_columnar=COO_weights_columnar_count)
	Rsp_rowsums_count
	
	
	
	
	# Make a plot of ancestral range size, versus number of possible descendent pairs
	ylims = c(0, max(max(Rsp_rowsums_count), max(rangesizes)))
	xlims = c(0, max(rangesizes))
	#plot(x=rangesizes, y=Rsp_rowsums, xlim=xlims, ylim=ylims, xlab="anc. range size", ylab="# desc. pairs with prob. > 0", main="Bias towards widespread ancestors")
	
	plot(x=rangesizes, y=Rsp_rowsums_count, xlim=xlims, ylim=ylims, main="", xaxt="n", yaxt="n", tck=0, ylab="", xlab="")
	axis(side=1, at=NULL, labels=NULL, tick=TRUE, padj=-1, las=0, cex.axis=1, tcl=-0.25, line=0, hadj=0.5)
	axis(side=2, at=NULL, labels=NULL, tick=TRUE, padj=1, las=3, cex.axis=1, tcl=-0.25, line=0, hadj=0.5)

	#mtext(text="param. est.", side=2, line=2.1, padj=0.5, adj=0.5, las=3, cex=0.8)#, adj=0.65)

	mtext("anc. range size", side=1, cex=0.8, line=2, font.main=1, cex.main=0.8)
	mtext("# desc. pairs with prob>0", side=2, line=2.1, padj=0.5, adj=0.5, las=3, cex=0.8)#, adj=0.65)
	title("Bias towards\nwidespread ancestors", font.main=2, cex.main=1)




	#######################################################
	# In the last plot, depict a bit of the cladogenesis matrix
	#######################################################
	plot(x=c(0,1), y=c(0,1), pch=".", col="white", xaxt="n", yaxt="n", xlab="", ylab="", bty="n")

	maxrows = 6
	if (length(rangesizes) < maxrows)
		{
		maxrows = length(rangesizes)
		}
	
	
	
	
	# Just take the 1st 8 rows of the COO_weights_columnar, and the last 8
	
	# First part of speciation matrix to print
	ancstates_1st = COO_weights_columnar[[1]][1:maxrows] + 1
	leftstates_1st = COO_weights_columnar[[2]][1:maxrows] + 1
	rightstates_1st = COO_weights_columnar[[3]][1:maxrows] + 1
	relprobs_1st = round(as.numeric(COO_weights_columnar[[4]][1:maxrows] / Rsp_rowsums[ancstates_1st]), 3)
	
	# Last part of speciation matrix to print
	tmpend = length(COO_weights_columnar[[1]])
	tmpstart = tmpend - maxrows
	ancstates_2nd = COO_weights_columnar[[1]][tmpstart:tmpend] + 1
	leftstates_2nd = COO_weights_columnar[[2]][tmpstart:tmpend] + 1
	rightstates_2nd = COO_weights_columnar[[3]][tmpstart:tmpend] + 1
	relprobs_2nd = round((COO_weights_columnar[[4]][tmpstart:tmpend] / Rsp_rowsums[ancstates_2nd]), 3)

	

	
	# Now make a big matrix, with a gap in the middle
	if (is.null(statenames) == TRUE)
		{
		# Fix include_null_range here to false, for plotting
		# cladogenesis probabilities
		statenames = areas_list_to_states_list_new(areas=areanames, maxareas=max_range_size, include_null_range=FALSE, split_ABC=FALSE)
		statenames
		}
	first_mat = rbind(statenames[ancstates_1st], statenames[leftstates_1st], statenames[rightstates_1st], relprobs_1st)
	second_mat = rbind(statenames[ancstates_2nd], statenames[leftstates_2nd], statenames[rightstates_2nd], relprobs_2nd)
	spacer = matrix(" ", nrow=4, ncol=1)
	printmat = cbind(first_mat, spacer, second_mat)
	printmat = adf2(printmat)
	names(printmat) = c(paste("", 1:maxrows, sep=""), "_", paste("", tmpstart:tmpend, sep=""))
	rownames(printmat) = c("Ancestor", "Left desc.", "Right desc.", "Rel. prob.")
	printmat
	
	
	addtable2plot(printmat, x=0, y=0.8, table=printmat, display.rownames=TRUE, display.colnames=TRUE, bty="o", hlines=TRUE, vlines=TRUE, xjust=0, yjust=0, cex=1.1)
	title("Conditional probabilities of example cladogenesis events", font.main=2, cex.main=1)
	
	} # END plot_cladogenesis_size_probabilities














#######################################################
# Convert stochastic map to state probabilities for 
# plotting with plot_BioGeoBEARS_results()
#######################################################


#' Convert stochastic map to state probabilities
#'
#' Converts a stochastic map to state probabilities for 
#' plotting with \code{\link{plot_BioGeoBEARS_results}}().
#'
#' See stochastic mapping example script at PhyloWiki.
#'
#' @param res A \code{BioGeoBEARS_results_object} that is the result of a 
#' \code{\link{bears_optim_run}}. (typically \code{res}, 
#' \code{resDEC}, \code{resDECj}, etc.)
#' @param master_table_cladogenetic_events The master table of cladogenetic
#' events (see stochastic mapping example script at PhyloWiki).
#' @param stratified Is the analysis time-stratified? Default is FALSE.
#' @return something
#' @export
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu}
#' @examples
#' test=1
#' 
stochastic_map_states_into_res <- function(res, master_table_cladogenetic_events, stratified=FALSE)
	{
	defaults='
	resmod = res
	master_table_cladogenetic_events = stochastic_mapping_results2$master_table_cladogenetic_events
	' # 
	
	# Load the tree
	tr = ape::read.tree(file=res$inputs$trfn)
	
	# Error checks -- is this a stratified analysis?
	strat_TF = ("SUBnode.type" %in% names(master_table_cladogenetic_events))
	if ( (strat_TF == TRUE) && (stratified == FALSE) )
		{
		errortxt = paste("\n\nError in stochastic_map_states_into_res():\n\nYou said 'stratified=FALSE' in your inputs (perhaps implicitly), but\nyour master_table_cladogenetic_events has e.g. 'SUBnode.type' in the column names,\nindicating a stratified analysis.\n\n", sep="")
		cat(errortxt)
		
		stop("Stopping on error.")
		}
		
	if ( (strat_TF == FALSE) && (stratified == TRUE) )
		{
		errortxt = paste("\n\nError in stochastic_map_states_into_res():\n\nYou said 'stratified=TRUE' in your inputs , but\nyour master_table_cladogenetic_events LACKS e.g. 'SUBnode.type' in the column names,\nindicating a non-stratified analysis.\n\n", sep="")
		cat(errortxt)
		
		stop("Stopping on error.")
		}	
	
	
	#######################################################
	# NOTE: YOU CANNOT REALLY PLOT THE 'BRANCH BOTTOM' STATES
	# AT BRANCH BOTTOMS, SINCE THEY ARE THE BOTTOM OF THE 
	# SUBTREE ROOT BRANCHES, NOT THE BOTTOM OF THE BRANCHES IN THE MASTER TREE
	# (although they will often, not always, be the same)
	#######################################################
	resmod = res

	# Input the sampled node states into the state probabilities
	resmod$ML_marginal_prob_each_state_at_branch_top_AT_node
	
	# Get the main nodes on the master tree from the master table
	# This has to be done somewhat differently
	if (stratified == TRUE)
		{
		# Get the main nodes on the master tree
		TF1 = master_table_cladogenetic_events$node.type != "tip"
		TF2 = master_table_cladogenetic_events$SUBnode.type != "tip"
		TF3 = master_table_cladogenetic_events$piececlass == "subtree"
		TF = ((TF1 + TF2 + TF3) == 3)
		sum(TF)
		} else {
		# Get the main nodes on the master tree
		TF = master_table_cladogenetic_events$node.type != "tip"
		sum(TF)
		} # END if (stratified == TRUE)
		
	# Look at the cladogenetic events table
	lastcol = ncol(master_table_cladogenetic_events)
	rownums_for_nodes_in_master_tree = (1:nrow(master_table_cladogenetic_events))[TF]
	internal_nodes = master_table_cladogenetic_events$node[rownums_for_nodes_in_master_tree]
	order_rows = order(internal_nodes)
	internal_nodes = sort(internal_nodes)
	master_table_cladogenetic_events[rownums_for_nodes_in_master_tree[order_rows],-lastcol]

	# Look at the sampled states
	master_table_cladogenetic_events$sampled_states_AT_nodes[rownums_for_nodes_in_master_tree[order_rows]]
	master_table_cladogenetic_events$samp_LEFT_dcorner[rownums_for_nodes_in_master_tree[order_rows]]
	master_table_cladogenetic_events$samp_RIGHT_dcorner[rownums_for_nodes_in_master_tree[order_rows]]

	# Zero out the old probabilities, and put in 1s for the sampled states

	# Nodes
	resmod$ML_marginal_prob_each_state_at_branch_top_AT_node[internal_nodes,] = 0
	states_1based_to_change = master_table_cladogenetic_events$sampled_states_AT_nodes[rownums_for_nodes_in_master_tree[order_rows]]
	resmod$ML_marginal_prob_each_state_at_branch_top_AT_node[internal_nodes,states_1based_to_change] = 1

	# Corners
#	if (strat_TF == FALSE)
#		{
		leftright_nodes_matrix = matrix(data=unlist(master_table_cladogenetic_events[rownums_for_nodes_in_master_tree[order_rows],]$daughter_nds), ncol=2, byrow=TRUE)
		Lcorners_above_nodenums = leftright_nodes_matrix[,2]
		Rcorners_above_nodenums = leftright_nodes_matrix[,1]
#		} else {
#		tr = read.tree(res$inputs$trfn)
#		tr_pruningwise = reorder(tr, "pruningwise")
#		leftright_nodes_matrix = get_leftright_nodes_matrix_from_results(tr_pruningwise)
#		Lcorners_above_nodenums = leftright_nodes_matrix[,2]
#		Rcorners_above_nodenums = leftright_nodes_matrix[,1]
#		}

	# Internal node states
	resmod$ML_marginal_prob_each_state_at_branch_top_AT_node[internal_nodes,] = 0
	states_1based_to_change = master_table_cladogenetic_events$sampled_states_AT_nodes[rownums_for_nodes_in_master_tree[order_rows]]
	for (i in 1:length(internal_nodes))
		{
		resmod$ML_marginal_prob_each_state_at_branch_top_AT_node[internal_nodes[i],states_1based_to_change[i]] = 1
		} # END for (i in 1:length(internal_nodes))


	# Get the row numbers in the master table corresponding to the 
	# standard tips and nodes in the unsliced master tree
	ntips = length(tr$tip.label)
	rownums_for_allnodes_in_master_tree = c(1:ntips, rownums_for_nodes_in_master_tree[order_rows])

	resmod$ML_marginal_prob_each_state_at_branch_bottom_below_node[,] = NA
	
	
	#######################################################
	# JUNK TO CUT
	#######################################################
	junk='
	# Branch bottoms OF SUBTREES 
	# (useful, but not for plotting states)
	for (i in 1:length(Lcorners_above_nodenums))
		{

		# Change left corners
		states_1based_to_change = master_table_cladogenetic_events$sampled_states_AT_brbots[rownums_for_allnodes_in_master_tree][Lcorners_above_nodenums[i]]
		resmod$ML_marginal_prob_each_state_at_branch_bottom_below_node[Lcorners_above_nodenums[i],states_1based_to_change] = 1
		resmod$ML_marginal_prob_each_state_at_branch_bottom_below_node[Lcorners_above_nodenums[i],-states_1based_to_change] = 0

		# Right corners
		states_1based_to_change = master_table_cladogenetic_events$sampled_states_AT_brbots[rownums_for_allnodes_in_master_tree][Rcorners_above_nodenums[i]]
		resmod$ML_marginal_prob_each_state_at_branch_bottom_below_node[Rcorners_above_nodenums[i],states_1based_to_change] = 1
		resmod$ML_marginal_prob_each_state_at_branch_bottom_below_node[Rcorners_above_nodenums[i],-states_1based_to_change] = 0
		}
	' #END junk
	#######################################################
	# END JUNK TO CUT
	#######################################################


	for (i in 1:length(internal_nodes))
		{

		# Change left corners
		states_1based_to_change = master_table_cladogenetic_events$samp_LEFT_dcorner[rownums_for_allnodes_in_master_tree][internal_nodes[i]]
		resmod$ML_marginal_prob_each_state_at_branch_bottom_below_node[Rcorners_above_nodenums[i],states_1based_to_change] = 1
		resmod$ML_marginal_prob_each_state_at_branch_bottom_below_node[Rcorners_above_nodenums[i],-states_1based_to_change] = 0

		# Right corners
		states_1based_to_change = master_table_cladogenetic_events$samp_RIGHT_dcorner[rownums_for_allnodes_in_master_tree][internal_nodes[i]]
		resmod$ML_marginal_prob_each_state_at_branch_bottom_below_node[Lcorners_above_nodenums[i],states_1based_to_change] = 1
		resmod$ML_marginal_prob_each_state_at_branch_bottom_below_node[Lcorners_above_nodenums[i],-states_1based_to_change] = 0
		}



	resmod$ML_marginal_prob_each_state_at_branch_bottom_below_node


	# States
	#analysis_titletxt = "Stochastic map"
	#scriptdir = np(system.file("extdata/a_scripts", package="BioGeoBEARS"))
	#res2 = plot_BioGeoBEARS_results(results_object=resmod, analysis_titletxt, addl_params=list("j"), plotwhat="text", label.offset=0.45, tipcex=0.7, statecex=0.7, splitcex=0.6, titlecex=0.8, plotsplits=plotsplits, cornercoords_loc=scriptdir, include_null_range=BioGeoBEARS_run_object$include_null_range, tr=tr, tipranges=tipranges)

	return(resmod)
	} # END plot_stochastic_map_states



#' Get default colors for a states_list
#'
#' Gets default BioGeoBEARS color scheme for a states_list.
#'
#' @param areanames The vector of area abbreviations (single letters)
#' @param states_list_0based A list of states, where each
#' list item is a vector of 0-based area numbers.
#' @param max_range_size The maximum allowed range size.
#' @param plot_null_range Should the null range be included in the list of colors? Default \code{FALSE}.
#' @return colors_list_for_states A vector of colors (hex color codes) for each input state.
#' @export
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu}
#' @examples
#' test=1
#' 
#' areanames = c("K", "O", "M", "H")
#' max_range_size = 4
#' plot_null_range = TRUE
#' states_list_0based=rcpp_areas_list_to_states_list(areas=areanames, maxareas=max_range_size, include_null_range=plot_null_range)
#' get_colors_for_states_list_0based(areanames, states_list_0based=NULL, max_range_size=NA, plot_null_range=FALSE)
#' 
get_colors_for_states_list_0based <- function(areanames, states_list_0based=NULL, max_range_size=NA, plot_null_range=FALSE)
	{
	defaults='
	areanames = c("K", "O", "M", "H")
	max_range_size = 4
	plot_null_range = TRUE
	states_list_0based=rcpp_areas_list_to_states_list(areas=areanames, maxareas=max_range_size, include_null_range=plot_null_range)
	get_colors_for_states_list_0based(areanames, states_list_0based=NULL, max_range_size=NA, plot_null_range=FALSE)
	'
	
	numareas = length(areanames)
	numareas
	
	if (is.na(max_range_size))
		{
		max_range_size = length(areanames)
		}
	max_range_size

	if (is.null(states_list_0based))
		{
		states_list_0based = rcpp_areas_list_to_states_list(areas=areanames, maxareas=max_range_size, include_null_range=plot_null_range)
		}

	# Set up colors
	colors_matrix = get_colors_for_numareas(length(areanames))
	states_list_0based_index = states_list_0based
	
	# Run it 
	# Fix plot_null_range
	colors_list_for_states = mix_colors_for_states(colors_matrix, states_list_0based_index=states_list_0based_index, plot_null_range=plot_null_range)
	colors_list_for_states
	
	return(colors_list_for_states)
	}




# tr is required input (since users could change it in a variety of ways)...

#' Paint branches according to stochastic map
#'
#' This function paints branches according to a BioGeoBEARS stochastic map. 
#'
#' @param res A \code{BioGeoBEARS_results_object} that is the result of a 
#' \code{\link{bears_optim_run}}. (typically \code{res}, 
#' \code{resDEC}, \code{resDECj}, etc.)
#' @param master_table_cladogenetic_events The master table of cladogenetic
#' events (see stochastic mapping example script at PhyloWiki).
#' @param tr An ape \code{phylo} object.
#' @param lwd line width, see \code{\link[graphics]{par}}
#' @param lty line type, see \code{\link[graphics]{par}}
#' @param root.edge An input for \code{node_coords}, which is passed to 
#'        \code{plot_phylo3_nodecoords}. Default \code{TRUE}.
#' @param cornercoords_loc The location of the extdata/scripts directory,
#'        necessary for complex reasons involving CRAN guidelines. Keep on 
#'        default.
#' @param stratified Is the analysis time-stratified? Default is \code{FALSE}.
#' @param plot_clado_1desc If \code{TRUE}, skip a branch if there are no anagenetic events.
#'        Default \code{FALSE}.
#' @param plot_clado_1desc_points If \code{TRUE}, plot symbols for each event, at the point
#'        at which they happened. Experimental. Default \code{FALSE}.
#' @param thickness_by_area When painting, make the thickness of the branch dependent
#'        on the number of areas occupied by the current range. This can help greatly
#'        in helping readers grasp the important issue of residence times in 
#'        different range sizes. Works well. Default \code{FALSE}.
#' @param states_list_0based A list of states, where each
#' list item is a vector of 0-based area numbers.
#' @param include_null_range Should the null range be included in the 
#'        state space? Default is \code{TRUE}.
#' @return times_in_each_state The total time spent in each state in the stochastic
#'         map. This can help greatly in helping readers grasp the important issue 
#'         of residence times in different range sizes.
#' @export
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu}
#' @examples
#' test=1
#' 
paint_stochastic_map_branches <- function(res, master_table_cladogenetic_events, colors_list_for_states, tr=NULL, lwd=5, lty=par("lty"), root.edge=TRUE, cornercoords_loc=np(system.file("extdata/a_scripts", package="BioGeoBEARS")), stratified=FALSE, plot_clado_1desc=FALSE, plot_clado_1desc_points=FALSE, thickness_by_area=FALSE, states_list_0based=NULL, include_null_range=TRUE)
	{
	defaults='
	res = resDEC
	master_table_cladogenetic_events = stochastic_mapping_results$master_table_cladogenetic_events

	plot_stochastic_map_states_stratified(res, tr, master_table_cladogenetic_events, plotsplits=TRUE)

	scriptdir = np(system.file("extdata/a_scripts", package="BioGeoBEARS"))
	#cornercoords_loc="manual"
	cornercoords_loc=scriptdir
	root.edge = TRUE
	lwd = 5
	lty=par("lty")
	
	
	# Get colors_list_for_states
	# Setup 
	#include_null_range = TRUE
	areanames = c("K", "O", "M", "H")
	areas = areanames
	max_range_size = 4
	states_list_0based = rcpp_areas_list_to_states_list(areas=areas, maxareas=max_range_size, include_null_range=include_null_range)

	plot_clado_1desc=FALSE
	plot_clado_1desc_points=FALSE
	thickness_by_area=FALSE

	# Get colors
	plot_null_range = FALSE
	colors_list_for_states = get_colors_for_states_list_0based(areanames=areanames, states_list_0based=states_list_0based, max_range_size=max_range_size, plot_null_range=plot_null_range)
	' # END defaults

	# Keep a 2-column list of the lengths in each state
	list_of_each_state = NULL
	lengths_in_each_state = NULL
	
	
	# Read a tree, if needed (but bad idea!)
	if (is.null(tr))
		{
		#tr = read.tree(res$inputs$trfn)
		tr = check_trfn(trfn=res$inputs$trfn)
		} # END if (is.null(tr))

	# If we want to paint by thickness, make a list of range areas
	if (thickness_by_area == FALSE)
		{
		rangesizes = rep(lwd, times=length(colors_list_for_states))
		} else {
		if (is.null(states_list_0based) == TRUE)
			{
			errortxt = paste("\n\nERROR in paint_stochastic_map_branches(): 'thickness_by_area' is TRUE,\nbut if you want to paint branches with their thickness, you need to supply 'states_list_0based', which you did not.\nHave a hoopy froody day.\n\n", sep="")
			cat(errortxt)
			rangesizes = rep(lwd, times=length(colors_list_for_states))
			#stop(errortxt)
			} else {
			# Calculate the number of areas
			
			# Plotted line with is lwd * numareas
			rangesizes = lwd*sapply(FUN=length, X=states_list_0based)
			rangesizes
			
			# Set the null range to size 0 areas
			if (is.na(states_list_0based[[1]]))
				{
				rangesizes[[1]] = 0
				} # if (is.na(states_list_0based[[1]]))
			} # if (is.null(states_list_0based))
		} # if (is.null(thickness_by_area) == FALSE)



	# Error checks -- is this a stratified analysis?
	strat_TF = ("SUBnode.type" %in% names(master_table_cladogenetic_events))
	if ( (strat_TF == TRUE) && (stratified == FALSE) )
		{
		errortxt = paste("\n\nError in paint_stochastic_map_branches():\n\nYou said 'stratified=FALSE' in your inputs (perhaps implicitly), but\nyour master_table_cladogenetic_events has e.g. 'SUBnode.type' in the column names,\nindicating a stratified analysis.\n\n", sep="")
		cat(errortxt)
		
		stop("Stopping on error.")
		}
		
	if ( (strat_TF == FALSE) && (stratified == TRUE) )
		{
		errortxt = paste("\n\nError in paint_stochastic_map_branches():\n\nYou said 'stratified=TRUE' in your inputs , but\nyour master_table_cladogenetic_events LACKS e.g. 'SUBnode.type' in the column names,\nindicating a non-stratified analysis.\n\n", sep="")
		cat(errortxt)
		
		stop("Stopping on error.")
		}	


	# Get the coordinates of the nodes
	nodes_xy = node_coords(tr, coords_fun="plot_phylo3_nodecoords", tmplocation=cornercoords_loc, root.edge=root.edge)
	nodes_xy

	# Get the maximum height / x value of the tree
	max_x = max(nodes_xy$x)

	# Plot branches with no change
	naTF = is.na(master_table_cladogenetic_events$anagenetic_events_txt_below_node) == TRUE
	master_table_cladogenetic_events$anagenetic_events_txt_below_node[naTF] = "none"
	nochange_TF = master_table_cladogenetic_events$anagenetic_events_txt_below_node == "none"
	sum(nochange_TF)

	# Plot branch histories with no change
	rownum=1
	for (rownum in 1:nrow(master_table_cladogenetic_events))
		{
		
		# Check for cladogenetic events at this node
		# (all nodes except tips)
		# Do the test differently depending on stratified or not
		if (stratified == TRUE)
			{
			test_TF = ((master_table_cladogenetic_events$SUBnode.type[rownum] == "internal") || (master_table_cladogenetic_events$SUBnode.type[rownum] == "root") )
			} else {
			test_TF = ((master_table_cladogenetic_events$node.type[rownum] == "internal") || (master_table_cladogenetic_events$node.type[rownum] == "root") )
			}
		
		
		# Plot the sampled cladogenesis events
		if (test_TF == TRUE )
			{
			# Get the nodenums of the descendants
			master_nodenum = master_table_cladogenetic_events$node[rownum]
			daughter_nds = master_table_cladogenetic_events$daughter_nds[rownum][[1]]
			Ldesc_nodenum = daughter_nds[1]
			Rdesc_nodenum = daughter_nds[2]
		
			# Draw the left descendant
			# Get the corresponding plot coordinates
			time_bp = master_table_cladogenetic_events$time_bp[rownum]
			start_x = max_x - time_bp
			end_x = max_x - time_bp
		
			# Y-coord at node
			start_y = nodes_xy$y[master_nodenum]
		
			# Left Node
			# Y-coord at corner (is just the y-coord of the desc. node)
			end_y = nodes_xy$y[Ldesc_nodenum]
		
			#######################################################
			# PAINT THE LEFT BRANCH (NODE TO CORNER)
			#######################################################
			# Plot the segment
			statenum_1based = as.numeric(master_table_cladogenetic_events$samp_LEFT_dcorner[rownum])
			tmpcolor = colors_list_for_states[statenum_1based]
			lwd_tmp = rangesizes[statenum_1based]
			segments(x0=start_x, x1=end_x, y0=start_y, y1=end_y, col=tmpcolor, lwd=lwd_tmp, lty=lty, lend="round")
			
			# Right Node
			# Y-coord at corner (is just the y-coord of the desc. node)
			end_y = nodes_xy$y[Rdesc_nodenum]
		
			#######################################################
			# PAINT THE RIGHT BRANCH (NODE TO CORNER)
			#######################################################
			# Plot the segment
			statenum_1based = as.numeric(master_table_cladogenetic_events$samp_RIGHT_dcorner[rownum])
			tmpcolor = colors_list_for_states[statenum_1based]
			lwd_tmp = rangesizes[statenum_1based]
			segments(x0=start_x, x1=end_x, y0=start_y, y1=end_y, col=tmpcolor, lwd=lwd_tmp, lty=lty, lend="round")
			} # END plot the sampled cladogenesis events
	
		# If there is no change on this piece of branch
		if (master_table_cladogenetic_events$anagenetic_events_txt_below_node[rownum] == "none")
			{
			starting_state_1based = master_table_cladogenetic_events$sampled_states_AT_brbots[rownum]
			ending_state_1based = master_table_cladogenetic_events$sampled_states_AT_nodes[rownum]
		
			# Get the relative start and stop of this "event" (non-event)
		
			# Get the master nodenum
			master_nodenum = master_table_cladogenetic_events$node[rownum]
			start_y = nodes_xy$y[master_nodenum]
			end_y = nodes_xy$y[master_nodenum]
		
			# Get the x of the top of the branch ( think...
			# The SUB times of this piece are relative to 
			# reltimept
			# This is absolute time before present
			if (stratified == TRUE)
				{
				time_top = master_table_cladogenetic_events$time_top[rownum]
				# This is absolute time before present for the branch piece we are
				# looking at
				branch_piece_TOP_time_bp = time_top + master_table_cladogenetic_events$SUBtime_bp[rownum]
				} else {
				time_top = master_table_cladogenetic_events$time_bp[rownum]
				branch_piece_TOP_time_bp = time_top
				} # END if (stratified == TRUE)

		

			# 2014-05-25_NJM: check if 
			# Is reltimept smaller?
			# Setup
			trtable = master_table_cladogenetic_events
			
			if (stratified == TRUE)
				{
				# Check for NA on branch length below node (e.g. at root)
				if (is.na(trtable$SUBedge.length[rownum]))
					{
					trtable$SUBedge.length[rownum] = 0
					}
		
				# Check if sub-edge length and edge length are the same:
				subedge_length_equals_edge_length_WORRY = FALSE
				if (trtable$SUBedge.length[rownum] > trtable$reltimept[rownum])
					{
					subedge_length_equals_edge_length_WORRY = TRUE

					# We can fix sub-branches easily:
					if (trtable$piececlass[rownum] == "subbranch")
						{
						# Fix to reltimept IF it's *NOT* a fossil:
						if ( (is.na(trtable$fossils[rownum])) || (trtable$fossils[rownum] == FALSE) )
							{
							brlen_in_section = trtable$reltimept[rownum]
							subedge_length_equals_edge_length_WORRY = FALSE
							} else {
							# It *IS* a fossil, it's brlen is based on time_bp
							brlen_in_section = trtable$time_bot[rownum] - trtable$time_bp[rownum]
							subedge_length_equals_edge_length_WORRY = FALSE
							}
						} # END check for fossils
					} else {
					brlen_in_section = trtable$SUBedge.length[rownum]
					}

				if (subedge_length_equals_edge_length_WORRY == TRUE)
					{
					errortxt = paste("\n\nError in plotting_stochastic_maps() (NO CHANGE ON BRANCH): your master_tree table, at row 'rownum'=", rownum, "\nhas an SUBedge.length > reltimept and was not corrected, as it's not a subbranch.\n\n", sep="")
					cat(errortxt)

					cat("\n\n")
					print("rownum:")
					cat("\n\n")
					print(rownum)
					cat("\n\n")
					print("trtable[rownum,]:")
					cat("\n\n")
					print(trtable[rownum,])
					cat("\n\n")
					}
				} else {
				# Stratified == FALSE
				# Check for NA on branch length below node (e.g. at root)
				if (is.na(trtable$edge.length[rownum]))
					{
					trtable$edge.length[rownum] = 0
					}
				# If stratified == FALSE, brlen is just the edge.length
				brlen_in_section = trtable$edge.length[rownum]
				} # END if (stratified == TRUE)

		
			branch_piece_BOT_time_bp = branch_piece_TOP_time_bp + brlen_in_section
		
			# Get the corresponding plot coordinates
			start_x = max_x - branch_piece_BOT_time_bp
			end_x = max_x - branch_piece_TOP_time_bp
			
			# 2017-04-07
			if (start_x > end_x)
				{
				warning("WARNING: start_x > end_x")
				start_x = end_x
				}
			
			
			#######################################################
			# PAINT THE BRANCH WITH NO ANAGENETIC CHANGE
			#######################################################
			# Plot the segment
			statenum_1based = master_table_cladogenetic_events$sampled_states_AT_nodes[rownum]
			tmpcolor = colors_list_for_states[statenum_1based]
			lwd_tmp = rangesizes[statenum_1based]
			# For a shift to NULL, don't plot line
			if (include_null_range==TRUE && statenum_1based==1)
				{
				segments(x0=end_x, x1=end_x, y0=start_y, y1=end_y, col=tmpcolor, lwd=lwd_tmp, lty=lty, lend="butt")
				} else {
				segments(x0=start_x, x1=end_x, y0=start_y, y1=end_y, col=tmpcolor, lwd=lwd_tmp, lty=lty, lend="butt")
				}
			
			list_of_each_state = c(list_of_each_state, statenum_1based)
			lengths_in_each_state = c(lengths_in_each_state, end_x - start_x)

			
			# END if no anagenetic events
			} else {
			# ANAGENETIC EVENTS ON BRANCH
			# Get the master nodenum
			master_nodenum = master_table_cladogenetic_events$node[rownum]
			start_y = nodes_xy$y[master_nodenum]
			end_y = nodes_xy$y[master_nodenum]

			# Get the branch text and extract the history
			branch_events_txt = master_table_cladogenetic_events$anagenetic_events_txt_below_node[rownum]
			events_table_for_branch = events_txt_into_events_table(branch_events_txt)
			
			
			if (plot_clado_1desc == FALSE)
				{
				events_table_for_branch_orig = events_table_for_branch
				
				TF1 = events_table_for_branch$event_type == "d"
				TF2 = events_table_for_branch$event_type == "e"
				TF3 = events_table_for_branch$event_type == "a"
				anagenetic_TF = (TF1 + TF2 + TF3) == 1
				
				events_table_for_branch = events_table_for_branch[anagenetic_TF, ]
				
				# Skip the loop if no anagenetic events
				if (sum(anagenetic_TF) == 0)
					{
					next()
					}
				} # END if (plot_clado_1desc == FALSE)
			
		
			if (stratified == TRUE)
				{
				time_top = master_table_cladogenetic_events$time_top[rownum]
				# This is absolute time before present for the branch piece we are
				# looking at
				branch_piece_TOP_time_bp = time_top + master_table_cladogenetic_events$SUBtime_bp[rownum]
				} else {
				time_top = master_table_cladogenetic_events$time_bp[rownum]
				branch_piece_TOP_time_bp = time_top
				} # END if (stratified == TRUE)

			# Just store the branch length and recall it here
			
			if (stratified == TRUE)
				{
				brlen_in_section = as.numeric(events_table_for_branch$brlen)
				branch_piece_BOT_time_bp = branch_piece_TOP_time_bp + brlen_in_section[1]
				} else {
				brlen_in_section = master_table_cladogenetic_events$edge.length[rownum]
				branch_piece_BOT_time_bp = branch_piece_TOP_time_bp + brlen_in_section
				}
		
		
			# Go through the events along the branch
			current_time_in_subbranch = 0
			
			# I think everything is off by 1 event
			old_painting_method = FALSE
			
			# Sort the branch events by time!!
			events_table_for_branch = events_table_for_branch[order(events_table_for_branch$event_time),]
			
			if (old_painting_method == FALSE)
				{
				# The first "event" is actually the corner state
				# up to the first event
				branch_event_BOT_time_bp = branch_piece_BOT_time_bp - current_time_in_subbranch
				branch_event_TOP_time_bp = as.numeric(events_table_for_branch$abs_event_time[1])
				
				# Need special edit if the ending state is "_" (null) - 2018-01-05
				if (events_table_for_branch$new_rangetxt[1] == "_")
					{
					branch_event_TOP_time_bp = branch_piece_TOP_time_bp
					}
				
				# Get the state at the bottom of the branch
				statenum_1based = master_table_cladogenetic_events$sampled_states_AT_brbots[rownum]
				tmpcolor = colors_list_for_states[statenum_1based]
				lwd_tmp = rangesizes[statenum_1based]
			
				#######################################################
				# PAINT THE FIRST EVENT ON THE BRANCH (new painting method)
				#######################################################
				# Get the corresponding plot coordinates
				start_x = max_x - branch_event_BOT_time_bp
				end_x = max_x - branch_event_TOP_time_bp
				# For a shift to NULL, don't plot line
				if (include_null_range==TRUE && statenum_1based==1)
					{
					segments(x0=end_x, x1=end_x, y0=start_y, y1=end_y, col=tmpcolor, lwd=lwd_tmp, lty=lty, lend="butt")
					} else {
					segments(x0=start_x, x1=end_x, y0=start_y, y1=end_y, col=tmpcolor, lwd=lwd_tmp, lty=lty, lend="butt")
					}

				current_time_in_subbranch = branch_event_TOP_time_bp
				
				list_of_each_state = c(list_of_each_state, statenum_1based)
				lengths_in_each_state = c(lengths_in_each_state, end_x - start_x)
				}
			
			# Loop through the events
			for (j in 1:nrow(events_table_for_branch))
				{
				# Time above the bottom of the branch
				#event_time = branch_piece_BOT_time_bp - #as.numeric(events_table_for_branch$event_time[j])
				#as.numeric(events_table_for_branch$event_time[j])
			
				# Calculate the time_bp of event in x
				if (old_painting_method == TRUE)
					{
					branch_event_BOT_time_bp = branch_piece_BOT_time_bp - current_time_in_subbranch
					branch_event_TOP_time_bp = as.numeric(events_table_for_branch$abs_event_time[j])
					} else {
					# Event bottom
					#branch_event_BOT_time_bp = branch_piece_BOT_time_bp - current_time_in_subbranch
					branch_event_BOT_time_bp = as.numeric(events_table_for_branch$abs_event_time[j])
					# Event top
					if (j < nrow(events_table_for_branch))
						{
						branch_event_TOP_time_bp = as.numeric(events_table_for_branch$abs_event_time[j+1])
						} else {
						branch_event_TOP_time_bp = as.numeric(branch_piece_TOP_time_bp)
						}
					}
			
				#cat("\n\n")
				#cat(events_table_for_branch$nodenum_at_top_of_branch[j], ": ", branch_event_BOT_time_bp, " -- ", branch_event_TOP_time_bp, sep="")
				#cat("\n\n")
			
				
			
				# Get the corresponding plot coordinates
				start_x = max_x - branch_event_BOT_time_bp
				end_x = max_x - branch_event_TOP_time_bp
			
				# Plot the segment
				# The old method used the current rangenum
				# The new method uses the new rangenum
				if (old_painting_method == TRUE)
					{
					statenum_1based = as.numeric(events_table_for_branch$current_rangenum_1based[j])
					} else {
					statenum_1based = as.numeric(events_table_for_branch$new_rangenum_1based[j])
					}
				
				#######################################################
				# PAINT THE NEXT EVENT ON THE BRANCH (new painting method)
				#######################################################
				# Paint the segment
				tmpcolor = colors_list_for_states[statenum_1based]
				lwd_tmp = rangesizes[statenum_1based]

				# For a shift to NULL, don't plot line
				if (include_null_range==TRUE && statenum_1based==1)
					{
					segments(x0=end_x, x1=end_x, y0=start_y, y1=end_y, col=tmpcolor, lwd=lwd_tmp, lty=lty, lend="butt")
					} else {
					segments(x0=start_x, x1=end_x, y0=start_y, y1=end_y, col=tmpcolor, lwd=lwd_tmp, lty=lty, lend="butt")
					}

				
				list_of_each_state = c(list_of_each_state, statenum_1based)
				lengths_in_each_state = c(lengths_in_each_state, end_x - start_x)


				if (plot_clado_1desc_points == TRUE)
					{
					# Default is NA (don't plot)
					pointchar = NA
					
					if ((events_table_for_branch$event_type[j] == "sympatry (y)") || (events_table_for_branch$event_type == "y") )
						{
						# Character to plot
						pointchar = "y\n|"
						}
					if ((events_table_for_branch$event_type[j] == "vicariance (v)") || (events_table_for_branch$event_type == "v") )
						{
						# Character to plot
						pointchar = "v\n|"
						}
					if ((events_table_for_branch$event_type[j] == "subset (s)") || (events_table_for_branch$event_type == "s") )
						{
						# Character to plot
						pointchar = "s\n|"
						}
					if ((events_table_for_branch$event_type[j] == "founder (j)") || (events_table_for_branch$event_type == "j") )
						{
						# Character to plot
						pointchar = "j\n|"
						}
					if (events_table_for_branch$event_type[j] == "d")
						{
						# Character to plot
						pointchar = "d\n|"
						}
					if (events_table_for_branch$event_type[j] == "e")
						{
						# Character to plot
						pointchar = "e\n|"
						}
					if (events_table_for_branch$event_type[j] == "a")
						{
						# Character to plot
						pointchar = "a\n|"
						}
					
					# Plot the point
					if (!is.na(pointchar))
						{
						# Plot points
						# points(x=start_x, y=start_y, pch="+", col="black", cex=1)
						
						# Plot some text
						text(x=start_x, y=start_y, labels=pointchar, adj=c(0.5,0.15), col=tmpcolor, cex=0.85)
						#text(x=start_x, y=start_y, labels=pointchar, adj=c(0.5,0.15), col="black", cex=0.85)
						} # END if (!is.na(pointchar))
					} # END if (plot_clado_1desc_points == TRUE)

				# Update the time of the base of the next event
				current_time_in_subbranch = branch_event_BOT_time_bp - branch_event_TOP_time_bp
				} # END for (j in 1:nrow(events_table_for_branch))
		

			# Calculate the time_bp of event in x
			branch_event_BOT_time_bp = branch_event_TOP_time_bp
			branch_event_TOP_time_bp = branch_piece_TOP_time_bp

			# This shouldn't be needed anymore
			if (old_painting_method == TRUE)
				{
				# Then, after the last change event, you still HAVE to plot the last chunk
				# of branch events
				j = nrow(events_table_for_branch)

				# Time above the bottom of the branch
				#event_time = as.numeric(events_table_for_branch$event_time[j])
				#event_time = branch_event_TOP_time_bp
		
				#cat("\n\n")
				#cat(events_table_for_branch$nodenum_at_top_of_branch[j], ": ", branch_event_BOT_time_bp, " -- ", branch_event_TOP_time_bp, sep="")
				#cat("\n\n")

		
		
				# Get the corresponding plot coordinates
				start_x = max_x - branch_event_BOT_time_bp
				end_x = max_x - branch_event_TOP_time_bp
		
				#######################################################
				# PAINT the remaining part of the branch (old painting method)
				#######################################################
				# Plot the segment
				statenum_1based = as.numeric(events_table_for_branch$new_rangenum_1based[j])
				tmpcolor = colors_list_for_states[statenum_1based]
				lwd_tmp = rangesizes[statenum_1based]
				
				# For a shift to NULL, don't plot line
				if (include_null_range==TRUE && statenum_1based==1)
					{
					segments(x0=end_x, x1=end_x, y0=start_y, y1=end_y, col=tmpcolor, lwd=lwd_tmp, lty=lty, lend="butt")
					} else {
					segments(x0=start_x, x1=end_x, y0=start_y, y1=end_y, col=tmpcolor, lwd=lwd_tmp, lty=lty, lend="butt")
					}

				list_of_each_state = c(list_of_each_state, statenum_1based)
				lengths_in_each_state = c(lengths_in_each_state, end_x - start_x)
				} # END if (old_painting_method == TRUE)
			} # END if anagenetic events
		} # END for (rownum in 1:nrow(master_table_cladogenetic_events))
		  # (ENDING loop through the branches and subbranches in
		  #  master_table_cladogenetic_events)
	
	times_in_each_state = NULL
	times_in_each_state$list_of_each_state = list_of_each_state
	times_in_each_state$lengths_in_each_state = lengths_in_each_state

	return(times_in_each_state)
	} # END function



#######################################################
# Plot stochastic maps
#######################################################


#' Plot stochastic maps
#' 
#' This function is copied and modified from \code{\link{plot_BioGeoBEARS_results}} to 
#' plot biogeographical stochastic maps. Many of the inputs are used only to 
#' plot the background tree. See PhyloWiki for example scripts.
#' 
#' @param results_object The results object from \code{\link{bears_optim_run}} (with ancestral states on).
#' @param clado_events_table The table of cladogenetic events (see stochastic mapping example script at PhyloWiki).
#' @param stratified Is the analysis time-stratified? Default is \code{FALSE}.
#' @param analysis_titletxt The main title of the plot. If NULL, \code{results_object$inputs$description} is checked.
#' @param addl_params The function will plot the log-likelihood (LnL) and the ML values of the free parameters. If you want additional parameters plotted, list them here.
#' @param plotwhat To plot the ML discrete states, "text".  To plot a piechart of the relative probability of all the states, "pie".
#' @param label.offset Offset for the tree tip labels. If \code{NULL}, program chooses 0.05 x tree height.
#' @param tipcex \code{cex} value for the tiplabels (scaling factor, i.e. 0.5 is half size)
#' @param statecex \code{cex} value for the states (scaling factor, i.e. 0.5 is half size). Used on piecharts if plotwhat="pie".
#' @param splitcex \code{cex} value for the splits (scaling factor, i.e. 0.5 is half size). Used on piecharts if plotwhat="pie".
#' @param titlecex \code{cex} value for the title (scaling factor, i.e. 0.5 is half size). 
#' @param plotsplits If \code{TRUE}, plot states on the corners -- text or pie charts, depending on \code{plotwhat}.
#' @param plotlegend If \code{TRUE}, make a (separate) plot with a legend giving the colors for each state/range, using \code{\link{colors_legend}}.
#' @param legend_ncol The number of columns in the legend.  If \code{NULL} (default), the function calculates \code{floor(sqrt(length(possible_ranges_list_txt) / 2))} 
#' when the number of states is <=64, and \code{sqrt(ceiling(length(possible_ranges_list_txt)))} when > 64. Note that when you have hundreds of states, there is probably 
#' no good way to have a readable legend, and it is easier to just rely upon printing the 
#' character codes for the ML states in the plots, with the colors, and users can then see and trace the common colors/states by eye.
#' @param legend_cex The cex (character expansion size) for the legend.  Defaults to 1, which means the \code{\link[graphics]{legend}} function determines the 
#' size.  The value 2.5 works well for 15 or 16 states/ranges.
#' @param cornercoords_loc The directory location containing the R script \code{plot_phylo3_nodecoords.R}. This function, modified from the APE function
#' \code{\link[ape]{plot.phylo}}, cannot be included directly in the R package as it contains C code that does not pass CRAN's R CMD check. The default, 
#' cornercoords_loc="manual", will not allow split states to be plot.  The R script \code{plot_phylo3_nodecoords.R} is located in the BioGeoBEARS extension data 
#' directory, \code{extdata/a_scripts}.  You should be able to get the full path with \code{list.files(system.file("extdata/a_scripts", package="BioGeoBEARS"), full.names=TRUE)}.
#' @param include_null_range If \code{TRUE} (default), the null range is included in calculation of colors. (Safest for now.)
#' @param tr Tree to plot on. Default \code{NULL}, which means the tree will be read from the file at \code{results_object$inputs$trfn}.
#' @param tipranges Tip geography data. Default \code{NULL}, which means the tree will be read from the file at \code{results_object$inputs$geogfn}.
#' @param if_ties What to do with ties in probability. Currently, 
#' the options are:
#' (1) "takefirst", which takes the first tied state in the 
#' probabilities column (The full probabilities of all states will be
#' shown in e.g. pie charts, of course); (2) "asterisk", which 
#' returns returns the first tied state, but marks it with an 
#' asterisk ("*").
#' @param pie_tip_statecex  \code{cex} value for pies at the tips (scaling factor, i.e. 0.5 is half size). 
#' @param juststats If \code{FALSE} (default), plots are plotted. If \code{TRUE}, 
#' no plots are done, the function just returns the summary statistics.
#' @param root.edge Should the root edge be plotted, if it exists in the tree?  Passed to
#' plot.phylo().  This can be handy if the root state display is getting cut off.
#' @param colors_list_for_states Default \code{NULL} auto-generates colors with 
#' \code{get_colors_for_states_list_0based}. Otherwise, users can specify colors for 
#' each state (remember that e.g. 4 areas can mean 2^4=16 states).
#' @param skiptree Skip the plotting of the tree -- useful for plotting the state labels
#' e.g. on top of a stochastic map. Default \code{FALSE}.
#' @param show.tip.label Same as for APE's \code{plot.phylo}.
#' @param tipcol The tip text color. Default "black".
#' @param dej_params_row Parameters used to generate an SSE simulation. dej_params_row can be 
#' obtained from \code{get_dej_params_row_from_SSEsim_results_processed}.
#' @param plot_max_age The maximum age of the plot (age below the tips). Default
#' is tree height
#' @param skiplabels If \code{TRUE}, skip the nodelabels command, resulting in plotting just the 
#' tree. (Yes, you could have just used \code{FALSE}).  Default \code{FALSE}.
#' @param plot_stratum_lines If \code{TRUE} (default), plot dotted lines at the stratum 
#' boundaries, *if* it's a time-stratified results object.
#' @param simplify_piecharts If \code{TRUE}, just plot one slice for the maximum probability, 
#' and white for "other". This should help with large plots with many ranges, 
#' which can overwhelm some graphics programs (imagine 1000 slices per piechart,
#' times 1000+nodes). Default \code{FALSE}.
#' @param include_null_range Should the null range be included in the 
#'        state space? Default is \code{NULL} (inferred from objects).
#' @param plot_null_range Should the null range be plotted? Default is \code{FALSE}.
#' @return NULL
#' @export
#' @seealso \code{\link{get_leftright_nodes_matrix_from_results}}, \code{\link{corner_coords}}, \code{\link[ape]{plot.phylo}}, \code{\link[ape]{plot.phylo}}, \code{\link[ape]{tiplabels}}, \code{\link[graphics]{legend}}, \code{\link[base]{floor}}, \code{\link[base]{ceiling}}, \code{\link[base]{floor}}, \code{\link[cladoRcpp]{numstates_from_numareas}}, \code{\link[base]{system.file}}, \code{\link[base]{list.files}}
#' @note Go (BioGeo)BEARS!
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu} 
#' @references
#' \url{https://code.google.com/p/lagrange/}
#' @examples
#' test=1
#' # See example scripts at PhyloWiki. 
#' 
plot_BSM <- function(results_object, clado_events_table, stratified, analysis_titletxt="Stochastic map", addl_params=list(), plotwhat="text", label.offset=NULL, tipcex=0.8, statecex=0.7, splitcex=0.6, titlecex=0.8, plotsplits=TRUE, plotlegend=FALSE, legend_ncol=NULL, legend_cex=1, cornercoords_loc="manual", tr=NULL, tipranges=NULL, if_ties="takefirst", pie_tip_statecex=0.7, juststats=FALSE, xlab="Millions of years ago", root.edge=TRUE, colors_list_for_states, lwd=5, skiptree=FALSE, show.tip.label=TRUE, tipcol="black", dej_params_row=NULL, plot_max_age=NULL, skiplabels=FALSE, plot_stratum_lines=TRUE, include_null_range=NULL, plot_null_range=FALSE)
	{
	scriptdir = np(system.file("extdata/a_scripts", package="BioGeoBEARS"))
	
	# Convert the BSM into a modified res object
	resmod = stochastic_map_states_into_res(res=results_object, master_table_cladogenetic_events=clado_events_table, stratified=stratified)

	# Plot the tree and states at nodes/corners
	# (copying everything from the inputs; mostly these should be kept on defaults)
	# (skiptree=FALSE the first time, TRUE the second time)
	plot_BioGeoBEARS_results(results_object=resmod, analysis_titletxt=analysis_titletxt, addl_params=addl_params, plotwhat=plotwhat, label.offset=label.offset, tipcex=tipcex, statecex=statecex, splitcex=splitcex, titlecex=titlecex, plotsplits=plotsplits, plotlegend=plotlegend, legend_ncol=legend_ncol, legend_cex=legend_cex, cornercoords_loc=cornercoords_loc, tr=tr, tipranges=tipranges, if_ties=if_ties, pie_tip_statecex=pie_tip_statecex, juststats=juststats, xlab=xlab, root.edge=root.edge, colors_list_for_states=colors_list_for_states, skiptree=FALSE, show.tip.label=show.tip.label, tipcol=tipcol, dej_params_row=dej_params_row, plot_max_age=plot_max_age, skiplabels=skiplabels, plot_stratum_lines=plot_stratum_lines, include_null_range=include_null_range, plot_null_range=plot_null_range)
	
	# Paint on the branch states
	paint_stochastic_map_branches(res=resmod, master_table_cladogenetic_events=clado_events_table, colors_list_for_states=colors_list_for_states, lwd=lwd, lty=par("lty"), root.edge=TRUE, stratified=stratified)

	# Re-plot the tree to get the states on top
	# (skiptree=TRUE this time)
	plot_BioGeoBEARS_results(results_object=resmod, analysis_titletxt=analysis_titletxt, addl_params=addl_params, plotwhat=plotwhat, label.offset=label.offset, tipcex=tipcex, statecex=statecex, splitcex=splitcex, titlecex=titlecex, plotsplits=plotsplits, plotlegend=plotlegend, legend_ncol=legend_ncol, legend_cex=legend_cex, cornercoords_loc=cornercoords_loc, tr=tr, tipranges=tipranges, if_ties=if_ties, pie_tip_statecex=pie_tip_statecex, juststats=juststats, xlab=xlab, root.edge=root.edge, colors_list_for_states=colors_list_for_states, skiptree=TRUE, show.tip.label=show.tip.label, tipcol=tipcol, dej_params_row=dej_params_row, plot_max_age=plot_max_age, skiplabels=skiplabels, plot_stratum_lines=plot_stratum_lines, include_null_range=include_null_range, plot_null_range=plot_null_range)	
	
	return(NULL)
	} # END plot_BSM


# Various problems emerge from "ladderize" in some versions
# https://www.mail-archive.com/r-sig-phylo@r-project.org/msg04176.html

#' Ladderize nodes, ensuring that resulting tree is appropriately ordered.
#'
#' This function avoids some potential problems with \code{\link[ape]{ladderize}}
#' by writing the ladderized tree to a newick string, and reading back
#' in.
#'
#' For discussion of the issues, see: 
#' https://www.mail-archive.com/r-sig-phylo@r-project.org/msg04176.html
#'
#' @param phy An ape phylo object
#' @param right TRUE for right-ordering, FALSE for left-ordering.
#' @return ltr A ladderized tree
#' @export
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu}
#' @examples
#' test=1
#' 
#' # You can see the issues here
#' tr = ape::read.tree(file="", text="((human:1,chimp:1):1,gorilla:2);")
#' ltr = ladderize(tr, right=TRUE)
#' ltr$edge
#' ltr = ladderize(tr, right=FALSE)
#' ltr$edge
#' ltr = ladderize_and_reorder(tr, right=TRUE)
#' ltr$edge
#' ltr = ladderize_and_reorder(tr, right=FALSE)
#' ltr$edge
#' 
ladderize_and_reorder <- function(phy, right=TRUE)
	{
	ex='
	# You can see the issues here
	tr = read.tree(file="", text="((human:1,chimp:1):1,gorilla:2);")
	ltr = ape::ladderize(tr, right=TRUE)
	ltr$edge
	ltr = ape::ladderize(tr, right=FALSE)
	ltr$edge
	ltr = ladderize_and_reorder(tr, right=TRUE)
	ltr$edge
	ltr = ladderize_and_reorder(tr, right=FALSE)
	ltr$edge
	'

	ltr = ape::ladderize(phy, right=right)
	# MAKE FREAKING SURE that this tree has the right node order etc.
	ltr = ape::read.tree(file="", text=write.tree(phy=ltr, file="") )
	return(ltr)
	} # END ladderize_and_reorder <- function(tr, right=TRUE)


# Returns:
# indexes_to_convert_tr2nodes_to_tr1
# NA for non-matches

# E.g. 
# tr1 = new tree
# tr2 = original tree, used the BioGeoBEARS analysis

#' Get the indexes to match tree2 nodes to tree1
#'
#' Given 2 trees with different rotations etc., this function
#' returns, for each tree2 node, the index of that tree2 node 
#' in the tree1 nodes list.
#'
#' (Basically, this is \code{\link[base]{match}} for tree nodes.)
#'
#' @param tr1 An ape \code{phylo} object
#' @param tr2 An ape \code{phylo} object
#' @return indexes_to_convert_tr2nodes_to_tr1 The vector of indices
#' @export
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu}
#' @examples
#' test=1
#' 
#' tr = ape::read.tree(file="", text="((human:1,chimp:1):1,gorilla:2);")
#' ltr1 = ladderize_and_reorder(tr, right=TRUE)
#' ltr2 = ladderize_and_reorder(tr, right=FALSE)
#' 
#' # Compare the tree tables
#' prt(ltr1, printflag=FALSE, get_tipnames=TRUE)
#' prt(ltr2, printflag=FALSE, get_tipnames=TRUE)
#' 
#' # Get the matching
#' indexes_to_convert_tr2nodes_to_tr1 = ordernodes(tr1=ltr1, tr2=ltr2)
#' indexes_to_convert_tr2nodes_to_tr1
#' 
#' # What if 1 tree is missing a taxon?
#' tr1 = ape::read.tree(file="", text="((human:1,chimp:1):1,gorilla:2);")
#' tr2 = ape::read.tree(file="", text="(((human:1,chimp:1):1,gorilla:2):1,orang:3);")
#' 
#' # Compare the tree tables
#' prt(tr1, printflag=FALSE, get_tipnames=TRUE)
#' prt(tr2, printflag=FALSE, get_tipnames=TRUE)
#' 
#' # Get the matching
#' indexes_to_convert_tr2nodes_to_tr1 = ordernodes(tr1=tr1, tr2=tr2)
#' indexes_to_convert_tr2nodes_to_tr1
#' 
#' # Reverse matching
#' \dontrun{
#' # (Works, but throws a warning)
#' indexes_to_convert_tr1nodes_to_tr2 = ordernodes(tr1=tr2, tr2=tr1)
#' indexes_to_convert_tr1nodes_to_tr2
#' }
#' 
ordernodes <- function(tr1, tr2)
	{
	ex='
	# Make 2 differently-rotated trees
	tr = ape::read.tree(file="", text="((human:1,chimp:1):1,gorilla:2);")
	ltr1 = ladderize_and_reorder(tr, right=TRUE)
	ltr2 = ladderize_and_reorder(tr, right=FALSE)

	# Compare the tree tables
	prt(ltr1, printflag=FALSE, get_tipnames=TRUE)
	prt(ltr2, printflag=FALSE, get_tipnames=TRUE)
	
	# Get the matching
	indexes_to_convert_tr2nodes_to_tr1 = ordernodes(tr1=ltr1, tr2=ltr2)
	indexes_to_convert_tr2nodes_to_tr1
	
	# What if 1 tree is missing a taxon?
	tr1 = ape::read.tree(file="", text="((human:1,chimp:1):1,gorilla:2);")
	tr2 = ape::read.tree(file="", text="(((human:1,chimp:1):1,gorilla:2):1,orang:3);")

	# Compare the tree tables
	prt(tr1, printflag=FALSE, get_tipnames=TRUE)
	prt(tr2, printflag=FALSE, get_tipnames=TRUE)
	
	# Get the matching
	indexes_to_convert_tr2nodes_to_tr1 = ordernodes(tr1=tr1, tr2=tr2)
	indexes_to_convert_tr2nodes_to_tr1
	
	# Reverse matching
	# (Works, but throws a warning)
	indexes_to_convert_tr1nodes_to_tr2 = ordernodes(tr1=tr2, tr2=tr1)
	indexes_to_convert_tr1nodes_to_tr2
	
	'

	tr1_table = prt(tr1, printflag=FALSE, get_tipnames=TRUE)
	tr2_table = prt(tr2, printflag=FALSE, get_tipnames=TRUE)
	indexes_to_convert_tr2nodes_to_tr1 = match(x=tr1_table$tipnames, table=tr2_table$tipnames)
	
	if (any(is.na(indexes_to_convert_tr2nodes_to_tr1)) == TRUE)
		{
		txt = "WARNING in ordernodes() or BioGeoBEARS_reorder() -- not all nodes match between tr1 and tr2. Any non-matching nodes will get NAs."
		cat("\n\n")
		cat(txt)
		cat("\n\n")
		
		warning(txt)
		} # END if (any(is.na(indexes_to_convert_tr2nodes_to_tr1)) == TRUE)
	
	return(indexes_to_convert_tr2nodes_to_tr1)
	} # END ordernodes <- function(tr1, tr2)


# Convert BioGeoBEARS object from one tree to another
# tr1 = new tree
# tr2 = original tree, used the BioGeoBEARS analysis

#' Convert BioGeoBEARS results object from one tree to another
#'
#' Converts a BioGeoBEARS results object from one tree to another. This can be 
#' useful if you ran a long, slow analysis on a particular tree, but then
#' decide that you want to plot the analysis on a different tree. If you just
#' switch the tree files, you will get nonsense, as the ancestral state/range
#' probabilities are tied to the node numbers of the original tree.
#' 
#' This function reorders the rows of all of the ancestral states/ranges tables
#' (and likelihoods etc.) to match a new input tree.
#'
#' @param res A \code{BioGeoBEARS_results_object} that is the result of a 
#' \code{\link{bears_optim_run}}. (typically \code{res}, 
#' \code{resDEC}, \code{resDECj}, etc.)
#' @param tr1 The new tree. An ape \code{phylo} object
#' @param tr2 The original tree used the BioGeoBEARS analysis. An ape \code{phylo} object
#' @param trfn_for_BGB_inputs The filename of the new tree, if you want to update the
#' the \code{res} object. (You probably should, as some other functions will need to 
#' read in the tree.)
#' @return res2 The results object with the probabilities reordered for the new tree.
#' @export
#' @author Nicholas J. Matzke \email{matzke@@berkeley.edu}
#' @examples
#' test=1
BioGeoBEARS_reorder <- function(res, tr1, tr2, trfn_for_BGB_inputs=NULL)
	{
	indexes_to_convert_tr2nodes_to_tr1 = ordernodes(tr1, tr2)
	res2 = res
	
	# Vector
	res2$computed_likelihoods_at_each_node = res$computed_likelihoods_at_each_node[indexes_to_convert_tr2nodes_to_tr1]
	
	# Matrices
	res2$relative_probs_of_each_state_at_branch_top_AT_node_DOWNPASS = res$relative_probs_of_each_state_at_branch_top_AT_node_DOWNPASS[indexes_to_convert_tr2nodes_to_tr1,]
	
	res2$condlikes_of_each_state = res$condlikes_of_each_state[indexes_to_convert_tr2nodes_to_tr1,]
	
	res2$relative_probs_of_each_state_at_branch_bottom_below_node_DOWNPASS = res$relative_probs_of_each_state_at_branch_bottom_below_node_DOWNPASS[indexes_to_convert_tr2nodes_to_tr1,]
	
	res2$relative_probs_of_each_state_at_branch_bottom_below_node_UPPASS = res$relative_probs_of_each_state_at_branch_bottom_below_node_UPPASS[indexes_to_convert_tr2nodes_to_tr1,]
	
	res2$relative_probs_of_each_state_at_branch_top_AT_node_UPPASS = res$relative_probs_of_each_state_at_branch_top_AT_node_UPPASS[indexes_to_convert_tr2nodes_to_tr1,]
	
	res2$ML_marginal_prob_each_state_at_branch_bottom_below_node = res$ML_marginal_prob_each_state_at_branch_bottom_below_node[indexes_to_convert_tr2nodes_to_tr1,]
	
	res2$ML_marginal_prob_each_state_at_branch_top_AT_node = res$ML_marginal_prob_each_state_at_branch_top_AT_node[indexes_to_convert_tr2nodes_to_tr1,]
	
	# Input the tree filename into BioGeoBEARS inputs, if desired
	if (is.null(trfn_for_BGB_inputs) == FALSE)
		{
		res2$inputs$trfn = trfn_for_BGB_inputs
		}
	
	return(res2)
	} # END BioGeoBEARS_reorder <- function(res, tr1, tr2)









# Not exported

#######################################################
# These functions fix graphics issues that arise in APE 5.0
# 
# They were pointed out by Liz, here:
#
# https://groups.google.com/forum/#!topic/biogeobears/gQ4bZ4U3FU8
# 
# BioGeoBEARS
# node_height error
# 1 post by 1 author  
# 
# Liz	
# 
# 12:25 PM (2 hours ago)
# 
# Hi all,
# 
# I'm a new user of BioGeoBEARS and I'm just trying to get the sample 
# script running.  I'm hoping the following is a simple issue....  
# 
# When I run the test script, everything runs fine until this command:
# 
# res1 = plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params=list("j"), plotwhat="text", label.offset=0.45, tipcex=0.7, statecex=0.7, splitcex=0.6, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges)
# 
# I get the following traceback error in R Studio:
# 
#  Error in .C("node_height", as.integer(Ntip), as.integer(Nnode), 
#               as.integer(edge[,  : 
#   Incorrect number of arguments (6), expecting 4 for 'node_height' 
# 6. .nodeHeight(Ntip, Nnode, z$edge, Nedge, yy) at plot_phylo3_nodecoords.R#156
# 5. plot_phylo3_nodecoords(tr, plot = FALSE, root.edge = root.edge) at <text>#1
# 4. eval(parse(text = cmdstr)) 
# 3. eval(parse(text = cmdstr)) at BioGeoBEARS_plots_v1.R#1871
# 2. node_coords(tr, tmplocation = cornercoords_loc, root.edge = root.edge) at BioGeoBEARS_plots_v1.R#457
# 1. plot_BioGeoBEARS_results(results_object, analysis_titletxt, addl_params = list("j"), 
#  plotwhat = "text", label.offset = 0.45, tipcex = 0.7, statecex = 0.7, 
#  splitcex = 0.6, titlecex = 0.8, plotsplits = TRUE, 
#  cornercoords_loc=scriptdir, 
#  include_null_range = TRUE, tr = tr, tipranges = tipranges) 
# 
# Please note I am using ape version 5.0.
# 
# Thanks for your help!
# Liz
#######################################################


# Access par (graphics parameters) without opening a &%$@ plot!
# https://stackoverflow.com/questions/20363266/how-can-i-suppress-the-creation-of-a-plot-while-calling-a-function-in-r
# Internal function
par_invisible <- function(parname)
	{
	setup='
	pin1 = par_invisible(parname="pin")
	'

	ff <- tempfile()
	png(filename=ff)
	tmp = do.call(par, args=list(parname))
	res = tmp[1]
	dev.off()
	unlink(ff)
	return(res)
	}


# Internal function
plot_phylo3_nodecoords_APE5 <- function (x, type = "phylogram", use.edge.length = TRUE, node.pos = NULL, 
    show.tip.label = TRUE, show.node.label = FALSE, edge.color = "black", 
    edge.width = 1, edge.lty = 1, font = 3, 
    adj = NULL, srt = 0, no.margin = FALSE, root.edge = FALSE, 
    label.offset = 0, underscore = FALSE, x.lim = NULL, y.lim = NULL, 
    direction = "rightwards", lab4ut = NULL, tip.color = "black", 
    plot = TRUE, rotate.tree = 0, open.angle = 0, node.depth = 1, 
    align.tip.label = FALSE, cex_original=FALSE, ...) 
{
    
    # Original behavior in APE 5.0 plot.phylo
    if (cex_original == TRUE)
    	{
    	cex = par("cex")
    	} else {
    	cex = par_invisible(parname="cex")
    	}
    
    Ntip <- length(x$tip.label)
    if (Ntip < 2) {
        warning("found less than 2 tips in the tree")
        return(NULL)
    }
    .nodeHeight <- function(edge, Nedge, yy) .C(node_height, 
        as.integer(edge[, 1]), as.integer(edge[, 2]), as.integer(Nedge), 
        as.double(yy))[[4]]
    .nodeDepth <- function(Ntip, Nnode, edge, Nedge, node.depth) .C(node_depth, 
        as.integer(Ntip), as.integer(edge[, 1]), as.integer(edge[, 
            2]), as.integer(Nedge), double(Ntip + Nnode), as.integer(node.depth))[[5]]
    .nodeDepthEdgelength <- function(Ntip, Nnode, edge, Nedge, 
        edge.length) .C(node_depth_edgelength, as.integer(edge[, 
        1]), as.integer(edge[, 2]), as.integer(Nedge), as.double(edge.length), 
        double(Ntip + Nnode))[[5]]
    Nedge <- dim(x$edge)[1]
    Nnode <- x$Nnode
    if (any(x$edge < 1) || any(x$edge > Ntip + Nnode)) 
        stop("tree badly conformed; cannot plot. Check the edge matrix.")
    ROOT <- Ntip + 1
    type <- match.arg(type, c("phylogram", "cladogram", "fan", 
        "unrooted", "radial"))
    direction <- match.arg(direction, c("rightwards", "leftwards", 
        "upwards", "downwards"))
    if (is.null(x$edge.length)) {
        use.edge.length <- FALSE
    } else {
        if (use.edge.length && type != "radial") {
            tmp <- sum(is.na(x$edge.length))
            if (tmp) {
                warning(paste(tmp, "branch length(s) NA(s): branch lengths ignored in the plot"))
                use.edge.length <- FALSE
            }
        }
    }
    if (is.numeric(align.tip.label)) {
        align.tip.label.lty <- align.tip.label
        align.tip.label <- TRUE
    } else {
        if (align.tip.label) 
            align.tip.label.lty <- 3
    }
    if (align.tip.label) {
        if (type %in% c("unrooted", "radial") || !use.edge.length || 
            is.ultrametric(x)) 
            align.tip.label <- FALSE
    }
    if (type %in% c("unrooted", "radial") || !use.edge.length || 
        is.null(x$root.edge) || !x$root.edge) 
        root.edge <- FALSE
    phyloORclado <- type %in% c("phylogram", "cladogram")
    horizontal <- direction %in% c("rightwards", "leftwards")
    xe <- x$edge
    if (phyloORclado) {
        phyOrder <- attr(x, "order")
        if (is.null(phyOrder) || phyOrder != "cladewise") {
            x <- reorder(x)
            if (!identical(x$edge, xe)) {
                ereorder <- match(x$edge[, 2], xe[, 2])
                if (length(edge.color) > 1) {
                  edge.color <- rep(edge.color, length.out = Nedge)
                  edge.color <- edge.color[ereorder]
                }
                if (length(edge.width) > 1) {
                  edge.width <- rep(edge.width, length.out = Nedge)
                  edge.width <- edge.width[ereorder]
                }
                if (length(edge.lty) > 1) {
                  edge.lty <- rep(edge.lty, length.out = Nedge)
                  edge.lty <- edge.lty[ereorder]
                }
            }
        }
        yy <- numeric(Ntip + Nnode)
        TIPS <- x$edge[x$edge[, 2] <= Ntip, 2]
        yy[TIPS] <- 1:Ntip
    }
    z <- reorder(x, order = "postorder")
    if (phyloORclado) {
        if (is.null(node.pos)) 
            node.pos <- if (type == "cladogram" && !use.edge.length) 
                2
            else 1
        if (node.pos == 1) 
            yy <- .nodeHeight(z$edge, Nedge, yy)
        else {
            ans <- .C(node_height_clado, as.integer(Ntip), as.integer(z$edge[, 
                1]), as.integer(z$edge[, 2]), as.integer(Nedge), 
                double(Ntip + Nnode), as.double(yy))
            xx <- ans[[5]] - 1
            yy <- ans[[6]]
        }
        if (!use.edge.length) {
            if (node.pos != 2) 
                xx <- .nodeDepth(Ntip, Nnode, z$edge, Nedge, 
                  node.depth) - 1
            xx <- max(xx) - xx
        } else {
            xx <- .nodeDepthEdgelength(Ntip, Nnode, z$edge, Nedge, 
                z$edge.length)
        }
    } else {
        twopi <- 2 * pi
        rotate.tree <- twopi * rotate.tree/360
        if (type != "unrooted") {
            TIPS <- x$edge[which(x$edge[, 2] <= Ntip), 2]
            xx <- seq(0, twopi * (1 - 1/Ntip) - twopi * open.angle/360, 
                length.out = Ntip)
            theta <- double(Ntip)
            theta[TIPS] <- xx
            theta <- c(theta, numeric(Nnode))
        }
        switch(type, fan = {
            theta <- .nodeHeight(z$edge, Nedge, theta)
            if (use.edge.length) {
                r <- .nodeDepthEdgelength(Ntip, Nnode, z$edge, 
                  Nedge, z$edge.length)
            } else {
                r <- .nodeDepth(Ntip, Nnode, z$edge, Nedge, node.depth)
                r <- 1/r
            }
            theta <- theta + rotate.tree
            if (root.edge) r <- r + x$root.edge
            xx <- r * cos(theta)
            yy <- r * sin(theta)
        }, unrooted = {
            nb.sp <- .nodeDepth(Ntip, Nnode, z$edge, Nedge, node.depth)
            XY <- if (use.edge.length) ape_unrooted_xy(Ntip, Nnode, 
                z$edge, z$edge.length, nb.sp, rotate.tree) else ape_unrooted_xy(Ntip, 
                Nnode, z$edge, rep(1, Nedge), nb.sp, rotate.tree)
            xx <- XY$M[, 1] - min(XY$M[, 1])
            yy <- XY$M[, 2] - min(XY$M[, 2])
        }, radial = {
            r <- .nodeDepth(Ntip, Nnode, z$edge, Nedge, node.depth)
            r[r == 1] <- 0
            r <- 1 - r/Ntip
            theta <- .nodeHeight(z$edge, Nedge, theta) + rotate.tree
            xx <- r * cos(theta)
            yy <- r * sin(theta)
        })
    }
    
    
    
    if (phyloORclado) {
        if (!horizontal) {
            tmp <- yy
            yy <- xx
            xx <- tmp - min(tmp) + 1
        }
        if (root.edge) {
            if (direction == "rightwards") 
                xx <- xx + x$root.edge
            if (direction == "upwards") 
                yy <- yy + x$root.edge
        }
    }
    if (no.margin)
    	{
			# Original behavior in APE 5.0 plot.phylo
			if (cex_original == TRUE)
				{
				par(mai = rep(0, 4))
				}
    	}

    if (show.tip.label) 
        nchar.tip.label <- nchar(x$tip.label)
    max.yy <- max(yy)
    getLimit <- function(x, lab, sin, cex, cex_original=TRUE) {
        if (cex_original == TRUE)
        	{
	        s <- strwidth(lab, "inches", cex = cex)
	        } else {
	        ff <- tempfile()
					png(filename=ff)
					s = strwidth(lab, "inches", cex = cex)
					dev.off()
					unlink(ff)
	        }
        if (any(s > sin)) 
            return(1.5 * max(x))
        Limit <- 0
        while (any(x > Limit)) {
            i <- which.max(x)
            alp <- x[i]/(sin - s[i])
            Limit <- x[i] + alp * s[i]
            x <- x + alp * s
        }
        Limit
    }
    
    
    if (is.null(x.lim)) {
        if (phyloORclado) {
            if (horizontal) {
                xx.tips <- xx[1:Ntip]
                if (show.tip.label) {
									# Original behavior in APE 5.0 plot.phylo
									if (cex_original == TRUE)
										{
	                  pin1 <- par("pin")[1]
	                  } else {
										res = par_invisible(parname="pin")
										pin1 = res[1]
	                  }
                  tmp <- getLimit(xx.tips, x$tip.label, pin1, 
                    cex, cex_original=cex_original)
                  tmp <- tmp + label.offset
                }
                else tmp <- max(xx.tips)
                x.lim <- c(0, tmp)
            }
            else x.lim <- c(1, Ntip)
        } else switch(type, fan = {
            if (show.tip.label) {
                offset <- max(nchar.tip.label * 0.018 * max.yy * 
                  cex)
                x.lim <- range(xx) + c(-offset, offset)
            } else x.lim <- range(xx)
        }, unrooted = {
            if (show.tip.label) {
                offset <- max(nchar.tip.label * 0.018 * max.yy * 
                  cex)
                x.lim <- c(0 - offset, max(xx) + offset)
            } else x.lim <- c(0, max(xx))
        }, radial = {
            if (show.tip.label) {
                offset <- max(nchar.tip.label * 0.03 * cex)
                x.lim <- c(-1 - offset, 1 + offset)
            } else x.lim <- c(-1, 1)
        })
    } else if (length(x.lim) == 1) {
        x.lim <- c(0, x.lim)
        if (phyloORclado && !horizontal) 
            x.lim[1] <- 1
        if (type %in% c("fan", "unrooted") && show.tip.label) 
            x.lim[1] <- -max(nchar.tip.label * 0.018 * max.yy * 
                cex)
        if (type == "radial") 
            x.lim[1] <- if (show.tip.label) 
                -1 - max(nchar.tip.label * 0.03 * cex)
            else -1
    }
    if (phyloORclado && direction == "leftwards") 
        xx <- x.lim[2] - xx
    if (is.null(y.lim)) {
        if (phyloORclado) {
            if (horizontal) 
                y.lim <- c(1, Ntip)
            else {
								# Original behavior in APE 5.0 plot.phylo
								if (cex_original == TRUE)
									{
									pin2 <- par("pin")[2]
									} else {
									res = par_invisible(parname="pin")
									pin2 = res[2]
									}

                yy.tips <- yy[1:Ntip]
                if (show.tip.label) {
                  tmp <- getLimit(yy.tips, x$tip.label, pin2, 
                    cex, cex_original=cex_original)
                  tmp <- tmp + label.offset
                }
                else tmp <- max(yy.tips)
                y.lim <- c(0, tmp)
            }
        } else switch(type, fan = {
            if (show.tip.label) {
                offset <- max(nchar.tip.label * 0.018 * max.yy * 
                  cex)
                y.lim <- c(min(yy) - offset, max.yy + offset)
            } else y.lim <- c(min(yy), max.yy)
        }, unrooted = {
            if (show.tip.label) {
                offset <- max(nchar.tip.label * 0.018 * max.yy * 
                  cex)
                y.lim <- c(0 - offset, max.yy + offset)
            } else y.lim <- c(0, max.yy)
        }, radial = {
            if (show.tip.label) {
                offset <- max(nchar.tip.label * 0.03 * cex)
                y.lim <- c(-1 - offset, 1 + offset)
            } else y.lim <- c(-1, 1)
        })
    } else if (length(y.lim) == 1) {
        y.lim <- c(0, y.lim)
        if (phyloORclado && horizontal) 
            y.lim[1] <- 1
        if (type %in% c("fan", "unrooted") && show.tip.label) 
            y.lim[1] <- -max(nchar.tip.label * 0.018 * max.yy * 
                cex)
        if (type == "radial") 
            y.lim[1] <- if (show.tip.label) 
                -1 - max(nchar.tip.label * 0.018 * max.yy * cex)
            else -1
    }
    if (phyloORclado && direction == "downwards") 
        yy <- y.lim[2] - yy
    if (phyloORclado && root.edge) {
        if (direction == "leftwards") 
            x.lim[2] <- x.lim[2] + x$root.edge
        if (direction == "downwards") 
            y.lim[2] <- y.lim[2] + x$root.edge
    }
    asp <- if (type %in% c("fan", "radial", "unrooted")) 
        1
    else NA
    
    # 2017-11-02_NJM add for plot_phylo3_nodecoords_APE5
    if (plot)
    	{
	    plot.default(0, type = "n", xlim = x.lim, ylim = y.lim, xlab = "", 
        ylab = "", axes = FALSE, asp = asp, ...)
      }
    if (plot) {
        if (is.null(adj)) 
            adj <- if (phyloORclado && direction == "leftwards") 
                1
            else 0
        if (phyloORclado && show.tip.label) {
            MAXSTRING <- max(strwidth(x$tip.label, cex = cex))
            loy <- 0
            if (direction == "rightwards") {
                lox <- label.offset + MAXSTRING * 1.05 * adj
            }
            if (direction == "leftwards") {
                lox <- -label.offset - MAXSTRING * 1.05 * (1 - 
                  adj)
            }
            if (!horizontal) {
                # Original behavior in APE 5.0 plot.phylo
								if (cex_original == TRUE)
									{
									psr <- par("usr")
									} else {
									psr = par_invisible(parname="usr")
									}
                MAXSTRING <- MAXSTRING * 1.09 * (psr[4] - psr[3])/(psr[2] - 
                  psr[1])
                loy <- label.offset + MAXSTRING * 1.05 * adj
                lox <- 0
                srt <- 90 + srt
                if (direction == "downwards") {
                  loy <- -loy
                  srt <- 180 + srt
                }
            }
        }
        if (type == "phylogram") {
            ape_phylogram_plot(x$edge, Ntip, Nnode, xx, yy, horizontal, 
                edge.color, edge.width, edge.lty)
        } else {
            if (type == "fan") {
                ereorder <- match(z$edge[, 2], x$edge[, 2])
                if (length(edge.color) > 1) {
                  edge.color <- rep(edge.color, length.out = Nedge)
                  edge.color <- edge.color[ereorder]
                }
                if (length(edge.width) > 1) {
                  edge.width <- rep(edge.width, length.out = Nedge)
                  edge.width <- edge.width[ereorder]
                }
                if (length(edge.lty) > 1) {
                  edge.lty <- rep(edge.lty, length.out = Nedge)
                  edge.lty <- edge.lty[ereorder]
                }
                ape_circular_plot(z$edge, Ntip, Nnode, xx, yy, theta, 
                  r, edge.color, edge.width, edge.lty)
            }
            else ape_cladogram_plot(x$edge, xx, yy, edge.color, edge.width, 
                edge.lty)
        }
        if (root.edge) {
            rootcol <- if (length(edge.color) == 1) 
                edge.color
            else "black"
            rootw <- if (length(edge.width) == 1) 
                edge.width
            else 1
            rootlty <- if (length(edge.lty) == 1) 
                edge.lty
            else 1
            if (type == "fan") {
                tmp <- ape_polar2rect(x$root.edge, theta[ROOT])
                segments(0, 0, tmp$x, tmp$y, col = rootcol, lwd = rootw, 
                  lty = rootlty)
            }
            else {
                switch(direction, rightwards = segments(0, yy[ROOT], 
                  x$root.edge, yy[ROOT], col = rootcol, lwd = rootw, 
                  lty = rootlty), leftwards = segments(xx[ROOT], 
                  yy[ROOT], xx[ROOT] + x$root.edge, yy[ROOT], 
                  col = rootcol, lwd = rootw, lty = rootlty), 
                  upwards = segments(xx[ROOT], 0, xx[ROOT], x$root.edge, 
                    col = rootcol, lwd = rootw, lty = rootlty), 
                  downwards = segments(xx[ROOT], yy[ROOT], xx[ROOT], 
                    yy[ROOT] + x$root.edge, col = rootcol, lwd = rootw, 
                    lty = rootlty))
            }
        }
        if (show.tip.label) {
            if (is.expression(x$tip.label)) 
                underscore <- TRUE
            if (!underscore) 
                x$tip.label <- gsub("_", " ", x$tip.label)
            if (phyloORclado) {
                if (align.tip.label) {
                  xx.tmp <- switch(direction, rightwards = max(xx[1:Ntip]), 
                    leftwards = min(xx[1:Ntip]), upwards = xx[1:Ntip], 
                    downwards = xx[1:Ntip])
                  yy.tmp <- switch(direction, rightwards = yy[1:Ntip], 
                    leftwards = yy[1:Ntip], upwards = max(yy[1:Ntip]), 
                    downwards = min(yy[1:Ntip]))
                  segments(xx[1:Ntip], yy[1:Ntip], xx.tmp, yy.tmp, 
                    lty = align.tip.label.lty)
                }
                else {
                  xx.tmp <- xx[1:Ntip]
                  yy.tmp <- yy[1:Ntip]
                }
                text(xx.tmp + lox, yy.tmp + loy, x$tip.label, 
                  adj = adj, font = font, srt = srt, cex = cex, 
                  col = tip.color)
            }
            else {
                angle <- if (type == "unrooted") 
                  XY$axe
                else atan2(yy[1:Ntip], xx[1:Ntip])
                lab4ut <- if (is.null(lab4ut)) {
                  if (type == "unrooted") 
                    "horizontal"
                  else "axial"
                }
                else match.arg(lab4ut, c("horizontal", "axial"))
                xx.tips <- xx[1:Ntip]
                yy.tips <- yy[1:Ntip]
                if (label.offset) {
                  xx.tips <- xx.tips + label.offset * cos(angle)
                  yy.tips <- yy.tips + label.offset * sin(angle)
                }
                if (lab4ut == "horizontal") {
                  y.adj <- x.adj <- numeric(Ntip)
                  sel <- abs(angle) > 0.75 * pi
                  x.adj[sel] <- -strwidth(x$tip.label)[sel] * 
                    1.05
                  sel <- abs(angle) > pi/4 & abs(angle) < 0.75 * 
                    pi
                  x.adj[sel] <- -strwidth(x$tip.label)[sel] * 
                    (2 * abs(angle)[sel]/pi - 0.5)
                  sel <- angle > pi/4 & angle < 0.75 * pi
                  y.adj[sel] <- strheight(x$tip.label)[sel]/2
                  sel <- angle < -pi/4 & angle > -0.75 * pi
                  y.adj[sel] <- -strheight(x$tip.label)[sel] * 
                    0.75
                  text(xx.tips + x.adj * cex, yy.tips + y.adj * 
                    cex, x$tip.label, adj = c(adj, 0), font = font, 
                    srt = srt, cex = cex, col = tip.color)
                }
                else {
                  if (align.tip.label) {
                    POL <- ape_rect2polar(xx.tips, yy.tips)
                    POL$r[] <- max(POL$r)
                    REC <- ape_polar2rect(POL$r, POL$angle)
                    xx.tips <- REC$x
                    yy.tips <- REC$y
                    segments(xx[1:Ntip], yy[1:Ntip], xx.tips, 
                      yy.tips, lty = align.tip.label.lty)
                  }
                  if (type == "unrooted") {
                    adj <- abs(angle) > pi/2
                    angle <- angle * 180/pi
                    angle[adj] <- angle[adj] - 180
                    adj <- as.numeric(adj)
                  }
                  else {
                    s <- xx.tips < 0
                    angle <- angle * 180/pi
                    angle[s] <- angle[s] + 180
                    adj <- as.numeric(s)
                  }
                  font <- rep(font, length.out = Ntip)
                  tip.color <- rep(tip.color, length.out = Ntip)
                  cex <- rep(cex, length.out = Ntip)
                  for (i in 1:Ntip) text(xx.tips[i], yy.tips[i], 
                    x$tip.label[i], font = font[i], cex = cex[i], 
                    srt = angle[i], adj = adj[i], col = tip.color[i])
                }
            }
        }
        if (show.node.label) 
            text(xx[ROOT:length(xx)] + label.offset, yy[ROOT:length(yy)], 
                x$node.label, adj = adj, font = font, srt = srt, 
                cex = cex)
    }
    # 2017-11-02_NJM add for plot_phylo3_nodecoords_APE5
    # added: , edge = xe, xx = xx, yy = yy
    L <- list(type = type, use.edge.length = use.edge.length, 
        node.pos = node.pos, node.depth = node.depth, show.tip.label = show.tip.label, 
        show.node.label = show.node.label, font = font, cex = cex, 
        adj = adj, srt = srt, no.margin = no.margin, label.offset = label.offset, 
        x.lim = x.lim, y.lim = y.lim, direction = direction, 
        tip.color = tip.color, Ntip = Ntip, Nnode = Nnode, root.time = x$root.time, 
        align.tip.label = align.tip.label, edge = xe, xx = xx, yy = yy)
    assign("last_plot.phylo", c(L, list(edge = xe, xx = xx, yy = yy)), 
        envir = .PlotPhyloEnv)
    invisible(L)
}





# From: ape:::floating.pie.asp
# Internal function
ape_floating_pie_asp <- function (xpos, ypos, x, edges = 200, radius = 1, col = NULL, 
    startpos = 0, ...) 
{
    
    # Correction in case x, i.e. the ancestral states table, is a data.frame
    x = unname(matrix(x))
    
    u <- par("usr")
    user.asp <- diff(u[3:4])/diff(u[1:2])
    p <- par("pin")
    inches.asp <- p[2]/p[1]
    asp <- user.asp/inches.asp
    if (!is.numeric(x) || any(is.na(x) | x < 0)) 
        stop("floating.pie: x values must be non-negative. You may get this error because 'pievals', i.e. the ancestral state probabilities, is a data.frame instead of a matrix.")
    x <- c(0, cumsum(x)/sum(x))
    dx <- diff(x)
    nx <- length(dx)
    col <- if (is.null(col)) 
        rainbow(nx)
    else rep_len(col, nx)
    if (length(i <- which(dx == 1))) {
        symbols(xpos, ypos, circles = radius, inches = FALSE, 
            add = TRUE, fg = par("fg"), bg = col[i])
    }
    else {
        bc <- 2 * pi * (x[1:nx] + dx/2) + startpos
        for (i in seq_len(nx)) {
            n <- max(2, floor(edges * dx[i]))
            t2p <- 2 * pi * seq(x[i], x[i + 1], length = n) + 
                startpos
            xc <- c(cos(t2p) * radius + xpos, xpos)
            yc <- c(sin(t2p) * radius * asp + ypos, ypos)
            polygon(xc, yc, col = col[i], ...)
        }
    }
}

# From ape:::unrooted.xy
# Internal function
ape_unrooted_xy <- function (Ntip, Nnode, edge, edge.length, nb.sp, rotate.tree) 
{
    foo <- function(node, ANGLE, AXIS) {
        ind <- which(edge[, 1] == node)
        sons <- edge[ind, 2]
        start <- AXIS - ANGLE/2
        for (i in 1:length(sons)) {
            h <- edge.length[ind[i]]
            angle[sons[i]] <<- alpha <- ANGLE * nb.sp[sons[i]]/nb.sp[node]
            axis[sons[i]] <<- beta <- start + alpha/2
            start <- start + alpha
            xx[sons[i]] <<- h * cos(beta) + xx[node]
            yy[sons[i]] <<- h * sin(beta) + yy[node]
        }
        for (i in sons) if (i > Ntip) 
            foo(i, angle[i], axis[i])
    }
    Nedge <- dim(edge)[1]
    yy <- xx <- numeric(Ntip + Nnode)
    axis <- angle <- numeric(Ntip + Nnode)
    foo(Ntip + 1L, 2 * pi, 0 + rotate.tree)
    M <- cbind(xx, yy)
    axe <- axis[1:Ntip]
    axeGTpi <- axe > pi
    axe[axeGTpi] <- axe[axeGTpi] - 2 * pi
    list(M = M, axe = axe)
}


# From ape:::phylogram.plot
# Internal function
ape_phylogram_plot <- function (edge, Ntip, Nnode, xx, yy, horizontal, edge.color, 
    edge.width, edge.lty) 
{
    nodes <- (Ntip + 1):(Ntip + Nnode)
    if (!horizontal) {
        tmp <- yy
        yy <- xx
        xx <- tmp
    }
    x0v <- xx[nodes]
    y0v <- y1v <- numeric(Nnode)
    NodeInEdge1 <- vector("list", Nnode)
    e1 <- edge[, 1]
    for (i in seq_along(e1)) {
        j <- e1[i] - Ntip
        NodeInEdge1[[j]] <- c(NodeInEdge1[[j]], i)
    }
    for (i in 1:Nnode) {
        j <- NodeInEdge1[[i]]
        tmp <- range(yy[edge[j, 2]])
        y0v[i] <- tmp[1]
        y1v[i] <- tmp[2]
    }
    x0h <- xx[edge[, 1]]
    x1h <- xx[edge[, 2]]
    y0h <- yy[edge[, 2]]
    nc <- length(edge.color)
    nw <- length(edge.width)
    nl <- length(edge.lty)
    if (nc + nw + nl == 3) {
        color.v <- edge.color
        width.v <- edge.width
        lty.v <- edge.lty
    }
    else {
        Nedge <- dim(edge)[1]
        edge.color <- rep(edge.color, length.out = Nedge)
        edge.width <- rep(edge.width, length.out = Nedge)
        edge.lty <- rep(edge.lty, length.out = Nedge)
        DF <- data.frame(edge.color, edge.width, edge.lty, stringsAsFactors = FALSE)
        color.v <- rep("black", Nnode)
        width.v <- rep(1, Nnode)
        lty.v <- rep(1, Nnode)
        for (i in 1:Nnode) {
            br <- NodeInEdge1[[i]]
            if (length(br) == 1) {
                A <- br[1]
                color.v[i] <- edge.color[A]
                width.v[i] <- edge.width[A]
                lty.v[i] <- edge.lty[A]
            }
            else if (length(br) > 2) {
                x <- unique(DF[br, 1])
                if (length(x) == 1) 
                  color.v[i] <- x
                x <- unique(DF[br, 2])
                if (length(x) == 1) 
                  width.v[i] <- x
                x <- unique(DF[br, 3])
                if (length(x) == 1) 
                  lty.v[i] <- x
            }
            else {
                A <- br[1]
                B <- br[2]
                if (any(DF[A, ] != DF[B, ])) {
                  color.v[i] <- edge.color[B]
                  width.v[i] <- edge.width[B]
                  lty.v[i] <- edge.lty[B]
                  y0v <- c(y0v, y0v[i])
                  y1v <- c(y1v, yy[i + Ntip])
                  x0v <- c(x0v, x0v[i])
                  color.v <- c(color.v, edge.color[A])
                  width.v <- c(width.v, edge.width[A])
                  lty.v <- c(lty.v, edge.lty[A])
                  y0v[i] <- yy[i + Ntip]
                }
                else {
                  color.v[i] <- edge.color[A]
                  width.v[i] <- edge.width[A]
                  lty.v[i] <- edge.lty[A]
                }
            }
        }
    }
    if (horizontal) {
        segments(x0h, y0h, x1h, y0h, col = edge.color, lwd = edge.width, 
            lty = edge.lty)
        segments(x0v, y0v, x0v, y1v, col = color.v, lwd = width.v, 
            lty = lty.v)
    }
    else {
        segments(y0h, x0h, y0h, x1h, col = edge.color, lwd = edge.width, 
            lty = edge.lty)
        segments(y0v, x0v, y1v, x0v, col = color.v, lwd = width.v, 
            lty = lty.v)
    }
}


# From ape:::circular.plot
# Internal function
ape_circular_plot <- function (edge, Ntip, Nnode, xx, yy, theta, r, edge.color, edge.width, 
    edge.lty) 
{
    r0 <- r[edge[, 1]]
    r1 <- r[edge[, 2]]
    theta0 <- theta[edge[, 2]]
    costheta0 <- cos(theta0)
    sintheta0 <- sin(theta0)
    x0 <- r0 * costheta0
    y0 <- r0 * sintheta0
    x1 <- r1 * costheta0
    y1 <- r1 * sintheta0
    segments(x0, y0, x1, y1, col = edge.color, lwd = edge.width, 
        lty = edge.lty)
    tmp <- which(diff(edge[, 1]) != 0)
    start <- c(1, tmp + 1)
    Nedge <- dim(edge)[1]
    end <- c(tmp, Nedge)
    foo <- function(edge.feat, default) {
        if (length(edge.feat) == 1) 
            return(as.list(rep(edge.feat, Nnode)))
        edge.feat <- rep(edge.feat, length.out = Nedge)
        feat.arc <- as.list(rep(default, Nnode))
        for (k in 1:Nnode) {
            tmp <- edge.feat[start[k]]
            if (tmp == edge.feat[end[k]]) {
                feat.arc[[k]] <- tmp
            }
            else {
                if (nodedegree[k] == 2) 
                  feat.arc[[k]] <- rep(c(tmp, edge.feat[end[k]]), 
                    each = 50)
            }
        }
        feat.arc
    }
    nodedegree <- tabulate(edge[, 1L])[-seq_len(Ntip)]
    co <- foo(edge.color, "black")
    lw <- foo(edge.width, 1)
    ly <- foo(edge.lty, 1)
    for (k in 1:Nnode) {
        i <- start[k]
        j <- end[k]
        X <- rep(r[edge[i, 1]], 100)
        Y <- seq(theta[edge[i, 2]], theta[edge[j, 2]], length.out = 100)
        x <- X * cos(Y)
        y <- X * sin(Y)
        x0 <- x[-100]
        y0 <- y[-100]
        x1 <- x[-1]
        y1 <- y[-1]
        segments(x0, y0, x1, y1, col = co[[k]], lwd = lw[[k]], 
            lty = ly[[k]])
    }
}



# From ape::cladogram.plot
# Internal function
ape_cladogram_plot <- function (edge, xx, yy, edge.color, edge.width, edge.lty) 
	{
	segments(xx[edge[, 1]], yy[edge[, 1]], xx[edge[, 2]], yy[edge[, 2]], col = edge.color, lwd = edge.width, lty = edge.lty)
	}

# From ape:::polar2rect
# Internal function
ape_polar2rect <- function (r, angle) 
	{
	list(x = r * cos(angle), y = r * sin(angle))
	}

# From ape:::rect2polar
# Internal function
ape_rect2polar <- function (x, y)
	{
	list(r = sqrt(x^2 + y^2), angle = atan2(y, x))
	}
nmatzke/BioGeoBEARS documentation built on April 11, 2024, 2:16 a.m.