R/computations.R

In gonzalezgouveia/CosmoBetti: Compute Betti Curves in R

# betti curves computation


#' Compute the Betti number
#'
#' Compute the Betti number of a persistence diagram
#' for a desired dimension and radius.
#' The Betti number of dimension \code{d} and radius \code{r}
#' of a diagram \code{diag}
#' refers to the number of \code{d}-dimensional holes of
#' a simplicial complex at radius \code{r} of
#' a \code{diag} generated by a point cloud.
#' When \code{d = 0} it should be interpreted as the
#' number of conected components of the simplicial complex
#' at radius \code{r}.
#'
#' @param diag is a persistence diagram generated by a point cloud
#' @param dim is the dimension of the homology group
#' @param radius is the radius for which the betti number is computed
#' @return the betti number (birth - death)
#' @export
compute_betti_number <- function(diag, dim, radius){
  # alphaComplexDiag comes with Inf in first row
  if(diag[1,3] == Inf){
    filtered_diag <- diag[which(diag[,1]==dim), 2:3][-1,]
  } else {
    filtered_diag <- diag[which(diag[,1]==dim), 2:3]
  }
  n_birth <- sum(filtered_diag[,1] <= radius)
  n_death <- sum(filtered_diag[,2] <= radius)
  return(n_birth-n_death)
}

#' Compute Betti curves
#'
#' Compute Betti curves of a persistence diagram
#' until a desired dimension, a maximum radius, and output vector lenght.
#' The Betti curve of dimension \code{d} of a diagram \code{diag}
#' refers to the sequence of Betti numbers of dimension \code{d} and
#' diamgram \code{diag} along radius from zero to infinity (maximun radius).
#'
#' if \code{maxdim = -1} then the dimension is infered from diag
#' if \code{maxradius = -1} the maximum radius is infered from diag
#'
#' @param diag is a persistence diagram from TDA o TDAstats
#' @param maxdim is the maximum dimension the betti curve is computed
#' @param maxradius is the maximum radius the betti curve is computed
#' @param len is the length of each betticurve
#' @return A matrix with the Betti curves
#' @examples
#' library(TDAstats)
#' X <- cbind(runif(100), runif(100))
#' diag <- calculate_homology(X)
#' compute_betti_curves(diag)
#' @export
compute_betti_curves <- function(diag, maxdim = -1, maxradius=-1, len = 100){
  # check if diag is a matrix
  # PENDIENTE: HACER FUNCION PARA CHECAR QUE ES UN DIAG VALIDO
  if (is.matrix(diag)) {
  } else if (is.list(diag)) {
    diag <- diag$diagram
  } else {
    print('diag not a valid persistence matrix')
    break
  }

  # assign maxradius automatically, info in description
  if (maxradius == -1){
    maxradius <- max(diag[-1,2:3]) * 1.05 # maximum plus 5%
  }

  # assign dimension automatically, info in description
  if (maxdim == -1){
    maxdim <- max(diag[,1]) # starts at 0
  }

  # create the x dimension for the bettimatrix
  radius_seq <- seq((maxradius/(2*len)), maxradius, length.out = len)

  bettimatrix <- matrix(NA, ncol = 3, nrow = 0)
  colnames(bettimatrix) <- c('dim', 'radius', 'value')

  for(dim in 0:maxdim){
    for(r in radius_seq){
      bettivalue <- compute_betti_number(diag, dim, r)
      bettimatrix <- rbind(bettimatrix,
                           c(dim, r, bettivalue))
    }
  }
  bettimatrix <- tibble::as_tibble(bettimatrix)
  bettimatrix$dim <- factor(bettimatrix$dim)
  return(bettimatrix)
}