R/gels.R

In ELAplus: Energy Landscape Analysis for Ecological Communities

#' hamming_distance
#' @description Compute Hamming distance between two equal-length vectors.
#' @param x Vectors.
#' @param y Vectors.
#' @return Integer count of positions where x and y differ.
#' @useDynLib ELAplus, .registration=TRUE
#' @noRd
hamming_distance <- function(x, y) {
  sum(x != y)
}

#' gen_concompo
#' @description Merge overlapping index groups into connected components.
#' @param mm A list of vectors representing groups of indices.
#' @return A list of merged groups (connected components).
#' @useDynLib ELAplus, .registration=TRUE
#' @noRd
gen_concompo <- function(mm){
  concompo <- list()
  lapply(mm, function(mmm){
    matching_indices <- sapply(seq_along(concompo), function(i) {
      any(sapply(mmm, function(x) x %in% unlist(concompo[[i]])))
    })
    if (!any(matching_indices)) {
      concompo <- c(concompo, list(mmm))
    } else {
      matching_elements <- concompo[matching_indices]
      concompo <- c(concompo[!matching_indices], list(union(unlist(matching_elements), mmm)))
    }})
  return(concompo)
}
  
#' con_split
#' @description Split \code{data} into consecutive runs defined by \code{factor_variable}.
#' @param data Object to split
#' @param factor_variable Vector defining run-length groups.
#' @return A list of grouped subsets of \code{data} in run order.
#' @useDynLib ELAplus, .registration=TRUE
#' @noRd
con_split <- function(data, factor_variable){group_flag <- rle(factor_variable)$lengths
                        grouped_data <- split(data, rep(seq_along(group_flag), group_flag))
                        return(grouped_data)}

#' ssextend
#' @description subfunction for GELAObj
#' @param spli List of entries.
#' @param de Numeric/vector of values
#' @return A data.frame with columns (de, group_id, stb).
#' @useDynLib ELAplus, .registration=TRUE
#' @noRd
ssextend <- function(spli, de) {
  aaa <- sapply(spli, function(x) x[[1]])
  stb <- ifelse(de %in% aaa, 1, 0)
  data.frame(de, rep(spli[[1]][[3]], length(de)), stb)
}

#' replace_bin
#' @description Flip 0/1 values of \code{vec} at given indices.
#' @param vec Numeric/binary vector
#' @param index Integer indices to flip.
#' @return The modified vector with selected bits flipped.
#' @useDynLib ELAplus, .registration=TRUE
#' @noRd
replace_bin <- function(vec, index) {
  vec[index] <- abs(vec[index] - 1)
  return(vec)
}

#' foldList
#' @description Iteratively apply \code{func} over \code{vec}, returning all intermediate states.
#' @param func Binary function
#' @param init Initial value
#' @param vec Vector of inputs to fold over.
#' @return A list of accumulated results including the initial state.
#' @useDynLib ELAplus, .registration=TRUE
#' @noRd
foldList <- function(func, init, vec) {
  result <- list(init)
  for (i in seq_along(vec)) {
    result <- c(result, list(func(result[[length(result)]], vec[i])))
  }
  return(result)
}

#' MinTippingPath
#' @description Compute a minimum “tipping” energy path between two states.
#' @param s1 Another integer state ID
#' @param s2 Another integer state ID
#' @param hg explicit/implicit environmental factors
#' @param j species interaction matrix
#' @param tmax Iteration limit.
#' @return Numeric vector of energies along the selected path (or scalar if s1==s2).
#' @useDynLib ELAplus, .registration=TRUE
#' @export
MinTippingPath <- function(s1, s2, hg, j, tmax){
    n <- length(hg)
    ss1 <- id2bin(s1, n)
    ss2 <- id2bin(s2, n)
    dif <- ss2 - ss1
    pd <- which(abs(dif)==1)
    seq <- sample(pd)
    #############

    if(s1 == s2){return(Energy(ss1, hg, j))}
    if(length(seq) <= 7){
      seq_perms <- gtools::permutations(length(seq), length(seq), seq)

      pathes <- apply(seq_perms, 1, function(x) {
          samplepath <- foldList(replace_bin, ss1, x)
          sapply(samplepath, function(y) Energy(y, hg, j))
          })

      tote <- apply(pathes, 2, function(x) {
          dife <-  x[-1] - x[-length(x)]
          sum(dife[which(dife > 0)])})
      pe <- pathes[,which(tote==min(tote))]
      if(is.null(ncol(pe))){minpath <- pe}else{minpath <- rowMeans(pe)}
    }
    else{
    tt <- 0
    pree <- Inf
    mintote <- Inf

    while (tt < tmax) {
      tt <- tt + 1
      rpl <- sample(length(seq), 2) 
      seqnew <- seq 
      seqnew[rpl] <- seqnew[c(rpl[2], rpl[1])]
      samplepath <- foldList(replace_bin, ss1, seqnew)
      pathe <- sapply(samplepath, function(y) Energy(y, hg, j))
      dife <- diff(pathe)
      tote <- sum(dife[dife > 0])
      if (tote < mintote) {
        minpath <- pathe
        mintote <- tote
        pree <- tote
      }
      if (runif(1) < min(1, exp(pree) / exp(tote))) {
        seq <- seqnew
        pree <- tote
      }
    }}
    return(minpath)}


#' MarginePath
#' @description sub function of GELAObj
#' @param s1 Integer state ID
#' @param l Path length
#' @param hg explicit/implicit environmental factors
#' @param j species interaction matrix
#' @return Numeric vector giving the mean energy profile along the sampled paths.
#' @useDynLib ELAplus, .registration=TRUE
#' @noRd
MarginePath <- function(s1, l, hg, j) {
  n <- length(hg)
  ss1 <- id2bin(s1, n)
  
  pathes <- replicate(min(1024, 2^l), {
    samplepath <- replicate(l, {
        red <-  Reduce(function(vec, i) replace_bin(vec, sample(length(vec), 1)), seq_len(l), init = ss1, accumulate = TRUE)
        sapply(red, function(x) Energy(x, hg, j))
    })
    samplepath
  })
  
  mean_path <- apply(pathes, 1, mean)
  return(mean_path)}

#' GELAObj
#' @description Construct smoothed energy landscape surface objects from GELA output and parameters.
#' @param gela GELA object/list
#' @param sa Parameter list generated by \link{runSA}
#' @param threads Number of parallel workers.
#' @return A list containing surface grid data and state-specific surface points.
#' @useDynLib ELAplus, .registration=TRUE
#' @importFrom utils combn
#' @importFrom parallel makeCluster stopCluster clusterCall
#' @importFrom stats smooth.spline predict median embed
#' @importFrom dplyr bind_rows
#' @importFrom purrr map
#' @importFrom foreach foreach %do% %dopar%
#' @importFrom doParallel registerDoParallel
#' @importFrom magrittr %>%
#' @export
GELAObj <- function(gela, sa, threads = 1){
  de <- gela[[2]]
  s <- length(sa[1][[1]][,1])
  refenv <- gela[[3]]
  #derange <- range((sa[[2]])[,is.na(refenv)])
  #de0 <- seq(derange[[1]], derange[[2]], length.out = length(de))

  #elsurf <- foreach(x=gela[[1]], .packages=c("tidyverse","ELAplus","gtools","igraph")) %dopar% {GraphObj(x[[1]])}
  elsurf <- foreach(x=gela[[1]]) %do% {GraphObj(x[[1]])}
  ssp <- lapply(elsurf, function(x) {
    x[which(sapply(x$ccc_str, function(y) length(eval(parse(text = y))) == 1)), ]
  })
  sssp <- foreach(x=ssp, y=de) %do% {
      mat <- matrix(c(rep(y, nrow(x)), x[,3], x[,5]), nrow(x))
      split(mat, seq_len(nrow(mat)))}
  sssp <- unlist(sssp, recursive = FALSE)
  sslist <- unique(unlist(lapply(ssp, function(df) df[, 5])))
  ###
  if(length(sslist)>1){
    ss_combi <- combn(sslist, 2)
    #neibor states in assembly graph
    mm <- apply(ss_combi, 2, function(pair) {
      dist <- hamming_distance(id2bin(pair[1], s), id2bin(pair[2], s))
      if (dist == 1) return(pair)})
    mm <- Filter(function(x) !is.null(x), mm)
  }else{
    mm <- NULL
  }
  ###
  concompo <- gen_concompo(mm)
  grp <- c(concompo, sapply(setdiff(sslist, unlist(concompo)), function(x) list(x)))
  ssspa <- lapply(grp, function(gg) {
    lapply(sssp, function(ss) {
      ss[grepl(paste(gg, collapse = "|"), ss[3])]
    })})
  ssspa <- lapply(ssspa, function(x) x[sapply(x, function(y) length(y) > 0)]) #Remove empty elements
  if(length(grp)==1){
    ssspa <- lapply(ssspa, function(x) x[order(as.numeric(sapply(x, '[[', 1)))]) #Sorting within each element
  }else{
    ssspa <- sapply(ssspa, function(x) x[order(as.numeric(sapply(x, '[[', 1)))]) #Sorting within each element
  }
  ssspa <- ssspa[order(sapply(ssspa, length))] #Sort by length
  ssspa <- ssspa[order(sapply(ssspa, function(x){as.numeric(x[[1]][[1]])}))] #Sort by starting point
  ssspa_unlist <- unlist(ssspa, recursive = FALSE)
  spl <- con_split(ssspa_unlist, sapply(ssspa_unlist, '[[', 3))
  tssspastb <- lapply(spl, function(x) ssextend(x, de))
  ssid <- sapply(tssspastb, function(x){x[[2]][[1]]})
  if(length(ssid) > 1){
    hh <- apply(matrix(stats::embed(ssid, 2)[, 2:1],ncol=2),
                1, function(x) hamming_distance(id2bin(x[1], s), id2bin(x[2], s)))
    ml <- max(7, max(hh))
    
    posit <- Reduce(function(x, y) x + y, hh, accumulate = TRUE, init = ml + 1)
  }else{
    ml <- 7
    posit <- 8
  }
  
  cluster <- parallel::makeCluster(threads, type = if (.Platform$OS.type == "windows") "PSOCK" else "FORK")
  parallel::clusterCall(cluster, function(x) .libPaths(x), .libPaths())
  parallel::clusterSetRNGStream(cluster)
  doParallel::registerDoParallel(cluster)
  on.exit(parallel::stopCluster(cluster), add = TRUE)

  filledsurf <- foreach(i=seq(length(de)), .packages=c("ELAplus", "tidyverse", "purrr", "foreach", "gtools"),
                        .export = c("tssspastb","refenv","sa","posit","ml", "MinTippingPath", "MarginePath", "foldList", "replace_bin")) %dopar% {
      v <- lapply(tssspastb, function(x) x[i,])
      v <- do.call(rbind, lapply(v, function(x) as.matrix(x[,])))
      xde <- as.numeric(v[1,1])
      basenv <- refenv
      basenv[is.na(basenv)] <- xde
      basenv <- t(basenv)
      sa2p <- sa2params(sa, c(basenv))
      he <- sa2p[[1]]
      je <- sa2p[[2]]
      ge <- sa2p[[3]]
      hge <- sa2p[[4]]
      ###
      pstb <- as.numeric(v[,3])
      sid <- ssid[as.numeric(v[,3])==1]
      zx <- posit[as.numeric(v[,3])==1]
      ###
      scape <- unlist(map(seq_along(sid)[-length(sid)], function(j) {
          MinTippingPath(sid[j], sid[j + 1], hge, je, 10000)}))
      if(length(scape) == 0){
          scaledscp <- data.frame(z = zx,  energy = MinTippingPath(sid[[1]], sid[[1]], hge, je, 10000))}
      else{
          scaledscp <- data.frame(z = seq(zx[1], zx[length(zx)], length.out = length(scape) - 1), energy = scape[-length(scape)])}
      leftmargine <- MarginePath(ssid[1], zx[1] - 1, hge, je)
      leftmargine <- data.frame(
        z = seq_along(leftmargine[-1]),
        energy = rev(leftmargine[-1]))
      rightmargine <- MarginePath(ssid[length(ssid)], ml + (posit[length(posit)] - zx[length(zx)]), hge, je)
      rightmargine <- data.frame(
        z = seq(zx[length(zx)] + 1, posit[length(posit)] + 1 + ml),
        energy = rightmargine)
      rbind(leftmargine, scaledscp, rightmargine)
  }

  cc <- 0
  intpol <- lapply(filledsurf, function(df) {
    cc <<- cc + 1
    spline_fit <- smooth.spline(df$z, df$energy, spar = 0.2)  # Set adjustable smoothing parameter spar
    pol <- predict(spline_fit, seq(df$z[1], df$z[length(df$z)], length.out = 100))$y
    dx <- seq(1, df$z[length(df$z)], length.out = 100)
    dex <- rep(de[cc], length(pol))#dex <- rep(de0[cc], length(pol))
    data.frame(dx = dx, dex = dex, pol = pol)
  })
  surf3d <- lapply(intpol, function(df) {
    spline_fit <- smooth.spline(df$dx, df$pol, spar = 0.2)
    smoothed_pol <- predict(spline_fit, df$dx)$y
    filtered_pol <- stats::filter(smoothed_pol, filter = rep(1/5, 5), sides = 2, circular = TRUE)
    data.frame(x = df$dx, y = df$dex, z = filtered_pol)
  })
  surf3dx <- dplyr::bind_rows(surf3d)

  l3d <- lapply(tssspastb, function(y){
    si <- which(y$stb==1)
    xe <- posit[which(ssid==y[1,2])]
    xval <- surf3d[[1]][,1]
    nearest_x <- which.min(abs(xval - xe))
    data.frame(t(sapply(si, function(x) surf3d[[x]][nearest_x,])))})
  l3dx <- dplyr::bind_rows(lapply(seq_along(l3d), function(i) {
    df <- l3d[[i]]
    df$Index <- i
    df$ssid <- sapply(tssspastb,function(x){unique(x[,2])})[i]
    df <- df[,c(4,5,1,2,3)]
    return(df)}))
  return(list(surf3dx, l3dx))
}