R/Renyi.R

In rasterdiv: Diversity Indices for Numerical Matrices

Renyi <- function(x, window=3, alpha=1, base=exp(1), rasterOut=TRUE, np=1, na.tolerance=1, cluster.type="SOCK", debugging=FALSE){

  # Initial checks
  if( !((is(x,"matrix") | is(x,"SpatialGridDataFrame") | is(x,"RasterLayer") | is(x,"list"))) ) {
    stop("\nNot a valid x object. Exiting...")
  }
  else if( is(x,"matrix") ) {
    rasterm <- x
  }
  else if( is(x,"SpatialGridDataFrame") ) {
    rasterm <- raster(x)
  }
  else if( is(x,"RasterLayer")) {
    rasterm <- matrix(getValues(x), ncol = ncol(x), nrow = nrow(x), byrow=TRUE)
  } 
  else if( is(x,"list") ) {
    message("x is a list, only first element will be taken.")
    if( !((is(x[[1]],"matrix") | is(x[[1]],"SpatialGridDataFrame") | is(x[[1]],"RasterLayer"))) ) {
      stop("The first element of list x is not a valid object. Exiting...")
    }
    rasterm<-x[[1]]
    if( is(rasterm,"RasterLayer") ) {
      rasterm <- matrix(getValues(rasterm), ncol = ncol(rasterm), nrow = nrow(rasterm), byrow=TRUE)
    }
  }
  if ( any(!is.numeric(alpha)) ){
    stop("alpha must be a numeric vector. Exiting...")
  }
  if ( any(alpha<0) ){
    stop("alpha must be only positive numbers. Exiting...")
  }

  # Assign mode according to the length of alpha
  if( length(alpha)==1 ) mode <- "single" else if( length(alpha)>1 ) mode <- "iterative"
  # One single alpha
  if (mode=="single"){
    # If alpha is ~1 then Shannon is calculated
    if ( abs(alpha-1)<.Machine$double.eps ) {
      Shannon <- TRUE
    } else{
      Shannon <- FALSE
    }
    # If alpha approaches positive infinity then Berger-Parker is calculated
    if ( alpha >= .Machine$integer.max ) {
      BergerParker <- TRUE
    } else {
      BergerParker <- FALSE
    }
  }
  # Print messages about output
  message(paste(c("\nObject x check OK: \nRenyi with alpha parameter value in ", alpha," will be returned."),collapse=" "))
  # Derive operational moving window
  if( window%%2==1 ){
    w <- (window-1)/2
  } else {
    stop("The size of the moving window must be an odd number. Exiting...")
  }

  # If one single process
  if (np == 1){
    if(mode == "single") {
      if( BergerParker ) {
        outS <- log(1/BergerParkerS(rasterm, w, na.tolerance, debugging))
      }
      else if( Shannon ) {
        outS <- ShannonS(rasterm, w, na.tolerance, debugging)
      }
      else{
        outS <- RenyiS(rasterm, w, alpha, base, na.tolerance, debugging)
      }
      message(("\nCalculation complete.\n"))
      if(rasterOut==T & class(x)[1]=="RasterLayer") {
        outR <- lapply(list(outS),raster, template=x)
        return(outR)
      }else{
        return(outS)
      }
    }
    else if(mode == "iterative"){
      outS <- list()
      for (ALPHA in alpha){
        message("\nProcessing alpha ",ALPHA,"\n")
        if((abs(ALPHA-1)<.Machine$double.eps)) {
          s <- "Shannon_Renyi_alpha_1"
          outS[[s]] <- ShannonS(rasterm, w, na.tolerance, debugging)
        }
        else if (ALPHA >= .Machine$integer.max) {
          s <- "Berger-Parker"
          outS[[s]] <- log(1/BergerParkerS(rasterm, w, na.tolerance, debugging))
        }
        else{
          s <- paste("Renyi_alpha_",as.character(ALPHA),sep="")
          outS[[s]] <- RenyiS(rasterm, w, ALPHA, base, na.tolerance, debugging)
        }
      }
      message(("\nCalculation complete.\n"))
      if(rasterOut==T & class(x)[1]=="RasterLayer") {
        outR <- lapply(outS, raster, template=x)
        return(outR)
      }else{
        return(outS)
      }
    }
  }
  else if (np>1){

    # If more than 1 process
    message("\n##################### Starting parallel calculation #######################")     
    # Opening the cluster
    if(debugging){cat("#check: Before parallel function.")}
    if( cluster.type=="SOCK" || cluster.type=="FORK" ) {
      cls <- makeCluster(np,type=cluster.type, outfile="",useXDR=FALSE,methods=FALSE,output="")
    } 
    else if( cluster.type=="MPI" ) {
      cls <- makeCluster(np,outfile="",useXDR=FALSE,methods=FALSE,output="")
    } 
    else {
      message("Wrong definition for cluster.type. Exiting...")
    }
    registerDoParallel(cls)
    # Close clusters on exit
    on.exit(stopCluster(cls))
    # Garbage collection
    gc()
    if(mode == "single") {
      if(BergerParker){
        outP <- 1/log(BergerParkerP(rasterm, w, na.tolerance, debugging))
      }
      if( Shannon ) {
        outP <- ShannonP(rasterm, w, na.tolerance, debugging)
      }
      else{
        outP <- RenyiP(rasterm, w, alpha, base, na.tolerance, debugging)
      }
      if(rasterOut==T & class(x)[1]=="RasterLayer") {
        outR <- lapply(list(do.call(cbind,outP)),raster, template=x)
        return(outR)
      }else{
        return(outP)
      }
    }
    else if(mode == "iterative"){
      outP <- list()
      for (ALPHA in alpha){
        if( (abs(ALPHA-1)<.Machine$double.eps) ) {
          s <- "Shannon_Renyi_alpha_1"
          out <- ShannonP(rasterm, w, na.tolerance, debugging)
          outP[[s]] <- do.call(cbind,out)
        } 
        else if ( ALPHA >= .Machine$integer.max ){
          s <- "Berger-Parker"
          out <- BergerParkerP(rasterm, w, na.tolerance, debugging)
          outP[[s]] <- 1/log(do.call(cbind,out))
        }
        else{
          s <- paste("Renyi_alpha_",as.character(ALPHA),sep="")
          out <- RenyiP(rasterm, w, ALPHA, base, na.tolerance, debugging)
          outP[[s]] <- do.call(cbind,out)
        }
      }
      if(rasterOut==T & class(x)[1]=="RasterLayer") {
        outR <- lapply(outP,raster, template=x)
        return(outR)
      }else{
        return(outP)
      }
      message("\nCalculation complete.\n")
    }
  }
}