## ----setup, include=FALSE----------------------------------------------------- library(knitr) opts_chunk$set( fig.align = "center", out.extra = 'style="max-width:100%; overflow-x:auto; white-space: nowrap;"', tidy = FALSE ) knitr::opts_chunk$set( collapse = TRUE, comment = "#>", message = FALSE, warning = FALSE ) ## ----------------------------------------------------------------------------- # load package library(ELAplus) # load test data set # species abundance table data("baseabtable") head(baseabtable) # metadata (environmental factors) data("basemetadata") head(basemetadata) ## ----------------------------------------------------------------------------- # load package library(ELAplus) # Generate random parameters hb_params <- hb.paramgen(24, ne = 2) # the number of environmental factors h.act <- hb_params[[1]] j.act <- hb_params[[2]] g.act <- hb_params[[3]] # Simulate community composition data rand_data <- HeatBath(100, 512, h.act, j.act, g = g.act) rand_mat <- rand_data[[1]] rand_enmat <- rand_data[[2]] ## ----------------------------------------------------------------------------- formatted <- Formatting(baseabtable = baseabtable, basemetadata = basemetadata) ## ----eval=FALSE--------------------------------------------------------------- # formatted <- Formatting(baseabtable = baseabtable, basemetadata = basemetadata, # normalize = 1, parameters = c(0.01, 0.01, 0.99)) # #> Processed 256 samples. # #> Relative abundance threshold = 0.01 # #> Occurrence threshold (lower) = 0.01 # #> Occurrence threshold (upper) = 0.99 # #> Selected 16 out of 16 species. # ocvecs <- formatted[[1]] # community composition binary matrix # abvecs <- formatted[[2]] # community abundance matrix (not binarized) # envecs <- formatted[[3]] # environmental variables (optional) ## ----echo=FALSE, out.width="80%"---------------------------------------------- knitr::include_graphics("figures/ParameterMatrix.png") ## ----echo=FALSE, out.width="95%"---------------------------------------------- knitr::include_graphics("figures/EPMEM_abstract2.png") ## ----eval=FALSE--------------------------------------------------------------- # bpresult <- Findbp(ocvecs, rep = 2, threads = 2, # cvmode="10fold",nrepeat = 10, we=c(0.001),totalit=4000, # lmd=c(1e-05,1e-03,0.01),maxlrs = c(0.005), # tol=0.03, runadamW=TRUE,sparse=TRUE,fastfitting=TRUE) # #Start hyper-parameter search: # #There are 300 tasks # #Finished 300 tasks # #SA: elapsed time 15.08 sec # # bp <- bpresult[[1]] # typically best bpresult # bpm <- as.numeric(unlist(strsplit(names(bp)[1], split = "_"))) # allresults <- bpresult[[2]] ## ----plotSAtest, eval=FALSE--------------------------------------------------- # plotSAtest(allresults) # # "id" corresponds to hyperparameter: lmd_we_maxlr ## ----echo=FALSE, out.width="95%"---------------------------------------------- knitr::include_graphics("figures/plotSAtest.png") ## ----eval=FALSE--------------------------------------------------------------- # sa <- runSA( # ocvecs, # enmat = NULL, # rep = 32, # getall = FALSE, # lambda = bpm[[1]], # we = bpm[[2]], # maxlr = bpm[[3]], # totalit = bpm[[4]] # ) # #> Start parameter fitting: # #> .SA: elapsed time 0.13 sec ## ----eval=FALSE--------------------------------------------------------------- # ela <- ELA( # sa, # env = NULL, # threads = 1, # SS.itr = 20000, # FindingTip.itr = 10000, # reporting = TRUE # ) # #> Start ELA: # #> 7 stable states were found. # #> Checking 21 tipping points. # #> converting... # #> ELA: elapsed time 2.19 sec ## ----eval=FALSE--------------------------------------------------------------- # ela[[1]] # stable state IDs # #> [1] "09x" "01t" "EWB" "EY8" "1uV" "1mL" "17s" # ela[[2]] # stable state energies # #> [1] -11.099418 -10.297931 -9.544050 -8.961634 -8.005040 -6.440385 -6.269140 ## ----eval=FALSE--------------------------------------------------------------- # elap <- ELPruning(ela, th=0.1, type="relative") # #> Start pruning: # #> *... # #> ELPruning: elapsed time 0.05 sec ## ----eval=FALSE--------------------------------------------------------------- # elap[[1]][[1]] # updated stable state IDs # #> [1] "09x" "EWB" "EY8" "1uV" # elap[[1]][[2]] # updated stable state energies # #> [1] -11.099418 -9.544050 -8.961634 -8.005040 ## ----eval=FALSE--------------------------------------------------------------- # elap[[2]] # #> ss.before.pruning ss.after.pruning # #> 1 09x 09x # #> 2 01t 09x # #> 3 EWB EWB # #> 4 EY8 EY8 # #> 5 1uV 1uV # #> 6 1mL 09x # #> 7 17s 09x ## ----echo=FALSE, out.width="95%"---------------------------------------------- knitr::include_graphics("figures/deltaE_scheme.png") ## ----eval=FALSE--------------------------------------------------------------- # elap <- ELPruning(ela, th=0.2, type="relative") # #> Start pruning: # #> *... # #> ELPruning: elapsed time 0.05 sec # elap[[1]][[1]] # updated stable state IDs # #> [1] "09x" "EWB" "1uV" ## ----showDG, eval=FALSE------------------------------------------------------- # # Before pruning # showDG(ela, ocvecs, "example") ## ----echo=FALSE, out.width="95%"---------------------------------------------- knitr::include_graphics("figures/DGexample.png") ## ----showDG_pruned, eval=FALSE------------------------------------------------ # # After pruning # showDG(elap[[1]], ocvecs, "example") ## ----echo=FALSE, out.width="95%"---------------------------------------------- knitr::include_graphics("figures/DGexample_pruned.png") ## ----echo=FALSE, out.width="95%"---------------------------------------------- knitr::include_graphics("figures/deltaE_hist.png") ## ----showDG_pvalue, eval=FALSE------------------------------------------------ # bs_params <- bootstrap_SA(ocvecs, rep = 16, threads=8,　bootstrap.it = 1000, # lambda = bpm[[1]],we = bpm[[2]], # maxlr = bpm[[3]], totalit = bpm[[4]]) # # bootstrap <- bootstrap_ELA(bs_params,ocvecs,ela=ela,threads = 1) # ela_rep <- bootstrap[[1]] # bootstrap_ene <- bootstrap[[2]] # showDG(ela_rep,ocvecs,bootstrap_ene = bootstrap_ene, min_pval = 0.01, # scale_labels = c(0.05,0.1,0.2,0.5)) ## ----echo=FALSE, out.width="95%"---------------------------------------------- knitr::include_graphics("figures/DGexample_pval.png") ## ----showDG_IQR_pruned, eval=FALSE-------------------------------------------- # elap_rep <- ELPruning(ela, type = "bootstrap",bootstrap_ene=bootstrap_ene,lower_qr = 0.25) # showDG(ela_rep,ocvecs,bootstrap_ene = bootstrap_ene, min_pval = 0.01, # scale_labels = c(0.05,0.1,0.2,0.5)) ## ----echo=FALSE, out.width="95%"---------------------------------------------- knitr::include_graphics("figures/DGexample_pval_pruned.png") ## ----showDG_manual, eval=FALSE------------------------------------------------ # grobj_to_plot <- showDG(ela_rep,ocvecs,bootstrap_ene = bootstrap_ene,getGraphObj = TRUE) # grobj_point <- grobj_to_plot[grobj_to_plot$point_str != "",] # g <- ggplot( # grobj_to_plot, # aes(x=.data$nodes2xposi, y=.data$energy, label=.data$point_str)) + # xlab("") + geom_text(hjust=0.75, vjust=2, aes(fontface=2)) + geom_path() + # geom_point( # data = grobj_point, # mapping = aes( # x = .data$nodes2xposi, # y = .data$energy, # color = .data$energy, # shape = .data$type, # ),size = 4 # ) + # scale_color_viridis_c(option = "plasma",alpha=0.8, # direction=-1,na.value = "grey") + # coord_cartesian(xlim=c(0.75,7.5), ylim=c(-5.8,-11.5)) # plot(g) ## ----echo=FALSE, out.width="95%"---------------------------------------------- knitr::include_graphics("figures/DG_manual.png") ## ----ShowIntrGraph, eval=FALSE------------------------------------------------ # showIntrGraph( # elap[[1]], # sa, # th = 0.01, # annot_adj = c(0.75, 1.00) # ) ## ----echo=FALSE, out.width="95%"---------------------------------------------- knitr::include_graphics("figures/IntrGraph.png") ## ----eval=FALSE--------------------------------------------------------------- # bpresult <- Findbp(ocvecs, enmat = envecs, rep = 2, threads = 2, # cvmode="10fold",nrepeat = 10, we=c(0.001),totalit=4000, # lmd=c(1e-05,1e-03,0.01),maxlrs = c(0.005), # tol=0.03, runadamW=TRUE,sparse=TRUE,fastfitting=TRUE) # bp <- bpresult[[1]] # bpm <- as.numeric(unlist(strsplit(names(bp)[1], split = "_"))) # # sa <- runSA( # ocvecs, # enmat = envecs, # rep = 16, # getall = FALSE, # lambda = bpm[[1]], # we = bpm[[2]], # maxlr = bpm[[3]], # totalit = bpm[[4]] # ) # ## ----eval=FALSE--------------------------------------------------------------- # # This process takes few minutes with a small number of threads, and here pre-computed results are shown. # gela <- GradELA( # sa = sa, # eid = "factor.1", # enmat = envecs, # env = NULL, # range = c(-1, 1), # steps = 32, # th = 0.05, # threads = 8, # pruning_type = "relative", # bs_params = bs_params, # ocmat = ocvecs) ## ----eval=FALSE--------------------------------------------------------------- # # This process takes few minutes with a small number of threads, and here pre-computed results are shown. # bs_params <- bootstrap_SA(ocvecs, enmat = envecs,rep = 16, threads=8,bootstrap.it = 1000, # lambda = bpm[[1]],we = bpm[[2]], # maxlr = bpm[[3]], totalit = bpm[[4]]) # # gela <- GradELA( # sa = sa, # eid = "factor.1", # enmat = envecs, # env = NULL, # range = c(-1, 1), # steps = 32, # lower_qr = 0.25, # threads = 8, # pruning_type = "bootstrap", # bs_params = bs_params, # ocmat = ocvecs) ## ----eval=FALSE--------------------------------------------------------------- # gela[[1]][[1]][[1]][[1]] # stable state IDs at 1st environmental parameter # # [1] "0Ht" "EWB" "1uV" # gela[[1]][[20]][[1]][[1]] # stable state IDs at 20th environmental parameter # # [1] "09x" "EWB" ## ----showSSD, eval=FALSE, fig.width=8, fig.height=6, message=FALSE------------ # # GradELA result (ELPruning with type="relative"). # showSSD(gela) ## ----echo=FALSE, out.width="95%"---------------------------------------------- knitr::include_graphics("figures/SSD.png") ## ----showSSD_bootstrap, eval=FALSE, fig.width=8, fig.height=6, message=FALSE---- # # GradELA result (ELPruning with type="bootstrap"). # showSSD(gela) ## ----echo=FALSE, out.width="95%"---------------------------------------------- knitr::include_graphics("figures/SSD_pruned.png") ## ----showSSD_ggplot, eval=FALSE, fig.width=8, fig.height=6, message=FALSE----- # ssd_df <- showSSD_ggplot(gela,plot = FALSE,getSSDobj = TRUE) # g <- ggplot(ssd_df, aes(x = env, y = Energy, color = SS_ID)) + # geom_point(shape = 19) + # geom_line(aes(group = SS_ID)) # ## ----echo=FALSE, out.width="95%"---------------------------------------------- knitr::include_graphics("figures/SSD_ggplot.png") ## ----GELA3D-demo, eval=FALSE, fig.width=8, fig.height=6, message=FALSE-------- # # This process takes few minutes, and here pre-computed result is drawn. # gelaobj <- GELAObj(gela, sa=sa,threads=2) # showGELA3D(gelaobj,theta = 30,phi = 30) ## ----echo=FALSE, out.width="95%"---------------------------------------------- knitr::include_graphics("figures/GELA3D.png")