R/runBSM_v1.R

Defines functions runBSM

# Runs biogeographic stochastic mapping 
# (non-stratified, or stratified, versions, read from the 
#  input res object)
# 
# The 'stochastic_mapping_inputs_list' object comes from:
# 
# Non-stratified:
# 
# 
# stochastic_mapping_inputs_list2 = get_inputs_for_stochastic_mapping_stratified(res=resDIVALIKEj)

# wait_before_save = number of iterations of a meaningless for-loop, to 
# space-out the saving operations to avoid memory faults in R.app

# 
runBSM <- function(res, stochastic_mapping_inputs_list, maxnum_maps_to_try=1, nummaps_goal=maxnum_maps_to_try, maxtries_per_branch=40000, save_after_every_try=TRUE, savedir=getwd(), seedval=as.numeric(Sys.time()), wait_before_save=0.1, master_nodenum_toPrint=0)
	{
	defaults='
	maxnum_maps_to_try=1
	nummaps_goal=maxnum_maps_to_try
	maxtries_per_branch=40000
	save_after_every_try=TRUE
	savedir=getwd()
	seedval=12345
	
	
	res; stochastic_mapping_inputs_list=stochastic_mapping_inputs_list; maxnum_maps_to_try=100; nummaps_goal=1; maxtries_per_branch=40000; save_after_every_try=TRUE; savedir=getwd(); seedval=12345; wait_before_save=0.01; master_nodenum_toPrint=0
	'
	
	
	# Necessary default setting to avoid numbers instead of states
	options(stringsAsFactors=FALSE)
	
	# Starter
	clado_events_tables = list()
	ana_events_tables = list()
	
	# 2016-05-05_bug fix
	lnum = 0
	
	# Reconcile # of maps to try with nummaps goal
	if (nummaps_goal > maxnum_maps_to_try)
		{
		maxnum_maps_to_try = 2 * nummaps_goal
		} # END if (nummaps_goal > maxnum_maps_to_try)
	
	# Set time-stratification
	stratified = FALSE
	if ((is.numeric(res$inputs$timeperiods))) #&& (length(inputs$timeperiods) > 1))
		{
		stratified = TRUE
		strat_TF = TRUE
		} else {
		stratified = FALSE
		strat_TF = FALSE		
		}

	

	# Loop through the number of attempts
	m=1
	for (m in 1:maxnum_maps_to_try)
		{
		# Initialize
		stochastic_mapping_results = NA
		clado_events_table = NA
		ana_events_table = NA
		
		# Skip loop once you get to nummaps_goal
		# (the desired number of successes)
		if (lnum >= nummaps_goal)
			{
			next()
			}
		
		startseed = seedval + m

		#print(paste0("runBSM startseed = ", startseed))
		
		# Run stratified or non-stratified stochastic mapping
		if (strat_TF == FALSE)
			{
			cmdstr = paste0("stochastic_mapping_results = stochastic_map_given_inputs(stochastic_mapping_inputs=stochastic_mapping_inputs_list, piecenum=NULL, maxtries=", maxtries_per_branch, ", seedval=", startseed, ", include_null_range=res$inputs$include_null_range, master_nodenum_toPrint=master_nodenum_toPrint, allow_null_tips=res$inputs$allow_null_tips)")
			} else {
			# Stratified
			#maxtries=40000; seedval=521; master_nodenum_toPrint=0
			
			# stochastic_mapping_results = stochastic_mapping_on_stratified(res=res, stochastic_mapping_inputs_list=stochastic_mapping_inputs_list, maxtries=40000, seedval=12346, master_nodenum_toPrint=master_nodenum_toPrint)
			
			# stochastic_mapping_inputs_list=stochastic_mapping_inputs_list; maxtries=40000; seedval=12346; master_nodenum_toPrint=9
			
			cmdstr = paste0("stochastic_mapping_results = stochastic_mapping_on_stratified(res=res, stochastic_mapping_inputs_list=stochastic_mapping_inputs_list, maxtries=", maxtries_per_branch, ", seedval=", startseed, ", master_nodenum_toPrint=master_nodenum_toPrint)")
			} # END if (strat_TF == FALSE)
		#print("AA_stochastic_mapping_on_stratified")
		#print(cmdstr)
		
		# Run either of the options
		#cat(" ", m, sep="")
		#print("try_result:")
		# maxtries=40000; seedval=12346; master_nodenum_toPrint=master_nodenum_toPrint
		# stochastic_mapping_on_stratified
		try_result = try(expr=eval(parse(text=cmdstr)))
		#eval(parse(text=cmdstr))
		#print(try_result)
		#cat("\n")
		#print(warnings())
		#print(class(try_result))
		#cat("\n")
		#cat(",", sep="")
		
		# 2016-05-05
		# Now that stochastic_mapping_results is initialized to NA each loop, should
		# have no weird blank results
		#TF2 = ( (length(stochastic_mapping_results)==1) && (is.na(stochastic_mapping_results) == TRUE) )
		success_TF = ( class(stochastic_mapping_results) == "data.frame" )
		
		#if ( (class(try_result) != "try-error") && TF2==FALSE )
		if ( (class(try_result) == "try-error") && (success_TF == FALSE) )
			{
			# FAILURE!
			cat("\nFailure on stochastic mapping attempt #", m, "/", maxnum_maps_to_try, " tries. Holding at success #", lnum, "/", nummaps_goal, ".", sep="")
			cat("\n")
			print(try_result)
			# END if (class(try_result) != "try-error")
			} else {
			# SUCCESS!
			# Non-stratified:
			if (strat_TF == FALSE)
				{
				##########################################
				# For NON-stratified stochastic mapping
				##########################################
				clado_events_table = stochastic_mapping_results
				
				#print(clado_events_table[20,])
				
				# Error check, for no anagenetic events
				TF1 = stochastic_mapping_results$anagenetic_events_txt_below_node == "none"
				TF2 = is.na(stochastic_mapping_results$anagenetic_events_txt_below_node)
				TF3 = stochastic_mapping_results$anagenetic_events_txt_below_node == "NA"
				TF4 = stochastic_mapping_results$anagenetic_events_txt_below_node == ""
				TFall = (TF1 + TF2 + TF3 + TF4) 
				TF = TFall == 0
				
				if (sum(TF, na.rm=TRUE) > 0)
					{
					# Some anagenetic events, make table
					#print("here1")
					ana_events_table = events_txt_list_into_events_table(events_txt_list=clado_events_table$anagenetic_events_txt_below_node, trtable=clado_events_table, recalc_abs_ages=TRUE)
					#print("here2")
					} else {
					# No anagenetic events
					ana_events_table = NA
					}
				} else {
				####################################daughter_nds
				# For STRATIFIED stochastic mapping
				####################################
				#clado_colnames = c("node", "node.type", "edge.length", "time_bp", "stratum", "piecenum", "piececlass", "SUBnode", "SUBnode.type", "SUBedge.length", "clado_event_type", "clado_event_txt", "clado_dispersal_to", "sampled_states_AT_nodes", "sampled_states_AT_brbots", "left_desc_nodes", "right_desc_nodes", "samp_LEFT_dcorner", "samp_RIGHT_dcorner", "anagenetic_events_txt_below_node")
				clado_colnames = c("node", "node.type", "parent_br", "edge.length", "ancestor", "daughter_nds", "time_bp", "fossils", "label", "stratum", "time_top", "time_bot", "reltimept", "piecenum", "piececlass", "SUBnode", "SUBnode.type", "SUBparent_br", "SUBedge.length", "SUBancestor", "SUBdaughter_nds" , "SUBtime_bp", "SUBfossils", "SUBlabel", "clado_event_type", "clado_event_txt", "clado_dispersal_to", "sampled_states_AT_nodes", "sampled_states_AT_brbots", "left_desc_nodes", "right_desc_nodes", "samp_LEFT_dcorner", "samp_RIGHT_dcorner", "anagenetic_events_txt_below_node")
				#rowsTF = stochastic_mapping_results$master_table_cladogenetic_events$node.type != "tip"
				clado_events_table = stochastic_mapping_results$master_table_cladogenetic_events[, clado_colnames]
				ana_events_table = stochastic_mapping_results$table_w_anagenetic_events
				} # END if (strat_TF == FALSE)
			
			# Print update
			lnum = lnum+1
			cat("\nSuccess #", lnum, "/", nummaps_goal, " on stochastic mapping attempt #", m, "/", maxnum_maps_to_try, " tries.", sep="")
			
			# Wait a second or two, every 100 successes
			# -- attempting to avoid segmentation faults during
			# rapid stochastic mapping
			if (FALSE)
				{
				Sys.sleep(wait_before_save)
				}
			
			if (TRUE)
				{
				if (lnum/50 == round(lnum/50))
					{
					Sys.sleep(wait_before_save * 3)
					} # END if (lnum/100 == round(lnum/100))
				}
			
			# Store results
			clado_events_tables[[lnum]] = clado_events_table
			# Error trap -- if no anagenetic events, put NA
			# 2016-05-05 bug: is.na(ana_events_table) is flagging yes when whole 
			# row is NAs
			#if (suppressWarnings((length(ana_events_table) > 1) && (is.na(ana_events_table)==TRUE)))
# 			if (check_for_ana_events_table(ana_events_table) == TRUE)
# 				{
# 				error_msg = "STOP ERROR in runBSM(): your anagenetic events table seems to have row(s) that have multiple NAs. This could be due to an out of date 'BioGeoBEARS_stochastic_mapping_v1.R' file, missing this bug fix: '# NJM 2016-05-05 bug fix: add 'as.numeric'   rownums_in_trtable = as.numeric(tmptable_rows$nodenum_at_top_of_branch)'"
# 				cat("\n\n")
# 				cat(error_msg)
# 				cat("\n\n")
# 				stop(error_msg)
# 				}
# 			
			
			if ( check_for_ana_events_table(ana_events_table) == TRUE )
				{
				ana_events_tables[[lnum]] = ana_events_table
				} else {
				cat("\nNOTE: 'ana_events_table' for this BSM was null, NA or length <1 or >1; this suggests no anagenetic events in this BSM. Replacing with 'NA'.")
				ana_events_table = NA
				ana_events_tables[[lnum]] = ana_events_table
				}
			
			
			if (save_after_every_try == TRUE)
				{
				# Sometimes it looks like saving happens too fast and R crashes
				if (is.null(wait_before_save) == FALSE)
					{
					# Wait this long before saving
					Sys.sleep(wait_before_save)
					} # END if (is.null(wait_before_save) == FALSE)
				
				# Archive every round
				RES_clado_events_tables = clado_events_tables
				RES_ana_events_tables = ana_events_tables

				RES_clado_events_tables_fn = slashslash(paste0(savedir, "/", "RES_clado_events_tables_PARTIAL.Rdata"))
				RES_ana_events_tables_fn = slashslash(paste0(savedir, "/", "RES_ana_events_tables_PARTIAL.Rdata"))
				
				# Message about saving
				cat("...saving to RES_clado_events_tables_PARTIAL.Rdata")
				cmdstr = "save(RES_clado_events_tables, file=RES_clado_events_tables_fn)"
				savetry = try(expr=cmdstr)
				if (class(savetry) == "try-error")
					{
					cat("...save FAILED for 'RES_clado_events_tables_PARTIAL.Rdata' for some reason, sometimes rapid saves seem to crash on some systems.\n")
					} else {
					cat("...saved 'RES_clado_events_tables_PARTIAL.Rdata'\n")
					}
				#print(RES_ana_events_tables)
				
				cmdstr = "save(RES_ana_events_tables, file=RES_ana_events_tables_fn)"
				savetry = try(expr=eval(parse(text=cmdstr)))
				if (class(savetry) == "try-error")
					{
					cat("...save FAILED for 'RES_ana_events_tables_PARTIAL.Rdata' for some reason, sometimes rapid saves seem to crash on some systems.\n")
					} else {
					cat("...saved 'RES_ana_events_tables_PARTIAL.Rdata'")
					}
				} # END if (save_after_every_try == TRUE)
			} # END if ( (class(try_result) == "try-error") && (success_TF == FALSE) )
		} # END for (m in 1:maxnum_maps_to_try)
	
	# Archive at the end
	RES_clado_events_tables = clado_events_tables
	RES_ana_events_tables = ana_events_tables
	RES_clado_events_tables_fn = slashslash(paste0(savedir, "/", "RES_clado_events_tables.Rdata"))
	RES_ana_events_tables_fn = slashslash(paste0(savedir, "/", "RES_ana_events_tables.Rdata"))
	
	# Print results
	cat("\n\n\n=============================================\nBiogeographic Stochastic Mapping is finished.\n\n")
	txt = paste0("Saving cladogenetic and anagenetic events tables. ")
	cat(txt)
	txt2 = paste0("The commands:\n\nload(file=", RES_clado_events_tables_fn, ")\n\nload(file=", RES_ana_events_tables_fn, ")\n\n...will load them to objects RES_clado_events_tables and RES_ana_events_tables.")
	cat(txt2)
	cat("\n\n")
		
	save(RES_clado_events_tables, file=RES_clado_events_tables_fn)
	save(RES_ana_events_tables, file=RES_ana_events_tables_fn)
	
	cat("Now returning to console 'BSM_output', which contains:\n")
	cat("BSM_output$RES_clado_events_tables")
	cat("\n")
	cat("BSM_output$RES_ana_events_tables")
	cat("\n=============================================")
	cat("\n")
	cat("\n")
	
	BSM_output = NULL
	BSM_output$RES_clado_events_tables = RES_clado_events_tables
	BSM_output$RES_ana_events_tables = RES_ana_events_tables
	
	extract='
	clado_events_tables = BSM_output$RES_clado_events_tables
	ana_events_tables = BSM_output$RES_ana_events_tables
	'
	
	return(BSM_output)
	} # END runBSM()
nmatzke/BioGeoBEARS documentation built on July 1, 2024, 12:45 p.m.