exactextractr coverage - 89.19%

Files
Source

# Copyright (c) 2018-2020 ISciences, LLC.
# All rights reserved.
#
# This software is licensed under the Apache License, Version 2.0 (the "License").
# You may not use this file except in compliance with the License. You may
# obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

if (!isGeneric("exact_extract")) {
	setGeneric("exact_extract", function(x, y, ...)
		standardGeneric("exact_extract"))
}

#' Extract or summarize values from Raster* objects
#'
#' Extracts the values of cells in a Raster* that are covered by a
#' simple feature collection containing polygonal geometries, as well as the
#' fraction of each cell that is covered by the polygon. Returns either
#' the result of a summary operation or function applied to the values
#' and coverage fractions (if \code{fun} is specified), or a data frame
#' containing the values and coverage fractions themselves (if \code{fun}
#' is \code{NULL}.)
#'
#' The value of \code{fun} may be set to a string (or vector of strings)
#' representing summary operations supported by the exactextract library.
#' If the input raster has a single layer and a single summary operation
#' is specified, \code{exact_extract} will return a vector with the result
#' of the summary operation for each feature in the input. If the input
#' raster has multiple layers, or if multiple summary operations are specified,
#' \code{exact_extract} will return a data frame with a row for each feature
#' and a column for each summary operation / layer combination. (The
#' \code{force_df} can be used to always return a data frame instead of a vector.)
#' In all of the summary operations, \code{NA} values in the raster are ignored
#' (i.e., \code{na.rm = TRUE}.)
#'
#' The following summary operations are supported:
#'
#' \itemize{
#'  \item{\code{min} - the minimum defined value in any raster cell wholly or
#'                     partially covered by the polygon}
#'  \item{\code{max} - the maximum defined value in any raster cell wholly or
#'                     partially covered by the polygon}
#'  \item{\code{count} - the sum of fractions of raster cells with defined values
#'                       covered by the polygon}
#'  \item{\code{sum}   - the sum of defined raster cell values, multiplied by
#'                       the fraction of the cell that is covered by the polygon}
#'  \item{\code{mean} - the mean cell value, weighted by the fraction of each cell
#'                      that is covered by the polygon}
#'  \item{\code{median} - the median cell value, weighted by the fraction of each
#'                        cell that is covered by the polygon}
#'  \item{\code{quantile} - arbitrary quantile(s) of cell values, specified in
#'                          \code{quantiles}, weighted by the fraction of each
#'                          cell that is covered by the polygon}
#'  \item{\code{mode} - the most common cell value, weighted by the fraction of
#'                      each cell that is covered by the polygon. Where multiple
#'                      values occupy the same maximum number of weighted cells,
#'                      the largest value will be returned.}
#'  \item{\code{majority} - synonym for \code{mode}}
#'  \item{\code{minority} - the least common cell value, weighted by the fraction
#'                          of each cell that is covered by the polygon. Where
#'                          multiple values occupy the same minimum number of
#'                          weighted cells, the smallest value will be returned.}
#'  \item{\code{variety} - the number of distinct values in cells that are wholly
#'                         or partially covered by the polygon.}
#'  \item{\code{variance} - the population variance of cell values, weighted by the
#'                          fraction of each cell that is covered by the polygon.}
#'  \item{\code{stdev} - the population standard deviation of cell values, weighted
#'                       by the fraction of each cell that is covered by the polygon.}
#'  \item{\code{coefficient_of_variation} - the population coefficient of variation of
#'                       cell values, weighted by the fraction of each cell that is
#'                       covered by the polygon.}
#'  \item{\code{weighted_mean} - the mean cell value, weighted by the product of
#'                               the fraction of each cell covered by the polygon
#'                               and the value of a second weighting raster provided
#'                               as \code{weights}}
#'  \item{\code{weighted_sum} - the sum of defined raster cell values, multiplied by
#'                              the fraction of each cell that is covered by the polygon
#'                               and the value of a second weighting raster provided
#'                               as \code{weights}}
#' }
#'
#' Alternatively, an R function may be provided as \code{fun}. The function will be
#' called for each feature with with vectors of cell values and weights as arguments.
#' \code{exact_extract} will then return a vector of the return values of \code{fun}.
#'
#' If \code{fun} is not specified, \code{exact_extract} will return a list with
#' one data frame for each feature in the input feature collection. The data
#' frame will contain a column with values from each layer in the input `Raster*`,
#' and a final column indicating the fraction of the cell that is covered by the
#' polygon.
#'
#' @param     x a \code{RasterLayer}, \code{RasterStack}, or \code{RasterBrick}
#' @param     y a sf object with polygonal geometries
#' @param     fun an optional function or character vector, as described below
#' @param     weights  a weighting raster to be used with the \code{weighted_mean}
#'                     and \code{weighted_sum} summary operations.
#' @param     quantiles   quantiles to be computed when \code{fun == 'quantile'}
#' @param     append_cols when \code{fun} is not \code{NULL}, an optional
#'                        character vector of columns from \code{y} to be
#'                        included in returned data frame.
#' @param     force_df always return a data frame instead of a vector, even if
#'                     \code{x} has only one layer and \code{fun} has length 1
#' @param     full_colnames include the names of \code{x} in the names of the
#'                          returned data frame, even if \code{x} has only one
#'                          layer. This is useful when the results of multiple
#'                          calls to \code{exact_extract} are combined with
#'                          \code{cbind}.
#' @param     include_cell if \code{TRUE}, and \code{fun} is \code{NULL}, augment
#'                       the returned data frame for each feature with a column
#'                       for the cell index (\code{cell}). If \code{TRUE} and
#'                       \code{fun} is not \code{NULL}, add \code{cell} to the
#'                       data frame passed to \code{fun} for each feature.
#' @param     include_cols an optional character vector of column names in
#'                         \code{y} to be added to the data frame for each
#'                         feature that is either returned (when \code{fun} is
#'                         \code{NULL}) or passed to \code{fun}.
#' @param     include_xy if \code{TRUE}, and \code{fun} is \code{NULL}, augment
#'                       the returned data frame for each feature with columns
#'                       for cell center coordinates (\code{x} and \code{y}). If
#'                       \code{TRUE} and \code{fun} is not \code{NULL}, add
#'                       \code{x} and {y} to the data frame passed to \code{fun}
#'                       for each feature.
#' @param     stack_apply   if \code{TRUE}, apply \code{fun} to each layer of
#'                          \code{x} independently. If \code{FALSE}, apply \code{fun}
#'                          to all layers of \code{x} simultaneously.
#' @param     max_cells_in_memory the maximum number of raster cells to load at
#'                                a given time when using a named summary operation
#'                                for \code{fun} (as opposed to a function defined using
#'                                R code). If a polygon covers more than \code{max_cells_in_memory}
#'                                raster cells, it will be processed in multiple chunks.
#' @param     progress if \code{TRUE}, display a progress bar during processing
#' @param     ... additional arguments to pass to \code{fun}
#' @return a vector or list of data frames, depending on the type of \code{x} and the
#'         value of \code{fun} (see Details)
#' @examples
#' rast <- raster::raster(matrix(1:100, ncol=10), xmn=0, ymn=0, xmx=10, ymx=10)
#' poly <- sf::st_as_sfc('POLYGON ((2 2, 7 6, 4 9, 2 2))')
#'
#' # named summary operation on RasterLayer, returns vector
#' exact_extract(rast, poly, 'mean')
#'
#' # two named summary operations on RasterLayer, returns data frame
#' exact_extract(rast, poly, c('min', 'max'))
#'
#' # named summary operation on RasterStack, returns data frame
#' stk <- raster::stack(list(a=rast, b=sqrt(rast)))
#' exact_extract(stk, poly, 'mean')
#'
#' # named weighted summary operation, returns vector
#' weights <- raster::raster(matrix(runif(100), ncol=10), xmn=0, ymn=0, xmx=10, ymx=10)
#' exact_extract(rast, poly, 'weighted_mean', weights=weights)
#'
#' # custom summary function, returns vector
#' exact_extract(rast, poly, function(value, cov_frac) length(value[cov_frac > 0.9]))
#'
#' @name exact_extract
NULL

#' @import sf
#' @import raster
#' @useDynLib exactextractr
#' @rdname exact_extract
#' @export
setMethod('exact_extract', signature(x='Raster', y='sf'),
          function(x, y, fun=NULL, ...,
                   include_xy=FALSE,
                   progress=TRUE,
                   max_cells_in_memory=30000000,
                   include_cell=FALSE,
                   force_df=FALSE,
                   full_colnames=FALSE,
                   stack_apply=FALSE,
                   append_cols=NULL,
                   include_cols=NULL,
                   quantiles=NULL) {
  .exact_extract(x, y, fun=fun, ...,
                include_xy=include_xy,
                progress=progress,
                max_cells_in_memory=max_cells_in_memory,
                include_cell=include_cell,
                force_df=force_df,
                full_colnames=full_colnames,
                stack_apply=stack_apply,
                append_cols=append_cols,
                include_cols=include_cols,
                quantiles=quantiles)
})

# Return the number of standard (non-...) arguments in a supplied function that
# do not have a default value. This is used to fail if the summary function
# provided by the user cannot accept arguments of values and weights.
.num_expected_args <- function(fun) {
  a <- formals(args(fun))
  a <- a[names(a) != '...']
  sum(sapply(a, nchar) == 0)
}

emptyVector <- function(rast) {
  switch(substr(raster::dataType(rast), 1, 3),
         LOG=logical(),
         INT=integer(),
         numeric())
}

.exact_extract <- function(x, y, fun=NULL, ...,
                           weights=NULL,
                           include_xy=FALSE,
                           progress=TRUE,
                           max_cells_in_memory=30000000,
                           include_cell=FALSE,
                           force_df=FALSE,
                           full_colnames=FALSE,
                           stack_apply=FALSE,
                           append_cols=NULL,
                           include_cols=NULL,
                           quantiles=NULL) {
  if(!is.null(append_cols)) {
    if (!inherits(y, 'sf')) {
      stop(sprintf('append_cols only supported for sf arguments (received %s)',
                   paste(class(y), collapse = ' ')))
    }

    force_df <- TRUE
  }

  if(sf::st_geometry_type(y, by_geometry = FALSE) == 'GEOMETRY') {
    if (!all(sf::st_dimension(y) == 2)) {
      stop("Features in sfc_GEOMETRY must be polygonal")
    }
    y <- sf::st_cast(y, 'MULTIPOLYGON')
  }

  if(!is.null(weights)) {
    if (!startsWith(class(weights), 'Raster')) {
      stop("Weights must be a Raster object.")
    }

    if (is.character(fun) && !any(startsWith(fun, "weighted"))) {
      warning("Weights provided but no requested operations use them.")
    }

    if (!is.na(sf::st_crs(x))) {
      if (is.na(sf::st_crs(weights))) {
        warning("No CRS specified for weighting raster; assuming it has the same CRS as the value raster.")
      } else if (sf::st_crs(x) != sf::st_crs(weights)) {
        stop("Weighting raster does not have the same CRS as value raster.")
      }
    }
  }

  if(is.na(sf::st_crs(x)) && !is.na(sf::st_crs(y))) {
    warning("No CRS specified for raster; assuming it has the same CRS as the polygons.")
  } else if(is.na(sf::st_crs(y)) && !is.na(sf::st_crs(x))) {
    warning("No CRS specified for polygons; assuming they have the same CRS as the raster.")
  } else if(sf::st_crs(x) != sf::st_crs(y)) {
    y <- sf::st_transform(y, sf::st_crs(x))
    warning("Polygons transformed to raster CRS (EPSG:", sf::st_crs(x)$epsg, ")")
  }

  if (!is.null(fun) && !is.character(fun) && .num_expected_args(fun) < 2) {
    stop("exact_extract was called with a function that does not appear to ",
         "be of the form `function(values, coverage_fractions, ...)`")
  }

  if (is.character(fun)) {
    if (length(list(...)) > 0) {
      stop("exact_extract was called with a named summary operation that ",
           "does not accept additional arguments ...")
    }

    if (include_xy) {
      stop("include_xy must be FALSE for named summary operations")
    }

    if (include_cell) {
      stop("include_cell must be FALSE for named summary operations")
    }

    if (!is.null(include_cols)) {
      stop("include_cols not supported for named_summary operations (see argument append_cols)")
    }
  }

  geoms <- sf::st_geometry(y)

  if (progress && length(geoms) > 1) {
    n <- length(geoms)
    pb <- utils::txtProgressBar(min = 0, max = n, initial=0, style=3)
    update_progress <- function() {
      i <- 1 + utils::getTxtProgressBar(pb)
      utils::setTxtProgressBar(pb, i)
      if (i == n) {
        close(pb)
      }
    }
  } else {
    update_progress <- function() {}
  }

  tryCatch({
    x <- raster::readStart(x)
    if (!is.null(weights)) {
      weights <- raster::readStart(weights)
    }

    if (is.character(fun)) {
      # Compute all stats in C++.
      # CPP_stats returns a matrix, which gets turned into a column by sapply
      # Results has one column per feature and one row per stat/raster layer
      if(is.null(weights) && any(.isWeighted(fun))) {
        stop("Weighted stat requested but no weights provided.")
      }

      results <- sapply(sf::st_as_binary(geoms, EWKB=TRUE), function(wkb) {
        ret <- CPP_stats(x, weights, wkb, fun, max_cells_in_memory, quantiles)
        update_progress()
        return(ret)
      })

      if (!is.matrix(results) && !force_df) {
        # Single stat? Return a vector unless asked otherwise via force_df.
        return(results)
      } else {
        # Return a data frame with a column for each stat
        colnames <- .resultColNames(names(x), names(weights), fun, full_colnames, quantiles)

        if (is.matrix(results)) {
          results <- t(results)
        } else {
          results <- matrix(results, nrow=length(results))
        }

        dimnames(results) <- list(NULL, colnames)
        ret <- as.data.frame(results)

        # drop duplicated columns (occurs when an unweighted stat is
        # requested alongside a weighted stat, with a stack of weights)
        ret <- ret[, unique(names(ret)), drop = FALSE]

        if (!is.null(append_cols)) {
          ret <- cbind(sf::st_drop_geometry(y[, append_cols]), ret)
        }

        return(ret)
      }
    } else {
      num_values <- raster::nlayers(x)
      num_weights <- ifelse(is.null(weights), 0, raster::nlayers(weights))
      value_names <- names(x)
      weight_names <- names(weights)

      if (stack_apply || (num_values == 1 && num_weights <= 1)) {
        apply_layerwise <- TRUE

        if (num_values > 1 && num_weights > 1 && num_values != num_weights) {
          stop(sprintf("Can't apply function layerwise with stacks of %d value layers and %d layers", num_values, num_weights))
        }

        result_names <- .resultColNames(value_names, weight_names, fun, full_colnames)
        num_results <- max(num_weights, num_values)
        ind <- .valueWeightIndexes(num_values, num_weights)
      } else {
        apply_layerwise <- FALSE
      }

      # For R summary functions and data frame output, we avoid using
      # input layer names if there is only one value/weight layer
      if (length(value_names) == 1) {
        value_names <- 'value'
      }
      if (length(weight_names) == 1) {
        weight_names <- 'weight'
      }

      ret <- lapply(seq_along(geoms), function(feature_num) {
        wkb <- sf::st_as_binary(geoms[[feature_num]], EWKB=TRUE)

        if (!is.null(include_cols)) {
          include_col_values <- sf::st_drop_geometry(y[feature_num, include_cols])
        } else {
          include_col_values <- NULL
        }

        # only raise a disaggregation warning for the first feature
        warn_on_disaggregate <- feature_num == 1

        col_list <- CPP_exact_extract(x, weights, wkb, include_xy, include_cell, include_col_values, value_names, weight_names, warn_on_disaggregate)
        if (!is.null(include_cols)) {
          # Replicate the include_cols vectors to be as long as the other columns,
          # so we can use quickDf
          nrow <- length(col_list$coverage_fraction)
          col_list[include_cols] <- lapply(col_list[include_cols], rep, nrow)
        }
        df <- .quickDf(col_list)

        update_progress()

        if (is.null(fun)) {
          return(df)
        }

        included_cols_df <- df[, !(names(df) %in% c(value_names, weight_names, 'coverage_fraction')), drop = FALSE]
        vals_df <- df[, value_names, drop = FALSE]
        weights_df <- df[, weight_names, drop = FALSE]
        cov_fracs <- df$coverage_fraction

        if (apply_layerwise) {
          result <- lapply(seq_len(num_results), function(i) {
            vx <- vals_df[, ind$values[i]]
            if (ncol(included_cols_df) > 0) {
              vx <- cbind(data.frame(value = vx), included_cols_df)
            }
            if (num_weights == 0) {
              fun(vx, cov_fracs, ...)
            } else {
              fun(vx, cov_fracs, weights_df[, ind$weights[i]], ...)
            }
          })

          if (num_results == 1) {
            return(result[[1]])
          }

          names(result) <- result_names

          .quickDf(result)
        } else {
          # Pass all layers to callback, to be handled together
          # Included columns (x/y/cell) are passed with the values.
          # Pass single-column data frames as vectors.
          vals_df <- .singleColumnToVector(cbind(vals_df, included_cols_df))
          weights_df <- .singleColumnToVector(weights_df)

          if (num_weights == 0) {
            return(fun(vals_df, cov_fracs, ...))
          } else {
            return(fun(vals_df, cov_fracs, weights_df, ...))
          }
        }
      })

      if (!is.null(fun)) {
        if (all(sapply(ret, is.data.frame))) {
        # function returned a data frame for each polygon? rbind them
          if (requireNamespace('dplyr', quietly = TRUE)) {
            ret <- dplyr::bind_rows(ret) # handle column name mismatches
          } else {
            ret <- do.call(rbind, ret)
          }
        } else {
          # function returned something else; combine the somethings into
          # an array
          ret <- simplify2array(ret)

          if (force_df) {
            ret <- data.frame(result = ret)
          }
        }
      }

      if (!is.null(append_cols)) {
        ret <- cbind(sf::st_drop_geometry(y[, append_cols]), ret)
      }

      return(ret)
    }
  }, finally={
    raster::readStop(x)
    if (!is.null(weights)) {
      raster::readStop(weights)
    }
  })
}

# faster replacement for as.data.frame when input is a named list
# with equal-length columns
# from Advanced R, sec. 24.4.2
.quickDf <- function(lst) {
  class(lst) <- 'data.frame'
  attr(lst, 'row.names') <- .set_row_names(length(lst[[1]]))
  lst
}

.singleColumnToVector <- function(df) {
  if (ncol(df) == 1) {
    df[, 1]
  } else {
    df
  }
}

.appendXY <- function(vals_df, rast, first_row, nrow, first_col, ncol) {
  if (nrow(vals_df) == 0) {
    vals_df$x <- numeric()
    vals_df$y <- numeric()
  } else {
    x_coords <- raster::xFromCol(rast, col=seq(first_col, first_col + ncol - 1))
    y_coords <- raster::yFromRow(rast, row=seq(first_row, first_row + nrow - 1))

    vals_df$x <- rep(x_coords, times=nrow)
    vals_df$y <- rep(y_coords, each=ncol)
  }

  return(vals_df)
}

.appendCell <- function(vals_df, rast, first_row, nrow, first_col, ncol) {
  if (nrow(vals_df) == 0) {
    vals_df$cell <- numeric()
  } else {
    rows <- rep(seq(first_row, first_row + nrow - 1), each = ncol)
    cols <- rep(seq(first_col, first_col + ncol - 1), times = nrow)

    vals_df$cell <- raster::cellFromRowCol(rast, row=rows, col=cols)
  }

  return(vals_df)
}

#' @useDynLib exactextractr
#' @rdname exact_extract
#' @export
setMethod('exact_extract', signature(x='Raster', y='SpatialPolygonsDataFrame'),
          function(x, y, ...) .exact_extract(x, sf::st_as_sf(y), ...))

#' @useDynLib exactextractr
#' @rdname exact_extract
#' @export
setMethod('exact_extract', signature(x='Raster', y='SpatialPolygons'),
          function(x, y, ...) .exact_extract(x, sf::st_as_sf(y), ...))

#' @useDynLib exactextractr
#' @rdname exact_extract
#' @export
setMethod('exact_extract', signature(x='Raster', y='sfc_MULTIPOLYGON'), .exact_extract)

#' @useDynLib exactextractr
#' @rdname exact_extract
#' @export
setMethod('exact_extract', signature(x='Raster', y='sfc_POLYGON'), .exact_extract)

#' @useDynLib exactextractr
#' @rdname exact_extract
#' @export
setMethod('exact_extract', signature(x='Raster', y='sfc_GEOMETRY'), .exact_extract)

# Copyright (c) 2018-2020 ISciences, LLC.
# All rights reserved.
#
# This software is licensed under the Apache License, Version 2.0 (the "License").
# You may not use this file except in compliance with the License. You may
# obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

#' Return column names to be used for summary operations
#'
#' @param value_names names of value raster layers
#' @param weight_names names of weighting raster layers
#' @param fun functions or names of summary operations
#' @param full_colnames return a complete column name even when there is no
#'                      ambiguity?
#' @param quantiles quantiles to use when \code{stat_names} contains \code{quantile}
#' @return character vector of column names
.resultColNames <- function(value_names, weight_names, fun, full_colnames, quantiles=numeric()) {
  if (inherits(fun, 'standardGeneric')) {
    stat_names <- fun@generic[1]
  } else if (is.function(fun)) {
    stat_names <- 'fun'
  } else {
    stat_names <- fun
  }

  quantile_index = which(stat_names == 'quantile')
  if (length(quantile_index) != 0) {
    stat_names <- c(stat_names[seq_along(stat_names) < quantile_index],
                    sprintf('q%02d', as.integer(100 * quantiles)),
                    stat_names[seq_along(stat_names) > quantile_index])
  }

  ind <- .valueWeightIndexes(length(value_names), length(weight_names))
  vn <- value_names[ind$values]
  wn <- weight_names[ind$weights]

  if (length(vn) > 1 || full_colnames) {
    # determine all combinations of index and stat
    z <- expand.grid(index=seq_along(vn),
                     stat=stat_names, stringsAsFactors=FALSE)
    z$value <- vn[z$index]
    if (is.null(wn)) {
      z$weights <- NA
    } else {
      z$weights <- wn[z$index]
    }

    # construct column names for each index, stat
    # add weight layer name only if layer is ambiguously weighted
    mapply(function(stat, value, weight) {
      ret <- stat
      if (full_colnames || length(value_names) > 1) {
        ret <- paste(ret, value, sep='.')
      }
      if (.includeWeightInColName(stat) && ((full_colnames & length(weight_names) > 0)
                                            || length(weight_names) > 1)) {
        ret <- paste(ret, weight, sep='.')
      }

      return(ret)
    }, z$stat, z$value, z$weight, USE.NAMES = FALSE)
  } else {
    stat_names
  }
}

.includeWeightInColName <- function(fun) {
  .isWeighted(fun) | fun == 'fun'
}

.isWeighted <- function(stat_name) {
  stat_name %in% c('weighted_mean', 'weighted_sum')
}

#' Compute indexes for the value and weight layers that should be
#' processed together
#'
#' @param num_values number of layers in value raster
#' @param num_weights number of layers in weighting raster
#' @return list with \code{values} and \code{weights} elements
#'         providing layer indexes
.valueWeightIndexes <- function(num_values, num_weights) {
  if (num_weights == 0) {
    vi <- seq_len(num_values)
    wi <- NA
  } else if (num_values == num_weights) {
    # process in parallel
    vi <- seq_len(num_values)
    wi <- seq_len(num_weights)
  } else if (num_values == 1 && num_weights > 1) {
    # recycle values
    vi <- rep.int(1, num_weights)
    wi <- seq_len(num_weights)
  } else if (num_values > 1 && num_weights == 1) {
    # recycle weights
    vi <- seq_len(num_values)
    wi <- rep.int(1, num_values)
  }

  list(values = vi, weights = wi)
}


# Copyright (c) 2020 ISciences, LLC.
# All rights reserved.
#
# This software is licensed under the Apache License, Version 2.0 (the "License").
# You may not use this file except in compliance with the License. You may
# obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

if (!isGeneric("exact_resample")) {
	setGeneric("exact_resample", function(x, y, ...)
		standardGeneric("exact_resample"))
}

#' Resample a raster to a new grid
#'
#' @param x a \code{RasterLayer} to be resampled
#' @param y a \code{RasterLayer} with a grid definition to which \code{x}
#'          should be resampled
#' @param fun a named summary operation to be used for the resampling
#' @return a resampled version of \code{x}
#'
#' @name exact_resample
NULL

#' @import sf
#' @import raster
#' @useDynLib exactextractr
#' @rdname exact_resample
#' @export
setMethod('exact_resample',
          signature(x='RasterLayer', y='RasterLayer'),
          function(x, y, fun) {
            if (sf::st_crs(x) != sf::st_crs(y)) {
              if (is.na(sf::st_crs(x))) {
                warning("No CRS specified for source raster; assuming it has the same CRS as destination raster.")
              } else if (is.na(sf::st_crs(y))) {
                warning("No CRS specified for destination raster; assuming it has the same CRS as source raster.")
              } else {
                stop('Destination raster must have same CRS as source.')
              }
            }

            x <- raster::readStart(x)
            tryCatch({
              CPP_resample(x, y, fun)
            }, finally={
              raster::readStop(x)
            })
          })

# Copyright (c) 2018-2020 ISciences, LLC.
# All rights reserved.
#
# This software is licensed under the Apache License, Version 2.0 (the "License").
# You may not use this file except in compliance with the License. You may
# obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

if (!isGeneric('coverage_fraction')) {
	setGeneric('coverage_fraction', function(x, y, crop=FALSE, ...)
		standardGeneric('coverage_fraction'))
}

.coverage_fraction <- function(x, y, crop) {
  if(is.na(sf::st_crs(x)) && !is.na(sf::st_crs(y))) {
    warning("No CRS specified for raster; assuming it has the same CRS as the polygons.")
  } else if(is.na(sf::st_crs(y)) && !is.na(sf::st_crs(x))) {
    warning("No CRS specified for polygons; assuming they have the same CRS as the raster.")
  } else if(sf::st_crs(x) != sf::st_crs(y)) {
    y <- sf::st_transform(y, sf::st_crs(x))
    warning("Polygons transformed to raster CRS (EPSG:", sf::st_crs(x)$epsg, ")")
  }

  lapply(sf::st_as_binary(y, EWKB=TRUE), function(wkb) {
    CPP_coverage_fraction(x, wkb, crop)
  })
}

#' Compute the fraction of raster cells covered by a polygon
#'
#' @param     x a (possibly empty) \code{RasterLayer} whose resolution and
#'            extent will be used for the generated \code{RasterLayer}.
#' @param     y a \code{sf} object with polygonal geometries
#' @param     crop if \code{TRUE}, each generated \code{RasterLayer} will be
#'                 cropped to the extent of its associated feature.
#' @return    a list with a \code{RasterLayer} for each feature in \code{y}.
#'            Values of the raster represent the fraction of each
#'            cell in \code{x} that is covered by \code{y}.
#' @examples
#' rast <- raster::raster(matrix(1:100, ncol=10), xmn=0, ymn=0, xmx=10, ymx=10)
#' poly <- sf::st_as_sfc('POLYGON ((2 2, 7 6, 4 9, 2 2))')
#'
#' cov_frac <- coverage_fraction(rast, poly)[[1]]
#' @name coverage_fraction
NULL

#' @import sf
#' @import raster
#' @useDynLib exactextractr
#' @rdname coverage_fraction
#' @export
setMethod('coverage_fraction', signature(x='RasterLayer', y='sf'), function(x, y, crop=FALSE) {
  coverage_fraction(x, sf::st_geometry(y), crop)
})

#' @rdname coverage_fraction
#' @export
setMethod('coverage_fraction', signature(x='RasterLayer', y='sfc_MULTIPOLYGON'), .coverage_fraction)

#' @rdname coverage_fraction
#' @export
setMethod('coverage_fraction', signature(x='RasterLayer', y='sfc_POLYGON'), .coverage_fraction)


// Copyright (c) 2018-2020 ISciences, LLC.
// All rights reserved.
//
// This software is licensed under the Apache License, Version 2.0 (the "License").
// You may not use this file except in compliance with the License. You may
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#include <Rcpp.h>

#include "geos_r.h"

#include "raster_utils.h"

#include "exactextract/src/raster.h"
#include "exactextract/src/raster_cell_intersection.h"

using exactextract::RasterView;
using exactextract::raster_cell_intersection;

// [[Rcpp::export]]
Rcpp::S4 CPP_coverage_fraction(Rcpp::S4 & rast, const Rcpp::RawVector & wkb, bool crop)
{
  GEOSAutoHandle geos;
  Rcpp::Environment raster = Rcpp::Environment::namespace_env("raster");
  Rcpp::Function rasterFn = raster["raster"];
  Rcpp::Function crsFn = raster["crs"];

  auto grid = make_grid(rast);
  auto coverage_fraction = raster_cell_intersection(grid, geos.handle, read_wkb(geos.handle, wkb).get());

  if (crop) {
    grid = coverage_fraction.grid();
  }
  RasterView<float> coverage_view(coverage_fraction, grid);

  Rcpp::NumericMatrix weights{static_cast<int>(grid.rows()),
                              static_cast<int>(grid.cols())};

  for (size_t i = 0; i < grid.rows(); i++) {
    for (size_t j = 0; j < grid.cols(); j++) {
      weights(i, j) = coverage_view(i, j);
      if (!crop && std::isnan(weights(i, j))) {
        weights(i, j) = 0;
      }
    }
  }

  return rasterFn(weights,
                  Rcpp::Named("xmn")=grid.xmin(),
                  Rcpp::Named("xmx")=grid.xmax(),
                  Rcpp::Named("ymn")=grid.ymin(),
                  Rcpp::Named("ymx")=grid.ymax(),
                  Rcpp::Named("crs")=crsFn(rast));
}

// Copyright (c) 2018-2020 ISciences, LLC.
// All rights reserved.
//
// This software is licensed under the Apache License, Version 2.0 (the "License").
// You may not use this file except in compliance with the License. You may
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#ifndef EXACTEXTRACT_MATRIX_H
#define EXACTEXTRACT_MATRIX_H

#include <iomanip>
#include <iterator>
#include <memory>
#include <cstring>
#include <vector>

namespace exactextract {

template<typename T>
class Matrix {

    public:
        using value_type = T;

        Matrix(size_t rows, size_t cols) :
            m_rows{rows},
            m_cols{cols}
        {
            if (m_rows > 0 && m_cols > 0) {
                // new T[]() initializes to zero
                m_data = std::unique_ptr<T[]>(new T[m_rows * m_cols]());
            }
        }

        Matrix(size_t rows, size_t cols, T value) :
                m_rows{rows},
                m_cols{cols}
        {
            if (m_rows > 0 && m_cols > 0) {
                // new T[] does not initialize
                m_data = std::unique_ptr<T[]>(new T[m_rows * m_cols]);
            }

            std::fill(m_data.get(), m_data.get() + m_rows*m_cols, value);
        }

        explicit Matrix(const std::vector<std::vector<T>> & data) :
            m_rows{data.size()},
            m_cols{data[0].size()}
        {
            m_data = std::unique_ptr<T[]>(new T[m_rows*m_cols]());

            auto lastpos = m_data.get();
            for (auto& row : data) {
                lastpos = std::copy(row.begin(), row.end(), lastpos);
            }
        }

        Matrix(Matrix<T>&& m) noexcept :
                m_rows{m.rows()},
                m_cols{m.cols()}
        {
            m_data = std::move(m.m_data);
        }

        T& operator()(size_t row, size_t col) {
            check(row, col);
            return m_data[row*m_cols + col];
        }

        T operator()(size_t row, size_t col) const {
            check(row, col);
            return m_data[row*m_cols + col];
        }

        bool operator==(const Matrix<T> & other) const {
            if (m_rows != other.m_rows) {
                return false;
            }
            if (m_cols != other.m_cols) {
                return false;
            }

            return 0 == memcmp(m_data.get(), other.m_data.get(), m_rows*m_cols*sizeof(T));
        }

        void increment(size_t row, size_t col, const T & val) {
            check(row, col);
            m_data[row*m_cols + col] += val;
        }

        size_t rows() const { return m_rows; }
        size_t cols() const { return m_cols; }

        T* row(size_t row) {
            return &(m_data[row*m_cols]);
        }

        T* data() {
            return m_data.get();
        }

#ifdef MATRIX_CHECK_BOUNDS
        void check(size_t row, size_t col) const {
                if (row + 1 > m_rows) {
                    throw std::out_of_range("Row " + std::to_string(row) + " is out of range.");
                }
                if (col + 1 > m_cols) {
                    throw std::out_of_range("Col " + std::to_string(col) + " is out of range.");
                }
        }
#else
        void check(size_t, size_t) const {}
#endif

    private:
        std::unique_ptr<T[]> m_data;

        size_t m_rows;
        size_t m_cols;

};

template<typename T>
std::ostream& operator<<(std::ostream & os, const Matrix<T> & m) {
    for (size_t i = 0; i < m.rows(); i++) {
        for (size_t j = 0; j < m.cols(); j++) {
            if (m(i, j) != 0) {
                os << std::right << std::fixed << std::setw(10) << std::setprecision(6) <<
                    m(i, j) << " ";
            } else {
                os << "           ";
            }
        }
        os << std::endl;
    }

    return os;
}

}

#endif

// Copyright (c) 2018-2020 ISciences, LLC.
// All rights reserved.
//
// This software is licensed under the Apache License, Version 2.0 (the "License").
// You may not use this file except in compliance with the License. You may
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#ifndef EXACTEXTRACT_RASTER_H
#define EXACTEXTRACT_RASTER_H

#include <cmath>
#include <limits>

#include "grid.h"
#include "matrix.h"

namespace exactextract {
    template<typename T>
    class AbstractRaster {
    public:
        using value_type = T;

        explicit AbstractRaster(const Grid<bounded_extent> & ex) :
            m_grid{ex},
            m_nodata{std::is_floating_point<T>::value ? std::numeric_limits<T>::quiet_NaN() : std::numeric_limits<T>::min()},
            m_has_nodata{false}
            {}

        AbstractRaster(const Grid<bounded_extent> & ex, const T& nodata_val) : m_grid{ex}, m_nodata{nodata_val}, m_has_nodata{true}  {}

        virtual ~AbstractRaster() = default;

        size_t rows() const {
            return m_grid.rows();
        }

        size_t cols() const {
            return m_grid.cols();
        }

        double xres() const {
            return m_grid.dx();
        }

        double yres() const {
            return m_grid.dy();
        }

        double xmin() const {
            return m_grid.xmin();
        }

        double ymin() const {
            return m_grid.ymin();
        }

        double xmax() const {
            return m_grid.xmax();
        }

        double ymax() const {
            return m_grid.ymax();
        }

        const Grid<bounded_extent>& grid() const {
            return m_grid;
        }

        virtual T operator()(size_t row, size_t col) const = 0;

        bool has_nodata() const { return m_has_nodata; }

        T nodata() const { return m_nodata; }

        bool get(size_t row, size_t col, T & val) {
            val = operator()(row, col);

            if (m_has_nodata && val == m_nodata) {
                return false;
            }
            if (std::is_floating_point<T>::value && std::isnan(val)) {
                return false;
            }

            return true;
        }

        void set_nodata(const T & val) {
            m_has_nodata = true;
            m_nodata = val;
        }

        bool operator==(const AbstractRaster<T> & other) const {
            if (rows() != other.rows())
                return false;

            if (cols() != other.cols())
                return false;

            if (xres() != other.xres())
                return false;

            if (yres() != other.yres())
                return false;

            if (xmin() != other.xmin())
                return false;

            if (ymin() != other.ymin())
                return false;

            // Do the rasters differ in their definition of NODATA? If so, we need to do some more detailed
            // checking below.
            bool nodata_differs = has_nodata() != other.has_nodata() ||
                    (has_nodata() && other.has_nodata() && nodata() != other.nodata());

            for (size_t i = 0; i < rows(); i++) {
                for (size_t j = 0; j < cols(); j++) {
                    if(operator()(i, j) != other(i, j)) {
                        // Override default behavior of NAN != NAN
                        if (!std::isnan(operator()(i, j)) || !std::isnan(other(i, j))) {
                            return false;
                        }
                    } else if (nodata_differs && (operator()(i, j) == nodata() || other(i, j) == other.nodata())) {
                        // For data types that do not have NAN, or even for floating point types where the user has
                        // selected some other value to represent NODATA, we need to reverse a positive equality test
                        // where the value is considered to be NODATA in one raster but not the other.
                        return false;
                    }
                }
            }

            return true;
        }

        bool operator!=(const AbstractRaster<T> & other) const {
            return !(operator==(other));
        }
    private:
        Grid<bounded_extent> m_grid;
        T m_nodata;
        bool m_has_nodata;
    };

    template<typename T>
    class Raster : public AbstractRaster<T> {
    public:
        Raster(Matrix<T>&& values, const Box & box) : AbstractRaster<T>(
                Grid<bounded_extent>(
                        box,
                        (box.xmax - box.xmin) / values.cols(),
                        (box.ymax - box.ymin) / values.rows())),
                m_values{std::move(values)} {}

        Raster(Matrix<T>&& values, const Grid<bounded_extent> & g) :
            AbstractRaster<T>(g),
            m_values{std::move(values)} {}

        Raster(const Box & box, size_t nrow, size_t ncol) :
                AbstractRaster<T>(Grid<bounded_extent>(box, (box.xmax-box.xmin) / ncol, (box.ymax-box.ymin) / nrow)),
                m_values{nrow, ncol}
                {}

        explicit Raster(const Grid<bounded_extent> & ex) :
            AbstractRaster<T>(ex),
            m_values{ex.rows(), ex.cols()}
            {}

        Matrix<T>& data() {
            return m_values;
        }

        T& operator()(size_t row, size_t col) {
            return m_values(row, col);
        }

        T operator()(size_t row, size_t col) const override {
            return m_values(row, col);
        }

    private:
        Matrix<T> m_values;
    };

    template<typename T>
    class RasterView : public AbstractRaster<T>{
    public:
        // Construct a view of a raster r at an extent ex that is larger
        // and/or of finer resolution than r
        RasterView(const AbstractRaster<T> & r, Grid<bounded_extent> ex) : AbstractRaster<T>(ex), m_raster{r} {
            double disaggregation_factor_x = r.xres() / ex.dx();
            double disaggregation_factor_y = r.yres() / ex.dy();

            if (std::abs(disaggregation_factor_x - std::floor(disaggregation_factor_x)) > 1e-6 ||
                std::abs(disaggregation_factor_y - std::floor(disaggregation_factor_y)) > 1e-6) {
                throw std::runtime_error("Must construct view at resolution that is an integer multiple of original.");
            }

            if (disaggregation_factor_x < 0 || disaggregation_factor_y < 0) {
                throw std::runtime_error("Must construct view at equal or higher resolution than original.");
            }

            m_x_off = static_cast<long>(std::round((ex.xmin() - r.xmin()) / ex.dx()));
            m_y_off = static_cast<long>(std::round((r.ymax() - ex.ymax()) / ex.dy()));
            m_rx = static_cast<size_t>(disaggregation_factor_x);
            m_ry = static_cast<size_t>(disaggregation_factor_y);

            if (r.has_nodata()) {
                this->set_nodata(r.nodata());
            }
        }

        T operator()(size_t row, size_t col) const override {
            if (m_x_off < 0 && static_cast<size_t>(-m_x_off) > col) {
                return this->nodata();
            }
            if (m_y_off < 0 && static_cast<size_t>(-m_y_off) > row) {
                return this->nodata();
            }

            size_t i0 = (row + m_y_off) / m_ry;
            size_t j0 = (col + m_x_off) / m_rx;

            if (i0 > m_raster.rows() - 1 || j0 > m_raster.cols() - 1) {
                return this->nodata();
            }

            return m_raster(i0, j0);
        }

    private:
        const AbstractRaster<T>& m_raster;

        long m_x_off;
        long m_y_off;
        size_t m_rx;
        size_t m_ry;
    };


    template<typename T>
    std::ostream& operator<<(std::ostream & os, const AbstractRaster<T> & m) {
        for (size_t i = 0; i < m.rows(); i++) {
            for (size_t j = 0; j < m.cols(); j++) {
                if (m(i, j) != 0) {
                    os << std::right << std::fixed << std::setw(10) << std::setprecision(6) <<
                       m(i, j) << " ";
                } else {
                    os << "           ";
                }
            }
            os << std::endl;
        }

        return os;
    }
}

#endif //EXACTEXTRACT_RASTER_H

// Copyright (c) 2018-2020 ISciences, LLC.
// All rights reserved.
//
// This software is licensed under the Apache License, Version 2.0 (the "License").
// You may not use this file except in compliance with the License. You may
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

// [[Rcpp::plugins("cpp14")]]

#pragma once

#include <memory>

#include <Rcpp.h>

#include <geos_c.h>

using geom_ptr = std::unique_ptr<GEOSGeometry, std::function<void(GEOSGeometry*)>>;
using wkb_reader_ptr = std::unique_ptr<GEOSWKBReader, std::function<void(GEOSWKBReader*)>>;

// GEOS warning handler
static void geos_warn(const char* fmt, ...) {
  char buf[BUFSIZ] = { '\0' };

  va_list msg;
  va_start(msg, fmt);
  vsnprintf(buf, BUFSIZ*sizeof(char), fmt, msg);
  va_end(msg);

  Rcpp::Function warning("warning");
  warning(buf);
}

// GEOS error handler
static void geos_error(const char* fmt, ...) {
  char buf[BUFSIZ] = { '\0' };

  va_list msg;
  va_start(msg, fmt);
  vsnprintf(buf, BUFSIZ*sizeof(char), fmt, msg);
  va_end(msg);

  Rcpp::stop(buf);
}

// GEOSContextHandle wrapper to ensure finishGEOS is called.
struct GEOSAutoHandle {
  GEOSAutoHandle() {
    handle = initGEOS_r(geos_warn, geos_error);
  }

  ~GEOSAutoHandle() {
    finishGEOS_r(handle);
  }

  GEOSContextHandle_t handle;
};

// Return a smart pointer to a Geometry, given WKB input
static inline geom_ptr read_wkb(const GEOSContextHandle_t & context, const Rcpp::RawVector & wkb) {
  wkb_reader_ptr wkb_reader{ GEOSWKBReader_create_r(context), [context](GEOSWKBReader* r) { GEOSWKBReader_destroy_r(context, r); } };

  geom_ptr geom{ GEOSWKBReader_read_r(context,
                                      wkb_reader.get(),
                                      std::addressof(wkb[0]),
                                      wkb.size()),
                 [context](GEOSGeometry* g) { GEOSGeom_destroy_r(context, g); } };

  if (geom.get() == nullptr) {
    Rcpp::stop("Failed to parse WKB geometry");
  }

  return geom;
}

// Copyright (c) 2018-2020 ISciences, LLC.
// All rights reserved.
//
// This software is licensed under the Apache License, Version 2.0 (the "License").
// You may not use this file except in compliance with the License. You may
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

// [[Rcpp::plugins("cpp14")]]
#include <memory>

#include <Rcpp.h>

#include "geos_r.h"
#include "raster_utils.h"

#include "s4_raster_source.h"

#include "exactextract/src/geos_utils.h"
#include "exactextract/src/grid.h"
#include "exactextract/src/matrix.h"
#include "exactextract/src/raster_cell_intersection.h"
#include "exactextract/src/raster_source.h"
#include "exactextract/src/raster_stats.h"

using exactextract::Box;
using exactextract::Matrix;
using exactextract::Grid;
using exactextract::subdivide;
using exactextract::bounded_extent;
using exactextract::Raster;
using exactextract::RasterView;
using exactextract::raster_cell_intersection;
using exactextract::RasterStats;
using exactextract::RasterSource;

#include <vector>
#include <string>

// [[Rcpp::export]]
Rcpp::List CPP_exact_extract(Rcpp::S4 & rast,
                             Rcpp::Nullable<Rcpp::S4> & weights,
                             const Rcpp::RawVector & wkb,
                             bool include_xy,
                             bool include_cell_number,
                             Rcpp::Nullable<Rcpp::List> & include_cols,
                             Rcpp::CharacterVector & src_names,
                             Rcpp::Nullable<Rcpp::CharacterVector> & p_weights_names,
                             bool warn_on_disaggregate) {
  GEOSAutoHandle geos;
  Rcpp::Function names("names");

  auto grid = make_grid(rast);
  auto weights_grid = exactextract::Grid<bounded_extent>::make_empty();
  auto common_grid = grid;

  S4RasterSource rsrc(rast);
  int src_nlayers = get_nlayers(rast);

  std::unique_ptr<S4RasterSource> rweights;
  int weights_nlayers = 0;
  Rcpp::CharacterVector weights_names;
  if (!weights.isNull()) {
    Rcpp::S4 weights_s4 = weights.get();
    weights_nlayers = get_nlayers(weights_s4);
    weights_grid = make_grid(weights_s4);
    common_grid = grid.common_grid(weights_grid);

    rweights = std::make_unique<S4RasterSource>(weights_s4);
    weights_names = p_weights_names.get();

    if (warn_on_disaggregate && (common_grid.dx() < grid.dx() || common_grid.dy() < grid.dy())) {
      Rcpp::warning("value raster implicitly disaggregated to match higher resolution of weights");
    }
  }

  auto geom = read_wkb(geos.handle, wkb);

  auto bbox = exactextract::geos_get_box(geos.handle, geom.get());

  common_grid = common_grid.crop(bbox);

  auto coverage_fractions = raster_cell_intersection(common_grid, geos.handle, geom.get());
  auto& cov_grid = coverage_fractions.grid();

  // We can't construct an Rcpp::DataFrame with a variable number of columns
  // because of a bug in Rcpp (see https://github.com/RcppCore/Rcpp/pull/1099)
  // Instead, we build up a list, and then create a data frame from the list.
  // Profiling shows that this ends up in as.data.frame, which is slow.
  // Once Rcpp 1.0.6 is released, this can be reworked to be simpler and faster.
  Rcpp::List cols;

  Rcpp::NumericVector coverage_fraction_vec = as_vector(coverage_fractions);
  Rcpp::LogicalVector covered = coverage_fraction_vec > 0;

  if (include_cols.isNotNull()) {
    Rcpp::List include_cols_list = include_cols.get();
    Rcpp::CharacterVector include_names = include_cols_list.attr("names");

    for (int i = 0; i < include_names.size(); i++) {
      std::string name(include_names[i]);
      cols[name] = include_cols_list[name];
    }
  }

  for (int i = 0; i < src_nlayers; i++) {
    auto values = rsrc.read_box(cov_grid.extent(), i);
    const NumericVectorRaster* r = static_cast<NumericVectorRaster*>(values.get());

    // TODO Perhaps extend this to preserve types (integer, logical.)
    // A bit challenging, since we don't know the type of the raster
    // until we call getValuesBlock, and even then we can't be sure
    // that the returned type is correct (as in a RasterStack with mixed types.)
    // Since R integers are only 32-bit, we are not going to lose data by
    // converting everything to numeric, although we pay a storage penalty.
    Rcpp::NumericVector value_vec = r->vec();
    if (grid.dx() != common_grid.dx() || grid.dy() != common_grid.dy() ||
        value_vec.size() != covered.size()) {
      // Transform values to common grid
      RasterView<double> rt(*r, common_grid);
      value_vec = as_vector(rt);
    }

    value_vec = value_vec[covered];
    cols[std::string(src_names[i])] = value_vec;
  }

  for (int i = 0; i < weights_nlayers; i++) {
    auto values = rweights->read_box(cov_grid.extent(), i);
    const NumericVectorRaster* r = static_cast<NumericVectorRaster*>(values.get());

    Rcpp::NumericVector weight_vec = r->vec();

    if (weights_grid.dx() != common_grid.dx() || weights_grid.dy() != common_grid.dy() ||
        weight_vec.size() != covered.size()) {
      // Transform weights to common grid
      RasterView<double> rt (*r, common_grid);

      weight_vec = as_vector(rt);
    }

    weight_vec = weight_vec[covered];

    std::string colname(weights_names[i]);
    if (cols.containsElementNamed(colname.c_str())) {
      // append ".1" to the column name, to match the behavior of
      // data.frame(). We're safe to just add ".1" (instead of incrementing
      // .1 to .2, for example) because duplicated names within the values
      // or weight stack will already have been made unique by raster::stack()
      colname = colname + ".1"; // append .1
    }
    cols[colname] = weight_vec;
  }

  if (include_xy) {
    // Include xy values from whichever input raster has the higher resolution
    if (weights_nlayers > 0 && (weights_grid.dx() < grid.dx() || weights_grid.dy() < grid.dy())) {
      Rcpp::S4 weights_s4 = weights.get();
      cols["x"] = get_x_values(weights_s4, cov_grid)[covered];
      cols["y"] = get_y_values(weights_s4, cov_grid)[covered];
    } else {
      cols["x"] = get_x_values(rast, cov_grid)[covered];
      cols["y"] = get_y_values(rast, cov_grid)[covered];
    }
  }

  if (include_cell_number) {
    cols["cell"] = get_cell_numbers(rast, cov_grid)[covered];
  }

  cols["coverage_fraction"] = coverage_fraction_vec[covered];

  return cols;
}

// Return a matrix with one row per stat and one row per raster layer
// [[Rcpp::export]]
Rcpp::NumericMatrix CPP_stats(Rcpp::S4 & rast,
                              Rcpp::Nullable<Rcpp::S4> weights,
                              const Rcpp::RawVector & wkb,
                              const Rcpp::StringVector & stats,
                              int max_cells_in_memory,
                              const Rcpp::Nullable<Rcpp::NumericVector> & quantiles) {
  try {
    GEOSAutoHandle geos;

    if (max_cells_in_memory < 1) {
      Rcpp::stop("Invalid value for max_cells_in_memory: %d", max_cells_in_memory);
    }

    int nlayers = get_nlayers(rast);

    S4RasterSource rsrc(rast);

    std::unique_ptr<S4RasterSource> rweights;
    bool weighted = false;
    int nweights = 0;
    if (weights.isNotNull()) {
      Rcpp::S4 weights_s4 = weights.get();
      nweights = get_nlayers(weights_s4);

      if (nlayers > 1 && nweights > 1 && nlayers != nweights) {
        Rcpp::stop("Incompatible number of layers in value and weighting rasters");
      }

      rweights = std::make_unique<S4RasterSource>(weights_s4);
      weighted = true;
    }

    auto geom = read_wkb(geos.handle, wkb);
    auto bbox = exactextract::geos_get_box(geos.handle, geom.get());

    auto grid = weighted ? rsrc.grid().common_grid(rweights->grid()) : rsrc.grid();

    bool disaggregated = (grid.dx() < rsrc.grid().dx() || grid.dy() < rsrc.grid().dy());

    bool store_values = false;
    int stat_result_size = 0;
    for (const auto & stat : stats) {
      store_values = store_values || requires_stored_values(stat);

      if (disaggregated && (stat == std::string("count") ||
                            stat == std::string("sum"))) {
        Rcpp::stop("Cannot compute 'count' or 'sum' when value raster is disaggregated to resolution of weights.");
      }

      if (stat == std::string("quantile")) {
        int num_quantiles = 0;

        if (quantiles.isNotNull()) {
          Rcpp::NumericVector qvec = quantiles.get();
          num_quantiles = qvec.size();
        }

        if (num_quantiles == 0) {
          Rcpp::stop("Quantiles not specified.");
        }

        stat_result_size += num_quantiles;
      } else {
        stat_result_size += 1;
      }
    }

    int nresults = std::max(nlayers, nweights);

    std::vector<RasterStats<double>> raster_stats;
    raster_stats.reserve(nresults);
    for (int i = 0; i < nresults; i++) {
      raster_stats.emplace_back(store_values);
    }

    Rcpp::NumericMatrix stat_results = Rcpp::no_init(nresults, stat_result_size);

    if (bbox.intersects(grid.extent())) {
      auto cropped_grid = grid.crop(bbox);

      for (const auto &subgrid : subdivide(cropped_grid, max_cells_in_memory)) {
        auto coverage_fraction = raster_cell_intersection(subgrid, geos.handle, geom.get());
        auto& cov_grid = coverage_fraction.grid();

        if (!cov_grid.empty()) {
          if (weighted) {
            if (nlayers > nweights) {
              // recycle weights
              auto weights = rweights->read_box(cov_grid.extent(), 0);

              for (int i = 0; i < nlayers; i++) {
                auto values = rsrc.read_box(cov_grid.extent(), i);
                raster_stats[i].process(coverage_fraction, *values, *weights);
              }
            } else if (nweights > nlayers) {
              // recycle values
              auto values = rsrc.read_box(cov_grid.extent(), 0);

              for (int i = 0; i < nweights; i++) {
                auto weights = rweights->read_box(cov_grid.extent(), i);
                raster_stats[i].process(coverage_fraction, *values, *weights);
              }
            } else {
              // process values and weights in parallel
              for (int i = 0; i < nlayers; i++) {
                auto values = rsrc.read_box(cov_grid.extent(), i);
                auto weights = rweights->read_box(cov_grid.extent(), i);

                raster_stats[i].process(coverage_fraction, *values, *weights);
              }
            }
          } else {
            for (int i = 0; i < nlayers; i++) {
              auto values = rsrc.read_box(cov_grid.extent(), i);
              raster_stats[i].process(coverage_fraction, *values);
            }
          }
        }
      }
    }

    for (int j = 0; j < nresults; j++) {
      const auto& rs = raster_stats[j];

      int i = 0;
      for(const auto& stat : stats) {
        if (stat == std::string("mean")) stat_results(j, i++) = rs.mean();

        else if (stat == std::string("sum")) stat_results(j, i++) = rs.sum();
        else if (stat == std::string("count")) stat_results(j, i++) = rs.count();

        else if (stat == std::string("min")) stat_results(j, i++) = rs.min().value_or(NA_REAL);
        else if (stat == std::string("max")) stat_results(j, i++) = rs.max().value_or(NA_REAL);

        else if (stat == std::string("median")) stat_results(j, i++) = rs.quantile(0.5).value_or(NA_REAL);

        else if (stat == std::string("mode")) stat_results(j, i++) = rs.mode().value_or(NA_REAL);
        else if (stat == std::string("majority")) stat_results(j, i++) = rs.mode().value_or(NA_REAL);
        else if (stat == std::string("minority")) stat_results(j, i++) = rs.minority().value_or(NA_REAL);

        else if (stat == std::string("variety")) stat_results(j, i++) = rs.variety();
        else if (stat == std::string("weighted_mean")) stat_results(j, i++) = rs.weighted_mean();
        else if (stat == std::string("weighted_sum")) stat_results(j, i++) = rs.weighted_sum();

        else if (stat == std::string("variance")) stat_results(j, i++) = rs.variance();
        else if (stat == std::string("stdev")) stat_results(j, i++) = rs.stdev();
        else if (stat == std::string("coefficient_of_variation")) stat_results(j, i++) = rs.coefficient_of_variation();

        else if (stat == std::string("quantile")) {
          Rcpp::NumericVector qvec = quantiles.get();
          for (double q : qvec) {
            stat_results(j, i++) = rs.quantile(q).value_or(NA_REAL);
          }
        }

        else Rcpp::stop("Unknown stat: " + stat);
      }
    }

    return stat_results;
  } catch (std::exception & e) {
    // throw predictable exception class
    Rcpp::stop(e.what());
  }
}

// Copyright (c) 2018-2019 ISciences, LLC.
// All rights reserved.
//
// This software is licensed under the Apache License, Version 2.0 (the "License").
// You may not use this file except in compliance with the License. You may
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#ifndef EXACTEXTRACT_BOX_H
#define EXACTEXTRACT_BOX_H

#include "coordinate.h"
#include "crossing.h"
#include "side.h"

#include <limits>

namespace exactextract {

    struct Box {
        double xmin;
        double ymin;
        double xmax;
        double ymax;

        Box(double p_xmin, double p_ymin, double p_xmax, double p_ymax) :
                xmin{p_xmin},
                ymin{p_ymin},
                xmax{p_xmax},
                ymax{p_ymax} {}

        static Box maximum_finite() {
            return {
                std::numeric_limits<double>::lowest(),
                std::numeric_limits<double>::lowest(),
                std::numeric_limits<double>::max(),
                std::numeric_limits<double>::max()
            };
        }

        static Box make_empty() {
            return {0, 0, 0, 0};
        }

        double width() const {
            return xmax - xmin;
        }

        double height() const {
            return ymax - ymin;
        }

        double area() const {
            return width() * height();
        }

        double perimeter() const {
            return 2 * width() + 2 * height();
        }

        bool intersects(const Box & other) const {
            if (other.ymin > ymax)
                return false;
            if (other.ymax < ymin)
                return false;
            if (other.xmin > xmax)
                return false;
            if (other.xmax < xmin)
                return false;

            return true;
        }

        Box intersection(const Box & other) const {
            return {
                std::max(xmin, other.xmin),
                std::max(ymin, other.ymin),
                std::min(xmax, other.xmax),
                std::min(ymax, other.ymax)
            };
        }

        Box translate(double dx, double dy) const {
            return {
                xmin + dx,
                ymin + dy,
                xmax + dx,
                ymax + dy
            };
        }

        Coordinate upper_left() const {
            return Coordinate{xmin, ymax};
        }

        Coordinate upper_right() const {
            return Coordinate{xmax, ymax};
        }

        Coordinate lower_left() const {
            return Coordinate{xmin, ymin};
        }

        Coordinate lower_right() const {
            return Coordinate{xmax, ymin};
        }

        Side side(const Coordinate &c) const;

        Crossing crossing(const Coordinate &c1, const Coordinate &c2) const;

        bool empty() const {
            return xmin >= xmax || ymin >= ymax;
        }

        Box expand_to_include(const Box & other) const {
            if (empty()) {
                return other;
            }

            if (other.empty()) {
                return *this;
            }

            return { std::min(xmin, other.xmin),
                     std::min(ymin, other.ymin),
                     std::max(xmax, other.xmax),
                     std::max(ymax, other.ymax) };
        }

        bool contains(const Box &b) const;

        bool contains(const Coordinate &c) const;

        bool strictly_contains(const Coordinate &c) const;

        bool operator==(const Box& other) const {
            return xmin == other.xmin && xmax == other.xmax && ymin == other.ymin && ymax == other.ymax;
        }

        friend std::ostream &operator<<(std::ostream &os, const Box &c);
    };

    std::ostream &operator<<(std::ostream &os, const Box &c);

}

#endif

// Copyright (c) 2018-2020 ISciences, LLC.
// All rights reserved.
//
// This software is licensed under the Apache License, Version 2.0 (the "License").
// You may not use this file except in compliance with the License. You may
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#ifndef EXACTEXTRACT_RASTER_STATS_H
#define EXACTEXTRACT_RASTER_STATS_H

#include <algorithm>
#include <iostream>
#include <limits>
#include <unordered_map>

#include "raster_cell_intersection.h"
#include "weighted_quantiles.h"
#include "variance.h"

#include "../vend/optional.hpp"

namespace exactextract {

    template<typename T>
    class RasterStats {

    public:
        /**
         * Compute raster statistics from a Raster representing intersection percentages,
         * a Raster representing data values, and (optionally) a Raster representing weights.
         * and a set of raster values.
         */
        explicit RasterStats(bool store_values = false) :
                m_min{std::numeric_limits<T>::max()},
                m_max{std::numeric_limits<T>::lowest()},
                m_sum_ciwi{0},
                m_sum_ci{0},
                m_sum_xici{0},
                m_sum_xiciwi{0},
                m_store_values{store_values} {}

        void process(const Raster<float> & intersection_percentages, const AbstractRaster<T> & rast) {
            RasterView<T> rv{rast, intersection_percentages.grid()};

            for (size_t i = 0; i < rv.rows(); i++) {
                for (size_t j = 0; j < rv.cols(); j++) {
                    float pct_cov = intersection_percentages(i, j);
                    T val;
                    if (pct_cov > 0 && rv.get(i, j, val)) {
                        process_value(val, pct_cov, 1.0);
                    }
                }
            }
        }

        void process(const Raster<float> & intersection_percentages, const AbstractRaster<T> & rast, const AbstractRaster<T> & weights) {
            // Process the entire intersection_percentages grid, even though it may
            // be outside the extent of the weighting raster. Although we've been
            // provided a weighting raster, we still need to calculate correct values
            // for unweighted stats.
            auto& common = intersection_percentages.grid();

            if (common.empty())
                return;

            RasterView<float> iv{intersection_percentages, common};
            RasterView<T> rv{rast,    common};
            RasterView<T> wv{weights, common};

            for (size_t i = 0; i < rv.rows(); i++) {
                for (size_t j = 0; j < rv.cols(); j++) {
                    float pct_cov = iv(i, j);
                    T weight;
                    T val;

                    if (pct_cov > 0 && rv.get(i, j, val)) {
                        if (wv.get(i, j, weight)) {
                            process_value(val, pct_cov, weight);
                        } else {
                            // Weight is NODATA, convert to NAN
                            process_value(val, pct_cov, std::numeric_limits<double>::quiet_NaN());
                        }
                    }
                }
            }
        }

        /**
         * The mean value of cells covered by this polygon, weighted
         * by the percent of the cell that is covered.
         */
        float mean() const {
            return sum() / count();
        }

        /**
         * The mean value of cells covered by this polygon, weighted
         * by the percent of the cell that is covered and a secondary
         * weighting raster.
         *
         * If any weights are undefined, will return NAN. If this is undesirable,
         * caller should replace undefined weights with a suitable default
         * before computing statistics.
         */
         float weighted_mean() const {
             return weighted_sum() / weighted_count();
         }

         /** The fraction of weighted cells to unweighted cells.
          *  Meaningful only when the values of the weighting
          *  raster are between 0 and 1.
          */
         float weighted_fraction() const {
             return weighted_sum() / sum();
         }

        /**
         * The raster value occupying the greatest number of cells
         * or partial cells within the polygon. When multiple values
         * cover the same number of cells, the greatest value will
         * be returned. Weights are not taken into account.
         */
        nonstd::optional<T> mode() const {
            if (variety() == 0) {
                return nonstd::nullopt;
            }

            return std::max_element(m_freq.cbegin(),
                                    m_freq.cend(),
                                    [](const auto &a, const auto &b) {
                                        return a.second < b.second || (a.second == b.second && a.first < b.first);
                                    })->first;
        }

        /**
         * The minimum value in any raster cell wholly or partially covered
         * by the polygon. Weights are not taken into account.
         */
        nonstd::optional<T> min() const {
            if (m_sum_ci == 0) {
                return nonstd::nullopt;
            }
            return m_min;
        }

        /**
         * The maximum value in any raster cell wholly or partially covered
         * by the polygon. Weights are not taken into account.
         */
        nonstd::optional<T> max() const {
            if (m_sum_ci == 0) {
                return nonstd::nullopt;
            }
            return m_max;
        }

        /**
         * The given quantile (0-1) of raster cell values. Coverage fractions
         * are taken into account but weights are not.
         */
        nonstd::optional<T> quantile(double q) const {
            if (m_sum_ci == 0) {
                return nonstd::nullopt;
            }

            // The weighted quantile computation is not processed incrementally.
            // Create it on demand and retain it in case we want multiple quantiles.
            if (!m_quantiles) {
                m_quantiles = std::make_unique<WeightedQuantiles>();

                for (const auto& entry : m_freq) {
                    m_quantiles->process(entry.first, entry.second);
                }
            }

            return m_quantiles->quantile(q);
        }

        /**
         * The sum of raster cells covered by the polygon, with each raster
         * value weighted by its coverage fraction.
         */
        float sum() const {
            return (float) m_sum_xici;
        }

        /**
         * The sum of raster cells covered by the polygon, with each raster
         * value weighted by its coverage fraction and weighting raster value.
         *
         * If any weights are undefined, will return NAN. If this is undesirable,
         * caller should replace undefined weights with a suitable default
         * before computing statistics.
         */
        float weighted_sum() const {
            return (float) m_sum_xiciwi;
        }

        /**
         * The number of raster cells with a defined value
         * covered by the polygon. Weights are not taken
         * into account.
         */
        float count() const {
            return (float) m_sum_ci;
        }

        /**
         * The population variance of raster cells touched
         * by the polygon. Cell coverage fractions are taken
         * into account; values of a weighting raster are not.
         */
        float variance() const {
            return static_cast<float>(m_variance.variance());
        }

        /**
         * The population standard deviation of raster cells
         * touched by the polygon. Cell coverage fractions
         * are taken into account; values of a weighting
         * raster are not.
         */
        float stdev() const {
            return static_cast<float>(m_variance.stdev());
        }

        /**
         * The population coefficient of variation of raster
         * cells touched by the polygon. Cell coverage fractions
         * are taken into account; values of a weighting
         * raster are not.
         */
        float coefficient_of_variation() const {
            return static_cast<float>(m_variance.coefficent_of_variation());
        }

        /**
         * The sum of weights for each cell covered by the
         * polygon, with each weight multiplied by the coverage
         * coverage fraction of each cell.
         *
         * If any weights are undefined, will return NAN. If this is undesirable,
         * caller should replace undefined weights with a suitable default
         * before computing statistics.
         */
        float weighted_count() const {
            return (float) m_sum_ciwi;
        }

        /**
         * The raster value occupying the least number of cells
         * or partial cells within the polygon. When multiple values
         * cover the same number of cells, the lowest value will
         * be returned.
         *
         * Cell weights are not taken into account.
         */
        nonstd::optional<T> minority() const {
            if (variety() == 0) {
                return nonstd::nullopt;
            }

            return std::min_element(m_freq.cbegin(),
                                    m_freq.cend(),
                                    [](const auto &a, const auto &b) {
                                        return a.second < b.second || (a.second == b.second && a.first < b.first);
                                    })->first;
        }

        /**
         * The number of distinct defined raster values in cells wholly
         * or partially covered by the polygon.
         */
        size_t variety() const {
            return m_freq.size();
        }

        bool stores_values() const {
            return m_store_values;
        }

    private:
        T m_min;
        T m_max;

        // ci: coverage fraction of pixel i
        // wi: weight of pixel i
        // xi: value of pixel i
        double m_sum_ciwi;
        double m_sum_ci;
        double m_sum_xici;
        double m_sum_xiciwi;
        WestVariance m_variance;

        mutable std::unique_ptr<WeightedQuantiles> m_quantiles;

        std::unordered_map<T, float> m_freq;

        bool m_store_values;

        void process_value(const T& val, float coverage, double weight) {
            m_sum_ci += static_cast<double>(coverage);
            m_sum_xici += val*static_cast<double>(coverage);

            m_variance.process(val, coverage);

            double ciwi = static_cast<double>(coverage)*weight;
            m_sum_ciwi += ciwi;
            m_sum_xiciwi += val * ciwi;

            if (val < m_min) {
                m_min = val;
            }

            if (val > m_max) {
                m_max = val;
            }

            // TODO should weights factor in here?
            if (m_store_values) {
                m_freq[val] += coverage;
                m_quantiles.reset();
            }
        }
    };

    template<typename T>
    std::ostream& operator<<(std::ostream& os, const RasterStats<T> & stats) {
        os << "{" << std::endl;
        os << "  \"count\" : " << stats.count() << "," << std::endl;

        os << "  \"min\" : ";
        if (stats.min().has_value()) {
            os << stats.min().value();
        } else {
            os << "null";
        }
        os << "," << std::endl;

        os << "  \"max\" : ";
        if (stats.max().has_value()) {
            os << stats.max().value();
        } else {
            os << "null";
        }
        os << "," << std::endl;

        os << "  \"mean\" : " << stats.mean() << "," << std::endl;
        os << "  \"sum\" : " << stats.sum() << "," << std::endl;
        os << "  \"weighted_mean\" : " << stats.weighted_mean() << "," << std::endl;
        os << "  \"weighted_sum\" : " << stats.weighted_sum();
        if (stats.stores_values()) {
            os << "," << std::endl;
            os << "  \"mode\" : ";
            if (stats.mode().has_value()) {
                os << stats.mode().value();
            } else {
                os << "null";
            }
            os << "," << std::endl;

            os << "  \"minority\" : ";
            if (stats.minority().has_value()) {
                os << stats.minority().value();
            } else {
                os << "null";
            }
            os << "," << std::endl;

            os << "  \"variety\" : " << stats.variety() << std::endl;
        } else {
            os << std::endl;
        }
        os << "}" << std::endl;
        return os;
    }


}

#endif

// Copyright (c) 2020 ISciences, LLC.
// All rights reserved.
//
// This software is licensed under the Apache License, Version 2.0 (the "License").
// You may not use this file except in compliance with the License. You may
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#ifndef EXACTEXTRACT_VARIANCE_H
#define EXACTEXTRACT_VARIANCE_H

#include <cmath>

namespace exactextract {

class WestVariance {
    /** \brief Implements an incremental algorithm for weighted standard
     * deviation, variance, and coefficient of variation, as described in
     * formula WV2 of West, D.H.D. (1979) "Updating Mean and Variance
     * Estimates: An Improved Method". Communications of the ACM 22(9).
     */

private:
    double sum_w = 0;
    double mean = 0;
    double t = 0;

public:
    /** \brief Update variance estimate with another value
     *
     * @param x value to add
     * @param w weight of `x`
     */
    void process(double x, double w) {
        double mean_old = mean;

        sum_w += w;
        mean += (w / sum_w) * (x - mean_old);
        t += w * (x - mean_old) * (x - mean);
    }

    /** \brief Return the population variance.
     */
    constexpr double variance() const {
        return t / sum_w;
    }

    /** \brief Return the population standard deviation
     */
    double stdev() const {
        return std::sqrt(variance());
    }

    /** \brief Return the population coefficient of variation
     */
    double coefficent_of_variation() const {
        return stdev() / mean;
    }

};

}

#endif //EXACTEXTRACT_VARIANCE_H

// Copyright (c) 2019-2020 ISciences, LLC.
// All rights reserved.
//
// This software is licensed under the Apache License, Version 2.0 (the "License").
// You may not use this file except in compliance with the License. You may
// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#ifndef EXACTEXTRACT_WEIGHTED_QUANTILES_H
#define EXACTEXTRACT_WEIGHTED_QUANTILES_H

#include <algorithm>
#include <cmath>
#include <stdexcept>
#include <vector>

namespace exactextract {

    class WeightedQuantiles {
        // Compute weighted quantiles based on https://stats.stackexchange.com/a/13223

    public:

        void process(double x, double w) {
            if (w < 0) {
                throw std::runtime_error("Weighted quantile calculation does not support negative weights.");
            }

            if (!std::isfinite(w)) {
                throw std::runtime_error("Weighted quantile does not support non-finite weights.");
            }

            m_ready_to_query = false;

            m_elems.emplace_back(x, w);
        }

        double quantile(double q) const;

    private:
        struct elem_t {
            elem_t(double _x, double _w) : x(_x), w(_w), cumsum(0), s(0) {}

            double x;
            double w;
            double cumsum;
            double s;
        };

        void prepare() const;

        mutable std::vector<elem_t> m_elems;
        mutable double m_sum_w;
        mutable bool m_ready_to_query;
    };

}

#endif //EXACTEXTRACT_WEIGHTED_QUANTILES_H

//
// Copyright (c) 2014-2018 Martin Moene
//
// https://github.com/martinmoene/optional-lite
//
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

#pragma once

#ifndef NONSTD_OPTIONAL_LITE_HPP
#define NONSTD_OPTIONAL_LITE_HPP

#define optional_lite_MAJOR  3
#define optional_lite_MINOR  2
#define optional_lite_PATCH  0

#define optional_lite_VERSION  optional_STRINGIFY(optional_lite_MAJOR) "." optional_STRINGIFY(optional_lite_MINOR) "." optional_STRINGIFY(optional_lite_PATCH)

#define optional_STRINGIFY(  x )  optional_STRINGIFY_( x )
#define optional_STRINGIFY_( x )  #x

// optional-lite configuration:

#define optional_OPTIONAL_DEFAULT  0
#define optional_OPTIONAL_NONSTD   1
#define optional_OPTIONAL_STD      2

#if !defined( optional_CONFIG_SELECT_OPTIONAL )
# define optional_CONFIG_SELECT_OPTIONAL  ( optional_HAVE_STD_OPTIONAL ? optional_OPTIONAL_STD : optional_OPTIONAL_NONSTD )
#endif

// Control presence of exception handling (try and auto discover):

#ifndef optional_CONFIG_NO_EXCEPTIONS
# if defined(__cpp_exceptions) || defined(__EXCEPTIONS) || defined(_CPPUNWIND)
#  define optional_CONFIG_NO_EXCEPTIONS  0
# else
#  define optional_CONFIG_NO_EXCEPTIONS  1
# endif
#endif

// C++ language version detection (C++20 is speculative):
// Note: VC14.0/1900 (VS2015) lacks too much from C++14.

#ifndef   optional_CPLUSPLUS
# if defined(_MSVC_LANG ) && !defined(__clang__)
#  define optional_CPLUSPLUS  (_MSC_VER == 1900 ? 201103L : _MSVC_LANG )
# else
#  define optional_CPLUSPLUS  __cplusplus
# endif
#endif

#define optional_CPP98_OR_GREATER  ( optional_CPLUSPLUS >= 199711L )
#define optional_CPP11_OR_GREATER  ( optional_CPLUSPLUS >= 201103L )
#define optional_CPP11_OR_GREATER_ ( optional_CPLUSPLUS >= 201103L )
#define optional_CPP14_OR_GREATER  ( optional_CPLUSPLUS >= 201402L )
#define optional_CPP17_OR_GREATER  ( optional_CPLUSPLUS >= 201703L )
#define optional_CPP20_OR_GREATER  ( optional_CPLUSPLUS >= 202000L )

// C++ language version (represent 98 as 3):

#define optional_CPLUSPLUS_V  ( optional_CPLUSPLUS / 100 - (optional_CPLUSPLUS > 200000 ? 2000 : 1994) )

// Use C++17 std::optional if available and requested:

#if optional_CPP17_OR_GREATER && defined(__has_include )
# if __has_include( <optional> )
#  define optional_HAVE_STD_OPTIONAL  1
# else
#  define optional_HAVE_STD_OPTIONAL  0
# endif
#else
# define  optional_HAVE_STD_OPTIONAL  0
#endif

#define optional_USES_STD_OPTIONAL  ( (optional_CONFIG_SELECT_OPTIONAL == optional_OPTIONAL_STD) || ((optional_CONFIG_SELECT_OPTIONAL == optional_OPTIONAL_DEFAULT) && optional_HAVE_STD_OPTIONAL) )

//
// in_place: code duplicated in any-lite, expected-lite, optional-lite, value-ptr-lite, variant-lite:
//

#ifndef nonstd_lite_HAVE_IN_PLACE_TYPES
#define nonstd_lite_HAVE_IN_PLACE_TYPES  1

// C++17 std::in_place in <utility>:

#if optional_CPP17_OR_GREATER

#include <utility>

namespace nonstd {

using std::in_place;
using std::in_place_type;
using std::in_place_index;
using std::in_place_t;
using std::in_place_type_t;
using std::in_place_index_t;

#define nonstd_lite_in_place_t(      T)  std::in_place_t
#define nonstd_lite_in_place_type_t( T)  std::in_place_type_t<T>
#define nonstd_lite_in_place_index_t(K)  std::in_place_index_t<K>

#define nonstd_lite_in_place(      T)    std::in_place_t{}
#define nonstd_lite_in_place_type( T)    std::in_place_type_t<T>{}
#define nonstd_lite_in_place_index(K)    std::in_place_index_t<K>{}

} // namespace nonstd

#else // optional_CPP17_OR_GREATER

#include <cstddef>

namespace nonstd {
namespace detail {

template< class T >
struct in_place_type_tag {};

template< std::size_t K >
struct in_place_index_tag {};

} // namespace detail

struct in_place_t {};

template< class T >
inline in_place_t in_place( detail::in_place_type_tag<T> /*unused*/ = detail::in_place_type_tag<T>() )
{
    return in_place_t();
}

template< std::size_t K >
inline in_place_t in_place( detail::in_place_index_tag<K> /*unused*/ = detail::in_place_index_tag<K>() )
{
    return in_place_t();
}

template< class T >
inline in_place_t in_place_type( detail::in_place_type_tag<T> /*unused*/ = detail::in_place_type_tag<T>() )
{
    return in_place_t();
}

template< std::size_t K >
inline in_place_t in_place_index( detail::in_place_index_tag<K> /*unused*/ = detail::in_place_index_tag<K>() )
{
    return in_place_t();
}

// mimic templated typedef:

#define nonstd_lite_in_place_t(      T)  nonstd::in_place_t(&)( nonstd::detail::in_place_type_tag<T>  )
#define nonstd_lite_in_place_type_t( T)  nonstd::in_place_t(&)( nonstd::detail::in_place_type_tag<T>  )
#define nonstd_lite_in_place_index_t(K)  nonstd::in_place_t(&)( nonstd::detail::in_place_index_tag<K> )

#define nonstd_lite_in_place(      T)    nonstd::in_place_type<T>
#define nonstd_lite_in_place_type( T)    nonstd::in_place_type<T>
#define nonstd_lite_in_place_index(K)    nonstd::in_place_index<K>

} // namespace nonstd

#endif // optional_CPP17_OR_GREATER
#endif // nonstd_lite_HAVE_IN_PLACE_TYPES

//
// Using std::optional:
//

#if optional_USES_STD_OPTIONAL

#include <optional>

namespace nonstd {

    using std::optional;
    using std::bad_optional_access;
    using std::hash;

    using std::nullopt;
    using std::nullopt_t;

    using std::operator==;
    using std::operator!=;
    using std::operator<;
    using std::operator<=;
    using std::operator>;
    using std::operator>=;
    using std::make_optional;
    using std::swap;
}

#else // optional_USES_STD_OPTIONAL

#include <cassert>
#include <utility>

// optional-lite alignment configuration:

#ifndef  optional_CONFIG_MAX_ALIGN_HACK
# define optional_CONFIG_MAX_ALIGN_HACK  0
#endif

#ifndef  optional_CONFIG_ALIGN_AS
// no default, used in #if defined()
#endif

#ifndef  optional_CONFIG_ALIGN_AS_FALLBACK
# define optional_CONFIG_ALIGN_AS_FALLBACK  double
#endif

// Compiler warning suppression:

#if defined(__clang__)
# pragma clang diagnostic push
# pragma clang diagnostic ignored "-Wundef"
#elif defined(__GNUC__)
# pragma GCC   diagnostic push
# pragma GCC   diagnostic ignored "-Wundef"
#elif defined(_MSC_VER )
# pragma warning( push )
#endif

// half-open range [lo..hi):
#define optional_BETWEEN( v, lo, hi ) ( (lo) <= (v) && (v) < (hi) )

// Compiler versions:
//
// MSVC++ 6.0  _MSC_VER == 1200 (Visual Studio 6.0)
// MSVC++ 7.0  _MSC_VER == 1300 (Visual Studio .NET 2002)
// MSVC++ 7.1  _MSC_VER == 1310 (Visual Studio .NET 2003)
// MSVC++ 8.0  _MSC_VER == 1400 (Visual Studio 2005)
// MSVC++ 9.0  _MSC_VER == 1500 (Visual Studio 2008)
// MSVC++ 10.0 _MSC_VER == 1600 (Visual Studio 2010)
// MSVC++ 11.0 _MSC_VER == 1700 (Visual Studio 2012)
// MSVC++ 12.0 _MSC_VER == 1800 (Visual Studio 2013)
// MSVC++ 14.0 _MSC_VER == 1900 (Visual Studio 2015)
// MSVC++ 14.1 _MSC_VER >= 1910 (Visual Studio 2017)

#if defined(_MSC_VER ) && !defined(__clang__)
# define optional_COMPILER_MSVC_VER      (_MSC_VER )
# define optional_COMPILER_MSVC_VERSION  (_MSC_VER / 10 - 10 * ( 5 + (_MSC_VER < 1900 ) ) )
#else
# define optional_COMPILER_MSVC_VER      0
# define optional_COMPILER_MSVC_VERSION  0
#endif

#define optional_COMPILER_VERSION( major, minor, patch )  ( 10 * (10 * (major) + (minor) ) + (patch) )

#if defined(__GNUC__) && !defined(__clang__)
# define optional_COMPILER_GNUC_VERSION   optional_COMPILER_VERSION(__GNUC__, __GNUC_MINOR__, __GNUC_PATCHLEVEL__)
#else
# define optional_COMPILER_GNUC_VERSION   0
#endif

#if defined(__clang__)
# define optional_COMPILER_CLANG_VERSION  optional_COMPILER_VERSION(__clang_major__, __clang_minor__, __clang_patchlevel__)
#else
# define optional_COMPILER_CLANG_VERSION  0
#endif

#if optional_BETWEEN(optional_COMPILER_MSVC_VERSION, 70, 140 )
# pragma warning( disable: 4345 )   // initialization behavior changed
#endif

#if optional_BETWEEN(optional_COMPILER_MSVC_VERSION, 70, 150 )
# pragma warning( disable: 4814 )   // in C++14 'constexpr' will not imply 'const'
#endif

// Presence of language and library features:

#define optional_HAVE(FEATURE) ( optional_HAVE_##FEATURE )

#ifdef _HAS_CPP0X
# define optional_HAS_CPP0X  _HAS_CPP0X
#else
# define optional_HAS_CPP0X  0
#endif

// Unless defined otherwise below, consider VC14 as C++11 for optional-lite:

#if optional_COMPILER_MSVC_VER >= 1900
# undef  optional_CPP11_OR_GREATER
# define optional_CPP11_OR_GREATER  1
#endif

#define optional_CPP11_90   (optional_CPP11_OR_GREATER_ || optional_COMPILER_MSVC_VER >= 1500)
#define optional_CPP11_100  (optional_CPP11_OR_GREATER_ || optional_COMPILER_MSVC_VER >= 1600)
#define optional_CPP11_110  (optional_CPP11_OR_GREATER_ || optional_COMPILER_MSVC_VER >= 1700)
#define optional_CPP11_120  (optional_CPP11_OR_GREATER_ || optional_COMPILER_MSVC_VER >= 1800)
#define optional_CPP11_140  (optional_CPP11_OR_GREATER_ || optional_COMPILER_MSVC_VER >= 1900)
#define optional_CPP11_141  (optional_CPP11_OR_GREATER_ || optional_COMPILER_MSVC_VER >= 1910)

#define optional_CPP14_000  (optional_CPP14_OR_GREATER)
#define optional_CPP17_000  (optional_CPP17_OR_GREATER)

// Presence of C++11 language features:

#define optional_HAVE_CONSTEXPR_11      optional_CPP11_140
#define optional_HAVE_IS_DEFAULT        optional_CPP11_140
#define optional_HAVE_NOEXCEPT          optional_CPP11_140
#define optional_HAVE_NULLPTR           optional_CPP11_100
#define optional_HAVE_REF_QUALIFIER     optional_CPP11_140

// Presence of C++14 language features:

#define optional_HAVE_CONSTEXPR_14      optional_CPP14_000

// Presence of C++17 language features:

#define optional_HAVE_NODISCARD         optional_CPP17_000

// Presence of C++ library features:

#define optional_HAVE_CONDITIONAL       optional_CPP11_120
#define optional_HAVE_REMOVE_CV         optional_CPP11_120
#define optional_HAVE_TYPE_TRAITS       optional_CPP11_90

#define optional_HAVE_TR1_TYPE_TRAITS   (!! optional_COMPILER_GNUC_VERSION )
#define optional_HAVE_TR1_ADD_POINTER   (!! optional_COMPILER_GNUC_VERSION )

// C++ feature usage:

#if optional_HAVE( CONSTEXPR_11 )
# define optional_constexpr  constexpr
#else
# define optional_constexpr  /*constexpr*/
#endif

#if optional_HAVE( IS_DEFAULT )
# define optional_is_default  = default;
#else
# define optional_is_default  {}
#endif

#if optional_HAVE( CONSTEXPR_14 )
# define optional_constexpr14  constexpr
#else
# define optional_constexpr14  /*constexpr*/
#endif

#if optional_HAVE( NODISCARD )
# define optional_nodiscard  [[nodiscard]]
#else
# define optional_nodiscard  /*[[nodiscard]]*/
#endif

#if optional_HAVE( NOEXCEPT )
# define optional_noexcept  noexcept
#else
# define optional_noexcept  /*noexcept*/
#endif

#if optional_HAVE( NULLPTR )
# define optional_nullptr  nullptr
#else
# define optional_nullptr  NULL
#endif

#if optional_HAVE( REF_QUALIFIER )
// NOLINTNEXTLINE( bugprone-macro-parentheses )
# define optional_ref_qual  &
# define optional_refref_qual  &&
#else
# define optional_ref_qual  /*&*/
# define optional_refref_qual  /*&&*/
#endif

// additional includes:

#if optional_CONFIG_NO_EXCEPTIONS
// already included: <cassert>
#else
# include <stdexcept>
#endif

#if optional_CPP11_OR_GREATER
# include <functional>
#endif

#if optional_HAVE( INITIALIZER_LIST )
# include <initializer_list>
#endif

#if optional_HAVE( TYPE_TRAITS )
# include <type_traits>
#elif optional_HAVE( TR1_TYPE_TRAITS )
# include <tr1/type_traits>
#endif

// Method enabling

#if optional_CPP11_OR_GREATER

#define optional_REQUIRES_0(...) \
    template< bool B = (__VA_ARGS__), typename std::enable_if<B, int>::type = 0 >

#define optional_REQUIRES_T(...) \
    , typename = typename std::enable_if< (__VA_ARGS__), nonstd::optional_lite::detail::enabler >::type

#define optional_REQUIRES_R(R, ...) \
    typename std::enable_if< (__VA_ARGS__), R>::type

#define optional_REQUIRES_A(...) \
    , typename std::enable_if< (__VA_ARGS__), void*>::type = nullptr

#endif

//
// optional:
//

namespace nonstd { namespace optional_lite {

namespace std11 {

#if optional_CPP11_OR_GREATER
    using std::move;
#else
    template< typename T > T & move( T & t ) { return t; }
#endif

#if optional_HAVE( CONDITIONAL )
    using std::conditional;
#else
    template< bool B, typename T, typename F > struct conditional              { typedef T type; };
    template<         typename T, typename F > struct conditional<false, T, F> { typedef F type; };
#endif // optional_HAVE_CONDITIONAL

} // namespace std11

#if optional_CPP11_OR_GREATER

/// type traits C++17:

namespace std17 {

#if optional_CPP17_OR_GREATER

using std::is_swappable;
using std::is_nothrow_swappable;

#elif optional_CPP11_OR_GREATER

namespace detail {

using std::swap;

struct is_swappable
{
    template< typename T, typename = decltype( swap( std::declval<T&>(), std::declval<T&>() ) ) >
    static std::true_type test( int /*unused*/ );

    template< typename >
    static std::false_type test(...);
};

struct is_nothrow_swappable
{
    // wrap noexcept(expr) in separate function as work-around for VC140 (VS2015):

    template< typename T >
    static constexpr bool satisfies()
    {
        return noexcept( swap( std::declval<T&>(), std::declval<T&>() ) );
    }

    template< typename T >
    static auto test( int /*unused*/ ) -> std::integral_constant<bool, satisfies<T>()>{}

    template< typename >
    static auto test(...) -> std::false_type;
};

} // namespace detail

// is [nothow] swappable:

template< typename T >
struct is_swappable : decltype( detail::is_swappable::test<T>(0) ){};

template< typename T >
struct is_nothrow_swappable : decltype( detail::is_nothrow_swappable::test<T>(0) ){};

#endif // optional_CPP17_OR_GREATER

} // namespace std17

/// type traits C++20:

namespace std20 {

template< typename T >
struct remove_cvref
{
    typedef typename std::remove_cv< typename std::remove_reference<T>::type >::type type;
};

} // namespace std20

#endif // optional_CPP11_OR_GREATER

/// class optional

template< typename T >
class optional;

namespace detail {

// for optional_REQUIRES_T

#if optional_CPP11_OR_GREATER
enum class enabler{};
#endif

// C++11 emulation:

struct nulltype{};

template< typename Head, typename Tail >
struct typelist
{
    typedef Head head;
    typedef Tail tail;
};

#if optional_CONFIG_MAX_ALIGN_HACK

// Max align, use most restricted type for alignment:

#define optional_UNIQUE(  name )       optional_UNIQUE2( name, __LINE__ )
#define optional_UNIQUE2( name, line ) optional_UNIQUE3( name, line )
#define optional_UNIQUE3( name, line ) name ## line

#define optional_ALIGN_TYPE( type ) \
    type optional_UNIQUE( _t ); struct_t< type > optional_UNIQUE( _st )

template< typename T >
struct struct_t { T _; };

union max_align_t
{
    optional_ALIGN_TYPE( char );
    optional_ALIGN_TYPE( short int );
    optional_ALIGN_TYPE( int );
    optional_ALIGN_TYPE( long int  );
    optional_ALIGN_TYPE( float  );
    optional_ALIGN_TYPE( double );
    optional_ALIGN_TYPE( long double );
    optional_ALIGN_TYPE( char * );
    optional_ALIGN_TYPE( short int * );
    optional_ALIGN_TYPE( int *  );
    optional_ALIGN_TYPE( long int * );
    optional_ALIGN_TYPE( float * );
    optional_ALIGN_TYPE( double * );
    optional_ALIGN_TYPE( long double * );
    optional_ALIGN_TYPE( void * );

#ifdef HAVE_LONG_LONG
    optional_ALIGN_TYPE( long long );
#endif

    struct Unknown;

    Unknown ( * optional_UNIQUE(_) )( Unknown );
    Unknown * Unknown::* optional_UNIQUE(_);
    Unknown ( Unknown::* optional_UNIQUE(_) )( Unknown );

    struct_t< Unknown ( * )( Unknown)         > optional_UNIQUE(_);
    struct_t< Unknown * Unknown::*            > optional_UNIQUE(_);
    struct_t< Unknown ( Unknown::* )(Unknown) > optional_UNIQUE(_);
};

#undef optional_UNIQUE
#undef optional_UNIQUE2
#undef optional_UNIQUE3

#undef optional_ALIGN_TYPE

#elif defined( optional_CONFIG_ALIGN_AS ) // optional_CONFIG_MAX_ALIGN_HACK

// Use user-specified type for alignment:

#define optional_ALIGN_AS( unused ) \
    optional_CONFIG_ALIGN_AS

#else // optional_CONFIG_MAX_ALIGN_HACK

// Determine POD type to use for alignment:

#define optional_ALIGN_AS( to_align ) \
    typename type_of_size< alignment_types, alignment_of< to_align >::value >::type

template< typename T >
struct alignment_of;

template< typename T >
struct alignment_of_hack
{
    char c;
    T t;
    alignment_of_hack();
};

template< size_t A, size_t S >
struct alignment_logic
{
    enum { value = A < S ? A : S };
};

template< typename T >
struct alignment_of

1		# Copyright (c) 2018-2020 ISciences, LLC.
2		# All rights reserved.
3		#
4		# This software is licensed under the Apache License, Version 2.0 (the "License").
5		# You may not use this file except in compliance with the License. You may
6		# obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
7		#
8		# Unless required by applicable law or agreed to in writing, software
9		# distributed under the License is distributed on an "AS IS" BASIS,
10		# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11		# See the License for the specific language governing permissions and
12		# limitations under the License.
13
14		#' Return column names to be used for summary operations
15		#'
16		#' @param value_names names of value raster layers
17		#' @param weight_names names of weighting raster layers
18		#' @param fun functions or names of summary operations
19		#' @param full_colnames return a complete column name even when there is no
20		#' ambiguity?
21		#' @param quantiles quantiles to use when \code{stat_names} contains \code{quantile}
22		#' @return character vector of column names
23		.resultColNames <- function(value_names, weight_names, fun, full_colnames, quantiles=numeric()) {
24	84x	if (inherits(fun, 'standardGeneric')) {
25	16x	stat_names <- fun@generic[1]
26	68x	} else if (is.function(fun)) {
27	24x	stat_names <- 'fun'
28		} else {
29	44x	stat_names <- fun
30		}
31
32	84x	quantile_index = which(stat_names == 'quantile')
33	84x	if (length(quantile_index) != 0) {
34	1x	stat_names <- c(stat_names[seq_along(stat_names) < quantile_index],
35	1x	sprintf('q%02d', as.integer(100 * quantiles)),
36	1x	stat_names[seq_along(stat_names) > quantile_index])
37		}
38
39	84x	ind <- .valueWeightIndexes(length(value_names), length(weight_names))
40	84x	vn <- value_names[ind$values]
41	84x	wn <- weight_names[ind$weights]
42
43	84x	if (length(vn) > 1 \|\| full_colnames) {
44		# determine all combinations of index and stat
45	33x	z <- expand.grid(index=seq_along(vn),
46	33x	stat=stat_names, stringsAsFactors=FALSE)
47	33x	z$value <- vn[z$index]
48	33x	if (is.null(wn)) {
49	14x	z$weights <- NA
50		} else {
51	19x	z$weights <- wn[z$index]
52		}
53
54		# construct column names for each index, stat
55		# add weight layer name only if layer is ambiguously weighted
56	33x	mapply(function(stat, value, weight) {
57	123x	ret <- stat
58	123x	if (full_colnames \|\| length(value_names) > 1) {
59	108x	ret <- paste(ret, value, sep='.')
60		}
61	123x	if (.includeWeightInColName(stat) && ((full_colnames & length(weight_names) > 0)
62	123x	\|\| length(weight_names) > 1)) {
63	45x	ret <- paste(ret, weight, sep='.')
64		}
65
66	123x	return(ret)
67	33x	}, z$stat, z$value, z$weight, USE.NAMES = FALSE)
68		} else {
69	51x	stat_names
70		}
71		}
72
73		.includeWeightInColName <- function(fun) {
74	123x	.isWeighted(fun) \| fun == 'fun'
75		}
76
77		.isWeighted <- function(stat_name) {
78	193x	stat_name %in% c('weighted_mean', 'weighted_sum')
79		}
80
81		#' Compute indexes for the value and weight layers that should be
82		#' processed together
83		#'
84		#' @param num_values number of layers in value raster
85		#' @param num_weights number of layers in weighting raster
86		#' @return list with \code{values} and \code{weights} elements
87		#' providing layer indexes
88		.valueWeightIndexes <- function(num_values, num_weights) {
89	128x	if (num_weights == 0) {
90	85x	vi <- seq_len(num_values)
91	85x	wi <- NA
92	43x	} else if (num_values == num_weights) {
93		# process in parallel
94	28x	vi <- seq_len(num_values)
95	28x	wi <- seq_len(num_weights)
96	15x	} else if (num_values == 1 && num_weights > 1) {
97		# recycle values
98	7x	vi <- rep.int(1, num_weights)
99	7x	wi <- seq_len(num_weights)
100	8x	} else if (num_values > 1 && num_weights == 1) {
101		# recycle weights
102	8x	vi <- seq_len(num_values)
103	8x	wi <- rep.int(1, num_values)
104		}
105
106	128x	list(values = vi, weights = wi)
107		}
108

1		# Copyright (c) 2020 ISciences, LLC.
2		# All rights reserved.
3		#
4		# This software is licensed under the Apache License, Version 2.0 (the "License").
5		# You may not use this file except in compliance with the License. You may
6		# obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
7		#
8		# Unless required by applicable law or agreed to in writing, software
9		# distributed under the License is distributed on an "AS IS" BASIS,
10		# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11		# See the License for the specific language governing permissions and
12		# limitations under the License.
13
14		if (!isGeneric("exact_resample")) {
15	7x	setGeneric("exact_resample", function(x, y, ...)
16	7x	standardGeneric("exact_resample"))
17		}
18
19		#' Resample a raster to a new grid
20		#'
21		#' @param x a \code{RasterLayer} to be resampled
22		#' @param y a \code{RasterLayer} with a grid definition to which \code{x}
23		#' should be resampled
24		#' @param fun a named summary operation to be used for the resampling
25		#' @return a resampled version of \code{x}
26		#'
27		#' @name exact_resample
28		NULL
29
30		#' @import sf
31		#' @import raster
32		#' @useDynLib exactextractr
33		#' @rdname exact_resample
34		#' @export
35		setMethod('exact_resample',
36		signature(x='RasterLayer', y='RasterLayer'),
37		function(x, y, fun) {
38	7x	if (sf::st_crs(x) != sf::st_crs(y)) {
39	3x	if (is.na(sf::st_crs(x))) {
40	1x	warning("No CRS specified for source raster; assuming it has the same CRS as destination raster.")
41	2x	} else if (is.na(sf::st_crs(y))) {
42	1x	warning("No CRS specified for destination raster; assuming it has the same CRS as source raster.")
43		} else {
44	1x	stop('Destination raster must have same CRS as source.')
45		}
46		}
47
48	6x	x <- raster::readStart(x)
49	6x	tryCatch({
50	6x	CPP_resample(x, y, fun)
51	6x	}, finally={
52	6x	raster::readStop(x)
53		})
54		})

1		// Copyright (c) 2018-2020 ISciences, LLC.
2		// All rights reserved.
3		//
4		// This software is licensed under the Apache License, Version 2.0 (the "License").
5		// You may not use this file except in compliance with the License. You may
6		// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
7		//
8		// Unless required by applicable law or agreed to in writing, software
9		// distributed under the License is distributed on an "AS IS" BASIS,
10		// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11		// See the License for the specific language governing permissions and
12		// limitations under the License.
13
14		#include <Rcpp.h>
15
16		#include "geos_r.h"
17
18		#include "raster_utils.h"
19
20		#include "exactextract/src/raster.h"
21		#include "exactextract/src/raster_cell_intersection.h"
22
23		using exactextract::RasterView;
24		using exactextract::raster_cell_intersection;
25
26		// [[Rcpp::export]]
27	12x	Rcpp::S4 CPP_coverage_fraction(Rcpp::S4 & rast, const Rcpp::RawVector & wkb, bool crop)
28		{
29	24x	GEOSAutoHandle geos;
30	36x	Rcpp::Environment raster = Rcpp::Environment::namespace_env("raster");
31	36x	Rcpp::Function rasterFn = raster["raster"];
32	36x	Rcpp::Function crsFn = raster["crs"];
33
34	12x	auto grid = make_grid(rast);
35	24x	auto coverage_fraction = raster_cell_intersection(grid, geos.handle, read_wkb(geos.handle, wkb).get());
36
37	12x	if (crop) {
38	2x	grid = coverage_fraction.grid();
39		}
40	24x	RasterView<float> coverage_view(coverage_fraction, grid);
41
42	12x	Rcpp::NumericMatrix weights{static_cast<int>(grid.rows()),
43	12x	static_cast<int>(grid.cols())};
44
45	629x	for (size_t i = 0; i < grid.rows(); i++) {
46	159018x	for (size_t j = 0; j < grid.cols(); j++) {
47	158401x	weights(i, j) = coverage_view(i, j);
48	158401x	if (!crop && std::isnan(weights(i, j))) {
49	139692x	weights(i, j) = 0;
50		}
51		}
52		}
53
54		return rasterFn(weights,
55	24x	Rcpp::Named("xmn")=grid.xmin(),
56	24x	Rcpp::Named("xmx")=grid.xmax(),
57	24x	Rcpp::Named("ymn")=grid.ymin(),
58	24x	Rcpp::Named("ymx")=grid.ymax(),
59	48x	Rcpp::Named("crs")=crsFn(rast));
60		}

1		// Copyright (c) 2018 ISciences, LLC.
2		// All rights reserved.
3		//
4		// This software is licensed under the Apache License, Version 2.0 (the "License").
5		// You may not use this file except in compliance with the License. You may
6		// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
7		//
8		// Unless required by applicable law or agreed to in writing, software
9		// distributed under the License is distributed on an "AS IS" BASIS,
10		// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11		// See the License for the specific language governing permissions and
12		// limitations under the License.
13
14		#ifndef EXACTEXTRACT_GRID_H
15		#define EXACTEXTRACT_GRID_H
16
17		#include <numeric>
18		#include <stdexcept>
19		#include <vector>
20
21		#include "box.h"
22
23		namespace exactextract {
24		struct infinite_extent {
25		static const size_t padding = 1;
26		};
27
28		struct bounded_extent {
29		static const size_t padding = 0;
30		};
31
32	239x	static inline bool is_integral(double d, double tol) {
33	239x	return std::abs(d - std::round(d)) <= tol;
34		}
35
36		template<typename extent_tag>
37		class Grid {
38
39		public:
40	1833890x	Grid(const Box & extent, double dx, double dy) :
41		m_extent{extent},
42		m_dx{dx},
43		m_dy{dy},
44	1833890x	m_num_rows{2*extent_tag::padding + (extent.ymax > extent.ymin ? static_cast<size_t>(std::round((extent.ymax - extent.ymin) / dy)) : 0)},
45	3667780x	m_num_cols{2*extent_tag::padding + (extent.xmax > extent.xmin ? static_cast<size_t>(std::round((extent.xmax - extent.xmin) / dx)) : 0)}
46		{}
47
48	118119x	static Grid make_empty() {
49	118119x	return Grid({0, 0, 0, 0}, 0, 0);
50		}
51
52	4789762x	size_t get_column(double x) const {
53		if (extent_tag::padding) {
54	6380700x	if (x < m_extent.xmin)
55	858x	return 0;
56	6379842x	if (x > m_extent.xmax)
57	1122x	return m_num_cols - 1;
58	6378720x	if (x == m_extent.xmax) // special case, returning the cell for which xmax is the right
59	3190154x	return m_num_cols - 2;
60		} else {
61	1599412x	if (x < m_extent.xmin \|\| x > m_extent.xmax)
62	!	throw std::out_of_range("x");
63
64	1599412x	if (x == m_extent.xmax)
65	800356x	return m_num_cols - 1;
66		}
67
68		// Since we've already range-checked our x value, make sure that
69		// the computed column index is no greater than the column index
70		// associated with xmax.
71	7180017x	return std::min(
72	2393339x	extent_tag::padding + static_cast<size_t>(std::floor((x - m_extent.xmin) / m_dx)),
73	4786678x	get_column(m_extent.xmax));
74		}
75
76	4788530x	size_t get_row(double y) const {
77		if (extent_tag::padding) {
78	6378008x	if (y > m_extent.ymax)
79	2206x	return 0;
80	6375802x	if (y < m_extent.ymin)
81	1848x	return m_num_rows - 1;
82	6373954x	if (y == m_extent.ymin) // special case, returning the cell for which ymin is the bottom
83	3188080x	return m_num_rows - 2;
84		} else {
85	1599526x	if (y < m_extent.ymin \|\| y > m_extent.ymax)
86	!	throw std::out_of_range("y");
87
88	1599526x	if (y == m_extent.ymin)
89	800356x	return m_num_rows - 1;
90		}
91
92		// Since we've already range-checked our y value, make sure that
93		// the computed row index is no greater than the column index
94		// associated with ymin.
95	7176321x	return std::min(
96	2392107x	extent_tag::padding + static_cast<size_t>(std::floor((m_extent.ymax - y) / m_dy)),
97	4784214x	get_row(m_extent.ymin));
98		}
99
100	1035458x	bool empty() const { return m_num_rows <= 2extent_tag::padding && m_num_cols <= 2extent_tag::padding; }
101
102	10621282x	size_t rows() const { return m_num_rows; }
103
104	12042395x	size_t cols() const { return m_num_cols; }
105
106	131x	size_t size() const { return rows()*cols(); }
107
108	5314174x	double xmin() const { return m_extent.xmin; }
109
110	3808835x	double xmax() const { return m_extent.xmax; }
111
112	3667155x	double ymin() const { return m_extent.ymin; }
113
114	5549438x	double ymax() const { return m_extent.ymax; }
115
116	5830049x	double dx() const { return m_dx; }
117
118	6067104x	double dy() const { return m_dy; }
119
120	2233705x	const Box& extent() const { return m_extent; }
121
122	399922x	size_t row_offset(const Grid & other) const { return static_cast<size_t>(std::round(std::abs(other.m_extent.ymax - m_extent.ymax) / m_dy)); }
123
124	399922x	size_t col_offset(const Grid & other) const { return static_cast<size_t>(std::round(std::abs(m_extent.xmin - other.m_extent.xmin) / m_dx)); }
125
126	283x	double x_for_col(size_t col) const { return m_extent.xmin + ((col - extent_tag::padding) + 0.5) * m_dx; }
127
128	898x	double y_for_row(size_t row) const { return m_extent.ymax - ((row - extent_tag::padding) + 0.5) * m_dy; }
129
130	1192x	Grid<extent_tag> crop(const Box & b) const {
131	1192x	if (extent().intersects(b)) {
132	1109x	return shrink_to_fit(b.intersection(extent()));
133		} else {
134	83x	return make_empty();
135		}
136		}
137
138	400178x	Grid<extent_tag> shrink_to_fit(const Box & b) const {
139	400178x	if (b.xmin < m_extent.xmin \|\| b.ymin < m_extent.ymin \|\| b.xmax > m_extent.xmax \|\| b.ymax > m_extent.ymax) {
140	!	throw std::range_error("Cannot shrink extent to bounds larger than original.");
141		}
142
143	400178x	size_t col0 = get_column(b.xmin);
144	400178x	size_t row1 = get_row(b.ymax);
145
146		// Shrink xmin and ymax to fit the upper-left corner of the supplied extent
147	400178x	double snapped_xmin = m_extent.xmin + (col0 - extent_tag::padding) * m_dx;
148	400178x	double snapped_ymax = m_extent.ymax - (row1 - extent_tag::padding) * m_dy;
149
150		// Make sure x0 and y1 are within the reduced extent. Because of
151		// floating point round-off errors, this is not always the case.
152	400178x	if (b.xmin < snapped_xmin) {
153	!	snapped_xmin -= m_dx;
154	!	col0--;
155		}
156
157	400178x	if (b.ymax > snapped_ymax) {
158	!	snapped_ymax += m_dy;
159	!	row1--;
160		}
161
162	400178x	size_t col1 = get_column(b.xmax);
163	400178x	size_t row0 = get_row(b.ymin);
164
165	400178x	size_t num_rows = 1 + (row0 - row1);
166	400178x	size_t num_cols = 1 + (col1 - col0);
167
168		// If xmax or ymin falls cleanly on a cell boundary, we don't
169		// need as many rows or columns as we otherwise would, because
170		// we assume that the rightmost cell of the grid is a closed
171		// interval in x, and the lowermost cell of the grid is a
172		// closed interval in y.
173	400178x	if (num_rows > 2 && (snapped_ymax - (num_rows-1)*m_dy <= b.ymin)) {
174	370x	num_rows--;
175		}
176	400178x	if (num_cols > 2 && (snapped_xmin + (num_cols-1)*m_dx >= b.xmax)) {
177	315x	num_cols--;
178		}
179
180		// Perform offsets relative to the new xmin, ymax origin
181		// points. If this is not done, then floating point roundoff
182		// error can cause progressive shrink() calls with the same
183		// inputs to produce different results.
184	400178x	Box reduced_box = {
185		snapped_xmin,
186	400178x	std::min(snapped_ymax - num_rows * m_dy, b.ymin),
187	400178x	std::max(snapped_xmin + num_cols * m_dx, b.xmax),
188		snapped_ymax
189		};
190
191		// Fudge computed xmax and ymin, if needed, to prevent extent
192		// from growing during a shrink operation.
193	400178x	if (reduced_box.xmax > m_extent.xmax) {
194	!	if (std::round((reduced_box.xmax - reduced_box.xmin)/m_dx) ==
195	!	std::round((m_extent.xmax - reduced_box.xmin)/m_dx)) {
196	!	reduced_box.xmax = m_extent.xmax;
197		} else {
198	!	throw std::runtime_error("Shrink operation failed.");
199		}
200		}
201	400178x	if (reduced_box.ymin < m_extent.ymin) {
202	!	if (std::round((reduced_box.ymax - reduced_box.ymin)/m_dy) ==
203	!	std::round((reduced_box.ymax - m_extent.ymin)/m_dy)) {
204	!	reduced_box.ymin = m_extent.ymin;
205		} else {
206	!	throw std::runtime_error("Shrink operation failed.");
207		}
208		}
209
210	400178x	Grid<extent_tag> reduced{reduced_box, m_dx, m_dy};
211
212	!	if (b.xmin < reduced.xmin() \|\| b.ymin < reduced.ymin() \|\| b.xmax > reduced.xmax() \|\| b.ymax > reduced.ymax()) {
213	!	throw std::runtime_error("Shrink operation failed.");
214		}
215
216	800356x	return reduced;
217		}
218
219		template<typename extent_tag2>
220	60x	bool compatible_with(const Grid<extent_tag2> &b) const {
221		// Define a tolerance for grid compatibility, to be used in the following ways:
222		// Grid resolutions must be integer multiples of each other within a factor of 'compatibility_tol'
223		// Grid origin points must differ by an integer multiple of the smaller of the two grid resolutions,
224		// within a factor of 'compatibility_tol'.
225		// Perhaps it's possible to remove this hardcoded tolerance and express something in terms of
226		// std::numeric_limits<double>::epsilon(), but I wasn't able to make that work.
227	60x	constexpr double compatability_tol = 1e-6;
228
229		if (empty() \|\| b.empty()) {
230	!	return true;
231		}
232
233		// Check x-resolution compatibility
234	60x	if (!is_integral(std::max(m_dx, b.m_dx) / std::min(m_dx, b.m_dx), std::min(m_dx, b.m_dx)*compatability_tol)) {
235	!	return false;
236		}
237
238		// Check y-resolution compatibility
239	60x	if (!is_integral(std::max(m_dy, b.m_dy) / std::min(m_dy, b.m_dy), std::min(m_dy, b.m_dy)*compatability_tol)) {
240	!	return false;
241		}
242
243		// Check left-hand boundary compatibility
244	60x	if (!is_integral(std::abs(b.m_extent.xmin - m_extent.xmin) / std::min(m_dx, b.m_dx), std::min(m_dx, b.m_dx)*compatability_tol)) {
245	1x	return false;
246		}
247
248		// Check upper boundary compatibility
249	59x	if (!is_integral(std::abs(b.m_extent.ymax - m_extent.ymax) / std::min(m_dy, b.m_dy), std::min(m_dy, b.m_dy)*compatability_tol)) {
250	!	return false;
251		}
252
253	59x	return true;
254		}
255
256		template<typename extent_tag2>
257	60x	Grid<extent_tag> common_grid(const Grid<extent_tag2> &b) const {
258	60x	if (!compatible_with(b)) {
259	1x	throw std::runtime_error("Incompatible extents.");
260		}
261
262	59x	if (b.empty()) {
263	!	return *this;
264		}
265
266	59x	const double common_dx = std::min(m_dx, b.m_dx);
267	59x	const double common_dy = std::min(m_dy, b.m_dy);
268
269	59x	const double common_xmin = std::min(m_extent.xmin, b.m_extent.xmin);
270	59x	const double common_ymax = std::max(m_extent.ymax, b.m_extent.ymax);
271
272	59x	double common_xmax = std::max(m_extent.xmax, b.m_extent.xmax);
273	59x	double common_ymin = std::min(m_extent.ymin, b.m_extent.ymin);
274
275	59x	const long nx = static_cast<long>(std::round((common_xmax - common_xmin) / common_dx));
276	59x	const long ny = static_cast<long>(std::round((common_ymax - common_ymin) / common_dy));
277
278	59x	common_xmax = std::max(common_xmax, common_xmin + nx*common_dx);
279	59x	common_ymin = std::min(common_ymin, common_ymax - ny*common_dy);
280
281	59x	return {{ common_xmin, common_ymin, common_xmax, common_ymax}, common_dx, common_dy };
282		}
283
284		template<typename extent_tag2>
285		Grid<extent_tag> overlapping_grid(const Grid<extent_tag2> &b) const {
286		if (!compatible_with(b)) {
287		throw std::runtime_error("Incompatible extents.");
288		}
289
290		if (empty() \|\| b.empty()) {
291		return make_empty();
292		}
293
294		const double common_dx = std::min(m_dx, b.m_dx);
295		const double common_dy = std::min(m_dy, b.m_dy);
296
297		const double common_xmin = std::max(m_extent.xmin, b.m_extent.xmin);
298		const double common_ymax = std::min(m_extent.ymax, b.m_extent.ymax);
299
300		double common_xmax = std::min(m_extent.xmax, b.m_extent.xmax);
301		double common_ymin = std::max(m_extent.ymin, b.m_extent.ymin);
302
303		const long nx = static_cast<long>(std::round((common_xmax - common_xmin) / common_dx));
304		const long ny = static_cast<long>(std::round((common_ymax - common_ymin) / common_dy));
305
306		common_xmax = std::max(common_xmax, common_xmin + nx*common_dx);
307		common_ymin = std::min(common_ymin, common_ymax - ny*common_dy);
308
309		return {{ common_xmin, common_ymin, common_xmax, common_ymax}, common_dx, common_dy };
310		}
311
312		bool operator==(const Grid<extent_tag> &b) const {
313		return
314		m_extent == b.m_extent &&
315		m_dx == b.m_dx &&
316		m_dy == b.m_dy;
317		}
318
319		bool operator!=(const Grid<extent_tag> &b) const {
320		return !(*this == b);
321		}
322
323		private:
324		Box m_extent;
325
326		double m_dx;
327		double m_dy;
328
329		size_t m_num_rows;
330		size_t m_num_cols;
331		};
332
333		Box grid_cell(const Grid<bounded_extent> & grid, size_t row, size_t col);
334		Box grid_cell(const Grid<infinite_extent> & grid, size_t row, size_t col);
335
336		Grid<infinite_extent> make_infinite(const Grid<bounded_extent> & grid);
337		Grid<bounded_extent> make_finite(const Grid<infinite_extent> & grid);
338
339		std::vector<Grid<bounded_extent>> subdivide(const Grid<bounded_extent> & grid, size_t max_size);
340
341		template<typename T>
342		Grid<bounded_extent> common_grid(T begin, T end) {
343		if (begin == end) {
344		return Grid<bounded_extent>::make_empty();
345		} else if (std::next(begin) == end) {
346		return begin->grid();
347		}
348		return std::accumulate(
349		std::next(begin),
350		end,
351		begin->grid(),
352		[](auto& acc, auto& op) {
353		return acc.common_grid(op.grid());
354		});
355		}
356
357		}
358
359		#endif //EXACTEXTRACT_INFINITEGRID_H

1		// Copyright (c) 2018-2019 ISciences, LLC.
2		// All rights reserved.
3		//
4		// This software is licensed under the Apache License, Version 2.0 (the "License").
5		// You may not use this file except in compliance with the License. You may
6		// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
7		//
8		// Unless required by applicable law or agreed to in writing, software
9		// distributed under the License is distributed on an "AS IS" BASIS,
10		// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11		// See the License for the specific language governing permissions and
12		// limitations under the License.
13
14		#ifndef EXACTEXTRACT_BOX_H
15		#define EXACTEXTRACT_BOX_H
16
17		#include "coordinate.h"
18		#include "crossing.h"
19		#include "side.h"
20
21		#include <limits>
22
23		namespace exactextract {
24
25		struct Box {
26		double xmin;
27		double ymin;
28		double xmax;
29		double ymax;
30
31	8721020x	Box(double p_xmin, double p_ymin, double p_xmax, double p_ymax) :
32		xmin{p_xmin},
33		ymin{p_ymin},
34		xmax{p_xmax},
35	8721020x	ymax{p_ymax} {}
36
37		static Box maximum_finite() {
38		return {
39		std::numeric_limits<double>::lowest(),
40		std::numeric_limits<double>::lowest(),
41		std::numeric_limits<double>::max(),
42		std::numeric_limits<double>::max()
43		};
44		}
45
46	267x	static Box make_empty() {
47	267x	return {0, 0, 0, 0};
48		}
49
50	6372234x	double width() const {
51	6372234x	return xmax - xmin;
52		}
53
54	6372234x	double height() const {
55	6372234x	return ymax - ymin;
56		}
57
58	6366768x	double area() const {
59	6366768x	return width() * height();
60		}
61
62	2733x	double perimeter() const {
63	2733x	return 2 * width() + 2 * height();
64		}
65
66	400715x	bool intersects(const Box & other) const {
67	400715x	if (other.ymin > ymax)
68	42x	return false;
69	400673x	if (other.ymax < ymin)
70	50x	return false;
71	400623x	if (other.xmin > xmax)
72	!	return false;
73	400623x	if (other.xmax < xmin)
74	!	return false;
75
76	400623x	return true;
77		}
78
79	4098885x	Box intersection(const Box & other) const {
80		return {
81	4098885x	std::max(xmin, other.xmin),
82	4098885x	std::max(ymin, other.ymin),
83	4098885x	std::min(xmax, other.xmax),
84	4098885x	std::min(ymax, other.ymax)
85		};
86		}
87
88		Box translate(double dx, double dy) const {
89		return {
90		xmin + dx,
91		ymin + dy,
92		xmax + dx,
93		ymax + dy
94		};
95		}
96
97		Coordinate upper_left() const {
98		return Coordinate{xmin, ymax};
99		}
100
101		Coordinate upper_right() const {
102		return Coordinate{xmax, ymax};
103		}
104
105		Coordinate lower_left() const {
106		return Coordinate{xmin, ymin};
107		}
108
109		Coordinate lower_right() const {
110		return Coordinate{xmax, ymin};
111		}
112
113		Side side(const Coordinate &c) const;
114
115		Crossing crossing(const Coordinate &c1, const Coordinate &c2) const;
116
117	516987x	bool empty() const {
118	516987x	return xmin >= xmax \|\| ymin >= ymax;
119		}
120
121	3x	Box expand_to_include(const Box & other) const {
122	3x	if (empty()) {
123	!	return other;
124		}
125
126	3x	if (other.empty()) {
127	!	return *this;
128		}
129
130	3x	return { std::min(xmin, other.xmin),
131	3x	std::min(ymin, other.ymin),
132	3x	std::max(xmax, other.xmax),
133	3x	std::max(ymax, other.ymax) };
134		}
135
136		bool contains(const Box &b) const;
137
138		bool contains(const Coordinate &c) const;
139
140		bool strictly_contains(const Coordinate &c) const;
141
142	1242x	bool operator==(const Box& other) const {
143	1242x	return xmin == other.xmin && xmax == other.xmax && ymin == other.ymin && ymax == other.ymax;
144		}
145
146		friend std::ostream &operator<<(std::ostream &os, const Box &c);
147		};
148
149		std::ostream &operator<<(std::ostream &os, const Box &c);
150
151		}
152
153		#endif

1		// Copyright (c) 2020 ISciences, LLC.
2		// All rights reserved.
3		//
4		// This software is licensed under the Apache License, Version 2.0 (the "License").
5		// You may not use this file except in compliance with the License. You may
6		// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
7		//
8		// Unless required by applicable law or agreed to in writing, software
9		// distributed under the License is distributed on an "AS IS" BASIS,
10		// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11		// See the License for the specific language governing permissions and
12		// limitations under the License.
13
14		#ifndef EXACTEXTRACT_VARIANCE_H
15		#define EXACTEXTRACT_VARIANCE_H
16
17		#include <cmath>
18
19		namespace exactextract {
20
21		class WestVariance {
22		/** \brief Implements an incremental algorithm for weighted standard
23		* deviation, variance, and coefficient of variation, as described in
24		* formula WV2 of West, D.H.D. (1979) "Updating Mean and Variance
25		* Estimates: An Improved Method". Communications of the ACM 22(9).
26		*/
27
28		private:
29		double sum_w = 0;
30		double mean = 0;
31		double t = 0;
32
33		public:
34		/** \brief Update variance estimate with another value
35		*
36		* @param x value to add
37		* @param w weight of `x`
38		*/
39	602617x	void process(double x, double w) {
40	602617x	double mean_old = mean;
41
42	602617x	sum_w += w;
43	602617x	mean += (w / sum_w) * (x - mean_old);
44	602617x	t += w * (x - mean_old) * (x - mean);
45		}
46
47		/** \brief Return the population variance.
48		*/
49	3x	constexpr double variance() const {
50	3x	return t / sum_w;
51		}
52
53		/** \brief Return the population standard deviation
54		*/
55	2x	double stdev() const {
56	2x	return std::sqrt(variance());
57		}
58
59		/** \brief Return the population coefficient of variation
60		*/
61	1x	double coefficent_of_variation() const {
62	1x	return stdev() / mean;
63		}
64
65		};
66
67		}
68
69		#endif //EXACTEXTRACT_VARIANCE_H

1		// Copyright (c) 2019-2020 ISciences, LLC.
2		// All rights reserved.
3		//
4		// This software is licensed under the Apache License, Version 2.0 (the "License").
5		// You may not use this file except in compliance with the License. You may
6		// obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0.
7		//
8		// Unless required by applicable law or agreed to in writing, software
9		// distributed under the License is distributed on an "AS IS" BASIS,
10		// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11		// See the License for the specific language governing permissions and
12		// limitations under the License.
13
14		#ifndef EXACTEXTRACT_WEIGHTED_QUANTILES_H
15		#define EXACTEXTRACT_WEIGHTED_QUANTILES_H
16
17		#include <algorithm>
18		#include <cmath>
19		#include <stdexcept>
20		#include <vector>
21
22		namespace exactextract {
23
24		class WeightedQuantiles {
25		// Compute weighted quantiles based on https://stats.stackexchange.com/a/13223
26
27		public:
28
29	55x	void process(double x, double w) {
30	55x	if (w < 0) {
31	!	throw std::runtime_error("Weighted quantile calculation does not support negative weights.");
32		}
33
34	55x	if (!std::isfinite(w)) {
35	!	throw std::runtime_error("Weighted quantile does not support non-finite weights.");
36		}
37
38	55x	m_ready_to_query = false;
39
40	55x	m_elems.emplace_back(x, w);
41		}
42
43		double quantile(double q) const;
44
45		private:
46		struct elem_t {
47	65x	elem_t(double _x, double _w) : x(_x), w(_w), cumsum(0), s(0) {}
48
49		double x;
50		double w;
51		double cumsum;
52		double s;
53		};
54
55		void prepare() const;
56
57		mutable std::vector<elem_t> m_elems;
58		mutable double m_sum_w;
59		mutable bool m_ready_to_query;
60		};
61
62		}
63
64		#endif //EXACTEXTRACT_WEIGHTED_QUANTILES_H

1		//
2		// Copyright (c) 2014-2018 Martin Moene
3		//
4		// https://github.com/martinmoene/optional-lite
5		//
6		// Distributed under the Boost Software License, Version 1.0.
7		// (See accompanying file LICENSE.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
8
9		#pragma once
10
11		#ifndef NONSTD_OPTIONAL_LITE_HPP
12		#define NONSTD_OPTIONAL_LITE_HPP
13
14		#define optional_lite_MAJOR 3
15		#define optional_lite_MINOR 2
16		#define optional_lite_PATCH 0
17
18		#define optional_lite_VERSION optional_STRINGIFY(optional_lite_MAJOR) "." optional_STRINGIFY(optional_lite_MINOR) "." optional_STRINGIFY(optional_lite_PATCH)
19
20		#define optional_STRINGIFY( x ) optional_STRINGIFY_( x )
21		#define optional_STRINGIFY_( x ) #x
22
23		// optional-lite configuration:
24
25		#define optional_OPTIONAL_DEFAULT 0
26		#define optional_OPTIONAL_NONSTD 1
27		#define optional_OPTIONAL_STD 2
28
29		#if !defined( optional_CONFIG_SELECT_OPTIONAL )
30		# define optional_CONFIG_SELECT_OPTIONAL ( optional_HAVE_STD_OPTIONAL ? optional_OPTIONAL_STD : optional_OPTIONAL_NONSTD )
31		#endif
32
33		// Control presence of exception handling (try and auto discover):
34
35		#ifndef optional_CONFIG_NO_EXCEPTIONS
36		# if defined(__cpp_exceptions) \|\| defined(__EXCEPTIONS) \|\| defined(_CPPUNWIND)
37		# define optional_CONFIG_NO_EXCEPTIONS 0
38		# else
39		# define optional_CONFIG_NO_EXCEPTIONS 1
40		# endif
41		#endif
42
43		// C++ language version detection (C++20 is speculative):
44		// Note: VC14.0/1900 (VS2015) lacks too much from C++14.
45
46		#ifndef optional_CPLUSPLUS
47		# if defined(_MSVC_LANG ) && !defined(__clang__)
48		# define optional_CPLUSPLUS (_MSC_VER == 1900 ? 201103L : _MSVC_LANG )
49		# else
50		# define optional_CPLUSPLUS __cplusplus
51		# endif
52		#endif
53
54		#define optional_CPP98_OR_GREATER ( optional_CPLUSPLUS >= 199711L )
55		#define optional_CPP11_OR_GREATER ( optional_CPLUSPLUS >= 201103L )
56		#define optional_CPP11_OR_GREATER_ ( optional_CPLUSPLUS >= 201103L )
57		#define optional_CPP14_OR_GREATER ( optional_CPLUSPLUS >= 201402L )
58		#define optional_CPP17_OR_GREATER ( optional_CPLUSPLUS >= 201703L )
59		#define optional_CPP20_OR_GREATER ( optional_CPLUSPLUS >= 202000L )
60
61		// C++ language version (represent 98 as 3):
62
63		#define optional_CPLUSPLUS_V ( optional_CPLUSPLUS / 100 - (optional_CPLUSPLUS > 200000 ? 2000 : 1994) )
64
65		// Use C++17 std::optional if available and requested:
66
67		#if optional_CPP17_OR_GREATER && defined(__has_include )
68		# if __has_include( <optional> )
69		# define optional_HAVE_STD_OPTIONAL 1
70		# else
71		# define optional_HAVE_STD_OPTIONAL 0
72		# endif
73		#else
74		# define optional_HAVE_STD_OPTIONAL 0
75		#endif
76
77		#define optional_USES_STD_OPTIONAL ( (optional_CONFIG_SELECT_OPTIONAL == optional_OPTIONAL_STD) \|\| ((optional_CONFIG_SELECT_OPTIONAL == optional_OPTIONAL_DEFAULT) && optional_HAVE_STD_OPTIONAL) )
78
79		//
80		// in_place: code duplicated in any-lite, expected-lite, optional-lite, value-ptr-lite, variant-lite:
81		//
82
83		#ifndef nonstd_lite_HAVE_IN_PLACE_TYPES
84		#define nonstd_lite_HAVE_IN_PLACE_TYPES 1
85
86		// C++17 std::in_place in <utility>:
87
88		#if optional_CPP17_OR_GREATER
89
90		#include <utility>
91
92		namespace nonstd {
93
94		using std::in_place;
95		using std::in_place_type;
96		using std::in_place_index;
97		using std::in_place_t;
98		using std::in_place_type_t;
99		using std::in_place_index_t;
100
101		#define nonstd_lite_in_place_t( T) std::in_place_t
102		#define nonstd_lite_in_place_type_t( T) std::in_place_type_t<T>
103		#define nonstd_lite_in_place_index_t(K) std::in_place_index_t<K>
104
105		#define nonstd_lite_in_place( T) std::in_place_t{}
106		#define nonstd_lite_in_place_type( T) std::in_place_type_t<T>{}
107		#define nonstd_lite_in_place_index(K) std::in_place_index_t<K>{}
108
109		} // namespace nonstd
110
111		#else // optional_CPP17_OR_GREATER
112
113		#include <cstddef>
114
115		namespace nonstd {
116		namespace detail {
117
118		template< class T >
119		struct in_place_type_tag {};
120
121		template< std::size_t K >
122		struct in_place_index_tag {};
123
124		} // namespace detail
125
126		struct in_place_t {};
127
128		template< class T >
129		inline in_place_t in_place( detail::in_place_type_tag<T> /unused/ = detail::in_place_type_tag<T>() )
130		{
131		return in_place_t();
132		}
133
134		template< std::size_t K >
135		inline in_place_t in_place( detail::in_place_index_tag<K> /unused/ = detail::in_place_index_tag<K>() )
136		{
137		return in_place_t();
138		}
139
140		template< class T >
141		inline in_place_t in_place_type( detail::in_place_type_tag<T> /unused/ = detail::in_place_type_tag<T>() )
142		{
143		return in_place_t();
144		}
145
146		template< std::size_t K >
147		inline in_place_t in_place_index( detail::in_place_index_tag<K> /unused/ = detail::in_place_index_tag<K>() )
148		{
149		return in_place_t();
150		}
151
152		// mimic templated typedef:
153
154		#define nonstd_lite_in_place_t( T) nonstd::in_place_t(&)( nonstd::detail::in_place_type_tag<T> )
155		#define nonstd_lite_in_place_type_t( T) nonstd::in_place_t(&)( nonstd::detail::in_place_type_tag<T> )
156		#define nonstd_lite_in_place_index_t(K) nonstd::in_place_t(&)( nonstd::detail::in_place_index_tag<K> )
157
158		#define nonstd_lite_in_place( T) nonstd::in_place_type<T>
159		#define nonstd_lite_in_place_type( T) nonstd::in_place_type<T>
160		#define nonstd_lite_in_place_index(K) nonstd::in_place_index<K>
161
162		} // namespace nonstd
163
164		#endif // optional_CPP17_OR_GREATER
165		#endif // nonstd_lite_HAVE_IN_PLACE_TYPES
166
167		//
168		// Using std::optional:
169		//
170
171		#if optional_USES_STD_OPTIONAL
172
173		#include <optional>
174
175		namespace nonstd {
176
177		using std::optional;
178		using std::bad_optional_access;
179		using std::hash;
180
181		using std::nullopt;
182		using std::nullopt_t;
183
184		using std::operator==;
185		using std::operator!=;
186		using std::operator<;
187		using std::operator<=;
188		using std::operator>;
189		using std::operator>=;
190		using std::make_optional;
191		using std::swap;
192		}
193
194		#else // optional_USES_STD_OPTIONAL
195
196		#include <cassert>
197		#include <utility>
198
199		// optional-lite alignment configuration:
200
201		#ifndef optional_CONFIG_MAX_ALIGN_HACK
202		# define optional_CONFIG_MAX_ALIGN_HACK 0
203		#endif
204
205		#ifndef optional_CONFIG_ALIGN_AS
206		// no default, used in #if defined()
207		#endif
208
209		#ifndef optional_CONFIG_ALIGN_AS_FALLBACK
210		# define optional_CONFIG_ALIGN_AS_FALLBACK double
211		#endif
212
213		// Compiler warning suppression:
214
215		#if defined(__clang__)
216		# pragma clang diagnostic push
217		# pragma clang diagnostic ignored "-Wundef"
218		#elif defined(__GNUC__)
219		# pragma GCC diagnostic push
220		# pragma GCC diagnostic ignored "-Wundef"
221		#elif defined(_MSC_VER )
222		# pragma warning( push )
223		#endif
224
225		// half-open range [lo..hi):
226		#define optional_BETWEEN( v, lo, hi ) ( (lo) <= (v) && (v) < (hi) )
227
228		// Compiler versions:
229		//
230		// MSVC++ 6.0 _MSC_VER == 1200 (Visual Studio 6.0)
231		// MSVC++ 7.0 _MSC_VER == 1300 (Visual Studio .NET 2002)
232		// MSVC++ 7.1 _MSC_VER == 1310 (Visual Studio .NET 2003)
233		// MSVC++ 8.0 _MSC_VER == 1400 (Visual Studio 2005)
234		// MSVC++ 9.0 _MSC_VER == 1500 (Visual Studio 2008)
235		// MSVC++ 10.0 _MSC_VER == 1600 (Visual Studio 2010)
236		// MSVC++ 11.0 _MSC_VER == 1700 (Visual Studio 2012)
237		// MSVC++ 12.0 _MSC_VER == 1800 (Visual Studio 2013)
238		// MSVC++ 14.0 _MSC_VER == 1900 (Visual Studio 2015)
239		// MSVC++ 14.1 _MSC_VER >= 1910 (Visual Studio 2017)
240
241		#if defined(_MSC_VER ) && !defined(__clang__)
242		# define optional_COMPILER_MSVC_VER (_MSC_VER )
243		# define optional_COMPILER_MSVC_VERSION (_MSC_VER / 10 - 10 * ( 5 + (_MSC_VER < 1900 ) ) )
244		#else
245		# define optional_COMPILER_MSVC_VER 0
246		# define optional_COMPILER_MSVC_VERSION 0
247		#endif
248
249		#define optional_COMPILER_VERSION( major, minor, patch ) ( 10 * (10 * (major) + (minor) ) + (patch) )
250
251		#if defined(__GNUC__) && !defined(__clang__)
252		# define optional_COMPILER_GNUC_VERSION optional_COMPILER_VERSION(__GNUC__, __GNUC_MINOR__, __GNUC_PATCHLEVEL__)
253		#else
254		# define optional_COMPILER_GNUC_VERSION 0
255		#endif
256
257		#if defined(__clang__)
258		# define optional_COMPILER_CLANG_VERSION optional_COMPILER_VERSION(__clang_major__, __clang_minor__, __clang_patchlevel__)
259		#else
260		# define optional_COMPILER_CLANG_VERSION 0
261		#endif
262
263		#if optional_BETWEEN(optional_COMPILER_MSVC_VERSION, 70, 140 )
264		# pragma warning( disable: 4345 ) // initialization behavior changed
265		#endif
266
267		#if optional_BETWEEN(optional_COMPILER_MSVC_VERSION, 70, 150 )
268		# pragma warning( disable: 4814 ) // in C++14 'constexpr' will not imply 'const'
269		#endif
270
271		// Presence of language and library features:
272
273		#define optional_HAVE(FEATURE) ( optional_HAVE_##FEATURE )
274
275		#ifdef _HAS_CPP0X
276		# define optional_HAS_CPP0X _HAS_CPP0X
277		#else
278		# define optional_HAS_CPP0X 0
279		#endif
280
281		// Unless defined otherwise below, consider VC14 as C++11 for optional-lite:
282
283		#if optional_COMPILER_MSVC_VER >= 1900
284		# undef optional_CPP11_OR_GREATER
285		# define optional_CPP11_OR_GREATER 1
286		#endif
287
288		#define optional_CPP11_90 (optional_CPP11_OR_GREATER_ \|\| optional_COMPILER_MSVC_VER >= 1500)
289		#define optional_CPP11_100 (optional_CPP11_OR_GREATER_ \|\| optional_COMPILER_MSVC_VER >= 1600)
290		#define optional_CPP11_110 (optional_CPP11_OR_GREATER_ \|\| optional_COMPILER_MSVC_VER >= 1700)
291		#define optional_CPP11_120 (optional_CPP11_OR_GREATER_ \|\| optional_COMPILER_MSVC_VER >= 1800)
292		#define optional_CPP11_140 (optional_CPP11_OR_GREATER_ \|\| optional_COMPILER_MSVC_VER >= 1900)
293		#define optional_CPP11_141 (optional_CPP11_OR_GREATER_ \|\| optional_COMPILER_MSVC_VER >= 1910)
294
295		#define optional_CPP14_000 (optional_CPP14_OR_GREATER)
296		#define optional_CPP17_000 (optional_CPP17_OR_GREATER)
297
298		// Presence of C++11 language features:
299
300		#define optional_HAVE_CONSTEXPR_11 optional_CPP11_140
301		#define optional_HAVE_IS_DEFAULT optional_CPP11_140
302		#define optional_HAVE_NOEXCEPT optional_CPP11_140
303		#define optional_HAVE_NULLPTR optional_CPP11_100
304		#define optional_HAVE_REF_QUALIFIER optional_CPP11_140
305
306		// Presence of C++14 language features:
307
308		#define optional_HAVE_CONSTEXPR_14 optional_CPP14_000
309
310		// Presence of C++17 language features:
311
312		#define optional_HAVE_NODISCARD optional_CPP17_000
313
314		// Presence of C++ library features:
315
316		#define optional_HAVE_CONDITIONAL optional_CPP11_120
317		#define optional_HAVE_REMOVE_CV optional_CPP11_120
318		#define optional_HAVE_TYPE_TRAITS optional_CPP11_90
319
320		#define optional_HAVE_TR1_TYPE_TRAITS (!! optional_COMPILER_GNUC_VERSION )
321		#define optional_HAVE_TR1_ADD_POINTER (!! optional_COMPILER_GNUC_VERSION )
322
323		// C++ feature usage:
324
325		#if optional_HAVE( CONSTEXPR_11 )
326		# define optional_constexpr constexpr
327		#else
328		# define optional_constexpr /constexpr/
329		#endif
330
331		#if optional_HAVE( IS_DEFAULT )
332		# define optional_is_default = default;
333		#else
334		# define optional_is_default {}
335		#endif
336
337		#if optional_HAVE( CONSTEXPR_14 )
338		# define optional_constexpr14 constexpr
339		#else
340		# define optional_constexpr14 /constexpr/
341		#endif
342
343		#if optional_HAVE( NODISCARD )
344		# define optional_nodiscard [[nodiscard]]
345		#else
346		# define optional_nodiscard /[[nodiscard]]/
347		#endif
348
349		#if optional_HAVE( NOEXCEPT )
350		# define optional_noexcept noexcept
351		#else
352		# define optional_noexcept /noexcept/
353		#endif
354
355		#if optional_HAVE( NULLPTR )
356		# define optional_nullptr nullptr
357		#else
358		# define optional_nullptr NULL
359		#endif
360
361		#if optional_HAVE( REF_QUALIFIER )
362		// NOLINTNEXTLINE( bugprone-macro-parentheses )
363		# define optional_ref_qual &
364		# define optional_refref_qual &&
365		#else
366		# define optional_ref_qual /&/
367		# define optional_refref_qual /&&/
368		#endif
369
370		// additional includes:
371
372		#if optional_CONFIG_NO_EXCEPTIONS
373		// already included: <cassert>
374		#else
375		# include <stdexcept>
376		#endif
377
378		#if optional_CPP11_OR_GREATER
379		# include <functional>
380		#endif
381
382		#if optional_HAVE( INITIALIZER_LIST )
383		# include <initializer_list>
384		#endif
385
386		#if optional_HAVE( TYPE_TRAITS )
387		# include <type_traits>
388		#elif optional_HAVE( TR1_TYPE_TRAITS )
389		# include <tr1/type_traits>
390		#endif
391
392		// Method enabling
393
394		#if optional_CPP11_OR_GREATER
395
396		#define optional_REQUIRES_0(...) \
397		template< bool B = (__VA_ARGS__), typename std::enable_if<B, int>::type = 0 >
398
399		#define optional_REQUIRES_T(...) \
400		, typename = typename std::enable_if< (__VA_ARGS__), nonstd::optional_lite::detail::enabler >::type
401
402		#define optional_REQUIRES_R(R, ...) \
403		typename std::enable_if< (__VA_ARGS__), R>::type
404
405		#define optional_REQUIRES_A(...) \
406		, typename std::enable_if< (__VA_ARGS__), void*>::type = nullptr
407
408		#endif
409
410		//
411		// optional:
412		//
413
414		namespace nonstd { namespace optional_lite {
415
416		namespace std11 {
417
418		#if optional_CPP11_OR_GREATER
419		using std::move;
420		#else
421		template< typename T > T & move( T & t ) { return t; }
422		#endif
423
424		#if optional_HAVE( CONDITIONAL )
425		using std::conditional;
426		#else
427		template< bool B, typename T, typename F > struct conditional { typedef T type; };
428		template< typename T, typename F > struct conditional<false, T, F> { typedef F type; };
429		#endif // optional_HAVE_CONDITIONAL
430
431		} // namespace std11
432
433		#if optional_CPP11_OR_GREATER
434
435		/// type traits C++17:
436
437		namespace std17 {
438
439		#if optional_CPP17_OR_GREATER
440
441		using std::is_swappable;
442		using std::is_nothrow_swappable;
443
444		#elif optional_CPP11_OR_GREATER
445
446		namespace detail {
447
448		using std::swap;
449
450		struct is_swappable
451		{
452		template< typename T, typename = decltype( swap( std::declval<T&>(), std::declval<T&>() ) ) >
453		static std::true_type test( int /unused/ );
454
455		template< typename >
456		static std::false_type test(...);
457		};
458
459		struct is_nothrow_swappable
460		{
461		// wrap noexcept(expr) in separate function as work-around for VC140 (VS2015):
462
463		template< typename T >
464		static constexpr bool satisfies()
465		{
466		return noexcept( swap( std::declval<T&>(), std::declval<T&>() ) );
467		}
468
469		template< typename T >
470		static auto test( int /unused/ ) -> std::integral_constant<bool, satisfies<T>()>{}
471
472		template< typename >
473		static auto test(...) -> std::false_type;
474		};
475
476		} // namespace detail
477
478		// is [nothow] swappable:
479
480		template< typename T >
481		struct is_swappable : decltype( detail::is_swappable::test<T>(0) ){};
482
483		template< typename T >
484		struct is_nothrow_swappable : decltype( detail::is_nothrow_swappable::test<T>(0) ){};
485
486		#endif // optional_CPP17_OR_GREATER
487
488		} // namespace std17
489
490		/// type traits C++20:
491
492		namespace std20 {
493
494		template< typename T >
495		struct remove_cvref
496		{
497		typedef typename std::remove_cv< typename std::remove_reference<T>::type >::type type;
498		};
499
500		} // namespace std20
501
502		#endif // optional_CPP11_OR_GREATER
503
504		/// class optional
505
506		template< typename T >
507		class optional;
508
509		namespace detail {
510
511		// for optional_REQUIRES_T
512
513		#if optional_CPP11_OR_GREATER
514		enum class enabler{};
515		#endif
516
517		// C++11 emulation:
518
519		struct nulltype{};
520
521		template< typename Head, typename Tail >
522		struct typelist
523		{
524		typedef Head head;
525		typedef Tail tail;
526		};
527
528		#if optional_CONFIG_MAX_ALIGN_HACK
529
530		// Max align, use most restricted type for alignment:
531
532		#define optional_UNIQUE( name ) optional_UNIQUE2( name, __LINE__ )
533		#define optional_UNIQUE2( name, line ) optional_UNIQUE3( name, line )
534		#define optional_UNIQUE3( name, line ) name ## line
535
536		#define optional_ALIGN_TYPE( type ) \
537		type optional_UNIQUE( _t ); struct_t< type > optional_UNIQUE( _st )
538
539		template< typename T >
540		struct struct_t { T _; };
541
542		union max_align_t
543		{
544		optional_ALIGN_TYPE( char );
545		optional_ALIGN_TYPE( short int );
546		optional_ALIGN_TYPE( int );
547		optional_ALIGN_TYPE( long int );
548		optional_ALIGN_TYPE( float );
549		optional_ALIGN_TYPE( double );
550		optional_ALIGN_TYPE( long double );
551		optional_ALIGN_TYPE( char * );
552		optional_ALIGN_TYPE( short int * );
553		optional_ALIGN_TYPE( int * );
554		optional_ALIGN_TYPE( long int * );
555		optional_ALIGN_TYPE( float * );
556		optional_ALIGN_TYPE( double * );
557		optional_ALIGN_TYPE( long double * );
558		optional_ALIGN_TYPE( void * );
559
560		#ifdef HAVE_LONG_LONG
561		optional_ALIGN_TYPE( long long );
562		#endif
563
564		struct Unknown;
565
566		Unknown ( * optional_UNIQUE(_) )( Unknown );
567		Unknown * Unknown::* optional_UNIQUE(_);
568		Unknown ( Unknown::* optional_UNIQUE(_) )( Unknown );
569
570		struct_t< Unknown ( * )( Unknown) > optional_UNIQUE(_);
571		struct_t< Unknown * Unknown::* > optional_UNIQUE(_);
572		struct_t< Unknown ( Unknown::* )(Unknown) > optional_UNIQUE(_);
573		};
574
575		#undef optional_UNIQUE
576		#undef optional_UNIQUE2
577		#undef optional_UNIQUE3
578
579		#undef optional_ALIGN_TYPE
580
581		#elif defined( optional_CONFIG_ALIGN_AS ) // optional_CONFIG_MAX_ALIGN_HACK
582
583		// Use user-specified type for alignment:
584
585		#define optional_ALIGN_AS( unused ) \
586		optional_CONFIG_ALIGN_AS
587
588		#else // optional_CONFIG_MAX_ALIGN_HACK
589
590		// Determine POD type to use for alignment:
591
592		#define optional_ALIGN_AS( to_align ) \
593		typename type_of_size< alignment_types, alignment_of< to_align >::value >::type
594
595		template< typename T >
596		struct alignment_of;
597
598		template< typename T >
599		struct alignment_of_hack
600		{
601		char c;
602		T t;
603		alignment_of_hack();
604		};
605
606		template< size_t A, size_t S >
607		struct alignment_logic
608		{
609		enum { value = A < S ? A : S };
610		};
611
612		template< typename T >
613		struct alignment_of