#' @title Plot Bipartite Network Metrics by Age

#' @description

#' Creates a faceted line-and-point plot of bipartite community

#' network metrics (e.g. connectance, nestedness, modularity)

#' across time slices (ages). Each metric is shown in its own

#' panel with a free y-axis scale.

#' @param data_network_metrics

#' A data frame or tibble with columns \code{age} (numeric),

#' \code{metric} (character), and \code{value} (numeric).

#' Typically the output of

#' \code{\link{compute_network_metrics}()} aggregated across

#' time slices and with the age column already parsed to numeric.

#' @param title

#' Optional character string for the plot title.

#' Defaults to \code{NULL} (no title).

#' @param subtitle

#' Optional character string for the plot subtitle.

#' Defaults to \code{NULL} (no subtitle).

#' @return

#' A \code{ggplot} object.

#' @details

#' The x-axis is reversed so that older ages appear on the left.

#' Each metric is shown in a separate facet panel with an

#' independent y-axis scale (\code{scales = "free_y"}) because

#' connectance, nestedness, and modularity occupy different

#' numeric ranges. Each panel shows a coloured line with points.

#' @seealso [compute_network_metrics()]

#' @export

plot_network_metrics_by_age <- function(

    data_network_metrics,

    title = NULL,

    subtitle = NULL) {

  )

    ),

    )

  )

  )

  )

      )

    ) +

    ) +

    ) +

    ) +

    )

}

#' @title Filter Rare Taxa

#' @description

#' Filters out rare taxa from community data based on a minimum proportion

#' threshold. Only taxa meeting or exceeding the threshold are retained.

#' @param data

#' A data frame containing taxon abundance data. Must include a column

#' named 'pollen_prop' with numeric proportions or abundances.

#' @param minimal_proportion

#' Numeric value between 0 and 1 specifying the minimum proportion

#' threshold for retaining taxa. Default is 0.01 (1%).

#' @return

#' A filtered data frame containing only taxa that meet or exceed the

#' minimal proportion threshold. Preserves all original columns.

#' @details

#' The function validates that minimal_proportion is a numeric value

#' between 0 and 1. After filtering, it checks that at least one taxon

#' remains in the dataset. If no taxa meet the threshold, an error is

#' raised suggesting the threshold may be too high.

#' @export

filter_rare_taxa <- function(

    data = NULL,

    minimal_proportion = 0.01) {

  )

  )

  )

  )

    )

  )

}

#' @title Get predictor collinearity

#' @description

#' Analyses collinearity among abiotic predictors in a long-format

#' data frame and returns a `collinear_output` object produced by

#' `collinear::collinear()`. The function pivots `data_source` from

#' long to wide format (one column per variable), removes any `age`

#' column, screens out zero-variance columns, and then performs the

#' collinearity analysis.

#' @param data_source

#' A data frame in long format containing at minimum the columns

#' `abiotic_variable_name` (character, predictor names) and

#' `abiotic_value` (numeric, predictor values).  An optional `age`

#' column is silently dropped before analysis.

#' @return

#' A `collinear_output` object as returned by

#' `collinear::collinear()`. The object contains a `result` element

#' with a `selection` character vector of the non-collinear predictor

#' names that were retained.

#' @details

#' The function validates inputs with `assertthat` and performs

#' post-hoc assertions on the output to guarantee structural

#' integrity before returning.  Missing values are filled with `NA`

#' when pivoting to wide format.  The `age` column is excluded

#' because it is a sampling dimension rather than a predictor.

#' Predictor names are captured before pivoting so that the

#' zero-variance check is scoped to predictor columns only —

#' ID or metadata columns that survive the pivot are never passed

#' to `collinear::collinear()`.  Any predictor whose standard

#' deviation is zero across all samples is dropped and reported via

#' `cli::cli_warn()`.  If no predictor with non-zero variance

#' remains, the function aborts via `cli::cli_abort()`.

#' @seealso

#' [collinear::collinear()] for the underlying collinearity method,

#' [get_abiotic_data()] for producing the expected input format.

#' @export

get_predictor_collinearity <- function(data_source) {

  )

    ),

    )

  )

  # Capture predictor names before pivoting so the variation check

  # is restricted to predictors only (not ID/metadata columns)

    ) |>

    )

    ) |>

    )

  if (

  ) {

        ),

      )

    )

  }

  if (

  ) {

        ),

        )

      )

    )

  }

    ) |>

    )

  )

    )

  )

    )

  )

  )

}

#' @title Compute Standard Errors for a Fitted sjSDM Model

#' @description

#' Post-hoc calculation of standard errors for a fitted sjSDM

#' model using `sjSDM::getSe()`. This is intentionally separated

#' from `fit_jsdm_model()` so that the GPU-based model fitting

#' step and the CPU-based SE computation step can be run

#' independently, allowing the SE step to exploit all available

#' CPU cores without being constrained by the GPU device setting.

#' @param mod_jsdm

#' An object of class `sjSDM` returned by `fit_jsdm_model()` or

#' `sjSDM::sjSDM()`. Must not be `NULL`.

#' @param parallel

#' Number of CPU cores to use for the data loader during SE

#' computation. Passed to the `parallel` argument of

#' `sjSDM::getSe()`. Default is `0L` (no parallelisation).

#' Use `config.model_fitting$n_cores` from the project

#' configuration to exploit all available cores.

#' @param step_size

#' Batch size for stochastic gradient descent used during SE

#' computation. Passed to `sjSDM::getSe()`. Default is `NULL`,

#' which lets sjSDM choose automatically.

#' @param verbose

#' Logical scalar. If `FALSE` (default), all console output

#' from `sjSDM::getSe()` — including Python/reticulate stdout

#' — is suppressed. Set to `TRUE` to see progress messages,

#' which is recommended when running inside a pipeline target

#' to monitor long-running SE computation.

#' @return

#' The input `mod_jsdm` object with its `$se` field populated

#' with the computed standard errors.

#' @details

#' `sjSDM::getSe()` uses CPU for SE computation regardless of the

#' device used for model fitting. Separating SE computation into

#' its own pipeline target therefore allows the caller to pass

#' the full number of available CPU cores via `parallel` without

#' conflict with the GPU device setting used during fitting.

#'

#' When `verbose = FALSE`, output is captured with

#' `reticulate::py_capture_output()` (Python stdout) and

#' `utils::capture.output()` (R stdout) so the console remains

#' clean. When `verbose = TRUE`, both output streams are printed

#' for monitoring progress.

#'

#' If SE computation fails, the function raises an error with

#' an informative message.

#' @seealso

#' `sjSDM::getSe`, `fit_jsdm_model`

#' @export

compute_jsdm_se <- function(

    mod_jsdm = NULL,

    parallel = 0L,

    step_size = NULL,

    verbose = FALSE) {

    )

  )

    )

  )

    )

  )

    )

  )

  if (

  ) {

      )

  } else {

    # Suppress both R stdout and Python/reticulate stdout.

    # Note: reticulate::py_capture_output() is needed because

    #   sjSDM::getSe() routes Python output through reticulate,

    #   which bypasses utils::capture.output(). The outer wrapper

    #   intercepts Python-side output; the inner wrapper

    #   intercepts any R-level stdout from the call.

        {

            )

        },

      )

    )

  }

}

#' @title Classify Taxonomic Resolution

#' @description

#' Classifies taxa in a data frame to a specified taxonomic resolution

#' using a classification table, and aggregates pollen proportions

#' accordingly. Supported resolutions are `kingdom`, `phylum`, `class`,

#' `order`, `family`, `genus`, and `species`.

#' @param data

#' A data frame containing taxon data with columns including 'taxon',

#' 'dataset_name', 'age', and 'pollen_prop'.

#' @param data_classification_table

#' A data frame mapping 'sel_name' to taxonomic levels. Must contain

#' at least one rank column at or below `taxonomic_resolution`

#' (e.g. 'family', 'genus', 'species').

#' @param taxonomic_resolution

#' A character string specifying the finest taxonomic level to use.

#' Must be one of `'kingdom'`, `'phylum'`, `'class'`, `'order'`,

#' `'family'`, `'genus'`, or `'species'`. Taxa will be classified at

#' this rank if possible, or at the coarsest available rank below it

#' if not (fallback behaviour).

#' @return

#' A data frame with taxa classified to the finest available rank at

#' or below `taxonomic_resolution` and pollen proportions aggregated

#' accordingly. The output preserves all dataset_name and age

#' combinations for true negatives.

#' @details

#' Performs a left join to map taxa to all available rank columns up

#' to and including `taxonomic_resolution`. The finest non-NA rank is

#' then selected via `dplyr::coalesce()` applied from finest to

#' coarsest. This means a taxon known only to family when genus is

#' requested will be assigned to its family name rather than dropped.

#' Taxa with no valid classification at any available rank are removed

#' with a `cli::cli_warn()` warning. Taxa that fall back to a coarser

#' rank are reported via `cli::cli_inform()`. Ranks finer than

#' `taxonomic_resolution` (e.g. species when genus is requested) are

#' never used, even when present in the classification table. The

#' NA-drop step prevents a column literally named NA appearing in the

#' community matrix produced by downstream `pivot_wider()` calls.

#' @seealso [filter_non_plantae_taxa()], [filter_rare_taxa()]

#' @export

classify_taxonomic_resolution <- function(

    data,

    data_classification_table,

    taxonomic_resolution) {

  )

    ),

    )

  )

  )

  )

  )

    )

  )

  # All ranks from kingdom up to and including the requested level.

  #   These are the only ranks eligible for fallback assignment.

      )

    ]

  # Subset to ranks actually present in the classification table.

    )

    )

  )

    )

  # Join all available rank columns, then coalesce from finest to

  #   coarsest so each taxon gets the most-specific non-NA label.

    ) |>

      )

    )

  # Report taxa that fell back to a coarser rank.

  #   These are kept in the data but flagged informally.

      ) |>

          )

        )

      )

    }

  }

  # Warn and drop taxa with no valid classification at any available

  #   rank. Without this filter, the NA taxon flows into pivot_wider()

  #   and creates a column literally named NA in the community matrix.

        ),

        )

      )

    )

  }

  # make dummy table with all dataset_name and age combinations

  #   this is needed to ensure that all combinations are present in the

  #   final output to preserve true negative values

    ) |>

    ) |>

    ) |>

    )

}

#' @title Fit an sjSDM Model

#' @description

#' Fits a joint Species Distribution Model (jSDM) using the sjSDM

#' package with specified abiotic, spatial, and error family

#' configurations.

#' @param data_to_fit

#' A list containing the data to fit the model. Must include:

#' - `data_community_to_fit`: matrix of community composition

#'   (constant-presence taxa already removed by

#'   `filter_constant_taxa()`)

#' - `data_abiotic_to_fit`: data frame of abiotic variables,

#'   already scaled by `scale_abiotic_for_fit()`

#' - `data_spatial_to_fit`: data frame of spatial predictors,

#'   already scaled by `scale_spatial_for_fit()`. Required

#'   only when `spatial_method` is not `"none"`.

#' @param sel_abiotic_formula

#' A formula object specifying the abiotic (environmental) predictors

#' @param abiotic_method

#' Method for modeling abiotic effects. One of "linear" (default) or

#' "DNN" (deep neural network)

#' @param sel_spatial_formula

#' A formula object specifying the spatial predictors. Defaults to an

#' interaction between longitude and latitude

#' (`~ 0 + coord_long:coord_lat`). Only used if

#' `spatial_method` is `"linear"` or `"DNN"`.

#' @param spatial_method

#' Method for modeling spatial effects. One of `"linear"` (default),

#' `"DNN"` (deep neural network), or `"none"` (no spatial structure)

#' @param error_family

#' Error family distribution. One of "gaussian" (default) or "binomial".

#' If "binomial", community data is converted to presence/absence and

#' a probit link is used

#' @param device

#' Computing device to use. One of "cpu" (default) or "gpu"

#' @param compute_se

#' Logical indicating whether to compute standard errors inline

#' during model fitting. Default is `FALSE`. Prefer the

#' post-hoc approach via `compute_jsdm_se()` in a separate

#' pipeline target, which allows CPU parallelisation

#' independent of the GPU device setting.

#' @param parallel

#' Number of CPU cores to use for parallel processing.

#' Only applicable if `device = "cpu"`. Default is `0L`

#' (no parallelisation).

#' @param ...

#' Additional arguments passed to `sjSDM::sjSDM()` (e.g.

#' `sampling`, `step_size`, `seed`). Do NOT pass `iter`,

#' `control` or `early_stopping` directly — use the dedicated

#' `iter` and `n_early_stopping` parameters instead.

#' @param iter

#' Positive integer. Number of training epochs. Default is `100L`.

#' @param n_early_stopping

#' Early stopping patience. Controls how many consecutive epochs

#' without improvement are tolerated before training is halted.

#' Three accepted values:

#' - `NULL` (default): auto-compute as `round(iter * 0.20)`, i.e.

#'   at least 20 \% of the epoch budget.

#' - `0L` or negative integer: disables early stopping entirely.

#' - Positive integer: uses `max(value, round(iter * 0.20))` to

#'   ensure patience is never set below 20 \% of `iter`. Passed

#'   as `early_stopping_training` in `sjSDM::sjSDMControl()`.

#' @return

#' An object of class sjSDM containing the fitted model

#' @details

#' This function prepares the data and fits a joint Species Distribution

#' Model using the sjSDM package. The spatial structure is modeled using

#' an interaction term between longitude and latitude coordinates. When

#' binomial error family is specified, the community data is converted

#' to binary presence/absence data.

#'

#' Standard error computation (`compute_se = TRUE`) may fail with certain

#' model configurations, particularly when using DNN methods or complex

#' spatial structures. If SE computation fails, the model will still be

#' returned with a warning.

#' @seealso sjSDM::sjSDM, sjSDM::linear, sjSDM::DNN, compute_jsdm_se

#' @export

fit_jsdm_model <- function(

    data_to_fit = NULL,

    sel_abiotic_formula = NULL,

    abiotic_method = c("linear", "DNN"),

    sel_spatial_formula = as.formula(~ 0 + coord_long:coord_lat),

    spatial_method = c("linear", "DNN", "none"),

    error_family = c("gaussian", "binomial"),

    device = c("cpu", "gpu"),

    parallel = 0L,

    compute_se = FALSE,

    ...,

    iter = 100L,

    n_early_stopping = NULL,

    verbose = FALSE) {

  # Validate `data_to_fit` structure

  )

  )

  )

  # Extract data components

  # Validate extracted data types

  )

  )

  # Validate formula arguments

  )

  )

  # Validate and match character arguments

  )

  )

  # Extract and validate spatial data when needed -----

  # Note: data_spatial_to_fit is pre-scaled by

  #   scale_spatial_for_fit(); no additional scaling applied.

  if (

  ) {

      )

    )

    )

  } else {

  }

  )

  )

  # Validate numeric and logical arguments

  )

  )

  )

    )

  )

    )

  )

  # Handle device/parallel conflict

  if (

  ) {

      )

    )

  }

  # Convert community data to presence/absence for binomial

  if (

  ) {

  } else {

  }

  # Build spatial structure

  # Note: sjSDM::linear/DNN use match.call() internally and re-evaluate

  #   formula symbols in parent.env(environment()) = namespace:sjSDM.

  #   Using do.call passes the formula as an already-evaluated object

  #   (class "formula"), so the bare-name eval branch is never triggered.

  if (

  ) {

        )

      )

        )

      )

  }

  # Build abiotic (environmental) structure

  # Note: data_abiotic is already scaled upstream by

  #   scale_abiotic_for_fit(); no additional scaling is applied.

  if (

  ) {

        )

      )

  } else {

        )

      )

  }

  # Three-tier early stopping patience:

  #  "NULL"  -> auto: round(iter * 0.20), ensuring >= 20% of budget

  #  <= 0  -> 0 (disabled, maps to sjSDMControl's "disabled" value)

  #  > 0   -> max(value, round(iter * 0.20)), floor at 20% of iter

    ) {

    ) {

    } else {

      )

    }

      ...

    )

}

#' @title Compute Spatiotemporal Moran Eigenvector Maps

#' @description

#' Computes Moran Eigenvector Maps (MEMs) from a 3-D

#' coordinate matrix `(x_km, y_km, age_kyr)` built at the

#' **sample level** (one row per site × time-slice).

#' All three dimensions are z-scored before eigenvector

#' computation so that spatial and temporal extents

#' contribute equally to the Euclidean distance structure.

#' Returns the first `n_mev` eigenvectors as a data frame

#' with row names in `"dataset_name__age"` format, ready

#' for `scale_spatial_for_fit()` and

#' `assemble_data_to_fit()`.

#' @param data_coords_projected

#' A data frame with `dataset_name` as row names and

#' columns `coord_x_km` and `coord_y_km`, as returned by

#' `project_coords_to_metric()`. One row per unique

#' site/core.

#' @param data_sample_ids

#' A data frame with columns `dataset_name` and `age`

#' giving the valid (site × time-slice) pairs to model,

#' as returned by `align_sample_ids()`.

#' @param n_mev

#' A positive integer giving the number of eigenvectors to

#' return. If it exceeds the number of positive Moran

#' eigenvectors produced by `sjSDM::generateSpatialEV()`

#' for the 3-D coordinate matrix, it is automatically

#' clamped down to that count and a `cli::cli_warn()`

#' message is emitted. Default is `20L`.

#' @return

#' A data frame with row names `"<dataset_name>__<age>"`

#' and `n_mev` columns named `mev_1`, `mev_2`, …,

#' `mev_{n_mev}`. The row order follows

#' `data_sample_ids` sorted by `dataset_name` then `age`,

#' matching the ordering produced by

#' `prepare_abiotic_for_fit()` and

#' `prepare_community_for_fit()`.

#' @details

#' Unlike `compute_spatial_mev()`, which operates on

#' unique site coordinates and must be expanded to sample

#' level via `prepare_spatial_predictors_for_fit()`, this

#' function builds the eigenvectors directly on the

#' sample-level 3-D coordinate matrix. This means each

#' observation (site × age) gets a unique row in the MEV

#' matrix, and within-core temporal autocorrelation is

#' captured alongside between-site spatial autocorrelation.

#'

#' The 3-D coordinate matrix is constructed as follows:

#' 1. Join `data_coords_projected` (x_km, y_km) onto

#'    `data_sample_ids` to get one row per sample.

#' 2. Convert age (years BP) to kiloyears: `age_kyr =

#'    age / 1000`.

#' 3. Z-score each dimension independently so that the

#'    spatial and temporal extents are on the same scale

#'    before Euclidean distances are computed.

#'

#' The z-scoring step is critical: if age is left in

#' kiloyears and coordinates in km, the time axis will

#' dominate or be negligible depending on the data range,

#' producing eigenvectors that are either purely temporal

#' or purely spatial.

#'

#' Because the eigenvectors are computed at the sample

#' level, **no further expansion step is required**.

#' Pass the returned data frame directly to

#' `scale_spatial_for_fit()`.

#' @seealso

#'   [compute_spatial_mev()],

#'   [project_coords_to_metric()],

#'   [align_sample_ids()],

#'   [scale_spatial_for_fit()],

#'   [assemble_data_to_fit()]

#' @export

compute_spatiotemporal_mev <- function(

    data_coords_projected = NULL,

    data_sample_ids = NULL,

    n_mev = 20L) {

  )

    ),

    )

  )

  )

    ),

    )

  )

  )

  # 1. Build sample-level 3-D coordinate data frame -----

    ) |>

    )

    )

  )

  # 2. Z-score each dimension separately -----

    )

  # 3. Compute Moran eigenvectors on 3-D matrix -----

    )

  # Force to matrix: sjSDM returns a vector when exactly

  # one positive eigenvalue is found (drops the dimension)

  # 4. Post-call validation: clamp n_mev if needed -----

  if (

  ) {

      )

    )

  }

  # 5. Select first n_mev columns -----

  # 6. Coerce to data frame with named columns -----

}

#' @title Save Progress Visualisation

#' @description

#' Generates a visualisation of project progress and saves it as HTML and PNG.

#' @param sel_script

#' The script file to be visualised.

#' @param sel_store

#' Path to the targets store directory. Defaults to the value from the active

#' configuration key "target_store".

#' @param output_file

#' The name of the output file (default: "project_status").

#' @param output_dir

#' Directory where the output files will be saved

#' (default: "Documentation/Progress").

#' @param background_color

#' Background color for the visualisation

#' (default: `"#141B22"` — BIODYNAMICS brand Console Panel; one step lighter

#' than the page background, used for cards and panels).

#' @param physics

#' Logical indicating whether to enable physics simulation in the network

#' graph (default: TRUE).

#' @param level_separation

#' Level separation for the visualisation graph (default: 250).

#' @return

#' No return value. Called for side effects: saves HTML and PNG files to a

#' store-specific subdirectory within `output_dir`.

#' @details

#' Uses `targets::tar_visnetwork` to create a network graph and saves two

#' HTML files (full and targets-only) using `visNetwork::visSave`, plus a

#' static PNG via `webshot2::webshot`. If the browser-backed PNG export fails,

#' the function emits a warning and keeps the HTML outputs so progress saving

#' does not abort the surrounding pipeline run. Files are written to

#' `output_dir/<store_name>/` where `<store_name>` is the last path segment

#' of `sel_store` after the `targets/` directory.

#'

#' The targets store is read only once. `targets_only = FALSE` retrieves the

#' full network; the targets-only variant for the static PNG graph is derived

#' by filtering non-target nodes (`type` values `"function"`, `"object"`,

#' `"value"`) from the same object, avoiding a second store read. If the

#' `type` column is absent (targets internal change), the function falls back

#' to a second `tar_visnetwork()` call so correctness is preserved.

#'

#' When `physics = TRUE`, the static graph for the PNG uses a

#' `stabilizationIterationsDone` vis.js event to disable physics as soon as

#' the silent background layout pass finishes. This ensures the screenshot

#' captures the fully settled layout without requiring a long wait or

#' rendering the graph twice.

#' @export

save_progress_visualisation <- function(

    sel_script,

    sel_store = get_active_config("target_store"),

    output_file = "project_status",

    output_dir = here::here(

      "Documentation/Progress"

    ),

    background_color = "#141B22",

    physics = TRUE,

    level_separation = 250) {

  # test to make sure {pandoc} is installed for webshot to work

  if (

  ) {

  }

  # Read the targets store once. targets_only = FALSE gives the full network;

  #   the targets-only static variant is derived below by filtering the same

  #   object, avoiding a second expensive store read.

    )

    # Apply BIODYNAMICS brand theme over the targets-generated colors.

    # Node fill colors (status: up-to-date / outdated / etc.) are set

    #   internally by targets and cannot be overridden here; only font

    #   and edge styling are applied.

      )

    ) |>

      ),

      )

    )

  # Derive the targets-only static base graph from the raw object already

  #   in memory. targets sets type = "function" / "object" / "value" for

  #   non-target nodes; filtering them mirrors targets_only = TRUE.

  #   The fallback re-calls tar_visnetwork() if targets ever drops the

  #   type column, so correctness is preserved across package changes.

  if (

  ) {

      ) |>

      )

      )

  } else {

    # Fallback: targets internals changed and type column is absent.

      )

  }

      )

    ) |>

      ),

      )

    ) |>

    # Freeze the graph the moment background stabilization finishes so the

    #   webshot captures the settled layout, not a mid-simulation frame.

    #   vis.js fires `stabilizationIterationsDone` after the silent pre-layout

    #   pass (fast, < 1 s), then starts the animated physics — disabling

    #   physics here prevents any further movement before the screenshot.

      }"

    )

    )

  # need to create the output directory if it doesn't exist

  if (

  ) {

  }

  )

  )

  # PNG export depends on launching a browser through chromote/webshot2,

  #   which can fail in interactive sessions even when the HTML outputs save.

  # `delay = 1` gives the browser time to complete the background

  #   stabilization pass before the screenshot is taken; the visEvents

  #   handler then keeps the graph frozen so the captured layout is stable.

    ),

          )

        )

      )

    }

  )

}

#' @title Interpolate 2-D Spatial MEVs to Prediction Grid

#' @description

#' Uses Inverse Distance Weighting (IDW, power = 2) to

#' approximate Moran Eigenvector Map (MEM) values from

#' training-site locations to arbitrary prediction locations,

#' then scales the result using training spatial scale

#' attributes so that the interpolated predictors are on the

#' same scale as those seen during model fitting.

#' @param data_coords_projected_train

#' A data frame with `coord_x_km` and `coord_y_km` columns

#' and `dataset_name` as row names, as returned by

#' `project_coords_to_metric()`. One row per unique training

#' site.

#' @param data_mev_core

#' A data frame with unscaled MEV columns (`mev_1`, `mev_2`,

#' …) and `dataset_name` as row names, as returned by

#' `compute_spatial_mev()`.

#' @param data_coords_projected_pred

#' A data frame with `coord_x_km` and `coord_y_km` columns

#' and arbitrary row names identifying prediction locations

#' (e.g. `"grid_1"`, `"grid_2"`), as returned by

#' `project_coords_to_metric()`.

#' @param spatial_scale_attributes

#' A named list of `"scaled:center"` and `"scaled:scale"`

#' attributes per MEV column, as returned by

#' `scale_spatial_for_fit()` in the `spatial_scale_attributes`

#' element. Used to bring interpolated MEV values onto the

#' same scale as the training spatial predictors.

#' @return

#' A data frame with the same row names as

#' `data_coords_projected_pred` and one column per MEV

#' (names matching `data_mev_core`). All columns are scaled

#' to match the training MEV distribution.

#' @details

#' MEMs are eigenvectors of the spatial connectivity matrix

#' at training sites and cannot be analytically evaluated at

#' new locations. IDW (power = 2) with a small epsilon

#' (1e-10) to prevent division-by-zero provides a smooth

#' spatial interpolation.

#'

#' This function handles the **2-D spatial case** only

#' (x_km, y_km). For models fitted with

#' `spatial_mode = "spatiotemporal"` use

#' `interpolate_st_mev_to_grid()` instead.

#' @seealso

#'   [compute_spatial_mev()],