diff --git a/DESCRIPTION b/DESCRIPTION index d385caa..b4fb66c 100644 --- a/DESCRIPTION +++ b/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' diff --git a/R/RcppExports.R b/R/RcppExports.R index debbb75..754905c 100644 --- a/R/RcppExports.R +++ b/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) } diff --git a/R/alloc_clumpy.R b/R/alloc_clumpy.R new file mode 100644 index 0000000..28cbdf6 --- /dev/null +++ b/R/alloc_clumpy.R @@ -0,0 +1,292 @@ +#' CLUMPY-style Allocation Methods +#' +#' @description +#' Methods for running CLUMPY-style LULC allocation. The algorithm works in +#' three stages per period: +#' +#' 1. **Prediction** – raw transition potentials are predicted and stored in +#' `trans_pot_t` via [predict_trans_pot()]. +#' 2. **Adjustment** – the adjusted view [adjusted_trans_pot_v()] rescales +#' potentials to match target rates and closes rows to \[0, 1\]. +#' 3. **Allocation** – the whole pivot-selection + patch-growth routine runs in +#' C++ ([allocate_clumpy_cpp()]). The method is chosen automatically from the +#' patch parameters: +#' * **uSAM** (Unbiased Simple Allocation Method, Mazy sec. 3.4.1) when every +#' transition is mono-pixel (`area_mean == 1` and `area_var == 0`): one MuST +#' (Multinomial Sampling Test, Mazy App. 3.B; the same test the reference +#' `clumpy` calls "GART") pass per anterior class, each selected pivot +#' allocated as a single cell. Quantity of change is enforced in +#' expectation. +#' * **uPAM** (Unbiased Patch Allocation Method, Mazy sec. 3.4.2, Fig. 3.2) +#' otherwise: iterative MuST with a per-transition pixel quota and sampling +#' without replacement. Affordable here because evoland's potentials come +#' from a fixed fitted model, so the marginal density does not need to be +#' re-estimated between patches. +#' +#' (Multi-pixel patches require uPAM; "uSAM with patches larger than one +#' pixel" is not a valid method, hence the automatic selection rather than a +#' user switch.) +#' +#' The per-cell pivot probability is divided by the mean patch area (the +#' 1/E(sigma) factor, Mazy Fig. 3.2) so the allocated quantity of change +#' matches the target transition rate; without it allocation over-shoots by +#' roughly the mean patch size. +#' +#' @references Mazy, 2022 (\url{https://theses.hal.science/tel-04382012v1}), Ch. 3. +#' +#' @name alloc_clumpy +#' @include trans_models_t.R alloc_params_t.R alloc_dinamica.R +NULL + +# Map the patch-area distribution name to the integer code used by +# allocate_clumpy_cpp (0 = log-normal, 1 = normal). +.clumpy_area_dist_code <- function(area_dist) { + area_dist <- match.arg(area_dist, c("lognormal", "normal")) + if (area_dist == "lognormal") 0L else 1L +} + +# --------------------------------------------------------------------------- +# Single-period CLUMPY allocation +# --------------------------------------------------------------------------- + +#' @describeIn alloc_clumpy +#' Allocate LULC changes for a single period using the CLUMPY algorithm. +#' +#' @param self An [evoland_db] instance. +#' @param id_period_ant Integer anterior period ID. +#' @param id_period_post Integer posterior period ID. +#' @param anterior_rast [terra::SpatRaster] of the anterior LULC state. +#' @param select_score Character; mlr3 measure ID for model selection. +#' @param select_maximize Logical; whether to maximise `select_score`. +#' @param area_dist Character; patch-area distribution, `"lognormal"` (default) +#' or `"normal"` (Gaussian with sd = `sqrt(area_var)`, clamped to >= 1). +#' @param avoid_aggregation Logical; if `TRUE` (default) uPAM patches that would +#' merge with an existing patch fail and allocate nothing (clumpy +#' `GaussianPatcher` semantics). Ignored for the mono-pixel uSAM path. +#' @param batch_size Integer; uPAM pivots attempted per MuST re-draw. `0` +#' (default) auto-scales to ~1% of each class's source pool, bounding the +#' number of MuST passes so large rasters stay tractable; `> 0` is an explicit +#' cap (1 = strict uPAM); `< 0` processes all candidates in a single pass. +#' @return An [lulc_data_t] with the simulated posterior LULC. +#' @keywords internal +alloc_clumpy_one_period <- function( + self, + id_period_ant, + id_period_post, + anterior_rast, + select_score, + select_maximize, + area_dist = "lognormal", + avoid_aggregation = TRUE, + batch_size = 0L +) { + # TODO as in dinamica, see if we can set potentials externally so we can manipulate + # them? or does this need something more elaborate like passing in a callback? + # 1. Predict and store raw transition potentials + self$predict_trans_pot( + id_period_post = id_period_post, + select_score = select_score, + select_maximize = select_maximize + ) + + # 2. Retrieve adjusted potentials, patch params and target rates + adj_pots <- self$adjusted_trans_pot_v(id_period_post) + clumpy_params <- self$alloc_params_clumpy_v() + rates <- self$trans_rates_t[ + id_period == id_period_post, + .(id_trans, rate) + ] + + # 3. Viable transitions (stable order: by anterior class, then transition id) + viable_trans <- self$trans_meta_t[is_viable == TRUE] + data.table::setorder(viable_trans, id_lulc_anterior, id_trans) + n_trans <- nrow(viable_trans) + if (n_trans == 0L) { + stop("No viable transitions found in trans_meta_t") + } + + # 4. Raster representation (row-major, 1-based cell indices) + nrow_r <- terra::nrow(anterior_rast) + ncol_r <- terra::ncol(anterior_rast) + n_cells <- nrow_r * ncol_r + ant_vec <- as.integer(terra::values(anterior_rast)) + + # 5. id_coord <-> raster cell mapping + coords_minimal <- self$coords_minimal + xy_mat <- as.matrix(coords_minimal[, .(lon, lat)]) + cell_idx <- terra::cellFromXY(anterior_rast, xy_mat) + coord_to_cell <- stats::setNames(cell_idx, coords_minimal$id_coord) + + # 6. Sparse per-transition potential columns: for each transition, the cells + # (1-based raster index) carrying a nonzero adjusted potential and the + # matching values. The adjusted-potential table is already sparse, so we + # pass it as per-transition lists and avoid materialising a dense + # n_cells x n_trans matrix (prohibitive at >1e7 cells). + prob_cell <- vector("list", n_trans) + prob_value <- vector("list", n_trans) + for (t in seq_len(n_trans)) { + id_trans_t <- viable_trans$id_trans[t] + pots_t <- adj_pots[id_trans == id_trans_t, .(id_coord, value)] + if (nrow(pots_t) == 0L) { + prob_cell[[t]] <- integer(0) + prob_value[[t]] <- numeric(0) + next + } + cells_t <- as.integer(coord_to_cell[as.character(pots_t$id_coord)]) + ok <- !is.na(cells_t) & cells_t >= 1L & cells_t <= n_cells + prob_cell[[t]] <- cells_t[ok] + prob_value[[t]] <- as.numeric(pots_t$value[ok]) + } + + # 7. Align patch params / rates to the viable-transition order + key <- data.table::data.table(id_trans = viable_trans$id_trans) + params_aligned <- clumpy_params[key, on = "id_trans"] + rates_aligned <- rates[key, on = "id_trans"] + + area_mean <- as.numeric(params_aligned$area_mean) + area_var <- as.numeric(params_aligned$area_var) + elongation <- as.numeric(params_aligned$elongation) + elongation[is.na(elongation)] <- 0.5 + target_rate <- as.numeric(rates_aligned$rate) + target_rate[is.na(target_rate)] <- 0.0 + + # 8. Select the method from the patch parameters: every transition mono-pixel + # (area_mean == 1 & area_var == 0) -> uSAM, otherwise uPAM. + # TODO see if this can be decided on a per-transition basis? + is_mono <- all(!is.na(area_mean) & area_mean == 1 & (is.na(area_var) | area_var == 0)) + method_code <- if (is_mono) 0L else 1L + method_name <- if (is_mono) "uSAM" else "uPAM" + + message(glue::glue( + "Running CLUMPY allocation ({method_name}): ", + "period {id_period_ant} -> {id_period_post}" + )) + + # 9. Run the full allocation routine in C++ + post_vec <- allocate_clumpy_cpp( + landscape = ant_vec, + nrow = nrow_r, + ncol = ncol_r, + trans_from = as.integer(viable_trans$id_lulc_anterior), + trans_to = as.integer(viable_trans$id_lulc_posterior), + prob_cell = prob_cell, + prob_value = prob_value, + area_mean = area_mean, + area_var = area_var, + elongation = elongation, + target_rate = target_rate, + method = method_code, + batch_size = as.integer(batch_size), + rarefy = TRUE, + shuffle = TRUE, + avoid_aggregation = avoid_aggregation, + area_dist = .clumpy_area_dist_code(area_dist) + ) + + # 10. Convert result vector back to lulc_data_t + message(" Converting posterior vector to lulc_data_t...") + coord_ids <- as.integer(names(coord_to_cell)) + cell_ids <- as.integer(coord_to_cell) + valid <- { + !is.na(cell_ids) & + cell_ids >= 1L & + cell_ids <= n_cells & + !is.na(post_vec[cell_ids]) + } + + lulc_result <- data.table::data.table( + id_run = self$id_run, + id_coord = coord_ids[valid], + id_lulc = as.integer(post_vec[cell_ids[valid]]), + id_period = id_period_post + ) |> + as_lulc_data_t() + + message(glue::glue(" Allocated {nrow(lulc_result)} cells")) + + lulc_result +} + +# --------------------------------------------------------------------------- +# Multi-period orchestration +# --------------------------------------------------------------------------- + +#' @describeIn alloc_clumpy +#' Run CLUMPY-style allocation over multiple periods. +#' +#' @param self An [evoland_db] instance. +#' @param id_periods Integer vector of posterior period IDs to simulate. +#' @param select_score Character; mlr3 measure ID for model selection. +#' @param select_maximize Logical; whether to maximise `select_score`. +#' @param area_dist Character; patch-area distribution, `"lognormal"` (default) +#' or `"normal"`. +#' @param avoid_aggregation Logical; uPAM merge avoidance (default `TRUE`). +#' @param batch_size Integer; uPAM pivots attempted per MuST re-draw. `0` +#' (default) auto-scales with the source pool; see [alloc_clumpy_one_period()]. +#' @param seed Optional integer random seed for reproducibility. +alloc_clumpy <- function( + self, + id_periods, + select_score, + select_maximize, + area_dist = "lognormal", + avoid_aggregation = TRUE, + batch_size = 0L, + seed = NULL +) { + stopifnot( + "id_periods must be a numeric vector" = is.numeric(id_periods), + "id_periods must be contiguous" = all(diff(id_periods) == 1L), + "id_run must be set" = !is.null(self$id_run) + ) + area_dist <- match.arg(area_dist, c("lognormal", "normal")) + + available_periods <- self$periods_t$id_period + missing_periods <- setdiff(id_periods, available_periods) + if (length(missing_periods) > 0L) { + stop(glue::glue( + "Periods not found in periods_t: {paste(missing_periods, collapse = ', ')}" + )) + } + + if (!is.null(seed)) { + set.seed(seed) + } + + message(glue::glue( + "Starting CLUMPY allocation simulation\n", + " Periods: {paste(id_periods, collapse = ' -> ')}\n", + " Run: {self$id_run}" + )) + + current_rast <- self$lulc_data_as_rast(id_period = id_periods[1L] - 1L) + + i <- 1L + for (id_period_post in id_periods) { + id_period_ant <- id_period_post - 1L + + message(glue::glue("\n=== Iteration {i}/{length(id_periods)} ===")) + + lulc_result <- alloc_clumpy_one_period( + self = self, + id_period_ant = id_period_ant, + id_period_post = id_period_post, + anterior_rast = current_rast, + select_score = select_score, + select_maximize = select_maximize, + area_dist = area_dist, + avoid_aggregation = avoid_aggregation, + batch_size = batch_size + ) + + self$commit(lulc_result, "lulc_data_t", method = "upsert") + self$upsert_new_neighbors(id_period_post) + current_rast <- self$lulc_data_as_rast(id_period = id_period_post) + + message(glue::glue("Iteration {i} complete")) + i <- i + 1L + } + + message("CLUMPY allocation complete!") + invisible(NULL) +} diff --git a/R/alloc_dinamica.R b/R/alloc_dinamica.R index ebdb8a3..037e74f 100644 --- a/R/alloc_dinamica.R +++ b/R/alloc_dinamica.R @@ -22,9 +22,7 @@ alloc_dinamica_setup_inputs <- function( id_period_ant, id_period_post, anterior_rast, - temp_dir, - select_score, - select_maximize + temp_dir ) { # Get metadata coords_meta <- self$get_table_metadata("coords_t") @@ -102,7 +100,7 @@ alloc_dinamica_setup_inputs <- function( message(glue::glue(" Wrote anterior LULC to {basename(anterior_path)}")) - # 6. Generate probability maps + # 6. Generate probability maps from the adjusted transition potentials prob_map_dir <- file.path(temp_dir, "probability_map_dir") |> ensure_dir() @@ -110,17 +108,13 @@ alloc_dinamica_setup_inputs <- function( message(" Writing probability maps...") coords_minimal <- self$coords_minimal - trans_pots_t <- self$predict_trans_pot( - id_period_post = id_period_post, - select_score = select_score, - select_maximize = select_maximize - ) + adj_trans_pots <- self$adjusted_trans_pot_v(id_period_post) # Iterate over viable transitions and write probability maps for (i in seq_len(nrow(viable_trans))) { id_trans_sel <- viable_trans$id_trans[i] prob_spatial <- coords_minimal[ - trans_pots_t[id_trans == id_trans_sel], + adj_trans_pots[id_trans == id_trans_sel], .(lon, lat, value), on = "id_coord" ] @@ -171,15 +165,21 @@ alloc_dinamica_one_period <- function( "Running Dinamica allocation: period {id_period_ant} -> {id_period_post}" )) + # Predict and store raw transition potentials in trans_pot_t + # TODO add argument to manually predict_trans_pot if manipulation is desired + self$predict_trans_pot( + id_period_post = id_period_post, + select_score = select_score, + select_maximize = select_maximize + ) + # Set up input files input_files <- alloc_dinamica_setup_inputs( self = self, id_period_ant = id_period_ant, id_period_post = id_period_post, anterior_rast = anterior_rast, - temp_dir = iteration_dir, - select_score = select_score, - select_maximize = select_maximize + temp_dir = iteration_dir ) gc() # just in case diff --git a/R/alloc_params_t.R b/R/alloc_params_t.R index 42fe43b..119d1dd 100644 --- a/R/alloc_params_t.R +++ b/R/alloc_params_t.R @@ -12,8 +12,11 @@ #' - `id_run`: Foreign key to runs_t #' - `id_trans`: Foreign key to trans_meta_t #' - `mean_patch_size`: Mean area of new patches (in cell units) -#' - `patch_size_variance`: Standard deviation of patch area -#' - `patch_isometry`: Measure of patch shape regularity +#' - `patch_size_variance`: Variance of patch area (in cell units) +#' - `patch_elongation`: Mean patch elongation (\eqn{e = 1 - \sqrt{\lambda_2 / \lambda_1}}, +#' range 0–1); the raw shape summary from `calculate_class_stats_cpp()` +#' - `patch_isometry`: Dinamica-specific isometry parameter derived from `patch_elongation` +#' via [isometry_from_elongation()] #' - `frac_expander`: Fraction of transition cells adjacent to existing patches #' - `frac_patcher`: Fraction of transition cells forming new patches #' - `similarity`: Similarity metric for allocation parameters, see @@ -27,6 +30,7 @@ as_alloc_params_t <- function(x) { id_trans = integer(0), mean_patch_size = numeric(0), patch_size_variance = numeric(0), + patch_elongation = numeric(0), patch_isometry = numeric(0), frac_expander = numeric(0), frac_patcher = numeric(0), @@ -39,6 +43,7 @@ as_alloc_params_t <- function(x) { cast_dt_col("id_trans", "int") |> cast_dt_col("mean_patch_size", "float") |> cast_dt_col("patch_size_variance", "float") |> + cast_dt_col("patch_elongation", "float") |> cast_dt_col("patch_isometry", "float") |> cast_dt_col("frac_expander", "float") |> cast_dt_col("frac_patcher", "float") |> @@ -62,6 +67,7 @@ validate.alloc_params_t <- function(x, ...) { "id_trans", "mean_patch_size", "patch_size_variance", + "patch_elongation", "patch_isometry", "frac_expander", "frac_patcher", @@ -148,9 +154,10 @@ isometry_from_elongation <- function( #' @param id_lulc_post Integer ID of the posterior LULC class #' #' @return A named list with allocation parameters: -#' - mean_patch_size: Mean area of patches (hectares) -#' - patch_size_variance: Standard deviation of patch area (hectares) -#' - patch_isometry: Measure of patch shape regularity (0-1) +#' - mean_patch_size: Mean area of patches (cell units) +#' - patch_size_variance: Variance of patch area (cell units) +#' - patch_elongation: Mean patch elongation (raw, range 0–1) +#' - patch_isometry: Dinamica isometry derived from elongation #' - frac_expander: Fraction of transition cells adjacent to old patches in \[0, 1\] #' - frac_patcher: Fraction of transition cells forming new patches in \[0, 1\] #' @@ -176,10 +183,11 @@ compute_alloc_params_single <- function( as.integer() if (is.na(n_trans_cells) || n_trans_cells == 0) { - # No transitions occurred - return NULL or default values + # No transitions occurred - return default values return(list( mean_patch_size = 0, patch_size_variance = 0, + patch_elongation = NA_real_, patch_isometry = 0, frac_expander = 0, frac_patcher = 0 @@ -225,12 +233,15 @@ compute_alloc_params_single <- function( # cellsize = 1 because we want the patch characteristics in cell edge units calculate_class_stats_cpp(cellsize = 1) + # Raw elongation from patch_stats.cpp; used as-is for CLUMPY (elongation) + # and converted to Dinamica isometry via isometry_from_elongation(). + raw_elongation <- trans_patch_stats$patch_elongation_mean[1] + list( - # Patch parameters for now intended for the Dinamica patcher - # https://dinamicaego.com/dokuwiki/doku.php?id=patcher mean_patch_size = trans_patch_stats$patch_area_mean[1], patch_size_variance = trans_patch_stats$patch_area_variance[1], - patch_isometry = isometry_from_elongation(trans_patch_stats$patch_elongation_mean[1]), + patch_elongation = raw_elongation, + patch_isometry = isometry_from_elongation(raw_elongation), frac_expander = frac_expander, frac_patcher = frac_patcher ) @@ -316,6 +327,7 @@ create_alloc_params_t <- function(self, n_perturbations = 5L, sd = 0.05) { id_period = period_post, mean_patch_size = alloc_params$mean_patch_size, patch_size_variance = alloc_params$patch_size_variance, + patch_elongation = alloc_params$patch_elongation, patch_isometry = alloc_params$patch_isometry, frac_expander = alloc_params$frac_expander, frac_patcher = alloc_params$frac_patcher @@ -341,6 +353,7 @@ create_alloc_params_t <- function(self, n_perturbations = 5L, sd = 0.05) { .( mean_patch_size = mean_na(mean_patch_size), patch_size_variance = mean_na(patch_size_variance), + patch_elongation = mean_na(patch_elongation), patch_isometry = mean_na(patch_isometry), frac_expander = mean_na(frac_expander), frac_patcher = mean_na(frac_patcher), diff --git a/R/evoland_db.R b/R/evoland_db.R index c6dd236..a8ddbc0 100644 --- a/R/evoland_db.R +++ b/R/evoland_db.R @@ -164,6 +164,34 @@ evoland_db <- R6::R6Class( create_method_binding(alloc_dinamica) }, + #' @description Runs CLUMPY-style LULC allocation, see [alloc_clumpy()]. + #' The method (uSAM vs uPAM) is selected automatically from the patch + #' parameters: mono-pixel patches (`area_mean == 1`, `area_var == 0`) use + #' uSAM, otherwise uPAM. + #' @param id_periods Integer vector of period IDs to include in the simulation. + #' @param select_score Character string; mlr3 measure ID (e.g. `"classif.auc"`) used + #' to select model for extrapolation. + #' @param select_maximize Logical; maximize (`TRUE`) or minimize (`FALSE`) the score. + #' @param area_dist Character; patch-area distribution, `"lognormal"` + #' (default) or `"normal"`. See [alloc_clumpy()]. + #' @param avoid_aggregation Logical; if `TRUE` (default) uPAM patches that + #' would merge fail and allocate nothing. Ignored for uSAM. + #' @param batch_size Integer; uPAM pivots attempted per MuST re-draw. `0` + #' (default) auto-scales with the source pool; `> 0` is an explicit cap + #' (`1` = strict uPAM); `< 0` = all candidates in one pass. Ignored for uSAM. + #' @param seed Optional integer random seed for reproducibility. + alloc_clumpy = function( + id_periods, + select_score, + select_maximize, + area_dist = "lognormal", + avoid_aggregation = TRUE, + batch_size = 0L, + seed = NULL + ) { + create_method_binding(alloc_clumpy) + }, + #' @description #' Evaluates allocation parameters in dinamica, see [eval_alloc_params_t()] #' @param select_score Character string; mlr3 measure ID (e.g. `"classif.auc"`) used @@ -266,7 +294,10 @@ evoland_db <- R6::R6Class( }, #' @description - #' Predict the transition potential for a given period, see [trans_pot_t()] + #' Predict the raw transition potential for a given period and store in + #' `trans_pot_t`, see [predict_trans_pot()]. Raw potentials are per-transition + #' model probabilities (not yet allocation-ready); use [adjusted_trans_pot_v()] + #' to obtain column-scaled, row-closed values. #' @param id_period_post Integerish, posterior period of the transition potential interval #' @param select_score Character string; mlr3 measure ID (e.g. `"classif.auc"`) used #' to select model for extrapolation @@ -278,6 +309,34 @@ evoland_db <- R6::R6Class( #' @description Get the transition rates that were observed, see [trans_rates_t] get_obs_trans_rates = function() { create_method_binding(get_obs_trans_rates) + }, + + # TODO the following pattern is different from the create_method_binding used elsewhere. + # should be fixed together with the other evoland_db$set calls in evoland_db_views.R + + #' @description + #' Return transition rates formatted for Dinamica export for a specific period, + #' see [evoland_db_views]. + #' @param id_period Integer period ID for which to export rates. + trans_rates_dinamica_v = function(id_period) { + stop("implemented by evoland_db_views.R via $set()") + }, + + #' @description + #' Return allocation-ready transition potentials for a given posterior period. + #' Raw potentials stored in `trans_pot_t` are column-scaled to match target + #' transition rates and row-closed so per-cell change probabilities sum to at + #' most 1. See [evoland_db_views] for details. + #' @param id_period_post Integer posterior period ID. + adjusted_trans_pot_v = function(id_period_post) { + stop("implemented by evoland_db_views.R via $set()") + }, + + #' @description + #' Return allocation parameters in CLUMPY-compatible format (area_mean, + #' area_var, elongation per transition). See [evoland_db_views]. + alloc_params_clumpy_v = function() { + stop("implemented by evoland_db_views.R via $set()") } ), @@ -309,6 +368,10 @@ evoland_db <- R6::R6Class( trans_models_t = create_table_binding("trans_models_t", "upsert"), #' @field alloc_params_t Get or upsert [alloc_params_t] alloc_params_t = create_table_binding("alloc_params_t", "upsert"), + #' @field trans_pot_t Get or upsert raw transition potentials [trans_pot_t]. + #' These are per-transition model probabilities stored by [predict_trans_pot()]. + #' Use [adjusted_trans_pot_v()] for allocation-ready values. + trans_pot_t = create_table_binding("trans_pot_t", "upsert"), #' @field neighbors_t Get or upsert [neighbors_t] neighbors_t = create_table_binding("neighbors_t", "write_once"), #' @field reporting_t Get or upsert [reporting_t] diff --git a/R/evoland_db_views.R b/R/evoland_db_views.R index 707c6b6..5bec05d 100644 --- a/R/evoland_db_views.R +++ b/R/evoland_db_views.R @@ -18,9 +18,14 @@ #' a specific transition. Used as input to covariance filtering. #' - `trans_rates_dinamica_v(id_period)` - Returns transition rates formatted for Dinamica export #' for a specific period. +#' - `adjusted_trans_pot_v(id_period_post)` - Returns allocation-ready transition potentials: +#' column-scaled to match target transition rates, then row-closed so per-cell change +#' probabilities sum to at most 1. +#' - `alloc_params_clumpy_v()` - Returns allocation parameters in CLUMPY format +#' (area_mean, area_var, elongation per transition). #' #' @name evoland_db_views -#' @aliases lulc_meta_long_v pred_sources_v trans_v coords_minimal trans_rates_dinamica_v +#' @aliases lulc_meta_long_v pred_sources_v trans_v coords_minimal trans_rates_dinamica_v adjusted_trans_pot_v alloc_params_clumpy_v #' @include evoland_db.R NULL @@ -111,6 +116,7 @@ evoland_db$set("active", "coords_minimal", function() { evoland_db$set( "public", "trans_rates_dinamica_v", + overwrite = TRUE, function(id_period) { stopifnot( "id_period must be a single integer" = { @@ -140,3 +146,108 @@ evoland_db$set( result } ) + +# Return allocation-ready transition potentials for a given posterior period. +# +# The raw potentials stored in trans_pot_t are per-transition MLR3 model +# probabilities that are NOT calibrated to the target transition demand. +# This view applies two adjustments (c.f. Mazy, 2022, section 2.5): +# +# 1. Column scaling: each transition's raw potentials are multiplied by +# rate / mean_potential so that the column mean matches the target +# transition rate from trans_rates_t. +# +# 2. Row closure: where the column-scaled probabilities for a cell sum to +# more than 1, all values for that cell are divided by the row sum. +# The implicit "no-change" probability is (1 - sum of stored values). +# +# id_period_post - integer posterior period ID +evoland_db$set( + "public", + "adjusted_trans_pot_v", + overwrite = TRUE, + function(id_period_post) { + stopifnot( + "id_period_post must be a single integer" = { + length(id_period_post) == 1L && id_period_post == as.integer(id_period_post) + } + ) + + pot_read_expr <- self$get_read_expr("trans_pot_t") + rates_read_expr <- self$get_read_expr("trans_rates_t") + + self$get_query(glue::glue( + r"{ + with raw as ( + select + t.id_trans, + t.id_coord, + t.id_period_post, + t.value, + r.rate, + avg(t.value) over (partition by t.id_trans) as mean_value + from {pot_read_expr} t + join {rates_read_expr} r + on t.id_trans = r.id_trans + and r.id_period = t.id_period_post + where t.id_period_post = {id_period_post} + ), + scaled as ( + select + id_trans, + id_coord, + id_period_post, + case + when mean_value > 0 then value * rate / mean_value + else 0.0 + end as scaled_value + from raw + ), + closed as ( + select + id_trans, + id_coord, + id_period_post, + case + when sum(scaled_value) over (partition by id_coord) > 1.0 + then scaled_value / sum(scaled_value) over (partition by id_coord) + else scaled_value + end as value + from scaled + ) + select id_trans, id_coord, id_period_post, value + from closed + }" + )) + } +) + +# Return allocation parameters in CLUMPY-compatible format. +# +# Maps the raw patch statistics stored in alloc_params_t to the three +# parameters consumed by the CLUMPY patcher: +# - area_mean <- mean_patch_size (mean cells per patch) +# - area_var <- patch_size_variance (variance, cell^2) +# - elongation <- patch_elongation (1 - sqrt(lambda_min/lambda_max)) +# +# Uses the active id_run / run lineage via get_read_expr. +evoland_db$set( + "public", + "alloc_params_clumpy_v", + overwrite = TRUE, + function() { + params_read_expr <- self$get_read_expr("alloc_params_t") + + self$get_query(glue::glue( + r"{ + select + id_run, + id_trans, + mean_patch_size as area_mean, + patch_size_variance as area_var, + patch_elongation as elongation + from {params_read_expr} + }" + )) + } +) diff --git a/R/trans_models_t.R b/R/trans_models_t.R index 5a6251f..aa0180b 100644 --- a/R/trans_models_t.R +++ b/R/trans_models_t.R @@ -219,7 +219,7 @@ fit_full_model_worker <- function(item, db, learner = NULL) { ) learner_params_val <- if (length(learner_params_val) == 0L) NULL else learner_params_val learner_spec_blob <- qs2::qs_serialize(trained_learner$clone(deep = TRUE)$reset()) - crossval_score_val <- list(list()) + crossval_score_val <- list(list(no.crossval = 1)) crossval_predictions_val <- list(NULL) } else { # Score-select mode: reconstruct from learner_spec; fall back to do.call diff --git a/R/trans_pot_t.R b/R/trans_pot_t.R index 41a81be..42b5010 100644 --- a/R/trans_pot_t.R +++ b/R/trans_pot_t.R @@ -87,13 +87,16 @@ print.trans_pot_t <- function(x, nrow = 10, ...) { } -#' @describeIn trans_pot_t For each viable transition, predict the transition potential -#' for a given period, with cumulative probabilities for a single id_coord capped to 1; -#' returns a `trans_pot_t` object +#' @describeIn trans_pot_t For each viable transition, predict the raw transition +#' potential for a given period and store it in `trans_pot_t` in the database. +#' Raw potentials are per-transition MLR3 model probabilities; they are **not** +#' yet allocation-ready (not column-scaled to target rates, not row-closed). +#' Use [adjusted_trans_pot_v()] to obtain allocation-ready values. #' @param self an [evoland_db] instance #' @param id_period_post scalar integerish, passed to `self$pred_data_wide_v()` #' @param select_score character scalar, name of score/measure to identify best fitting model #' @param select_maximize logical scalar, whether to maximize or minimize `select_score` +#' @return A `trans_pot_t` object (invisibly); the same data are committed to the DB. predict_trans_pot <- function( self, id_period_post, @@ -103,6 +106,24 @@ predict_trans_pot <- function( # TODO parallelize viable_trans <- self$trans_meta_t[is_viable == TRUE] + # Fail early with an actionable message naming missing models + modeled_ids <- self$get_query(glue::glue( + r"[ + select distinct id_trans + from {self$get_read_expr("trans_models_t")} + where learner_full is not null + ]" + ))$id_trans + + missing_models <- setdiff(viable_trans$id_trans, modeled_ids) + if (length(missing_models) > 0L) { + stop(glue::glue( + "No fitted model for viable transition(s): {toString(sort(missing_models))}. ", + "Every transition with is_viable == TRUE must have a non-null learner_full in ", + "trans_models_t." + )) + } + gather <- list() message(glue::glue("Predicting transition potential for {nrow(viable_trans)} transitions")) @@ -151,7 +172,6 @@ predict_trans_pot <- function( # Ensure probabilities are in [0, 1] probs <- pmax(0, pmin(1, probs)) - # Create a data.table with id_coord and probability gather[[id_trans]] <- data.table::data.table( id_trans = id_trans, id_coord = pred_data_post$id_coord, @@ -159,20 +179,13 @@ predict_trans_pot <- function( ) } - # normalize probabilities if they exceed 1 per id_coord - normalized <- - data.table::rbindlist(gather)[, - tot_pot := sum(value), - by = id_coord - ][ - tot_pot > 1, - value := value / tot_pot - ][, - `:=`( - tot_pot = NULL, - id_period_post = id_period_post - ) - ] - - as_trans_pot_t(normalized) + result <- data.table::rbindlist(gather)[, id_period_post := id_period_post] + + trans_pot <- as_trans_pot_t(result) + + # Store raw potentials in the DB so that adjusted_trans_pot_v() and allocation + # backends can retrieve them without re-running the models. + self$commit(trans_pot, "trans_pot_t", method = "upsert") + + invisible(trans_pot) } diff --git a/R/trans_preds_t.R b/R/trans_preds_t.R index 39de2c9..efcfb35 100644 --- a/R/trans_preds_t.R +++ b/R/trans_preds_t.R @@ -159,7 +159,7 @@ pred_filter_worker <- function(item, db, filter, ordered_pred_data = FALSE) { return(scores_dt) }, error = function(e) { - warning(glue::glue("Error processing id_trans?={id_trans}: {e$message}")) + warning(glue::glue("Error processing id_trans={id_trans}: {e$message}")) item[[filter_id]] <- NA_real_ return(item) } @@ -206,6 +206,6 @@ get_pred_filter_score <- function( filter = filter, ordered_pred_data = ordered_pred_data ) |> - data.table::rbindlist() |> + data.table::rbindlist(use.names = TRUE) |> as_trans_preds_t() } diff --git a/inst/tinytest/test_alloc_clumpy.R b/inst/tinytest/test_alloc_clumpy.R new file mode 100644 index 0000000..c7ecb1b --- /dev/null +++ b/inst/tinytest/test_alloc_clumpy.R @@ -0,0 +1,377 @@ +library(tinytest) + +# -------------------------------------------------------------------------- +# Unit tests for the CLUMPY allocation backend (all in C++) +# -------------------------------------------------------------------------- + +# --- evoland:::must_cpp() (Multinomial Sampling Test) --------------------------------------------------- + +# Simple case: 2 states, equal probability +set.seed(42L) +P <- matrix(c(0.5, 0.5), nrow = 1L, ncol = 2L) +result <- evoland:::must_cpp(P, c(1L, 2L)) +expect_equal(result, 2L) + +# With 100 cells, each gets exactly one state +set.seed(1L) +P100 <- matrix(rep(c(0.3, 0.7), each = 100L), nrow = 100L, ncol = 2L) +must_results_many <- evoland:::must_cpp(P100, c(10L, 20L)) +expect_equal( + must_results_many, + # fmt: skip + c( + 10, 20, 20, 20, 10, 20, 20, 20, 20, 10, 10, 10, 20, 20, 20, 20, 20, 20, 20, + 20, 20, 10, 20, 10, 10, 20, 10, 20, 20, 20, 20, 20, 20, 10, 20, 20, 20, 10, + 20, 20, 20, 20, 20, 20, 20, 20, 10, 20, 20, 20, 20, 20, 20, 10, 10, 10, 20, + 20, 20, 20, 20, 10, 20, 20, 20, 10, 20, 20, 10, 20, 20, 20, 20, 20, 20, 20, + 20, 20, 20, 20, 20, 20, 20, 20, 20, 10, 20, 10, 10, 10, 10, 10, 20, 20, 20, + 20, 20, 20, 20, 20 + ) +) + + +# "Stay" column: assign from_class when u > cumsum of all change probs +P_stay <- matrix(c(0.0, 0.0, 1.0), nrow = 1L, ncol = 3L) # all stay +res_stay <- evoland:::must_cpp(P_stay, c(1L, 2L, 3L)) +expect_equal(res_stay, 3L) + +# Negative / NaN probabilities are clamped to 0 (matches reference clumpy) +P_neg <- matrix(c(-1.0, 1.0), nrow = 1L, ncol = 2L) +expect_equal(evoland:::must_cpp(P_neg, c(1L, 2L)), 2L) + +# --- area samplers ---------------------------------------------------------- + +set.seed(1L) +expect_equal(evoland:::sample_lognorm_area_cpp(area_mean = 4, area_var = 2), 3) +expect_equal(evoland:::sample_lognorm_area_cpp(0, 1), 1L) +expect_equal(evoland:::sample_lognorm_area_cpp(NA_real_, 1), 1L) + +set.seed(666L) +expect_equal(evoland:::sample_normal_area_cpp(area_mean = 4, area_var = 2), 5) +expect_equal(evoland:::sample_normal_area_cpp(0, 1), 1L) +# Normal with zero variance returns the (rounded) mean +expect_equal(evoland:::sample_normal_area_cpp(5, 0), 5L) + +# --- evoland:::raster_neighbors_cpp() --------------------------------------- + +expect_equal( + evoland:::raster_neighbors_cpp(3L, 4L), # 3-row, 4-col raster (12 cells) + list( + above = c(0L, 0L, 0L, 0L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L), + below = c(5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 0L, 0L, 0L, 0L), + left = c(0L, 1L, 2L, 3L, 0L, 5L, 6L, 7L, 0L, 9L, 10L, 11L), + right = c(2L, 3L, 4L, 0L, 6L, 7L, 8L, 0L, 10L, 11L, 12L, 0L) + ) +) + +# --- evoland:::grow_patch_cpp() --------------------------------------------- + +# 4x4 raster, all class 1, no obstacles +n <- 16L +landscape <- as.integer(rep(1L, n)) +ant_land <- as.integer(rep(1L, n)) +probs <- rep(0.8, n) +nbrs <- evoland:::raster_neighbors_cpp(4L, 4L) + +patch <- evoland:::grow_patch_cpp( + landscape = landscape, # changes landscape by reference + ant_landscape = ant_land, + probs = probs, + nbr_above = nbrs$above, + nbr_below = nbrs$below, + nbr_left = nbrs$left, + nbr_right = nbrs$right, + pivot = 1L, + target_area = 4L, + from_class = 1L, + to_class = 2L, + elongation = 0.5, + ncol = 4L +) +expect_equal(patch, c(1L, 2L, 5L, 3L)) +expect_equal( + landscape, + c(2L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L) +) + +# Pivot with wrong class -> empty result +landscape_wrong <- as.integer(rep(2L, n)) +patch_empty <- evoland:::grow_patch_cpp( + landscape = landscape_wrong, + ant_landscape = ant_land, + probs = probs, + nbr_above = nbrs$above, + nbr_below = nbrs$below, + nbr_left = nbrs$left, + nbr_right = nbrs$right, + pivot = 1L, + target_area = 4L, + from_class = 1L, + to_class = 2L, + elongation = 0.5, + ncol = 4L +) +expect_length(patch_empty, 0L) + +# avoid_aggregation: a patch that would touch an existing patch fails entirely. +# 1x4 row [1,1,2,1] (cell 3 is an existing class-2 patch, originally class 1). +nbr14 <- evoland:::raster_neighbors_cpp(1L, 4L) +land_adj <- as.integer(c(1L, 1L, 2L, 1L)) +ant_adj <- as.integer(c(1L, 1L, 1L, 1L)) +probs_adj <- rep(1.0, 4L) + +# With avoid_aggregation = TRUE, growing from cell 1 toward area 3 hits cell 3 +# (a foreign patch) and fails -> nothing allocated. +patch_agg <- evoland:::grow_patch_cpp( + landscape = land_adj, + ant_landscape = ant_adj, + probs = probs_adj, + nbr_above = nbr14$above, + nbr_below = nbr14$below, + nbr_left = nbr14$left, + nbr_right = nbr14$right, + pivot = 1L, + target_area = 3L, + from_class = 1L, + to_class = 2L, + elongation = 0.0, + ncol = 4L, + avoid_aggregation = TRUE +) +expect_length(patch_agg, 0L) + +# With avoid_aggregation = FALSE, it grows as far as it can (cells 1 and 2). +land_adj2 <- as.integer(c(1L, 1L, 2L, 1L)) +patch_noagg <- evoland:::grow_patch_cpp( + landscape = land_adj2, + ant_landscape = ant_adj, + probs = probs_adj, + nbr_above = nbr14$above, + nbr_below = nbr14$below, + nbr_left = nbr14$left, + nbr_right = nbr14$right, + pivot = 1L, + target_area = 3L, + from_class = 1L, + to_class = 2L, + elongation = 0.0, + ncol = 4L, + avoid_aggregation = FALSE +) +expect_equal(sort(patch_noagg), c(1L, 2L)) + +# --- evoland:::allocate_clumpy_cpp() ---------------------------------------- + +# Sparse potentials are passed as per-transition lists. Helper: one transition +# over all `ncell` cells with constant potential `p`. +sparse_const <- function(ncell, p) { + list(cell = list(seq_len(ncell)), value = list(rep(p, ncell))) +} + +nr <- 5L +nc <- 5L +ncell <- nr * nc +ant <- as.integer(rep(1L, ncell)) # all class 1 +sp <- sparse_const(ncell, 0.5) # one transition 1 -> 2, potential 0.5 + +# uSAM (method 0): mono-pixel single pass +set.seed(1L) +res_usam <- evoland:::allocate_clumpy_cpp( + landscape = ant, + nrow = nr, + ncol = nc, + trans_from = 1L, + trans_to = 2L, + prob_cell = sp$cell, + prob_value = sp$value, + area_mean = 1.0, + area_var = 0.0, + elongation = 0.0, + target_rate = 0.3, + method = 0L, + batch_size = 1L, + rarefy = TRUE, + shuffle = TRUE, + avoid_aggregation = FALSE, + area_dist = 0L +) +expect_equal( + res_usam, + # fmt: skip + c( + 2L, 2L, 1L, 1L, 2L, + 1L, 1L, 1L, 1L, 2L, + 2L, 2L, 1L, 2L, 1L, + 2L, 1L, 1L, 2L, 1L, + 1L, 2L, 1L, 2L, 2L + ) +) + +# uPAM (method 1): iterative, multi-pixel, quota +set.seed(1L) +res_upam <- evoland:::allocate_clumpy_cpp( + landscape = ant, + nrow = nr, + ncol = nc, + trans_from = 1L, + trans_to = 2L, + prob_cell = sp$cell, + prob_value = sp$value, + area_mean = 2.0, + area_var = 1.0, + elongation = 0.0, + target_rate = 0.3, + method = 1L, + batch_size = 1L, + rarefy = TRUE, + shuffle = TRUE, + avoid_aggregation = TRUE, + area_dist = 0L +) +expect_equal( + res_upam, + # fmt: skip + c( + 1L, 1L, 1L, 2L, 1L, + 1L, 2L, 1L, 2L, 1L, + 1L, 2L, 1L, 1L, 1L, + 1L, 1L, 1L, 1L, 2L, + 1L, 1L, 1L, 1L, 2L + ) +) + +# Deterministic forcing: potential 1 + uSAM => every source cell transitions, +# regardless of the RNG draw (the MuST rejection step is bypassed). +sp1 <- sparse_const(ncell, 1.0) +set.seed(123L) +res_forced <- evoland:::allocate_clumpy_cpp( + landscape = ant, + nrow = nr, + ncol = nc, + trans_from = 1L, + trans_to = 2L, + prob_cell = sp1$cell, + prob_value = sp1$value, + area_mean = 1.0, + area_var = 0.0, + elongation = 0.0, + target_rate = 1.0, + method = 0L, + batch_size = 1L, + rarefy = FALSE, + shuffle = TRUE, + avoid_aggregation = FALSE, + area_dist = 0L +) +expect_true(all(res_forced == 2L)) + +# Empty sparse potentials (no entries) => nothing changes +set.seed(123L) +res_zero <- evoland:::allocate_clumpy_cpp( + landscape = ant, + nrow = nr, + ncol = nc, + trans_from = 1L, + trans_to = 2L, + prob_cell = list(integer(0)), + prob_value = list(numeric(0)), + area_mean = 2.0, + area_var = 1.0, + elongation = 0.0, + target_rate = 0.3, + method = 1L, + batch_size = 1L, + rarefy = TRUE, + shuffle = TRUE, + avoid_aggregation = TRUE, + area_dist = 0L +) +expect_true(all(res_zero == 1L)) + +# avoid_aggregation: rejecting merging patches leaves fewer cells allocated. +# Deterministic setup so the comparison is a true invariant (not RNG/quota/FP +# dependent): a 1x5 row, forced potential 1 (deterministic MuST), area exactly 2 +# (normal with variance 0), no shuffle, quota = everything (target_rate 1). +# - without avoidance, patches tile the row -> all 5 cells change; +# - with avoidance, the middle patch (cell 3) would touch cell 2's patch, so it +# is rejected and cell 3 stays -> 4 cells change. +ant_row <- as.integer(rep(1L, 5L)) +row_cell <- list(1:5L) +row_val <- list(rep(1.0, 5L)) +run_row <- function(avoid_agg) { + set.seed(1L) + evoland:::allocate_clumpy_cpp( + landscape = ant_row, + nrow = 1L, + ncol = 5L, + trans_from = 1L, + trans_to = 2L, + prob_cell = row_cell, + prob_value = row_val, + area_mean = 2.0, + area_var = 0.0, + elongation = 0.0, + target_rate = 1.0, + method = 1L, + batch_size = 1L, + rarefy = FALSE, + shuffle = FALSE, + avoid_aggregation = avoid_agg, + area_dist = 1L + ) +} +expect_equal(run_row(avoid_agg = FALSE), c(2, 2, 2, 2, 2)) +expect_equal(run_row(avoid_agg = TRUE), c(2, 2, 1, 2, 2)) + +# A sparse subset: only some cells carry a potential; only those can change. +set.seed(9L) +some_cells <- c(1L, 7L, 13L, 19L, 25L) +res_subset <- evoland:::allocate_clumpy_cpp( + landscape = ant, + nrow = nr, + ncol = nc, + trans_from = 1L, + trans_to = 2L, + prob_cell = list(some_cells), + prob_value = list(rep(1.0, length(some_cells))), + area_mean = 1.0, + area_var = 0.0, + elongation = 0.0, + target_rate = 1.0, + method = 0L, + batch_size = 1L, + rarefy = FALSE, + shuffle = TRUE, + avoid_aggregation = FALSE, + area_dist = 0L +) +# forced potential 1 on exactly those cells (mono-pixel) => exactly they change +expect_equal(which(res_subset == 2L), some_cells) + +# Auto batch (batch_size = 0): scales with the pool, runs and stays valid. +big2 <- 40L +antc <- as.integer(rep(1L, big2 * big2)) +spc <- sparse_const(big2 * big2, 0.4) +set.seed(11L) +res_auto <- evoland:::allocate_clumpy_cpp( + landscape = antc, + nrow = big2, + ncol = big2, + trans_from = 1L, + trans_to = 2L, + prob_cell = spc$cell, + prob_value = spc$value, + area_mean = 4.0, + area_var = 2.0, + elongation = 0.0, + target_rate = 0.3, + method = 1L, + batch_size = 0L, + rarefy = TRUE, + shuffle = TRUE, + avoid_aggregation = TRUE, + area_dist = 0L +) +expect_length(res_auto, big2 * big2) +# quota-bounded: changed count should be near rate * pool, not wildly over +expect_true(sum(res_auto == 2L) <= ceiling(0.3 * big2 * big2) + 10L) diff --git a/inst/tinytest/test_alloc_params_t.R b/inst/tinytest/test_alloc_params_t.R index a58250a..4083ee3 100644 --- a/inst/tinytest/test_alloc_params_t.R +++ b/inst/tinytest/test_alloc_params_t.R @@ -7,10 +7,11 @@ alloc_params_t <- as_alloc_params_t(list( id_trans = 1L, mean_patch_size = 1.3, patch_size_variance = 1.4, + patch_elongation = 0.15, patch_isometry = 0.2, frac_expander = 0.8, - gof_window_size = 11, - gof_fuzzy_similarity = 0.8 + frac_patcher = 0.2, + similarity = NA_real_ )) expect_silent(alloc_params_t) @@ -67,8 +68,20 @@ params <- evoland:::compute_alloc_params_single( id_lulc_ant = 1L, id_lulc_post = 2L ) + +# Result now includes patch_elongation (raw) alongside patch_isometry (Dinamica) +expect_true("patch_elongation" %in% names(params)) +expect_true("patch_isometry" %in% names(params)) +# Patch isometry is derived from elongation via isometry_from_elongation() +expect_true(is.numeric(params$patch_isometry)) expect_equal( - params, + params[c( + "mean_patch_size", + "patch_size_variance", + "patch_isometry", + "frac_expander", + "frac_patcher" + )], list( mean_patch_size = 3, patch_size_variance = NA_real_, @@ -93,3 +106,4 @@ params_empty <- evoland:::compute_alloc_params_single( expect_equal(params_empty$mean_patch_size, 0) expect_equal(params_empty$frac_expander, 0) expect_equal(params_empty$frac_patcher, 0) +expect_true(is.na(params_empty$patch_elongation)) diff --git a/inst/tinytest/test_integ_allocation.R b/inst/tinytest/test_integ_allocation.R index f70e3f8..664a370 100644 --- a/inst/tinytest/test_integ_allocation.R +++ b/inst/tinytest/test_integ_allocation.R @@ -36,12 +36,12 @@ expect_message( "Fitting full models for" ) +# Test Dinamica backend # switching run to no 1, which is the base estimate for allocation parameters db$id_run <- 1L # no data for period 4 yet expect_equal(nrow(db$fetch("lulc_data_t", where = "id_period = 4")), 0L) -# Test alloc_dinamica with a simple two-period simulation if (Sys.which("DinamicaConsole") == "") { expect_warning( db$alloc_dinamica( @@ -69,7 +69,84 @@ if (Sys.which("DinamicaConsole") == "") { # 900 coords in period 4 (same as period 3) expect_equal(nrow(db$fetch("lulc_data_t", cols = "id_coord", where = "id_period = 4")), 900L) -# Test eval_alloc_params_t +# trans_pot_t should now be populated +expect_equal(db$row_count("trans_pot_t"), 900L) + +# adjusted_trans_pot_v should return values for period 4 +adj_pots <- db$adjusted_trans_pot_v(4L) +expect_equal(nrow(adj_pots), 900L) +expect_true(all(adj_pots$value >= 0 & adj_pots$value <= 1)) + +# Test CLUMPY with uSAM +# alloc_params_clumpy_v should return CLUMPY-format params +clumpy_params <- db$alloc_params_clumpy_v() +expect_true(nrow(clumpy_params) > 0L) +expect_equal( + c("id_run", "id_trans", "area_mean", "area_var", "elongation"), + names(clumpy_params) +) + +# CLUMPY: the method (uSAM vs uPAM) is auto-selected from the patch params. +# The test DB has multi-pixel patches, so both runs exercise uPAM; the +# mono-pixel uSAM path is covered by the unit tests in test_alloc_clumpy.R. +db$id_run <- 2L +expect_equal(nrow(db$fetch("lulc_data_t", where = "id_period = 4")), 0L) + +expect_message( + db$alloc_clumpy( + id_periods = 4L, + select_score = "classif.auc", + select_maximize = TRUE, + seed = 42L + ), + "CLUMPY allocation" +) + +# Period 4 should now be populated +expect_equal(nrow(db$fetch("lulc_data_t", cols = "id_coord", where = "id_period = 4")), 900L) + +# CLUMPY again, this time with aggregation avoidance disabled. +db$id_run <- 3L +expect_equal(nrow(db$fetch("lulc_data_t", where = "id_period = 4")), 0L) + +expect_message( + db$alloc_clumpy( + id_periods = 4L, + select_score = "classif.auc", + select_maximize = TRUE, + avoid_aggregation = FALSE, + seed = 42L + ), + "CLUMPY allocation" +) + +# Period 4 should now be populated +expect_equal(nrow(db$fetch("lulc_data_t", cols = "id_coord", where = "id_period = 4")), 900L) + +# -------------------------------------------------------------------------- +# Test error handling +# -------------------------------------------------------------------------- +# alloc() with unknown method +expect_error( + db$alloc( + method = "unknown_backend", + id_periods = 4L, + select_score = "classif.auc", + select_maximize = TRUE + ) +) + +# Non-contiguous periods +expect_error( + db$alloc_dinamica( + id_periods = c(1L, 3L) + ), + "id_periods must be contiguous" +) + +# -------------------------------------------------------------------------- +# Test eval_alloc_params_t (Dinamica-only) +# -------------------------------------------------------------------------- expect_message( db$alloc_params_t <- evaluated_params <- @@ -91,11 +168,3 @@ expect_true(all( evaluated_params$similarity <= 1 & evaluated_params$similarity >= 0 )) - -# Test error handling - invalid id_periods -expect_error( - db$alloc_dinamica( - id_periods = c(1L, 3L) # Non-contiguous - ), - "id_periods must be contiguous" -) diff --git a/inst/tinytest/test_integ_trans_models_t.R b/inst/tinytest/test_integ_trans_models_t.R index 1999fe1..5d40993 100644 --- a/inst/tinytest/test_integ_trans_models_t.R +++ b/inst/tinytest/test_integ_trans_models_t.R @@ -188,11 +188,14 @@ expect_error( ) # Test fit function that throws an error: overwrite trans_preds_t -db$trans_preds_t <- as_trans_preds_t(data.table::data.table( - id_run = 0L, - id_pred = 99999L, # non-existent predictor - id_trans = 1L -)) +db$commit( + as_trans_preds_t(data.table::data.table( + id_run = 0L, + id_pred = 99999L, # non-existent predictor + id_trans = 1L + )), + "trans_preds_t" +) expect_warning( partial_models_error <- @@ -270,3 +273,41 @@ expect_equal( rep("FALSE", 220) ) ) + +# test prediction and trans rate based potential adjustment +expect_error( + db$predict_trans_pot(id_period_post = 4, select_score = "classif.auc", select_maximize = TRUE), + "No fitted model for viable transition" +) +expect_message( + db$trans_models_t <- db$fit_full_models(learner = mlr3::lrn("classif.rpart")), + "Fitting full" +) +expect_message( + db$predict_trans_pot(id_period_post = 4, select_score = "no.crossval", select_maximize = TRUE), + "Predicting transition potential" +) +db$trans_rates_t <- + db$get_obs_trans_rates() |> + extrapolate_trans_rates(db$periods_t[is_extrapolated == TRUE]) +expect_equal( + db$row_count("trans_pot_t"), + 900L +) + +# the mean potential must scale to the overall prescribed transition rate +expect_equal( + data.table::as.data.table( + db$adjusted_trans_pot_v(4)[ + order(id_trans), + .(rate = mean(value)), + by = id_trans + ] + ), + data.table::as.data.table( + db$trans_rates_t[ + order(id_trans), + .(id_trans, rate) + ] + ) +) diff --git a/inst/tinytest/test_trans_pot_t.R b/inst/tinytest/test_trans_pot_t.R new file mode 100644 index 0000000..04f4483 --- /dev/null +++ b/inst/tinytest/test_trans_pot_t.R @@ -0,0 +1,25 @@ +library(tinytest) + +# -------------------------------------------------------------------------- +# Unit tests for trans_pot_t schema and adjusted_trans_pot_v logic +# -------------------------------------------------------------------------- + +# as_trans_pot_t: basic construction +tp <- as_trans_pot_t(data.frame( + id_trans = 1L, + id_period_post = 2L, + id_coord = 1L, + value = 0.3 +)) +expect_inherits(tp, "trans_pot_t") +expect_equal(nrow(tp), 1L) +expect_true(all(c("id_trans", "id_period_post", "id_coord", "value") %in% names(tp))) + +# Values must remain in [0, 1] – validate() should not modify valid inputs +tp_valid <- as_trans_pot_t(data.frame( + id_trans = c(1L, 2L, 1L, 2L), + id_period_post = 2L, + id_coord = c(1L, 1L, 2L, 2L), + value = c(0.4, 0.4, 0.6, 0.3) +)) +expect_true(all(tp_valid$value >= 0 & tp_valid$value <= 1)) diff --git a/man/alloc_clumpy.Rd b/man/alloc_clumpy.Rd new file mode 100644 index 0000000..7cf4e78 --- /dev/null +++ b/man/alloc_clumpy.Rd @@ -0,0 +1,104 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/alloc_clumpy.R +\name{alloc_clumpy} +\alias{alloc_clumpy} +\alias{alloc_clumpy_one_period} +\title{CLUMPY-style Allocation Methods} +\usage{ +alloc_clumpy_one_period( + self, + id_period_ant, + id_period_post, + anterior_rast, + select_score, + select_maximize, + area_dist = "lognormal", + avoid_aggregation = TRUE, + batch_size = 0L +) + +alloc_clumpy( + self, + id_periods, + select_score, + select_maximize, + area_dist = "lognormal", + avoid_aggregation = TRUE, + batch_size = 0L, + seed = NULL +) +} +\arguments{ +\item{self}{An \link{evoland_db} instance.} + +\item{id_period_ant}{Integer anterior period ID.} + +\item{id_period_post}{Integer posterior period ID.} + +\item{anterior_rast}{\link[terra:SpatRaster]{terra::SpatRaster} of the anterior LULC state.} + +\item{select_score}{Character; mlr3 measure ID for model selection.} + +\item{select_maximize}{Logical; whether to maximise \code{select_score}.} + +\item{area_dist}{Character; patch-area distribution, \code{"lognormal"} (default) +or \code{"normal"}.} + +\item{avoid_aggregation}{Logical; uPAM merge avoidance (default \code{TRUE}).} + +\item{batch_size}{Integer; uPAM pivots attempted per MuST re-draw. \code{0} +(default) auto-scales with the source pool; see \code{\link[=alloc_clumpy_one_period]{alloc_clumpy_one_period()}}.} + +\item{id_periods}{Integer vector of posterior period IDs to simulate.} + +\item{seed}{Optional integer random seed for reproducibility.} +} +\value{ +An \link{lulc_data_t} with the simulated posterior LULC. +} +\description{ +Methods for running CLUMPY-style LULC allocation. The algorithm works in +three stages per period: +\enumerate{ +\item \strong{Prediction} – raw transition potentials are predicted and stored in +\code{trans_pot_t} via \code{\link[=predict_trans_pot]{predict_trans_pot()}}. +\item \strong{Adjustment} – the adjusted view \code{\link[=adjusted_trans_pot_v]{adjusted_trans_pot_v()}} rescales +potentials to match target rates and closes rows to [0, 1]. +\item \strong{Allocation} – the whole pivot-selection + patch-growth routine runs in +C++ (\code{\link[=allocate_clumpy_cpp]{allocate_clumpy_cpp()}}). The method is chosen automatically from the +patch parameters: +\itemize{ +\item \strong{uSAM} (Unbiased Simple Allocation Method, Mazy sec. 3.4.1) when every +transition is mono-pixel (\code{area_mean == 1} and \code{area_var == 0}): one MuST +(Multinomial Sampling Test, Mazy App. 3.B; the same test the reference +\code{clumpy} calls "GART") pass per anterior class, each selected pivot +allocated as a single cell. Quantity of change is enforced in +expectation. +\item \strong{uPAM} (Unbiased Patch Allocation Method, Mazy sec. 3.4.2, Fig. 3.2) +otherwise: iterative MuST with a per-transition pixel quota and sampling +without replacement. Affordable here because evoland's potentials come +from a fixed fitted model, so the marginal density does not need to be +re-estimated between patches. +} + +(Multi-pixel patches require uPAM; "uSAM with patches larger than one +pixel" is not a valid method, hence the automatic selection rather than a +user switch.) + +The per-cell pivot probability is divided by the mean patch area (the +1/E(sigma) factor, Mazy Fig. 3.2) so the allocated quantity of change +matches the target transition rate; without it allocation over-shoots by +roughly the mean patch size. +} +} +\section{Functions}{ +\itemize{ +\item \code{alloc_clumpy_one_period()}: Allocate LULC changes for a single period using the CLUMPY algorithm. + +\item \code{alloc_clumpy()}: Run CLUMPY-style allocation over multiple periods. + +}} +\references{ +Mazy, 2022 (\url{https://theses.hal.science/tel-04382012v1}), Ch. 3. +} +\keyword{internal} diff --git a/man/alloc_dinamica.Rd b/man/alloc_dinamica.Rd index 72a9438..f122c79 100644 --- a/man/alloc_dinamica.Rd +++ b/man/alloc_dinamica.Rd @@ -12,9 +12,7 @@ alloc_dinamica_setup_inputs( id_period_ant, id_period_post, anterior_rast, - temp_dir, - select_score, - select_maximize + temp_dir ) alloc_dinamica_one_period( diff --git a/man/alloc_params_t.Rd b/man/alloc_params_t.Rd index 3c5b049..e0e3b54 100644 --- a/man/alloc_params_t.Rd +++ b/man/alloc_params_t.Rd @@ -58,8 +58,11 @@ A data.table of class "alloc_params_t" with columns: \item \code{id_run}: Foreign key to runs_t \item \code{id_trans}: Foreign key to trans_meta_t \item \code{mean_patch_size}: Mean area of new patches (in cell units) -\item \code{patch_size_variance}: Standard deviation of patch area -\item \code{patch_isometry}: Measure of patch shape regularity +\item \code{patch_size_variance}: Variance of patch area (in cell units) +\item \code{patch_elongation}: Mean patch elongation (\eqn{e = 1 - \sqrt{\lambda_2 / \lambda_1}}, +range 0–1); the raw shape summary from \code{calculate_class_stats_cpp()} +\item \code{patch_isometry}: Dinamica-specific isometry parameter derived from \code{patch_elongation} +via \code{\link[=isometry_from_elongation]{isometry_from_elongation()}} \item \code{frac_expander}: Fraction of transition cells adjacent to existing patches \item \code{frac_patcher}: Fraction of transition cells forming new patches \item \code{similarity}: Similarity metric for allocation parameters, see @@ -68,9 +71,10 @@ A data.table of class "alloc_params_t" with columns: A named list with allocation parameters: \itemize{ -\item mean_patch_size: Mean area of patches (hectares) -\item patch_size_variance: Standard deviation of patch area (hectares) -\item patch_isometry: Measure of patch shape regularity (0-1) +\item mean_patch_size: Mean area of patches (cell units) +\item patch_size_variance: Variance of patch area (cell units) +\item patch_elongation: Mean patch elongation (raw, range 0–1) +\item patch_isometry: Dinamica isometry derived from elongation \item frac_expander: Fraction of transition cells adjacent to old patches in [0, 1] \item frac_patcher: Fraction of transition cells forming new patches in [0, 1] } diff --git a/man/allocate_clumpy_cpp.Rd b/man/allocate_clumpy_cpp.Rd new file mode 100644 index 0000000..32ad3a1 --- /dev/null +++ b/man/allocate_clumpy_cpp.Rd @@ -0,0 +1,77 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{allocate_clumpy_cpp} +\alias{allocate_clumpy_cpp} +\title{Run the CLUMPY allocation routine (C++)} +\usage{ +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 +) +} +\arguments{ +\item{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.} + +\item{nrow, ncol}{Raster dimensions.} + +\item{trans_from, trans_to}{IntegerVectors (length T) of the source/target +class for each transition. The set of anterior classes is derived from +\code{trans_from}.} + +\item{prob_cell, prob_value}{Lists of length T (one element per transition) +giving the SPARSE adjusted potentials: \code{prob_cell[[t]]} is an integer +vector of 1-based cell indices and \code{prob_value[[t]]} the matching numeric +potentials for transition t. Cells absent from a transition read as 0.} + +\item{area_mean, area_var, elongation}{NumericVectors (length T) of patch +parameters per transition.} + +\item{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.} + +\item{method}{0 = uSAM (mono-pixel single pass), 1 = uPAM (iterative, quota).} + +\item{batch_size}{uPAM only: pivots attempted per MuST re-draw. \verb{> 0} is an +explicit cap (1 = strict uPAM); \verb{< 0} processes all candidates in one pass; +\code{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).} + +\item{rarefy}{If TRUE, divide pivot probabilities by \code{area_mean} (the +1/E(sigma) factor) so the allocated quantity of change matches the target.} + +\item{shuffle}{If TRUE, randomise pivot processing order.} + +\item{avoid_aggregation}{uPAM only: if TRUE, patches that would merge fail +and allocate nothing (clumpy GaussianPatcher semantics).} + +\item{area_dist}{Patch-area distribution: 0 = log-normal, 1 = normal.} +} +\value{ +IntegerVector (length n_cells) of the posterior LULC state. +} +\description{ +Allocates LULC change for a single period. See the file header for the +uSAM vs uPAM methods and the meaning of \code{rarefy} / \code{avoid_aggregation} / +\code{area_dist}. The anterior reference is snapshotted internally from +\code{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). +} +\keyword{internal} diff --git a/man/evoland_db.Rd b/man/evoland_db.Rd index 141e570..2c997a9 100644 --- a/man/evoland_db.Rd +++ b/man/evoland_db.Rd @@ -49,6 +49,10 @@ Additional methods and active bindings are added to this class in separate files \item{\code{alloc_params_t}}{Get or upsert \link{alloc_params_t}} + \item{\code{trans_pot_t}}{Get or upsert raw transition potentials \link{trans_pot_t}. +These are per-transition model probabilities stored by \code{\link[=predict_trans_pot]{predict_trans_pot()}}. +Use \code{\link[=adjusted_trans_pot_v]{adjusted_trans_pot_v()}} for allocation-ready values.} + \item{\code{neighbors_t}}{Get or upsert \link{neighbors_t}} \item{\code{reporting_t}}{Get or upsert \link{reporting_t}} @@ -65,6 +69,8 @@ Additional methods and active bindings are added to this class in separate files \subsection{Public methods}{ \itemize{ \item \href{#method-evoland_db-trans_rates_dinamica_v}{\code{evoland_db$trans_rates_dinamica_v()}} + \item \href{#method-evoland_db-adjusted_trans_pot_v}{\code{evoland_db$adjusted_trans_pot_v()}} + \item \href{#method-evoland_db-alloc_params_clumpy_v}{\code{evoland_db$alloc_params_clumpy_v()}} \item \href{#method-evoland_db-initialize}{\code{evoland_db$new()}} \item \href{#method-evoland_db-get_read_expr}{\code{evoland_db$get_read_expr()}} \item \href{#method-evoland_db-print}{\code{evoland_db$print()}} @@ -76,6 +82,7 @@ Additional methods and active bindings are added to this class in separate files \item \href{#method-evoland_db-trans_pred_data_v}{\code{evoland_db$trans_pred_data_v()}} \item \href{#method-evoland_db-pred_data_wide_v}{\code{evoland_db$pred_data_wide_v()}} \item \href{#method-evoland_db-alloc_dinamica}{\code{evoland_db$alloc_dinamica()}} + \item \href{#method-evoland_db-alloc_clumpy}{\code{evoland_db$alloc_clumpy()}} \item \href{#method-evoland_db-eval_alloc_params_t}{\code{evoland_db$eval_alloc_params_t()}} \item \href{#method-evoland_db-create_alloc_params_t}{\code{evoland_db$create_alloc_params_t()}} \item \href{#method-evoland_db-lulc_data_as_rast}{\code{evoland_db$lulc_data_as_rast()}} @@ -114,6 +121,28 @@ Additional methods and active bindings are added to this class in separate files } } +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-evoland_db-adjusted_trans_pot_v}{}}} +\subsection{\code{evoland_db$adjusted_trans_pot_v()}}{ + \subsection{Usage}{ + \if{html}{\out{
}} + \preformatted{evoland_db$adjusted_trans_pot_v(id_period_post)} + \if{html}{\out{
}} + } +} + +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-evoland_db-alloc_params_clumpy_v}{}}} +\subsection{\code{evoland_db$alloc_params_clumpy_v()}}{ + \subsection{Usage}{ + \if{html}{\out{
}} + \preformatted{evoland_db$alloc_params_clumpy_v()} + \if{html}{\out{
}} + } +} + \if{html}{\out{
}} \if{html}{\out{}} \if{latex}{\out{\hypertarget{method-evoland_db-initialize}{}}} @@ -377,6 +406,47 @@ to select model for extrapolation} } } +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-evoland_db-alloc_clumpy}{}}} +\subsection{\code{evoland_db$alloc_clumpy()}}{ + Runs CLUMPY-style LULC allocation, see \code{\link[=alloc_clumpy]{alloc_clumpy()}}. +The method (uSAM vs uPAM) is selected automatically from the patch +parameters: mono-pixel patches (\code{area_mean == 1}, \code{area_var == 0}) use +uSAM, otherwise uPAM. + \subsection{Usage}{ + \if{html}{\out{
}} + \preformatted{evoland_db$alloc_clumpy( + id_periods, + select_score, + select_maximize, + area_dist = "lognormal", + avoid_aggregation = TRUE, + batch_size = 0L, + seed = NULL +)} + \if{html}{\out{
}} + } + \subsection{Arguments}{ + \if{html}{\out{
}} + \describe{ + \item{\code{id_periods}}{Integer vector of period IDs to include in the simulation.} + \item{\code{select_score}}{Character string; mlr3 measure ID (e.g. \code{"classif.auc"}) used +to select model for extrapolation.} + \item{\code{select_maximize}}{Logical; maximize (\code{TRUE}) or minimize (\code{FALSE}) the score.} + \item{\code{area_dist}}{Character; patch-area distribution, \code{"lognormal"} +(default) or \code{"normal"}. See \code{\link[=alloc_clumpy]{alloc_clumpy()}}.} + \item{\code{avoid_aggregation}}{Logical; if \code{TRUE} (default) uPAM patches that +would merge fail and allocate nothing. Ignored for uSAM.} + \item{\code{batch_size}}{Integer; uPAM pivots attempted per MuST re-draw. \code{0} +(default) auto-scales with the source pool; \verb{> 0} is an explicit cap +(\code{1} = strict uPAM); \verb{< 0} = all candidates in one pass. Ignored for uSAM.} + \item{\code{seed}}{Optional integer random seed for reproducibility.} + } + \if{html}{\out{
}} + } +} + \if{html}{\out{
}} \if{html}{\out{}} \if{latex}{\out{\hypertarget{method-evoland_db-eval_alloc_params_t}{}}} @@ -577,7 +647,10 @@ feature selection.} \if{html}{\out{}} \if{latex}{\out{\hypertarget{method-evoland_db-predict_trans_pot}{}}} \subsection{\code{evoland_db$predict_trans_pot()}}{ - Predict the transition potential for a given period, see \code{\link[=trans_pot_t]{trans_pot_t()}} + Predict the raw transition potential for a given period and store in +\code{trans_pot_t}, see \code{\link[=predict_trans_pot]{predict_trans_pot()}}. Raw potentials are per-transition +model probabilities (not yet allocation-ready); use \code{\link[=adjusted_trans_pot_v]{adjusted_trans_pot_v()}} +to obtain column-scaled, row-closed values. \subsection{Usage}{ \if{html}{\out{
}} \preformatted{evoland_db$predict_trans_pot(id_period_post, select_score, select_maximize)} diff --git a/man/evoland_db_views.Rd b/man/evoland_db_views.Rd index 9e54dff..6aed0d6 100644 --- a/man/evoland_db_views.Rd +++ b/man/evoland_db_views.Rd @@ -7,6 +7,8 @@ \alias{trans_v} \alias{coords_minimal} \alias{trans_rates_dinamica_v} +\alias{adjusted_trans_pot_v} +\alias{alloc_params_clumpy_v} \title{Views on the evoland-plus data model} \description{ This file adds view active bindings and methods to the \code{evoland_db} class using R6's \verb{$set()} @@ -30,6 +32,11 @@ method. These provide computed views on the database without storing additional a specific transition. Used as input to covariance filtering. \item \code{trans_rates_dinamica_v(id_period)} - Returns transition rates formatted for Dinamica export for a specific period. +\item \code{adjusted_trans_pot_v(id_period_post)} - Returns allocation-ready transition potentials: +column-scaled to match target transition rates, then row-closed so per-cell change +probabilities sum to at most 1. +\item \code{alloc_params_clumpy_v()} - Returns allocation parameters in CLUMPY format +(area_mean, area_var, elongation per transition). } } diff --git a/man/grow_patch_cpp.Rd b/man/grow_patch_cpp.Rd new file mode 100644 index 0000000..506e048 --- /dev/null +++ b/man/grow_patch_cpp.Rd @@ -0,0 +1,58 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{grow_patch_cpp} +\alias{grow_patch_cpp} +\title{Grow a single land-use patch from a pivot cell (C++)} +\usage{ +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 = FALSE +) +} +\arguments{ +\item{landscape}{IntegerVector of current LULC values (NA_INTEGER = no-data).} + +\item{ant_landscape}{IntegerVector of anterior (immutable) LULC values.} + +\item{probs}{NumericVector of transition probabilities (length == landscape).} + +\item{nbr_above, nbr_below, nbr_left, nbr_right}{Neighbour index vectors.} + +\item{pivot}{1-based pivot cell index.} + +\item{target_area}{Target patch size (cells).} + +\item{from_class, to_class}{Source/target LULC classes.} + +\item{elongation}{Target elongation in [0, 1] (0 = isometric).} + +\item{ncol}{Raster column count.} + +\item{avoid_aggregation}{If TRUE, the patch is all-or-nothing and fails if it +would merge with another patch or cannot reach \code{target_area}.} +} +\value{ +1-based integer vector of allocated cell indices (incl. pivot), or +empty if the patch failed / the pivot is not an available \code{from_class} cell. +} +\description{ +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 +\code{landscape} (set to \code{to_class}); on failure nothing is committed. Neighbour +vectors are 1-based with 0 == no neighbour, as produced by +\code{\link[=raster_neighbors_cpp]{raster_neighbors_cpp()}}. +} +\keyword{internal} diff --git a/man/must_cpp.Rd b/man/must_cpp.Rd new file mode 100644 index 0000000..0706a6f --- /dev/null +++ b/man/must_cpp.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{must_cpp} +\alias{must_cpp} +\title{Multinomial Sampling Test (MuST) in C++} +\usage{ +must_cpp(P, states, u = NULL) +} +\arguments{ +\item{P}{Numeric matrix (n_cells x n_states); each row should sum to ~1 +(include the "stay" column).} + +\item{states}{Integer vector of length \code{ncol(P)} giving the state id of each +column.} + +\item{u}{Optional NumericVector of length \code{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 \code{clumpy}) be replayed for an exact +cross-tool comparison of the pivot test.} +} +\value{ +Integer vector of length \code{nrow(P)} with the sampled state per cell. +} +\description{ +Inverse-CDF multinomial draw of a final state per cell (Mazy 2022, +Appendix 3.B). This is the same test the reference \code{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 \code{clumpy} Python). +} +\keyword{internal} diff --git a/man/raster_neighbors_cpp.Rd b/man/raster_neighbors_cpp.Rd new file mode 100644 index 0000000..266a841 --- /dev/null +++ b/man/raster_neighbors_cpp.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{raster_neighbors_cpp} +\alias{raster_neighbors_cpp} +\title{Rook-adjacency neighbour indices for a raster (C++)} +\usage{ +raster_neighbors_cpp(nrow, ncol) +} +\arguments{ +\item{nrow, ncol}{Raster dimensions.} +} +\value{ +Named list \code{above}/\code{below}/\code{left}/\code{right}, each a 1-based cell index +per cell (row-major) with 0 meaning "no neighbour" (edge). +} +\description{ +Rook-adjacency neighbour indices for a raster (C++) +} +\keyword{internal} diff --git a/man/sample_lognorm_area_cpp.Rd b/man/sample_lognorm_area_cpp.Rd new file mode 100644 index 0000000..73cca54 --- /dev/null +++ b/man/sample_lognorm_area_cpp.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{sample_lognorm_area_cpp} +\alias{sample_lognorm_area_cpp} +\title{Log-normal patch-area sampler (C++)} +\usage{ +sample_lognorm_area_cpp(area_mean, area_var) +} +\arguments{ +\item{area_mean}{Mean patch area (cells); NA / <= 0 returns 1.} + +\item{area_var}{Patch-area variance (cells^2); NA / <= 0 treated as 1.} +} +\value{ +Integer >= 1. +} +\description{ +Log-normal patch-area sampler (C++) +} +\keyword{internal} diff --git a/man/sample_normal_area_cpp.Rd b/man/sample_normal_area_cpp.Rd new file mode 100644 index 0000000..a1e8060 --- /dev/null +++ b/man/sample_normal_area_cpp.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{sample_normal_area_cpp} +\alias{sample_normal_area_cpp} +\title{Normal patch-area sampler (C++)} +\usage{ +sample_normal_area_cpp(area_mean, area_var) +} +\arguments{ +\item{area_mean}{Mean patch area (cells); NA / <= 0 returns 1.} + +\item{area_var}{Patch-area variance (cells^2); sd = sqrt(area_var); the +draw is clamped to >= 1.} +} +\value{ +Integer >= 1. +} +\description{ +Normal patch-area sampler (C++) +} +\keyword{internal} diff --git a/man/trans_pot_t.Rd b/man/trans_pot_t.Rd index 06aee92..b86fe7a 100644 --- a/man/trans_pot_t.Rd +++ b/man/trans_pot_t.Rd @@ -36,6 +36,8 @@ A data.table of class "trans_pot_t" with columns: \item \code{id_coord}: Foreign key to \code{\link[=coords_t]{coords_t()}} \item \code{value}: Map of model (hyper) parameters } + +A \code{trans_pot_t} object (invisibly); the same data are committed to the DB. } \description{ Estimate transition potential at \code{id_period_post}. Based on the LULC at \code{id_period_anterior} @@ -47,8 +49,10 @@ Estimate transition potential at \code{id_period_post}. Based on the LULC at \co }} \section{Functions}{ \itemize{ -\item \code{predict_trans_pot()}: For each viable transition, predict the transition potential -for a given period, with cumulative probabilities for a single id_coord capped to 1; -returns a \code{trans_pot_t} object +\item \code{predict_trans_pot()}: For each viable transition, predict the raw transition +potential for a given period and store it in \code{trans_pot_t} in the database. +Raw potentials are per-transition MLR3 model probabilities; they are \strong{not} +yet allocation-ready (not column-scaled to target rates, not row-closed). +Use \code{\link[=adjusted_trans_pot_v]{adjusted_trans_pot_v()}} to obtain allocation-ready values. }} diff --git a/rproject.toml b/rproject.toml index 915f5fd..3315767 100644 --- a/rproject.toml +++ b/rproject.toml @@ -1,6 +1,6 @@ [project] name = "evoland-plus" -r_version = "4.6" +r_version = "4.5" repositories = [ { alias = "CRAN", url = "https://stat.ethz.ch/CRAN/" }, @@ -13,4 +13,5 @@ dependencies = [ "mirai", "httpgd", "languageserver", + "profvis", ] diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index e15d123..d57b644 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -10,6 +10,106 @@ Rcpp::Rostream& Rcpp::Rcout = Rcpp::Rcpp_cout_get(); Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); #endif +// raster_neighbors_cpp +List raster_neighbors_cpp(int nrow, int ncol); +RcppExport SEXP _evoland_raster_neighbors_cpp(SEXP nrowSEXP, SEXP ncolSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< int >::type nrow(nrowSEXP); + Rcpp::traits::input_parameter< int >::type ncol(ncolSEXP); + rcpp_result_gen = Rcpp::wrap(raster_neighbors_cpp(nrow, ncol)); + return rcpp_result_gen; +END_RCPP +} +// must_cpp +IntegerVector must_cpp(NumericMatrix P, IntegerVector states, Nullable u); +RcppExport SEXP _evoland_must_cpp(SEXP PSEXP, SEXP statesSEXP, SEXP uSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericMatrix >::type P(PSEXP); + Rcpp::traits::input_parameter< IntegerVector >::type states(statesSEXP); + Rcpp::traits::input_parameter< Nullable >::type u(uSEXP); + rcpp_result_gen = Rcpp::wrap(must_cpp(P, states, u)); + return rcpp_result_gen; +END_RCPP +} +// sample_lognorm_area_cpp +int sample_lognorm_area_cpp(double area_mean, double area_var); +RcppExport SEXP _evoland_sample_lognorm_area_cpp(SEXP area_meanSEXP, SEXP area_varSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< double >::type area_mean(area_meanSEXP); + Rcpp::traits::input_parameter< double >::type area_var(area_varSEXP); + rcpp_result_gen = Rcpp::wrap(sample_lognorm_area_cpp(area_mean, area_var)); + return rcpp_result_gen; +END_RCPP +} +// sample_normal_area_cpp +int sample_normal_area_cpp(double area_mean, double area_var); +RcppExport SEXP _evoland_sample_normal_area_cpp(SEXP area_meanSEXP, SEXP area_varSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< double >::type area_mean(area_meanSEXP); + Rcpp::traits::input_parameter< double >::type area_var(area_varSEXP); + rcpp_result_gen = Rcpp::wrap(sample_normal_area_cpp(area_mean, area_var)); + return rcpp_result_gen; +END_RCPP +} +// grow_patch_cpp +IntegerVector grow_patch_cpp(IntegerVector landscape, const IntegerVector& ant_landscape, const NumericVector& probs, const IntegerVector& nbr_above, const IntegerVector& nbr_below, const IntegerVector& nbr_left, const IntegerVector& nbr_right, int pivot, int target_area, int from_class, int to_class, double elongation, int ncol, bool avoid_aggregation); +RcppExport SEXP _evoland_grow_patch_cpp(SEXP landscapeSEXP, SEXP ant_landscapeSEXP, SEXP probsSEXP, SEXP nbr_aboveSEXP, SEXP nbr_belowSEXP, SEXP nbr_leftSEXP, SEXP nbr_rightSEXP, SEXP pivotSEXP, SEXP target_areaSEXP, SEXP from_classSEXP, SEXP to_classSEXP, SEXP elongationSEXP, SEXP ncolSEXP, SEXP avoid_aggregationSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< IntegerVector >::type landscape(landscapeSEXP); + Rcpp::traits::input_parameter< const IntegerVector& >::type ant_landscape(ant_landscapeSEXP); + Rcpp::traits::input_parameter< const NumericVector& >::type probs(probsSEXP); + Rcpp::traits::input_parameter< const IntegerVector& >::type nbr_above(nbr_aboveSEXP); + Rcpp::traits::input_parameter< const IntegerVector& >::type nbr_below(nbr_belowSEXP); + Rcpp::traits::input_parameter< const IntegerVector& >::type nbr_left(nbr_leftSEXP); + Rcpp::traits::input_parameter< const IntegerVector& >::type nbr_right(nbr_rightSEXP); + Rcpp::traits::input_parameter< int >::type pivot(pivotSEXP); + Rcpp::traits::input_parameter< int >::type target_area(target_areaSEXP); + Rcpp::traits::input_parameter< int >::type from_class(from_classSEXP); + Rcpp::traits::input_parameter< int >::type to_class(to_classSEXP); + Rcpp::traits::input_parameter< double >::type elongation(elongationSEXP); + Rcpp::traits::input_parameter< int >::type ncol(ncolSEXP); + Rcpp::traits::input_parameter< bool >::type avoid_aggregation(avoid_aggregationSEXP); + rcpp_result_gen = Rcpp::wrap(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)); + return rcpp_result_gen; +END_RCPP +} +// allocate_clumpy_cpp +IntegerVector allocate_clumpy_cpp(IntegerVector landscape, int nrow, int ncol, IntegerVector trans_from, IntegerVector trans_to, List prob_cell, List prob_value, NumericVector area_mean, NumericVector area_var, NumericVector elongation, NumericVector target_rate, int method, int batch_size, bool rarefy, bool shuffle, bool avoid_aggregation, int area_dist); +RcppExport SEXP _evoland_allocate_clumpy_cpp(SEXP landscapeSEXP, SEXP nrowSEXP, SEXP ncolSEXP, SEXP trans_fromSEXP, SEXP trans_toSEXP, SEXP prob_cellSEXP, SEXP prob_valueSEXP, SEXP area_meanSEXP, SEXP area_varSEXP, SEXP elongationSEXP, SEXP target_rateSEXP, SEXP methodSEXP, SEXP batch_sizeSEXP, SEXP rarefySEXP, SEXP shuffleSEXP, SEXP avoid_aggregationSEXP, SEXP area_distSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< IntegerVector >::type landscape(landscapeSEXP); + Rcpp::traits::input_parameter< int >::type nrow(nrowSEXP); + Rcpp::traits::input_parameter< int >::type ncol(ncolSEXP); + Rcpp::traits::input_parameter< IntegerVector >::type trans_from(trans_fromSEXP); + Rcpp::traits::input_parameter< IntegerVector >::type trans_to(trans_toSEXP); + Rcpp::traits::input_parameter< List >::type prob_cell(prob_cellSEXP); + Rcpp::traits::input_parameter< List >::type prob_value(prob_valueSEXP); + Rcpp::traits::input_parameter< NumericVector >::type area_mean(area_meanSEXP); + Rcpp::traits::input_parameter< NumericVector >::type area_var(area_varSEXP); + Rcpp::traits::input_parameter< NumericVector >::type elongation(elongationSEXP); + Rcpp::traits::input_parameter< NumericVector >::type target_rate(target_rateSEXP); + Rcpp::traits::input_parameter< int >::type method(methodSEXP); + Rcpp::traits::input_parameter< int >::type batch_size(batch_sizeSEXP); + Rcpp::traits::input_parameter< bool >::type rarefy(rarefySEXP); + Rcpp::traits::input_parameter< bool >::type shuffle(shuffleSEXP); + Rcpp::traits::input_parameter< bool >::type avoid_aggregation(avoid_aggregationSEXP); + Rcpp::traits::input_parameter< int >::type area_dist(area_distSEXP); + rcpp_result_gen = Rcpp::wrap(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)); + return rcpp_result_gen; +END_RCPP +} // distance_neighbors_cpp List distance_neighbors_cpp(DataFrame coords_t, double max_distance, bool quiet); RcppExport SEXP _evoland_distance_neighbors_cpp(SEXP coords_tSEXP, SEXP max_distanceSEXP, SEXP quietSEXP) { @@ -37,6 +137,12 @@ END_RCPP } static const R_CallMethodDef CallEntries[] = { + {"_evoland_raster_neighbors_cpp", (DL_FUNC) &_evoland_raster_neighbors_cpp, 2}, + {"_evoland_must_cpp", (DL_FUNC) &_evoland_must_cpp, 3}, + {"_evoland_sample_lognorm_area_cpp", (DL_FUNC) &_evoland_sample_lognorm_area_cpp, 2}, + {"_evoland_sample_normal_area_cpp", (DL_FUNC) &_evoland_sample_normal_area_cpp, 2}, + {"_evoland_grow_patch_cpp", (DL_FUNC) &_evoland_grow_patch_cpp, 14}, + {"_evoland_allocate_clumpy_cpp", (DL_FUNC) &_evoland_allocate_clumpy_cpp, 17}, {"_evoland_distance_neighbors_cpp", (DL_FUNC) &_evoland_distance_neighbors_cpp, 3}, {"_evoland_calculate_class_stats_cpp", (DL_FUNC) &_evoland_calculate_class_stats_cpp, 2}, {NULL, NULL, 0} diff --git a/src/alloc_clumpy.cpp b/src/alloc_clumpy.cpp new file mode 100644 index 0000000..9be98b1 --- /dev/null +++ b/src/alloc_clumpy.cpp @@ -0,0 +1,729 @@ +// CLUMPY allocation backend. +// +// This file implements the full CLUMPY-style allocation routine in C++, so that +// the hot pivot-selection + patch-growth loop runs without crossing into R per +// patch. The R layer (R/alloc_clumpy.R) prepares numeric inputs (the landscape +// vector, per-transition probability columns and patch parameters) and calls +// `allocate_clumpy_cpp()` once per period. +// +// Two allocation methods are provided, selected by `method`: +// +// method = 0 (uSAM, Unbiased Simple Allocation Method, thesis sec. 3.4.1): +// a single MuST pass per anterior class; every pixel that draws a +// change is allocated as a *mono-pixel* patch. This is the +// bias-free simple method and is only meaningful for mono-pixel +// patches (patch area == 1); quantity of change is enforced in +// expectation. Patch parameters (area, elongation, +// avoid_aggregation) are ignored. +// +// method = 1 (uPAM, Unbiased Patch Allocation Method, thesis sec. 3.4.2, +// Fig. 3.2): iterative. MuST is run over the remaining pool, a +// batch of pivots is drawn and grown into patches, allocated (or +// failed) pixels are removed from the pool (sampling without +// replacement) and the per-transition quota is decremented; the +// loop repeats until the quota is met or the pool is exhausted. +// This is the general multi-pixel method (and also covers +// mono-pixel patches, with a quota). +// +// The R layer auto-selects the method from the patch parameters (all mono-pixel +// => uSAM, otherwise uPAM); the explicit `method` flag here is kept so the +// backend can be exercised directly (unit tests / cross-tool comparison). +// +// Pivot rarefaction (`rarefy`): the per-cell pivot probability fed to MuST is +// P(v|u,z) / E(sigma) (thesis Fig. 3.2, factor 1/E(sigma)), because each pivot +// grows into a patch of mean area E(sigma). Without it the allocated quantity +// of change is inflated by ~mean patch size. +// +// Aggregation avoidance (`avoid_aggregation`, uPAM only): replicates clumpy's +// GaussianPatcher. Patch growth is all-or-nothing and writes nothing until the +// whole patch is built: if at any growth step (or at the final step) a border +// cell already belongs to another patch of the same transition, or the patch +// cannot reach its sampled area, the patch fails and allocates nothing. Failed +// (attempted) cells are removed from the pool for the rest of the time step +// (sampling without replacement), so patches never merge. +// +// Patch area distribution (`area_dist`): 0 = log-normal, 1 = normal. Both are +// parameterised by (area_mean, area_var); the normal uses sd = sqrt(area_var) +// and clamps to >= 1 (matching clumpy's GaussianPatcher, whose `area_cov` was +// the standard deviation). +// +// All RNG goes through R's generator (R::unif_rand / R::rlnorm / R::rnorm), so +// set.seed() in R makes a run reproducible. + +#include "clumpy_geometry.h" +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace Rcpp; + +// NOTE: cell indices are held as `int` for compatibility with R's 32-bit +// integers (terra cell numbers, IntegerVector). This caps the raster at +// INT_MAX (~2.1e9) cells; larger rasters would need a 64-bit index type here +// and on the R side. + +// --------------------------------------------------------------------------- +// Internal helpers +// --------------------------------------------------------------------------- + +// Sparse per-transition probability column: the adjusted transition potentials +// are naturally sparse (only the transition's source cells carry a value), so +// instead of a dense n x T matrix we keep, per transition, the nonzero cells and +// their values sorted by cell, with an O(log nnz) lookup. Built once from the +// R-side per-transition lists; missing cells read as 0. +struct SparseColumn { + std::vector idx; // 0-based cell indices, sorted ascending + std::vector val; + + double value_at(int cell) const { + auto it = std::lower_bound(idx.begin(), idx.end(), cell); + if (it != idx.end() && *it == cell) { + return val[(size_t)(it - idx.begin())]; + } + return 0.0; + } +}; + +// Build a SparseColumn from R vectors of 1-based cell indices and values. +// Out-of-range cells are dropped; entries are sorted by cell. +static SparseColumn build_sparse_column(const IntegerVector &cells_1based, + const NumericVector &values, int n) { + const int m = cells_1based.size(); + std::vector> tmp; + tmp.reserve(m); + for (int i = 0; i < m; ++i) { + const int c = cells_1based[i] - 1; // -> 0-based + if (c < 0 || c >= n) continue; + tmp.emplace_back(c, values[i]); + } + std::sort(tmp.begin(), tmp.end(), + [](const std::pair &a, const std::pair &b) { + return a.first < b.first; + }); + SparseColumn col; + col.idx.reserve(tmp.size()); + col.val.reserve(tmp.size()); + for (const auto &p : tmp) { + col.idx.push_back(p.first); + col.val.push_back(p.second); + } + return col; +} + +// TODO build_neighbors: Can we also look at queen neighbors, does that make sense for +// the elongation calculus? Should we expand to triangle or hexagon tesselations (Mazy +// 3.I.2)? Can 60 or 45 degree relationships be easily summarised to primary and +// secondary moments, which are assumed to be orthogonal here? Just take minimum and +// maximum of the 60deg moments? CLUMPY implements the moment accumulation in +// Patcher.allocate + +// Rook-adjacency neighbour indices (0-based, -1 == no neighbour / edge) for a +// row-major (nrow x ncol) raster. +static void build_neighbors(int nrow, int ncol, std::vector &up, + std::vector &down, std::vector &left, + std::vector &right) { + const int n = nrow * ncol; + up.assign(n, -1); + down.assign(n, -1); + left.assign(n, -1); + right.assign(n, -1); + for (int idx = 0; idx < n; ++idx) { + const int r = idx / ncol; + const int c = idx % ncol; + if (r > 0) up[idx] = idx - ncol; + if (r < nrow - 1) down[idx] = idx + ncol; + if (c > 0) left[idx] = idx - 1; + if (c < ncol - 1) right[idx] = idx + 1; + } +} + +// Log-normal patch-area draw, parameterised by the area mean E and variance V +// (cell units). Returns an integer >= 1. NA / non-positive mean -> 1. +static int draw_lognorm_area(double area_mean, double area_var) { + if (ISNAN(area_mean) || area_mean <= 0.0) return 1; + const double E = area_mean; + const double V = (ISNAN(area_var) || area_var <= 0.0) ? 1.0 : area_var; + const double mu = std::log(E * E / std::sqrt(V + E * E)); + const double sigma = std::sqrt(std::log(V / (E * E) + 1.0)); + const double draw = R::rlnorm(mu, sigma); + long a = std::lround(draw); + return a < 1 ? 1 : (int)a; +} + +// Normal patch-area draw: N(area_mean, sd = sqrt(area_var)), clamped to >= 1 +// (matches clumpy's GaussianPatcher). NA / non-positive mean -> 1. +static int draw_normal_area(double area_mean, double area_var) { + if (ISNAN(area_mean) || area_mean <= 0.0) return 1; + const double sd = (ISNAN(area_var) || area_var <= 0.0) ? 0.0 : std::sqrt(area_var); + double draw = R::rnorm(area_mean, sd); + if (draw < 1.0) draw = 1.0; + long a = std::lround(draw); + return a < 1 ? 1 : (int)a; +} + +// Dispatch on area_dist: 0 = log-normal, 1 = normal. +static int draw_area(double area_mean, double area_var, int area_dist) { + return area_dist == 1 ? draw_normal_area(area_mean, area_var) + : draw_lognorm_area(area_mean, area_var); +} + +// Shuffle two parallel vectors with the same permutation, using R's RNG (so +// set.seed() applies). +static void shuffle_pair(std::vector &a, std::vector &b) { + for (int i = (int)a.size() - 1; i > 0; --i) { + int j = (int)std::floor(R::unif_rand() * (i + 1)); + if (j > i) j = i; + std::swap(a[i], a[j]); + std::swap(b[i], b[j]); + } +} + +// Grow a single patch from `pivot0` (0-based). This is the patch-construction +// step of Mazy (2022) Appendix 3.I: starting from the pivot, neighbours are +// added one at a time until the sampled area sigma is reached, steering the +// shape towards a target elongation. +// +// Patch cells are accumulated in a local list and only committed to `land` (set +// to `to_class`) on success, so a failed patch leaves the landscape untouched +// (deferred-write rollback, like clumpy's GaussianPatcher). Greedy growth: at +// each step pick the eligible border cell maximising +// prob / (|elongation_if_added - target| + eps) +// i.e. weight by the transition potential P(v|u,z) of the candidate while +// pulling the patch's elongation (Appendix 3.I, eq. 3.I.12; computed via the +// shared clumpy::elongation_from_raw_moments) towards the target. Running +// spatial moments make each candidate's elongation an O(1) update. +// +// `avoid_aggregation` enforces the no-patch-merging requirement (Mazy sec. 3.2.4 +// "patch merging" / sec. 3.4.2), which is necessary for the post-allocation +// patch-size distribution to stay unbiased: +// * if a border cell already belongs to another patch of this transition +// (ant == from_class && land == to_class), the patch FAILS (returns 0); +// * if no eligible cell remains before reaching the target area, the patch +// FAILS (all-or-nothing); +// * a final check is applied to the last added cell's neighbours. +// With avoid_aggregation == false, the patch grows greedily up to the target +// and a short patch is accepted (no merge checks). +// +// Returns the number of cells committed (>= 1 on success; 0 on failure / invalid +// pivot). `out_cells` always receives the cells that were *attempted* (0-based), +// so the caller can remove them from the pool on failure (without replacement). +template +static int grow_one_patch(std::vector &land, const std::vector &ant, + ProbFn prob_at, const std::vector &up, + const std::vector &down, + const std::vector &left, + const std::vector &right, int pivot0, + int target_area, int from_class, int to_class, + double elong_target, int ncol, bool avoid_aggregation, + std::vector &out_cells) { + out_cells.clear(); + const int n = (int)land.size(); + if (pivot0 < 0 || pivot0 >= n) return 0; + if (land[pivot0] == NA_INTEGER) return 0; + if (land[pivot0] != from_class || ant[pivot0] != from_class) return 0; + + std::vector patch; + patch.reserve(target_area > 0 ? target_area : 1); + std::unordered_set in_patch; + patch.push_back(pivot0); + in_patch.insert(pivot0); + + double m00 = 0, s_r = 0, s_c = 0, s_rr = 0, s_cc = 0, s_rc = 0; + auto add_moments = [&](int idx) { + const double r = idx / ncol; + const double c = idx % ncol; + m00 += 1.0; + s_r += r; + s_c += c; + s_rr += r * r; + s_cc += c * c; + s_rc += r * c; + }; + add_moments(pivot0); + + auto neighbours = [&](int c, int nb[4]) { + nb[0] = up[c]; + nb[1] = down[c]; + nb[2] = left[c]; + nb[3] = right[c]; + }; + auto is_foreign = [&](int c) { + return ant[c] == from_class && land[c] == to_class; // another committed patch + }; + + while ((int)patch.size() < target_area) { + // Border of the current patch (ordered set => deterministic tie-breaks). + std::set border; + for (int pc : patch) { + int nb[4]; + neighbours(pc, nb); + for (int x : nb) { + if (x < 0 || in_patch.count(x)) continue; + border.insert(x); + } + } + + if (avoid_aggregation) { + for (int b : border) { + if (is_foreign(b)) { // would merge with another patch -> fail + out_cells.assign(patch.begin(), patch.end()); + return 0; + } + } + } + + int best = -1; + double best_score = -1.0; + for (int b : border) { + if (ant[b] != from_class || land[b] != from_class) continue; // not available + double prob = prob_at(b); + if (ISNAN(prob) || prob < 0.0) prob = 0.0; + const double r = b / ncol; + const double c = b % ncol; + const double e = clumpy::elongation_from_raw_moments( + m00 + 1.0, s_r + r, s_c + c, s_rr + r * r, s_cc + c * c, s_rc + r * c); + const double score = prob / (std::abs(e - elong_target) + 1e-6); + if (score > best_score) { + best_score = score; + best = b; + } + } + + if (best < 0) { + if (avoid_aggregation) { // could not reach target area -> all-or-nothing + out_cells.assign(patch.begin(), patch.end()); + return 0; + } + break; // partial patch accepted when not avoiding aggregation + } + + patch.push_back(best); + in_patch.insert(best); + add_moments(best); + } + + if (avoid_aggregation) { + int nb[4]; + neighbours(patch.back(), nb); + for (int x : nb) { + if (x >= 0 && is_foreign(x)) { + out_cells.assign(patch.begin(), patch.end()); + return 0; + } + } + } + + for (int c : patch) land[c] = to_class; // commit + out_cells.assign(patch.begin(), patch.end()); + return (int)patch.size(); +} + +// Inverse-CDF Multinomial Sampling Test (MuST) draw for one cell over the active +// transitions. `cum_prob(q)` yields the (non-negative) probability of the q-th +// active transition. Returns the selected active-transition index, or -1 for +// "stay". +// +// Mazy (2022) Appendix 3.B.1: with the cumulative sum eta_w = sum_{v<=w} P(v|u,z) +// and a single uniform draw xi in [0,1), the final state is the w with +// eta_{w-1} <= xi < eta_w (one draw tests all final states at once). Here the +// "stay" mass is the remainder 1 - sum_q cum_prob(q), so xi falling past the last +// cumulative bound returns -1 (no change). +// +// `must_pick` takes the uniform explicitly (so an externally supplied stream can +// be replayed for deterministic cross-tool comparison); `must_draw_one` draws it +// from R's RNG. +template static int must_pick(int k, double u, F cum_prob) { + double cs = 0.0; + for (int q = 0; q < k; ++q) { + cs += cum_prob(q); + if (u < cs) return q; // first q with u < eta_q (eta_{q-1} <= u < eta_q) + } + return -1; // stay +} +template static int must_draw_one(int k, F cum_prob) { + return must_pick(k, R::unif_rand(), cum_prob); +} + +// --------------------------------------------------------------------------- +// Exporting small building blocks for unit-testing in R/tinytest. Not actually +// used in R codebase. +// --------------------------------------------------------------------------- + +//' 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 +// [[Rcpp::export]] +List raster_neighbors_cpp(int nrow, int ncol) { + std::vector up, down, left, right; + build_neighbors(nrow, ncol, up, down, left, right); + const int n = nrow * ncol; + IntegerVector above(n), below(n), lft(n), rgt(n); + for (int i = 0; i < n; ++i) { + above[i] = up[i] < 0 ? 0 : up[i] + 1; + below[i] = down[i] < 0 ? 0 : down[i] + 1; + lft[i] = left[i] < 0 ? 0 : left[i] + 1; + rgt[i] = right[i] < 0 ? 0 : right[i] + 1; + } + return List::create(_["above"] = above, _["below"] = below, _["left"] = lft, + _["right"] = rgt); +} + +//' 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 +// [[Rcpp::export]] +IntegerVector must_cpp(NumericMatrix P, IntegerVector states, + Nullable u = R_NilValue) { + const int n = P.nrow(); + const int k = P.ncol(); + if ((int)states.size() != k) { + stop("length(states) must equal ncol(P)"); + } + const bool have_u = u.isNotNull(); + NumericVector uu; + if (have_u) { + uu = u.get(); + if ((int)uu.size() != n) stop("length(u) must equal nrow(P)"); + } + IntegerVector y(n); + auto cum = [&](int i) { + return [&, i](int q) { + double p = P(i, q); + if (ISNAN(p) || p < 0.0) p = 0.0; + return p; + }; + }; + for (int i = 0; i < n; ++i) { + const int sel = + have_u ? must_pick(k, uu[i], cum(i)) : must_draw_one(k, cum(i)); + y[i] = states[sel < 0 ? k - 1 : sel]; // sel<0 only if row sums < u + } + return y; +} + +//' 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 +// [[Rcpp::export]] +int sample_lognorm_area_cpp(double area_mean, double area_var) { + return draw_lognorm_area(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 +// [[Rcpp::export]] +int sample_normal_area_cpp(double area_mean, double area_var) { + return draw_normal_area(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 +// [[Rcpp::export]] +IntegerVector grow_patch_cpp(IntegerVector landscape, + const IntegerVector &ant_landscape, + const NumericVector &probs, + const IntegerVector &nbr_above, + const IntegerVector &nbr_below, + const IntegerVector &nbr_left, + const IntegerVector &nbr_right, int pivot, + int target_area, int from_class, int to_class, + double elongation, int ncol, + bool avoid_aggregation = false) { + const int n = landscape.size(); + std::vector land(landscape.begin(), landscape.end()); + std::vector ant(ant_landscape.begin(), ant_landscape.end()); + std::vector up(n), down(n), left(n), right(n); + for (int i = 0; i < n; ++i) { + up[i] = nbr_above[i] - 1; // 1-based (0 = none) -> 0-based (-1 = none) + down[i] = nbr_below[i] - 1; + left[i] = nbr_left[i] - 1; + right[i] = nbr_right[i] - 1; + } + std::vector pr(probs.begin(), probs.end()); + auto prob_at = [&](int c) { + return (c >= 0 && c < (int)pr.size()) ? pr[c] : 0.0; + }; + + std::vector out; + int s = grow_one_patch(land, ant, prob_at, up, down, left, right, pivot - 1, + target_area, from_class, to_class, elongation, ncol, + avoid_aggregation, out); + + if (s > 0) { + for (int idx : out) landscape[idx] = to_class; // reflect commit in caller + IntegerVector res(out.size()); + for (size_t i = 0; i < out.size(); ++i) res[i] = out[i] + 1; + return res; + } + return IntegerVector(0); // failed / invalid pivot +} + +// --------------------------------------------------------------------------- +// Exported: full allocation routine +// --------------------------------------------------------------------------- + +//' 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 +// [[Rcpp::export]] +IntegerVector allocate_clumpy_cpp( + IntegerVector landscape, int nrow, int ncol, IntegerVector trans_from, + IntegerVector trans_to, List prob_cell, List prob_value, + NumericVector area_mean, NumericVector area_var, NumericVector elongation, + NumericVector target_rate, int method, int batch_size, bool rarefy, + bool shuffle, bool avoid_aggregation, int area_dist) { + const int n = landscape.size(); + const int T = trans_from.size(); + + std::vector land(landscape.begin(), landscape.end()); + const std::vector ant(land); // immutable anterior snapshot + + std::vector up, down, left, right; + build_neighbors(nrow, ncol, up, down, left, right); + + // Build the sparse per-transition potential columns from the R lists. + if (prob_cell.size() != T || prob_value.size() != T) { + stop("prob_cell and prob_value must each have length(trans_from) elements"); + } + std::vector cols; + cols.reserve(T); + for (int t = 0; t < T; ++t) { + cols.push_back(build_sparse_column(as(prob_cell[t]), + as(prob_value[t]), n)); + } + + // Anterior classes to process, derived from trans_from. + std::set from_set(trans_from.begin(), trans_from.end()); + + // Effective (clamped, optionally rarefied) pivot-selection probability. + // Mazy Fig. 3.2: the pivot-cell probability fed to MuST is the per-pixel + // transition probability divided by the mean patch area E(sigma). Each pivot + // grows into a patch of mean size E(sigma), so to allocate the target quantity + // of change P(v|u)#J (eq. 3.11: E(#Jv^c) E(sigma) = P(v|u)#J) the pivots must + // be rarefied by 1/E(sigma); otherwise allocation overshoots by ~E(sigma). + auto pivot_prob = [&](int cell, int t) -> double { + double p = cols[t].value_at(cell); + if (ISNAN(p) || p < 0.0) p = 0.0; + if (rarefy) { + const double am = area_mean[t]; + if (!ISNAN(am) && am > 1.0) p /= am; // 1/E(sigma); mono-pixel needs none + } + return p; + }; + + for (int fc : from_set) { + std::vector at; // transition indices with trans_from == fc + for (int t = 0; t < T; ++t) { + if (trans_from[t] == fc) at.push_back(t); + } + if (at.empty()) continue; + const int k = (int)at.size(); + + // Pool of available source cells (originally fc and not yet changed). + std::vector pool; + for (int idx = 0; idx < n; ++idx) { + if (land[idx] == fc && ant[idx] == fc) pool.push_back(idx); + } + if (pool.empty()) continue; + + if (method == 0) { + // ---- uSAM: single MuST pass, mono-pixel allocation ---------------- + const int m = (int)pool.size(); + std::vector sampled(m); + for (int i = 0; i < m; ++i) { + const int cell = pool[i]; + const int sel = + must_draw_one(k, [&](int q) { return pivot_prob(cell, at[q]); }); + sampled[i] = (sel < 0) ? -1 : at[sel]; + } + for (int i = 0; i < m; ++i) { + if (sampled[i] < 0) continue; + const int pv = pool[i]; + if (land[pv] == fc && ant[pv] == fc) land[pv] = trans_to[sampled[i]]; + } + } else { + // ---- uPAM: iterative MuST + quota + without-replacement ------------ + // Mazy sec. 3.4.2 (Fig. 3.2): per-transition quantity-of-change quota in + // pixels, N_{u->v} = P(v|u) #J (App. 3.E.2 Algorithm 3, line 3; eq. 3.11); + // decremented by each successful patch's size sigma until exhausted. The + // pool is sampled without replacement (allocated and failed-attempt cells + // are removed), which is what keeps the post-allocation distributions + // unbiased. + std::vector remaining(k); + const double m0 = (double)pool.size(); + for (int q = 0; q < k; ++q) { + double rt = target_rate[at[q]]; + if (ISNAN(rt) || rt < 0.0) rt = 0.0; + remaining[q] = rt * m0; // N_{u->v} = P(v|u) * #J + } + std::vector blocked(n, 0); // attempted-but-failed cells (removed) + // batch_size controls how many pivots are grown per MuST re-draw: + // > 0 explicit cap; + // < 0 all candidates in a single pass; + // 0 AUTO: ~1% of this class's pool, so the number of MuST passes is + // bounded (~100) and the per-pass pool re-scan does not blow up to + // the O(#patches x pool) cost of strict batch=1 on large rasters. + int bs; + if (batch_size > 0) { + bs = batch_size; + } else if (batch_size < 0) { + bs = INT_MAX; + } else { + bs = std::max(1, (int)(pool.size() / 100)); + } + + int guard = 0; + const int guard_max = n + 16; + while (guard++ < guard_max) { + bool any_quota = false; + for (int q = 0; q < k; ++q) { + if (remaining[q] > 0.0) { + any_quota = true; + break; + } + } + if (!any_quota) break; + + // Compact the pool to cells still available and not blocked. + pool.erase(std::remove_if(pool.begin(), pool.end(), + [&](int idx) { + return blocked[idx] || + !(land[idx] == fc && ant[idx] == fc); + }), + pool.end()); + if (pool.empty()) break; + + // MuST over the current pool, restricted to in-quota transitions. + std::vector cand_cell, cand_q; + for (int idx : pool) { + const int sel = must_draw_one(k, [&](int q) { + return remaining[q] > 0.0 ? pivot_prob(idx, at[q]) : 0.0; + }); + if (sel >= 0) { + cand_cell.push_back(idx); + cand_q.push_back(sel); + } + } + if (cand_cell.empty()) break; + if (shuffle) shuffle_pair(cand_cell, cand_q); + + int processed = 0; + std::vector out; + for (size_t c = 0; c < cand_cell.size() && processed < bs; ++c) { + const int q = cand_q[c]; + if (remaining[q] <= 0.0) continue; + const int pv = cand_cell[c]; + if (blocked[pv] || land[pv] != fc || ant[pv] != fc) continue; + const int t = at[q]; + const int to = trans_to[t]; + const int area = draw_area(area_mean[t], area_var[t], area_dist); + const SparseColumn &col = cols[t]; + const int s = grow_one_patch( + land, ant, [&](int c) { return col.value_at(c); }, up, down, left, + right, pv, area, fc, to, elongation[t], ncol, avoid_aggregation, + out); + ++processed; // every attempt counts toward the batch + if (s > 0) { + remaining[q] -= (double)s; + } else { + for (int b : out) blocked[b] = 1; // remove attempted cells + } + } + // Each iteration removes >= 1 cell from the pool (a success commits + // to_class, a failure blocks the attempted cells), so the loop is + // guaranteed to terminate. + } + } + } + + IntegerVector res(n); + for (int i = 0; i < n; ++i) res[i] = land[i]; + return res; +} diff --git a/src/clumpy_geometry.h b/src/clumpy_geometry.h new file mode 100644 index 0000000..7cde20c --- /dev/null +++ b/src/clumpy_geometry.h @@ -0,0 +1,77 @@ +#ifndef EVOLAND_CLUMPY_GEOMETRY_H +#define EVOLAND_CLUMPY_GEOMETRY_H + +#include + +// Shared patch-shape geometry for the CLUMPY allocation backend. +// +// `elongation_from_raw_moments()` is the SINGLE definition of the patch shape +// metric used both during calibration (`calculate_class_stats_cpp` in +// patch_stats.cpp, which measures the elongation of observed patches) and +// during allocation (`grow_patch_cpp` / `allocate_clumpy_cpp` in +// alloc_clumpy.cpp, which grows patches towards a target elongation). Keeping +// one definition guarantees the shape that allocation optimises is exactly the +// shape calibration reports -- previously these lived in two functions +// (`calculate_elongation` and `patch_eccentricity`) that computed the same +// quantity with different variable names and indexing conventions. +// +// Definition (Mazy 2022, https://theses.hal.science/tel-04382012, eq. 3.I.12): +// +// e = 1 - sqrt(lambda_min / lambda_max) +// +// where lambda_min <= lambda_max are the eigenvalues of the 2x2 matrix of +// second central moments of the cell coordinates. e = 0 is circular / +// isometric; e -> 1 is increasingly elongated (linear). +// +// The eigenvalues of the (symmetric) second-moment matrix are invariant under +// the choice of axes, so callers may pass row/col or x/y coordinates, and it +// does not matter whether the raster is stored row- or column-major. +// +// Arguments are raw (un-centred) moment accumulators over the cell set: +// m00 = number of cells +// s10 = sum of first-axis coordinates (e.g. sum of row indices) +// s01 = sum of second-axis coordinates (e.g. sum of col indices) +// s20 = sum of first-axis coordinates^2 +// s02 = sum of second-axis coordinates^2 +// s11 = sum of (first-axis * second-axis) +namespace clumpy { + +inline double elongation_from_raw_moments(double m00, double s10, double s01, + double s20, double s02, double s11) { + if (m00 <= 1.0) { + return 0.0; // single cell or empty: isometric by convention + } + + // Second central moments about the centre of mass: + // mu_pq = (sum x^p y^q)/n - mean_x^p mean_y^q + const double mean_1 = s10 / m00; + const double mean_2 = s01 / m00; + + const double mu20 = s20 / m00 - mean_1 * mean_1; + const double mu02 = s02 / m00 - mean_2 * mean_2; + const double mu11 = s11 / m00 - mean_1 * mean_2; + + const double trace = mu20 + mu02; + if (trace <= 0.0) { + return 0.0; + } + + const double delta = + std::sqrt((mu20 - mu02) * (mu20 - mu02) + 4.0 * mu11 * mu11); + + const double lambda_max = 0.5 * (trace + delta); + double lambda_min = 0.5 * (trace - delta); + + if (lambda_max <= 0.0) { + return 0.0; + } + if (lambda_min < 0.0) { + lambda_min = 0.0; // guard against floating-point round-off + } + + return 1.0 - std::sqrt(lambda_min / lambda_max); +} + +} // namespace clumpy + +#endif // EVOLAND_CLUMPY_GEOMETRY_H diff --git a/src/patch_stats.cpp b/src/patch_stats.cpp index 688a125..eb23a10 100644 --- a/src/patch_stats.cpp +++ b/src/patch_stats.cpp @@ -1,6 +1,6 @@ +#include "clumpy_geometry.h" #include #include -#include #include #include @@ -95,40 +95,10 @@ struct UnionFind { */ double calculate_elongation(double m00, double m10, double m01, double s20, double s02, double s11) { - if (m00 <= 1.0) - return 0.0; // Single pixel or empty - - // Normalized second order moments about the center of mass - double mean_x = m10 / m00; - double mean_y = m01 / m00; - - // eqs 3.I.6, 3.I.7, 3.I.8 - // mu20 = sum((x - mean_x)^2) / n becomes - // mu20 = (sum(x^2) / n) - mean_x^2 because - // sum(x) = n * mean_x - double mu20 = (s20 / m00) - (mean_x * mean_x); - double mu02 = (s02 / m00) - (mean_y * mean_y); - double mu11 = (s11 / m00) - (mean_x * mean_y); - - // Eigenvalues of the inertia tensor - // D = sqrt((mu20 - mu02)^2 + 4 * mu11^2) - double term1 = mu20 - mu02; - double D = std::sqrt(term1 * term1 + 4.0 * mu11 * mu11); - - double lambda1 = 0.5 * (mu20 + mu02 - D); - double lambda2 = 0.5 * (mu20 + mu02 + D); - - // lambda2 is the larger eigenvalue (major axis related variance) - // lambda1 is the smaller eigenvalue (minor axis related variance) - - if (lambda2 <= 1e-9) - return 0.0; // Should not happen for >1 pixel, but check for safety - - // Avoid negative values inside sqrt due to floating point precision - if (lambda1 < 0) - lambda1 = 0; - - return 1.0 - std::sqrt(lambda1 / lambda2); + // Single shared definition lives in clumpy_geometry.h so that the shape + // metric measured here during calibration is identical to the one the patch + // grower (alloc_clumpy.cpp) targets during allocation. + return clumpy::elongation_from_raw_moments(m00, m10, m01, s20, s02, s11); } /** @@ -242,6 +212,7 @@ DataFrame calculate_class_stats_cpp(IntegerMatrix mat, double cellsize) { // 3. Compute Final Statistics per Class std::vector out_class; + std::vector out_count; std::vector out_mean_area; std::vector out_variance_area; std::vector out_mean_elongation; @@ -275,10 +246,12 @@ DataFrame calculate_class_stats_cpp(IntegerMatrix mat, double cellsize) { } double mean_elong = sum_elong / n; out_mean_elongation.push_back(mean_elong); + out_count.push_back((int)areas.size()); } return DataFrame::create( - Named("class") = out_class, Named("patch_area_mean") = out_mean_area, + Named("class") = out_class, Named("patch_count") = out_count, + Named("patch_area_mean") = out_mean_area, Named("patch_area_variance") = out_variance_area, Named("patch_elongation_mean") = out_mean_elongation); } diff --git a/vignettes/evoland.qmd b/vignettes/evoland.qmd index 5ce924f..8787c2e 100644 --- a/vignettes/evoland.qmd +++ b/vignettes/evoland.qmd @@ -12,7 +12,8 @@ number-sections: true ```{r} #| label: cleanup #| include: false -# TODO render tables into nice HTML interactive ones?? +# TODO render tables into nice HTML interactive ones +# TODO use the nicer synthetic data routine from stochastic tutorial # make sure we're starting clean unlink("firstmodel.evolanddb", recursive = TRUE) @@ -263,16 +264,20 @@ instead, we generate a synthetic LULC raster with 3 layers (one per period) by a # autoregressive noise with skellam distribution n_cells <- dim(template_rast)[1] * dim(template_rast)[2] noise1 <- runif(n_cells, min = 0, max = 10) -noise2 <- noise1 + stats::rpois(n_cells, 1) - stats::rpois(n_cells, 1) -noise3 <- noise2 + stats::rpois(n_cells, 1) - stats::rpois(n_cells, 1) +noise2 <- noise1 + stats::rpois(n_cells, 0.2) - stats::rpois(n_cells, 0.2) +noise3 <- noise2 + stats::rpois(n_cells, 0.2) - stats::rpois(n_cells, 0.2) synthetic_lulc <- rast(template_rast, nlyrs = 3, vals = c(noise1, noise2, noise3)) |> focal(w = 3, fun = mean, na.rm = TRUE) |> clamp(lower = 0, upper = 10) |> - classify(rcl = data.frame( - from = 0:9, to = 1:10, becomes = c(3, 7, 1, 10, 5, 8, 2, 9, 4, 6) - )) + classify( + rcl = data.frame( + from = 0:9, + to = 1:10, + becomes = c(3, 7, 1, 10, 5, 8, 2, 9, 4, 6) + ) + ) plot(synthetic_lulc, nc = 3) ``` @@ -435,7 +440,7 @@ db$trans_rates_t <- ) ``` -We estimate allocation parameters for DinamicaEGO, which determine the shape and size of new patches, respectively which fraction of converted land use is in new versus expanded patches. +We estimate allocation parameters, which determine the shape and size of new patches, respectively which fraction of converted land use is in new versus expanded patches. This estimation procedure is not unbiased and hence a single estimate may not be enough: normally, we would perturb the estimate and use multiple `id_run`s to identify the best parametrization. For simplicity, we now just take the estimates for granted and assign `id_run=0`, i.e. the base run ID. @@ -449,17 +454,16 @@ db$alloc_params_t <- alloc_for_eval # Prediction + Allocation -With all components in place, we run the allocation step. -If Dinamica is not installed, you'll get a warning that the anterior LULC map is returned as the posterior. -For installation instructions, see the [Installing Dinamica EGO](install-dinamica.html) vignette. +For this tutorial, we will use the CLUMPY backend for a self-contained stochastic allocation that +does not require the presence of [DinamicaEGO](install-dinamica.html) as an external solver. ```{r} #| label: allocation -#| results: hold -db$alloc_dinamica( - id_period = db$periods_t[is_extrapolated == TRUE, id_period], +db$alloc_clumpy( + id_period = db$periods_t[is_extrapolated == TRUE, id_period], # select all extrapolation periods select_score = "classif.auc", - select_maximize = TRUE + select_maximize = TRUE, + seed = 42L # optional: reproducibility ) ``` @@ -468,15 +472,38 @@ db$alloc_dinamica( Finally, we can extract the simulated LULC maps into `SpatRaster` objects to visualize them. ```{r} -#| label: visualization +#| label: visualization-cats labels <- db$periods_t[ id_period != 0, paste0(year(start_date), " to ", year(end_date)) ] -plot_maps <- db$lulc_data_as_rast() |> setNames(labels) - -plot(plot_maps, type = "classes", levels = db$lulc_meta_t$pretty_name) +plot_maps <- + db$lulc_data_as_rast() |> + categories( + layer = 0, # set for all layers + value = data.frame(id = 1:4, name = db$lulc_meta_t$pretty_name) + ) |> + setNames(labels) + +plot(plot_maps) ``` -Of the four maps, only the first three show changes. -Due to Dinamica not being available on our github runner, the extrapolated step should look exactly like the step before. +Due to the extrapolated transition rates, the last step shows a similar difference as exists between the first three periods. +We can also show the cumulative difference like so: + +```{r} +#| label: visualization-changes +changes <- + c( + create_change_map(plot_maps[[1]], plot_maps[[2]]), + create_change_map(plot_maps[[1]], plot_maps[[3]]), + create_change_map(plot_maps[[1]], plot_maps[[4]]) + ) |> + categories( + layer = 0, + value = data.frame(id = 1:4, name = paste("Posterior:", db$lulc_meta_t$pretty_name)) + ) |> + setNames(paste("Period 1 to", 2:4)) + +plot(changes) +``` diff --git a/vignettes/stochastic-allocation-sensitivity.qmd b/vignettes/stochastic-allocation-sensitivity.qmd new file mode 100644 index 0000000..e94859a --- /dev/null +++ b/vignettes/stochastic-allocation-sensitivity.qmd @@ -0,0 +1,629 @@ +--- +title: "Stochastic allocation sensitivity with `runs_t`" +author: Jan Hartman +date: last-modified +vignette: > + %\VignetteIndexEntry{stochastic-allocation-sensitivity} + %\VignetteEngine{quarto::html} + %\VignetteEncoding{UTF-8} +number-sections: true +--- + +```{r} +#| label: cleanup +#| include: false +unlink("stochastic-sensitivity.evolanddb", recursive = TRUE) + +set.seed(666) +knitr::knit_hooks$set(seed = function(before, options, envir) { + if (before && !is.null(options$seed)) { + set.seed(options$seed) + } +}) + +local({ + chunk_messages <- character() + + knitr::knit_hooks$set(message = function(x, options) { + chunk_messages <<- c(chunk_messages, x) + return(character(0)) + }) + + default_chunk_hook <- knitr::knit_hooks$get("chunk") + knitr::knit_hooks$set(chunk = function(x, options) { + out <- default_chunk_hook(x, options) + if (length(chunk_messages) > 0) { + msg_text <- paste(chunk_messages, collapse = "") + callout <- paste0( + "\n\n
Messages
\n",
+        msg_text,
+        "
\n" + ) + chunk_messages <<- character() + out <- paste0(out, callout) + } + return(out) + }) +}) +``` + +*Note: this tutorial assumes that you already know the basic `evoland-plus` workflow and the concepts of transition potentials, constrained demand, and patch-based allocation. The preceding tutorial introduced those ideas; here we focus on some advantages of evoland-style sensitivity analyses.* + +## Why another tutorial? + +The default tutorial walks through a single calibration-and-allocation workflow, ending with one extrapolated realization. That is a good first contact, but it hides an important property of patch-based allocators: **once transition demand and transition potentials are fixed, the final pattern can still vary from run to run because allocation is stochastic**. + +For sensitivity analysis, a single map is therefore not enough. In this vignette we will: + +1. establish a simple synthetic model, +2. register 30 allocation realizations in `runs_t`, +3. allocate exactly one extrapolated period 30 times, +4. compare each extrapolated realization to the last observed map, and +5. summarize the results as a **change-frequency heatmap** showing where change happens consistently. + +This is useful whenever you want to distinguish **structurally robust change**, i.e. cells that nearly always switch class. + + +## Setup + +```{r} +#| label: setup +#| output: false +library(evoland) +library(data.table) +library(terra) +``` + +We create a fresh on-disk database with some synthetic land use categories; +we use a spatial domain of 30x30 cells on a square grid defined in meters and retain a template raster for later use; +and we set up three "observed" periods to be filled with synthetic data, plus a single extrapolated period. + +```{r} +#| label: create-db +#| collapse: true +db <- evoland_db$new(path = "stochastic-sensitivity.evolanddb") + +db$lulc_meta_t <- create_lulc_meta_t( + list( + forest = list(pretty_name = "Forest"), + arable = list(pretty_name = "Arable Land"), + urban = list(pretty_name = "Urban Areas"), + static = list(pretty_name = "Immutable") + ) +) + +template_rast <- terra::rast( + crs = "EPSG:2056", + extent = terra::ext(c( + xmin = 2697000, + xmax = 2697000 + 30 * 100, + ymin = 1252000, + ymax = 1252000 + 30 * 100 + )), + resolution = 100 +) + +db$coords_t <- create_coords_t_square( + epsg = terra::crs(template_rast, describe = TRUE)$code |> as.integer(), + extent = terra::ext(template_rast), + resolution = terra::res(template_rast)[1] +) + +db$periods_t <- create_periods_t( + period_length_str = "P10Y", + start_observed = "1995-01-01", + end_observed = "2020-01-01", + end_extrapolated = "2030-01-01" +) +``` + +## Synthetic observed data + +As in the introductory vignette, we create a small autoregressive synthetic landscape with three observed raster layers. +We then transform the synthetic raster stack into `lulc_data_t` records and attach `id_run = 0L`, which is the reserved base run for observed data. + +```{r} +#| seed: 123 +#| label: synthesize-lulc +#| code-fold: true +#| fig-asp: 0.3 +n_cells <- dim(template_rast)[1] * dim(template_rast)[2] +noise1 <- runif(n_cells, min = -0.5, max = 4.2) +noise2 <- noise1 + stats::rpois(n_cells, 0.1) - stats::rpois(n_cells, 0.1) +noise3 <- noise2 + stats::rpois(n_cells, 0.1) - stats::rpois(n_cells, 0.1) + +synthetic_lulc <- + rast(template_rast, nlyrs = 3, vals = c(noise1, noise2, noise3)) |> + focal(w = 3, fun = mean, na.rm = TRUE) |> + clamp(lower = 0, upper = 4) |> + classify( + rcl = data.frame(from = 0:3, to = 1:4, becomes = c(2, 1, 3, 4)), + right = FALSE + ) |> + setNames(paste0("synthetic LULC id_period=", 1:3)) + +plot( + synthetic_lulc, + nc = 3, + col = data.frame( + value = 1:4, + color = c("#91B690", "#EB9486", "#F3DE8A", "#CBC6D2") + ) +) + +db$lulc_data_t <- + extract_using_coords_t(synthetic_lulc, db$coords_t)[, + .( + id_run = 0L, + id_coord, + id_period = substr(layer, 26, 26) |> as.integer(), + id_lulc = value + ) + ] |> + as_lulc_data_t() +``` + +## Predictors and neighborhoods + +Here, we confabulate some predictors based on gradient based "latent" influences then used to construct an "accessibility" and "site quality" quantifier. + +```{r} +#| label: predictors +#| code-fold: true +#| seed: 321 +#| fig-asp: 0.3 +# helper: scale raster to [0, 1] +scale01 <- function(x) { + rng <- terra::global(x, c("min", "max"), na.rm = TRUE) + (x - rng[1, 1]) / (rng[1, 2] - rng[1, 1]) +} + +# helper: smooth random field +smooth_field <- function(template, sd = 1, w = 7) { + terra::setValues(template, rnorm(terra::ncell(template), sd = sd)) |> + terra::focal(w = w, fun = mean, na.rm = TRUE) |> + scale01() +} + +# coordinate-based gradients +xy <- terra::crds(template_rast, df = TRUE) +x_grad <- terra::setValues(template_rast, (xy$x - min(xy$x)) / (max(xy$x) - min(xy$x))) +y_grad <- terra::setValues(template_rast, (xy$y - min(xy$y)) / (max(xy$y) - min(xy$y))) + +latent1 <- smooth_field(template_rast, sd = 1, w = 9) +latent2 <- smooth_field(template_rast, sd = 1, w = 5) + +accessibility <- scale01(0.55 * (1 - x_grad) + 0.25 * (1 - y_grad) + 0.20 * latent1) +site_quality <- scale01(0.50 * y_grad + 0.35 * latent2 + 0.15 * x_grad) + +c(accessibility, site_quality) |> + setNames(c("accessibility", "site_quality")) |> + plot() +``` + +These "drivers" are used to construct predictor variables that are somewhat correlated with the synthetic land use map we constructed earlier - we combine the drivers, noise, and land class occurrence density and normalize to $[0,1]$. + +```{r} +#| label: predictors2 +#| code-fold: true +#| seed: 321 +#| fig-asp: 1 +as_indicator <- function(x, class) { + terra::app(x, fun = function(v) as.numeric(v %in% class)) +} + +forest_share <- + as_indicator(synthetic_lulc, class = 1) |> + terra::focal(w = 5, fun = mean, na.rm = TRUE) |> + terra::mean() +arable_share <- + as_indicator(synthetic_lulc, class = 2) |> + terra::focal(w = 5, fun = mean, na.rm = TRUE) |> + terra::mean() +urban_share <- + as_indicator(synthetic_lulc, class = 3) |> + terra::focal(w = 5, fun = mean, na.rm = TRUE) + +make_dynamic <- function(base_context, static1, static2, noise_weight = 0.10) { + out <- vector("list", terra::nlyr(base_context)) + for (i in seq_len(terra::nlyr(base_context))) { + noise_i <- smooth_field(template_rast, sd = 1, w = 5) + # fmt: skip + out[[i]] <- scale01( + 0.55 * base_context[[i]] + + 0.25 * static1 + + 0.10 * static2 + noise_weight * noise_i + ) + } + terra::rast(out) +} + +urban_pressure <- + make_dynamic( + base_context = urban_share, + static1 = accessibility, + static2 = 1 - site_quality + ) +add(urban_pressure) <- + # fmt: skip + scale01( + 0.60 * urban_pressure[[3]] + + 0.25 * urban_share[[3]] + + 0.10 * accessibility + + 0.05 * smooth_field(accessibility) + ) +names(urban_pressure) <- + paste0("urban_pressure id_period=", 1:4) + +forest_suitability <- + make_dynamic( + base_context = forest_share, + static1 = site_quality, + static2 = 1 - accessibility + ) |> + setNames("forest_suitability id_period=0") + +arable_yield <- + make_dynamic( + base_context = arable_share, + static1 = 1 - site_quality, + static2 = accessibility + ) |> + setNames("arable_yield id_period=0") + +# optional nuisance predictor that should usually score poorly +random_nuisance <- + smooth_field(template_rast, sd = 1, w = 3) |> + setNames("random_nuisance id_period=0") + +plot(c(urban_pressure, forest_suitability, arable_yield, random_nuisance)) +``` + +Now we can declare minimal predictor metadata and transform it into tabular form. + +```{r} +#| label: predictors3 +#| code-fold: true +#| seed: 321 +#| fig-asp: 0.3 +db$pred_meta_t <- + create_pred_meta_t(list( + urban_pressure = list( + description = "synthetic confabulation", + data_type = "float" + ), + forest_suitability = list( + description = "synthetic confabulation", + data_type = "float" + ), + arable_yield = list( + description = "synthetic confabulation", + data_type = "float" + ), + random_nuisance = list( + description = "synthetic confabulation", + data_type = "float" + ) + )) + +db$pred_data_t <- + extract_using_coords_t( + c(urban_pressure, forest_suitability, arable_yield, random_nuisance), + db$coords_minimal + )[, layer := as.character(layer)][, + .( + id_coord, + id_run = 0L, + id_period = substr(layer, nchar(layer), nchar(layer)) |> as.integer(), + id_pred = fcase( + grepl("urban_pressure", layer) , 1L , + grepl("forest_suitability", layer) , 2L , + grepl("arable_yield", layer) , 3L , + grepl("random_nuisance", layer) , 4L + ), + value + ) + ] |> + as_pred_data_t() +``` + +Now we add a simple set of neighbor relations as predictors. + +```{r} +#| label: predictors4 +#| code-fold: true +db$set_neighbors( + max_distance = 1000, + distance_breaks = c(0, 300, 1000), + quiet = TRUE +) +db$generate_neighbor_predictors() +``` + +## Calibration + +### Eligible transitions and predictor pruning + +Here we subset the number of eligible transitions according to number of observed transitions. +Then we construct a full set of transition/predictor permutations and use it to estimate each predictor's utility for each transition. +We then retain the best 3 predictors per transition. +These will then be used to train models. + +```{r} +#| label: transitions-and-predictors +#| seed: 666 +db$trans_meta_t <- create_trans_meta_t( + db$trans_v, + min_cardinality_abs = 20, + exclude_anterior = 4 # exclude immutable class +) + +db$set_full_trans_preds() + +trans_pred_scored <- db$get_pred_filter_score( + filter = mlr3filters::FilterImportance$new( + learner = mlr3::lrn("classif.rpart") + ) +) + +db$commit( + trans_pred_scored[order(-importance)][, head(.SD, 3), by = id_trans], # top three per id_trans + "trans_preds_t", + method = "overwrite" +) +``` + +### Transition models + +Here, we train one single model per transition; see the main tutorial for an illustration on how to train multiple partials on split samples to assess their goodness of fit. + +```{r} +#| label: fit-models +#| seed: 666 +trans_models <- db$fit_full_models(learner = mlr3::lrn("classif.rpart")) + +# which transitions produced a usable model (non-null learner_full)? read this +# from the returned object before committing it +modeled_trans <- unique( + trans_models$id_trans[!vapply(trans_models$learner_full, is.null, logical(1L))] +) + +db$trans_models_t <- trans_models +``` + +One subtlety is worth calling out, because it is easy to trip over in a stochastic +setting. Model fitting does not necessarily produce a usable model for *every* viable +transition: a rare transition can fail to train -- for example a `classif.rpart` tree +that never splits, or too few positive cases in the random draw of the synthetic data -- +in which case its `learner_full` is left empty (`NULL`). Such a transition would stay +`is_viable == TRUE` yet have no model, and allocation (which predicts a potential for +every viable transition) would then abort with `No fitted model for viable +transition(s): ...`. We keep the viable set and the fitted models in lockstep by +demoting any transition whose full model did not train, reading the result returned by +`fit_full_models()` directly: + +```{r} +#| label: reconcile-viable-transitions +trans_meta_reconciled <- db$trans_meta_t +trans_meta_reconciled[, is_viable := is_viable & id_trans %in% modeled_trans] +db$trans_meta_t <- trans_meta_reconciled +``` + +### Transition rates and baseline allocation parameters + +We prescribe a set of transition rates (i.e. share of cells to transition from one land use to another) that are to be used in period 4. + +```{r} +#| label: transition-rates +#| seed: 666 +db$trans_rates_t <- + # fmt: skip + rowwiseDT( + id_trans = , rate = , + 1, 0.3, + 2, 0.2, + 3, 0.1, + 4, 0, + 5, 0, + 6, 0.2 + )[, `:=`( + id_run = 0L, + id_period = 4L, + count = n_cells * rate + )] |> + as_trans_rates_t() + +db$alloc_params_t <- + db$create_alloc_params_t(n_perturbations = 0)[, id_run := 0L] +``` + +## Registering stochastic realizations in `runs_t` + +The central idea of this tutorial is that **each stochastic allocation is a separate run**. We therefore create 30 run records and propagate the same baseline allocation parameters to each run. + +As with all evoland tables, `runs_t` only ensures a minimum set of columns, so we can add the `kind` and `seed` columns when overwriting. + +```{r} +#| label: register-runs +runs_new <- + data.table( + id_run = 0:30, # base run 0 & stochastic runs 1:30 + parent_id_run = NA_integer_, # NA for the scenario root + description = "base", + kind = "synthetic 'historical' case", + seed = 666L + )[ + id_run > 0, + `:=`( + parent_id_run = 0L, + description = paste("stochastic allocation", id_run), + kind = "allocation sensitivity", + seed = id_run + 100L + ) + ] + +db$commit(as_runs_t(runs_new), "runs_t", method = "overwrite") + +runs_out <- db$runs_t[id_run > 0] +run_ids <- runs_out$id_run +``` + +Note that we do **not** copy the allocation parameters into each run. Every stochastic run is registered with `parent_id_run = 0`, so the run lineage lets each one inherit the baseline `alloc_params_t` (and the rest of the calibration) from run `0` automatically. This guarantees that **all differences among runs arise from the allocator's stochasticity -- driven by the per-run `seed` -- and not from differing parameters**. + +## Running one extrapolated period 30 times + +We only allocate the first extrapolated period. This makes the comparison easier to interpret: every realization starts from the same observed state. +We loop over the run records: for each run we set the active `db$id_run` -- which controls both the lineage that allocation reads from and the `id_run` stamped on the results it writes -- and then allocate the single extrapolated period (`id_period = 4`) using that run's `seed`. Because only the seed changes between iterations, the resulting maps differ purely through allocation stochasticity. + +```{r} +#| label: allocate-realizations +#| output: false +for (i in seq_len(nrow(runs_out))) { + db$id_run <- runs_out$id_run[i] + db$alloc_clumpy( + id_periods = 4, # first extrapolated period + select_score = "classif.auc", + select_maximize = TRUE, + seed = runs_out$seed[i] + ) +} +``` + +## Summarizing stochastic sensitivity + +A convenient first summary is binary: for each cell we ask whether its simulated class differs from the last observed class. Averaging that binary indicator over 30 realizations gives a **change frequency** between 0 and 1. + +- `0` means the cell never changes. +- `1` means the cell always changes. +- intermediate values indicate uncertainty induced by the allocator. + +```{r} +#| label: summarize-change-frequency +db$id_run <- NULL # setting to NULL to access all records + +# 900 locations x 1 observation = 900 rows +observed_last <- db$fetch("lulc_data_t", where = "id_run = 0 and id_period = 3")[, + .(id_coord, id_lulc_observed = id_lulc) +] + +# 900 locations * 30 runs = 27000 rows +simulated_extrap <- db$fetch("lulc_data_t", where = "id_run > 0 and id_period = 4")[, + .(id_run, id_coord, id_lulc_simulated = id_lulc) +] + +change_maps <- + simulated_extrap[ + observed_last, + on = .(id_coord) + ][, + .( + n_changed = sum(id_lulc_simulated != id_lulc_observed) |> as.double(), + p_changed = mean(id_lulc_simulated != id_lulc_observed), + observed = id_lulc_observed[1] |> as.double(), + mode_simulated = id_lulc_simulated |> tabulate() |> which.max() |> as.double() + ), + by = .(id_coord) + ] |> + data.table::melt(id.vars = "id_coord") |> + tabular_to_raster(coords = db$coords_minimal, value_col = "value") +``` + +The chunk above melts these four per-cell summaries into long form and rasterizes them with `tabular_to_raster()`, producing one map layer per statistic. Below we plot the change-frequency layer as a heatmap. + +## Heatmap of stochastic change consistency + +```{r} +#| label: plot-change-frequency +#| fig-width: 4 +#| fig-height: 4 +plot( + change_maps$variable_p_variable_changed, + main = "Share of runs in which a cell changes class" +) +``` + +The interpretation is straightforward: + +- **cooler / lighter cells** are stable under repeated allocation, +- **warmer / darker cells** change consistently and may reflect strong transition potential and/or limited competition, and +- **mid-range cells** reveal locations where several cells compete for a limited transition demand and small stochastic differences decide which ones actually convert. + +In other words, this map is not a transition-potential surface. It is a **realization-frequency surface**. + +## Inspecting individual realizations + +The aggregate heatmap compresses 30 runs into a single surface. It is also worth looking at the individual realizations to get an intuition for how much the allocator actually moves things around from one run to the next. +Rather than printing 30 static panels, we stitch the realizations into a short looping animation -- one frame per run. + +```{r} +#| label: animate-realisations +#| animation.hook: gifski +#| interval: 0.4 +#| fig-width: 3 +#| fig-height: 3 +#| fig-cap: "Each frame is one stochastic allocation of the single extrapolated period." +realisation_rasters <- + simulated_extrap |> + tabular_to_raster(coords = db$coords_minimal, value_col = "id_lulc_simulated") + +lulc_palette <- data.frame( + value = 1:4, + color = c("#91B690", "#EB9486", "#F3DE8A", "#CBC6D2") +) + +# one plot per layer -> one GIF frame per realization +for (i in seq_len(terra::nlyr(realisation_rasters))) { + plot( + realisation_rasters[[i]], + legend = FALSE, + axes = FALSE, + mar = c(0.1, 0.1, 0.1, 0.1), + main = paste0("Realisation ", i, " / ", terra::nlyr(realisation_rasters)), + col = lulc_palette + ) +} +``` + +Even with identical calibration, transition rates, and allocation parameters, local differences appear across runs. The heatmap above compresses those differences into a robust summary. + +## Optional: class-specific consistency maps + +The previous heatmap asked only whether a cell changes at all. You can build more targeted summaries in the same way. For example, the next chunk computes the frequency with which a cell becomes `urban`. + +```{r} +#| label: urban-frequency +urban_frequency <- + simulated_extrap[, + .( + p_urban = mean(id_lulc_simulated == 3L) #urban id + ), + by = .(id_coord) + ] |> + tabular_to_raster(db$coords_minimal, value_col = "p_urban") + +plot( + urban_frequency, + main = "Share of runs in which a cell ends up urban", +) +``` + +This is often the more policy-relevant view if stakeholders care about a specific target class rather than generic change. + +## What `runs_t` buys you + +In a small tutorial like this, we could have stored realizations in ad hoc objects. Using `runs_t` is better because it scales to more serious experiments: + +- each realization is explicitly registered, +- run-level metadata such as seeds, scenario labels, or parameter sets can be stored alongside the simulation, +- downstream summaries can be grouped by run families, and +- the same database structure can later support sensitivity analysis, parameter sweeps, or scenario ensembles. + +A natural next step would be to combine **between-run stochasticity** (shown here) with **between-parameter variability** by creating multiple allocation parameter sets per scenario and multiple stochastic replicates per parameter set. + +Repeated allocation changes the interpretation of the model output: + +- one extrapolated map is a **draw**, +- the stack of 30 realizations is an **ensemble**, and +- the heatmap of change frequency is a compact way to communicate **where the ensemble is robust and where it is uncertain**. + +That makes `runs_t` more than bookkeeping: it becomes the backbone for uncertainty-aware land use change analysis.