Feature: CLUMPY style allocation

DESCRIPTION

@@ -27,11 +27,13 @@ Imports:
     terra (>= 1.8-42)
 Suggests:
     base64enc,
+    gifski,
     knitr,
     mlr3viz,
     processx,
     quarto,
     ranger,
+    rpart,
     tinytest
 VignetteBuilder: quarto
 Config/testthat/edition: 3
@@ -42,6 +44,7 @@ Collate:
     'alloc_params_t.R'
     'trans_models_t.R'
     'alloc_dinamica.R'
+    'alloc_clumpy.R'
     'coords_t.R'
     'parquet_db_utils.R'
     'parquet_db.R'

R/RcppExports.R

@@ -1,6 +1,128 @@
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
+#' Rook-adjacency neighbour indices for a raster (C++)
+#'
+#' @param nrow,ncol Raster dimensions.
+#' @return Named list `above`/`below`/`left`/`right`, each a 1-based cell index
+#'   per cell (row-major) with 0 meaning "no neighbour" (edge).
+#' @keywords internal
+raster_neighbors_cpp <- function(nrow, ncol) {
+    .Call(`_evoland_raster_neighbors_cpp`, nrow, ncol)
+}
+
+#' Multinomial Sampling Test (MuST) in C++
+#'
+#' Inverse-CDF multinomial draw of a final state per cell (Mazy 2022,
+#' Appendix 3.B).  This is the same test the reference `clumpy` implementation
+#' calls the "generalized allocation rejection test" (GART); the thesis itself
+#' only uses the MuST name, so we follow that here.  NaN and negative
+#' probabilities are clamped to 0 (matching the reference `clumpy` Python).
+#'
+#' @param P Numeric matrix (n_cells x n_states); each row should sum to ~1
+#'   (include the "stay" column).
+#' @param states Integer vector of length `ncol(P)` giving the state id of each
+#'   column.
+#' @param u Optional NumericVector of length `nrow(P)` of uniform draws in
+#'   [0, 1) to replay instead of drawing from R's RNG.  Lets an external uniform
+#'   stream (e.g. numpy's, from the reference `clumpy`) be replayed for an exact
+#'   cross-tool comparison of the pivot test.
+#' @return Integer vector of length `nrow(P)` with the sampled state per cell.
+#' @keywords internal
+must_cpp <- function(P, states, u = NULL) {
+    .Call(`_evoland_must_cpp`, P, states, u)
+}
+
+#' Log-normal patch-area sampler (C++)
+#'
+#' @param area_mean Mean patch area (cells); NA / <= 0 returns 1.
+#' @param area_var  Patch-area variance (cells^2); NA / <= 0 treated as 1.
+#' @return Integer >= 1.
+#' @keywords internal
+sample_lognorm_area_cpp <- function(area_mean, area_var) {
+    .Call(`_evoland_sample_lognorm_area_cpp`, area_mean, area_var)
+}
+
+#' Normal patch-area sampler (C++)
+#'
+#' @param area_mean Mean patch area (cells); NA / <= 0 returns 1.
+#' @param area_var  Patch-area variance (cells^2); sd = sqrt(area_var); the
+#'   draw is clamped to >= 1.
+#' @return Integer >= 1.
+#' @keywords internal
+sample_normal_area_cpp <- function(area_mean, area_var) {
+    .Call(`_evoland_sample_normal_area_cpp`, area_mean, area_var)
+}
+
+#' Grow a single land-use patch from a pivot cell (C++)
+#'
+#' Low-level patch grower (the building block where the mutable working layer
+#' and the immutable anterior reference are distinct), kept for direct use /
+#' unit testing.  On success the allocated cells are written back into
+#' `landscape` (set to `to_class`); on failure nothing is committed.  Neighbour
+#' vectors are 1-based with 0 == no neighbour, as produced by
+#' [raster_neighbors_cpp()].
+#'
+#' @param landscape IntegerVector of current LULC values (NA_INTEGER = no-data).
+#' @param ant_landscape IntegerVector of anterior (immutable) LULC values.
+#' @param probs NumericVector of transition probabilities (length == landscape).
+#' @param nbr_above,nbr_below,nbr_left,nbr_right Neighbour index vectors.
+#' @param pivot 1-based pivot cell index.
+#' @param target_area Target patch size (cells).
+#' @param from_class,to_class Source/target LULC classes.
+#' @param elongation Target elongation in \[0, 1\] (0 = isometric).
+#' @param ncol Raster column count.
+#' @param avoid_aggregation If TRUE, the patch is all-or-nothing and fails if it
+#'   would merge with another patch or cannot reach `target_area`.
+#' @return 1-based integer vector of allocated cell indices (incl. pivot), or
+#'   empty if the patch failed / the pivot is not an available `from_class` cell.
+#' @keywords internal
+grow_patch_cpp <- function(landscape, ant_landscape, probs, nbr_above, nbr_below, nbr_left, nbr_right, pivot, target_area, from_class, to_class, elongation, ncol, avoid_aggregation = FALSE) {
+    .Call(`_evoland_grow_patch_cpp`, landscape, ant_landscape, probs, nbr_above, nbr_below, nbr_left, nbr_right, pivot, target_area, from_class, to_class, elongation, ncol, avoid_aggregation)
+}
+
+#' Run the CLUMPY allocation routine (C++)
+#'
+#' @description
+#' Allocates LULC change for a single period.  See the file header for the
+#' uSAM vs uPAM methods and the meaning of `rarefy` / `avoid_aggregation` /
+#' `area_dist`.  The anterior reference is snapshotted internally from
+#' `landscape`, so a cell is eligible as a pivot only while it still equals its
+#' original source class (prevents a cell changing twice in one time step).
+#'
+#' @param landscape IntegerVector of the anterior LULC state (row-major,
+#'   1-based class ids, NA_INTEGER for no-data).  Not modified; a copy is
+#'   returned with the allocated changes applied.
+#' @param nrow,ncol Raster dimensions.
+#' @param trans_from,trans_to IntegerVectors (length T) of the source/target
+#'   class for each transition.  The set of anterior classes is derived from
+#'   `trans_from`.
+#' @param prob_cell,prob_value Lists of length T (one element per transition)
+#'   giving the SPARSE adjusted potentials: `prob_cell[[t]]` is an integer
+#'   vector of 1-based cell indices and `prob_value[[t]]` the matching numeric
+#'   potentials for transition t.  Cells absent from a transition read as 0.
+#' @param area_mean,area_var,elongation NumericVectors (length T) of patch
+#'   parameters per transition.
+#' @param target_rate NumericVector (length T) of the target transition rate
+#'   P(v|u) per transition (fraction of source pixels that change).  Used only
+#'   by uPAM to set the per-transition pixel quota.
+#' @param method 0 = uSAM (mono-pixel single pass), 1 = uPAM (iterative, quota).
+#' @param batch_size uPAM only: pivots attempted per MuST re-draw. `> 0` is an
+#'   explicit cap (1 = strict uPAM); `< 0` processes all candidates in one pass;
+#'   `0` auto-scales to ~1% of each class's pool (bounds MuST passes so large
+#'   rasters avoid the O(#patches x pool) cost of strict batch=1).
+#' @param rarefy If TRUE, divide pivot probabilities by `area_mean` (the
+#'   1/E(sigma) factor) so the allocated quantity of change matches the target.
+#' @param shuffle If TRUE, randomise pivot processing order.
+#' @param avoid_aggregation uPAM only: if TRUE, patches that would merge fail
+#'   and allocate nothing (clumpy GaussianPatcher semantics).
+#' @param area_dist Patch-area distribution: 0 = log-normal, 1 = normal.
+#' @return IntegerVector (length n_cells) of the posterior LULC state.
+#' @keywords internal
+allocate_clumpy_cpp <- function(landscape, nrow, ncol, trans_from, trans_to, prob_cell, prob_value, area_mean, area_var, elongation, target_rate, method, batch_size, rarefy, shuffle, avoid_aggregation, area_dist) {
+    .Call(`_evoland_allocate_clumpy_cpp`, landscape, nrow, ncol, trans_from, trans_to, prob_cell, prob_value, area_mean, area_var, elongation, target_rate, method, batch_size, rarefy, shuffle, avoid_aggregation, area_dist)
+}
+
 distance_neighbors_cpp <- function(coords_t, max_distance, quiet = FALSE) {
     .Call(`_evoland_distance_neighbors_cpp`, coords_t, max_distance, quiet)
 }