moved specific code from package.R to files

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: divent
 Type: Package
 Title: Entropy Partitioning to Measure Diversity
-Version: 0.4-99.9013
+Version: 0.4-99.9014
 Authors@R: c(
     person("Eric", "Marcon", email="eric.marcon@agroparistech.fr", role=c("aut", "cre"), comment=c(ORCID="0000-0002-5249-321X"))
     )

NEWS.md

@@ -1,4 +1,4 @@
-# divent 0.4-99.9013
+# divent 0.4-99.9014
 
 ## Features
 

R/accum_sp.R

@@ -333,3 +333,84 @@ plot_map <- function(
   # Return the image for further processing
   return(invisible(the_image))
 }
+
+
+#' Helper to prepare parameters for `accum_sp` plot and autoplot.
+#'
+#' @param x the `accum_sp` object to plot.
+#' @param q the order of diversity.
+#' @param main the title of the plot.
+#' @param xlab X-axis label.
+#' @param ylab Y-axis label.
+#' @param ylim Y-axis limits
+#'
+#' @returns a vector of parameters for the plots
+#' @noRd
+#'
+accum_sp_plot_helper <- function(x, q, main, xlab, ylab, ylim) {
+
+  # Find the row in the accumulation table
+  q_row <- which(dimnames(x$accumulation)$q == q)
+  if (length(q_row) != 1) {
+    stop("The value of q does not correspond to any accumulation curve.")
+  }
+
+  if (is.null(ylim)) {
+    # Evaluate ylim if not set by an argument
+    ymin <- min(x$accumulation[q_row, , ])
+    ymax <- max(x$accumulation[q_row, , ])
+  } else {
+    ymin <- ylim[1]
+    ymax <- ylim[2]
+  }
+
+  if (main == "Accumulation of ...") {
+    # Prepare the main title
+    if (inherits(x, "accum_sp_entropy")) {
+      main <- paste("Accumulation of Entropy of order", q)
+    }
+    if (inherits(x, "accum_sp_diversity")) {
+      if (q == 0) {
+        main <- "Species Accumulation Curve"
+      } else {
+        main <- paste("Accumulation of Diversity of order", q)
+      }
+    }
+    if (inherits(x, "accum_sp_mixing")) main <- paste("Mixing index of order", q)
+  }
+
+  if (xlab == "Sample size...") {
+    if (names(dimnames(x$accumulation)[2]) == "n") {
+      xlab <- "Number of individuals"
+    } else {
+      xlab <- "Distance from the central individual"
+    }
+  }
+
+  if (ylab == "Diversity...") {
+    # Prepare Y-axis
+    if (inherits(x, "accum_sp_entropy")) {
+      ylab <- "Diversity"
+    }
+    if (inherits(x, "accum_sp_diversity")) {
+      if (q == 0) {
+        ylab <- "Richness"
+      } else {
+        ylab <- "Diversity"
+      }
+      if (inherits(x, "accum_sp_mixing")) {
+        ylab <- "Mixing index"
+      }
+    }
+    return(
+      list(
+        q_row = q_row, 
+        ymin = ymin, 
+        ymax = ymax, 
+        main = main, 
+        xlab = xlab, 
+        ylab = ylab
+      )
+    )
+  }
+}

R/div_hurlbert.R

@@ -117,3 +117,35 @@ div_hurlbert.species_distribution <- function(
 
   return(the_diversity)
 }
+
+
+#' Compute Hurlbert's diversity from its entropy
+#' 
+#' Find the effective number of species numerically
+#'
+#' @param hurlbert_entropy The entropy.
+#' @param k The order of entropy.
+#'
+#' @returns Hurlbert's effective number of species.
+#' @noRd
+#' 
+hurlbert_ent2div <- function(hurlbert_entropy, k) {
+  # Relation between diversity and entropy
+  # (D for diversity, S for entropy, k is the parameter)
+  f <- function(D, S, k) {D * (1 - (1 - 1 / D)^k) - S}
+  # Minimize it
+  return(
+    vapply(
+      hurlbert_entropy, 
+      FUN = function(S) {
+        stats::uniroot(
+          f = f, 
+          interval = c(1, 1E+7), 
+          S = S, 
+          k = k
+        )$root
+      },
+      FUN.VALUE = 0
+    )
+  )
+}

R/ent_phylo.R

@@ -205,3 +205,109 @@ ent_phylo.species_distribution <- function(
     }
   }
 }
+
+
+#' Phylogenetic entropies
+#' 
+#' Calculate entropies of a list of phylogenetic abundances (obtained by 
+#' [phylo_abd]).
+#' Each item of the list corresponds to a phylogenetic group, i.e. an interval
+#' of the tree (where the species do not change).
+#' 
+#' @param phylo_abd A list of matrices of abundance (caution: lines are species,
+#' columns are communities).
+#'
+#' @returns A vector. Each item is the entropy of a community.
+#'  
+#' @noRd
+#'
+phylo_entropy.phylo_abd <- function(
+    # Allow lapply along q
+  q,
+  phylo_abd,
+  tree,
+  normalize,
+  # Other arguments for ent_tsallis.numeric
+  estimator,
+  level, 
+  probability_estimator,
+  unveiling,
+  richness_estimator,
+  jack_alpha, 
+  jack_max,
+  coverage_estimator,
+  gamma) {
+
+  if (gamma) {
+    # Calculate gamma entropy of each group.
+    # simplify2array() makes a vector with the list of numbers.
+    phylo_entropies <- simplify2array(
+      lapply(
+        # Calculate entropy in each item of the list, i.e. group.
+        # Obtain a list.
+        phylo_abd,
+        FUN = function(group) {
+          ent_tsallis.species_distribution(
+            as_abundances(t(group)),
+            q = q,
+            estimator = estimator,
+            level = level, 
+            probability_estimator = probability_estimator,
+            unveiling = unveiling,
+            richness_estimator = richness_estimator,
+            jack_alpha  = jack_alpha, 
+            jack_max = jack_max,
+            coverage_estimator = coverage_estimator,
+            gamma = TRUE,
+            # Obtain a vector.
+            as_numeric = TRUE,
+            check_arguments = FALSE
+          )
+        }
+      )
+    )
+
+  } else {
+    # Calculate entropy of each community in each group.
+    # simplify2array() makes a matrix with the list of vectors.
+    phylo_entropies <- simplify2array(
+      lapply(
+        # Calculate entropy in each item of the list, i.e. group.
+        # Obtain a list.
+        phylo_abd,
+        FUN = function(group) {
+          apply(
+            group,
+            # Calculate entropy of each column of the matrix, i.e. community.
+            MARGIN = 2,
+            FUN = ent_tsallis.numeric,
+            # Arguments
+            q = q,
+            estimator = estimator,
+            level = level, 
+            probability_estimator = probability_estimator,
+            unveiling = unveiling,
+            richness_estimator = richness_estimator,
+            jack_alpha  = jack_alpha, 
+            jack_max = jack_max,
+            coverage_estimator = coverage_estimator,
+            # Obtain a vector.
+            as_numeric = TRUE,
+            check_arguments = FALSE
+          )
+        }
+      )
+    )
+  }
+  # Should be a matrix, but simplify2array() makes a vector instead of a 1-col 
+  # matrix and gamma entropy is a vector. Force a matrix.
+  if (is.vector(phylo_entropies)) {
+    phylo_entropies <- t(phylo_entropies)
+  }
+
+  # Calculate the weighted mean of entropy and normalize
+  the_entropy <- as.numeric(tree$intervals %*% t(phylo_entropies))
+  if (normalize) the_entropy <- the_entropy / sum(tree$intervals)
+
+  return(the_entropy)
+}

R/ent_similarity.R

@@ -442,3 +442,81 @@ S_v <- function(
   v_used <- seq_len(sample_size - abd[species_index])
   return(sum(w_v[v_used] * p_V_Ns[v_used, species_index]))
 }
+
+
+#' Similarity-Based Gamma entropy of a metacommunity
+#' 
+#' Build the metacommunity and check that abundances are integers.
+#' 
+#' See [ent_gamma_tsallis] for details.
+#' @param species_distribution An object of class [species_distribution].
+#' 
+#' @returns A tibble with the estimator used and the estimated entropy.
+#' @noRd
+#' 
+ent_gamma_similarity <- function(
+    species_distribution,
+    similarities,
+    q,
+    estimator,
+    probability_estimator,
+    unveiling,
+    jack_alpha,
+    jack_max,
+    coverage_estimator,
+    as_numeric) {
+
+  # Build the metacommunity
+  abd <- metacommunity.abundances(
+    species_distribution, 
+    as_numeric = TRUE, 
+    check_arguments = FALSE
+  )
+  if (is_integer_values(abd)) {
+    # Sample coverage is useless
+    sample_coverage <- NULL
+  } else {
+    # Non-integer values in the metacommunity. 
+    # Calculate the sample coverage and change the estimator.
+    sample_coverage <- coverage.numeric(
+      colSums(
+        species_distribution[
+          , !colnames(species_distribution) %in% non_species_columns
+        ]
+      ),
+      estimator = coverage_estimator,
+      as_numeric = TRUE,
+      check_arguments = FALSE
+    )
+    if (!estimator %in% c("Marcon", "ChaoShen")) {
+      estimator <- "Marcon"
+    }
+  }
+
+  # Compute the entropy.
+  the_entropy <- ent_similarity.numeric(
+    abd,
+    similarities = similarities,
+    q = q,
+    estimator = estimator,
+    probability_estimator = probability_estimator,
+    unveiling = unveiling,
+    jack_alpha  = jack_alpha,
+    jack_max = jack_max,
+    sample_coverage = sample_coverage,
+    as_numeric = as_numeric,
+    check_arguments = FALSE
+  )
+
+  # Return
+  if (as_numeric) {
+    return(the_entropy)
+  } else {
+    return(
+      dplyr::bind_cols(
+        site = "Metacommunity",
+        the_entropy
+      )
+    )
+  }
+}

R/ent_sp_simpson.R

@@ -189,3 +189,23 @@ ent_sp_simpsonEnvelope <- function(
   # Return the envelope
   return(the_envelope)
 }
+
+
+#' Extract a column from an fv object
+#' according to an edge-effect correction
+#'
+#' @param fv the function value object, see [spatstat.explore::fv.object].
+#' @param correction the edge-effect correction: 
+#' "isotropic", "translate" or "none"
+#'
+#' @returns a vector with the function values
+#' @noRd
+#'
+correction_fv <- function(fv, correction) {
+  switch(
+    correction,
+    "isotropic" = fv$iso,
+    "translate" = fv$trans,
+    "none" = fv$un
+  )
+}