## ---- echo=FALSE--------------------------------------------------------------
knitr::opts_chunk$set(
  cache = FALSE,
  fig.width = 9,
  message = FALSE,
  eval=FALSE,  
  warning = FALSE)

## ----install-bioc, eval=FALSE-------------------------------------------------
#  if (!requireNamespace("BiocManager", quietly = TRUE))
#      install.packages("BiocManager")

## ----load, eval=TRUE----------------------------------------------------------
library(miaSim)

## ----deps---------------------------------------------------------------------
#  library(ggplot2)
#  library(miaSim)
#  library(vegan)
#  library(ggplot2)
#  library(reshape2)

## ----seed---------------------------------------------------------------------
#  set.seed(42)
#  # n_rep <- 50 # Full (and slow) example
#  n_rep <- 5 # Speed up demo example

## ----initialize---------------------------------------------------------------
#  result_df <- data.frame(
#      n_species = integer(),
#      theta = numeric(),
#      i = integer(),
#      n_resources = integer(),
#      value = numeric()
#  )
#  result_df2 <- data.frame(
#    matrix(NA, nrow = 0, ncol = 13, dimnames = list(NULL, paste0("sp", seq_len(13))))
#  )
#  sorensen_df <- data.frame(
#      n_species = integer(),
#      theta = numeric(),
#      rho_mean = numeric(),
#      rho_sd = numeric()
#  )

## ----funcs--------------------------------------------------------------------
#  gradient_df_generator <- function(n_row, n_col, density_row, max_gradient, error_interval){
#    list_initial <- list()
#    dissimilarity.gradient <- seq(from = 0, to = max_gradient, length.out = n_row)
#    for (i in seq_len(n_row)){
#      print(i)
#      if (i == 1){
#        row_temp <- rbeta(n_col, 1, n_col)
#        col_to_remove <- sample(x = seq_len(n_col), size = n_col-n_col*density_row)
#        row_temp[col_to_remove] <- 0
#        list_initial[[i]] <- row_temp
#      } else {
#        while (length(list_initial) < i) {
#          row_temp <- rbeta(n_col, 1, n_col)
#          col_to_remove <- sample(x = seq_len(n_col), size = n_col-n_col*density_row)
#          row_temp[col_to_remove] <- 0
#          diff_temp <- abs(vegdist(rbind(list_initial[[1]], row_temp), method = "bray") - dissimilarity.gradient[i])
#          if (diff_temp < error_interval) {
#            list_initial[[i]] <- row_temp
#          }
#        }
#      }
#    }
#  
#    # Return
#    as.data.frame(t(matrix(unlist(list_initial), ncol = n_row)))
#  
#  }

## ----initialize2--------------------------------------------------------------
#  # Full example
#  #n_species_types <- c(13, 3, 4)
#  #theta_types <- c(1, 0.75, 0.5, 0.25, 0.1, 0.05)
#  #n_resources_types <- c(1,2,4,8,16,32)
#  
#  # Speeding up the example by reducing the number of combinations
#  n_species_types <- c(3, 4)
#  theta_types <- c(1, 0.5)
#  n_resources_types <- c(1,2)
#  
#  community.initial.df <- as.list(
#    lapply(n_species_types,
#           gradient_df_generator,
#           n_row = n_rep,
#           density_row = 1,
#           max_gradient = 0.7,
#           error_interval = 0.15)
#  )

## ----loop---------------------------------------------------------------------
#  n_species_types <- c(3, 4)
#  for (n_species in n_species_types){
#      for (theta in theta_types) {
#          sorensen <- c()
#          for (i in seq_len(n_rep)){
#              for (n_resources in n_resources_types) {
#                  ### generate E ####
#                  Etest <- randomE(n_species = n_species,
#                                   n_resources = n_resources,
#                                   mean_consumption = theta*n_resources,
#                                   exact = TRUE)
#  
#                  ### calculate rho ####
#                  if (n_resources == max(n_resources_types)){
#                      Etest_pos <- Etest
#                      Etest_pos[Etest_pos<0] <- 0
#                      for (j in seq_len(n_species - 1)){
#                          for (k in 2:n_species){
#                              sorensen <- c(sorensen,
#                                            sum(apply(Etest_pos[c(j,k),], 2, min)))
#                          }
#                      }
#                  }
#  
#                  if (n_resources > 1){
#                      Priority <- t(apply(matrix(sample(n_species * n_resources), nrow = n_species), 1, order))
#                      Priority <- (Etest > 0) * Priority
#                  } else {
#                      Priority <- NULL
#                  }
#  
#                  print(paste0("n_species=",n_species, " theta=",theta, " i=", i, " n_resources=", n_resources))
#                  x0temp <- as.numeric(community.initial.df[[match(n_species, n_species_types)]][i,])
#                  x0temp <- 10*x0temp/sum(x0temp)
#                  CRMtest <- simulateConsumerResource(n_species = n_species,
#                                                      n_resources = n_resources,
#                                                      x0 = x0temp, #rep(10, n_species),
#                                                      resources = rep(100, n_resources),
#                                                      E = Etest,
#                                                      # trophic_priority = Priority,
#                                                      stochastic = TRUE,
#                                                      t_end = 1000,
#                                                      t_step = 1,
#                                                      t_store = 1000)
#                  CRMspecies <- assay(CRMtest, "counts")[, ncol(CRMtest)] # last time point; was getCommunity(CRMtest)
#                  CRMspeciesTotal <- sum(CRMspecies)
#                  result_df[nrow(result_df)+1,] <- c(n_species, theta, i, n_resources, CRMspeciesTotal)
#                  result_df2[nrow(result_df2)+1,] <- c(CRMspecies, rep(NA, 13-length(CRMspecies)))
#  
#              }
#          }
#          rho_mean <- mean(sorensen)
#          rho_sd <- var(sorensen)
#          sorensen_df[nrow(sorensen_df)+1, ] <- c(n_species, theta, rho_mean, rho_sd)
#      }
#  }
#  

## ----visualize----------------------------------------------------------------
#  n_species_types <- c(3, 4)
#  p <- ggplot(result_df, aes(x = n_resources, y = value, group = n_resources)) +
#      geom_boxplot(outlier.shape = NA) +
#      geom_jitter(alpha = 0.2, width = 0.2) +
#      facet_grid(. ~ factor(n_species, levels = n_species_types)) +
#      theme_bw() +
#      scale_x_continuous(trans = "log2", breaks = n_resources_types) +
#      labs(x = "number of resources",
#           y = "growth yield (number of individuals)")
#  print(p)

