R/tabmachine.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

# ============================================================================
# R4VN tabmachine() - comprehensive machine-learning workflow
# Package-ready implementation
# ============================================================================

#' Comprehensive machine-learning analysis with publication-ready reporting
#'
#' @description
#' `tabmachine()` is the comprehensive machine-learning command in R4VN. It is
#' designed for health and biomedical research where the user needs a complete,
#' reproducible workflow rather than only a fitted prediction model.
#'
#' The function can detect the prediction task, split development data,
#' preprocess predictors, handle missing values, encode categorical predictors,
#' detect problematic predictors, standardize predictors when required,
#' transform skewed numeric variables when explicitly requested, address class
#' imbalance, select predictors, tune candidate algorithms, perform cross-
#' validation, compare models, select a final model, determine a classification
#' threshold using training data only, evaluate the untouched test set, calculate
#' confidence intervals for performance measures whenever a defensible interval
#' is implemented, assess calibration and clinical utility, calculate variable
#' importance, and generate prediction-ready model objects.
#'
#' The central R4VN principle is that automation must remain transparent.
#' `"auto"` may choose an analysis action, but every action is stored in the
#' returned object and shown in the report. Preprocessing, feature selection,
#' class balancing, tuning, and threshold optimization are learned from training
#' data only. The held-out test data are not used to make those decisions.
#'
#' @param outcome Outcome variable. Supply an unquoted variable name or a
#'   one-element character name. Binary, multiclass, and continuous outcomes are
#'   supported.
#' @param x Candidate predictors. Use `vars(age, sex, bmi)`, a character vector,
#'   a single bare variable, or `.` to use every variable except `outcome` and
#'   `exclude`. If omitted, all other variables are used.
#' @param data Optional data frame. If omitted, the active R4VN data selected by
#'   `usedf()` or `opendata(..., active = TRUE)` is used.
#' @param exclude Optional predictors to exclude, supplied as `vars(...)`, a bare
#'   variable, or a character vector. Typical examples are identifiers, names,
#'   dates that would leak future information, or variables unavailable at the
#'   intended prediction time.
#' @param task Prediction task: `"auto"` (default), `"binary"`, `"multiclass"`,
#'   or `"regression"`. In auto mode, an outcome with two observed levels is
#'   binary; a factor/character outcome or a low-cardinality integer outcome is
#'   multiclass; otherwise a numeric outcome is regression.
#' @param event Positive/event level for binary classification. If omitted,
#'   R4VN recognizes common positive encodings such as 1, TRUE, Yes, Positive,
#'   Case, or Co; otherwise the second observed level is used. The chosen event
#'   is always reported.
#'
#' @param preprocess Preprocessing policy. `"auto"` performs safe structural
#'   checks, imputation, factor-level alignment, dummy encoding, zero/near-zero
#'   variance removal, and optional correlation filtering. `"none"` disables
#'   automatic structural filtering but still creates a usable model matrix.
#' @param missing Missing-value handling for predictors: `"auto"`, `"median"`,
#'   `"mean"`, `"mode"`, or `"complete"`. `"auto"` uses median for numeric and
#'   mode for categorical predictors. Imputation values are estimated on each
#'   training fold and then applied to its validation fold. Outcome missingness
#'   is never imputed.
#' @param missing_max Maximum allowed proportion missing in a candidate predictor
#'   before automatic structural filtering removes it. Default 0.50.
#' @param encode Encoding of categorical predictors. Currently `"auto"` and
#'   `"dummy"` use treatment/dummy coding through a training-derived
#'   `model.matrix` blueprint. New/unseen levels in validation data are mapped
#'   to the training reference fallback and reported.
#' @param standardize Standardization policy: `"auto"`, `"none"`, `"z"`,
#'   `"minmax"`, or `"robust"`. In auto mode, z-standardization is applied only
#'   inside algorithms that materially benefit from scaling (penalized models,
#'   SVM, KNN, and multinomial neural optimization). Tree-based models receive
#'   the unscaled encoded matrix.
#' @param transform Numeric transformation: `"none"` (default), `"auto"`,
#'   `"log"`, or `"yeojohnson"`. `"auto"` is intentionally conservative and
#'   currently performs no transformation unless a future R4VN rule explicitly
#'   justifies one. `"log"` uses a training-derived shift when necessary.
#'   `"yeojohnson"` estimates lambda on the training data by profile likelihood.
#' @param outlier Outlier policy: `"none"` (default), `"detect"`, `"winsor"`, or
#'   `"robust"`. Detection uses training-data IQR rules and never deletes rows.
#'   `"winsor"` caps numeric values at the 1st and 99th training percentiles;
#'   `"robust"` caps at median +/- 5 MAD. Cut points are then reused unchanged
#'   for validation/test data.
#' @param corr Correlation filtering for numeric candidate predictors. `"auto"`
#'   uses 0.95, a numeric value in (0,1) supplies a custom absolute-correlation
#'   threshold, and `FALSE`/`"none"` disables correlation filtering. Filtering
#'   is learned separately inside each training fold.
#'
#' @param feature Feature engineering: `"none"` (default), `"auto"`,
#'   `"polynomial"`, `"interaction"`, or `"all"`. `"polynomial"` adds powers
#'   of numeric predictors up to `degree`; `"interaction"` adds pairwise
#'   products among numeric predictors; `"all"` adds both. `"auto"` is
#'   intentionally conservative and currently behaves as `"none"`. All
#'   engineered features are created from a fold-specific training blueprint.
#' @param degree Highest polynomial degree for numeric feature engineering.
#'   Default 2; values 2 or 3 are supported.
#' @param reduce Dimensionality reduction: `"none"` (default), `"auto"`, or
#'   `"pca"`. PCA is learned only on the training fold after encoding/feature
#'   engineering. `"auto"` uses PCA only in strongly high-dimensional settings
#'   (more encoded features than training observations and at least 50 encoded
#'   features); otherwise it remains off to preserve clinical interpretability.
#' @param variance Target cumulative variance retained by PCA. Default 0.95.
#'
#' @param select Feature-selection strategy: `"auto"`, `"none"`, `"filter"`,
#'   `"lasso"`, `"stepwise"`, `"importance"`, `"rfe"`, `"boruta"`, or
#'   `"compare"`. `"auto"` deliberately resolves to the native R4VN `"filter"`
#'   strategy so the default workflow never depends on an optional selection
#'   package. LASSO uses `glmnet` when explicitly requested. Boruta uses the optional `Boruta`
#'   package. `"compare"` compares available strategies by cross-validation on
#'   the development training set and does not inspect the held-out test set.
#' @param nfeatures Number of predictors to retain for importance/RFE selection,
#'   or `"auto"`. With `"auto"`, R4VN chooses a compact candidate size based on
#'   the training sample size and number of available encoded predictors.
#' @param simplify Logical; after choosing the best algorithm, search for a
#'   smaller predictor set whose development cross-validated performance is
#'   within `simplify_tol` of the larger selected model. This step is training-
#'   only. Default TRUE.
#' @param simplify_tol Maximum acceptable loss in the primary metric when
#'   preferring a smaller model. For metrics where larger is better this is an
#'   absolute decrease; for RMSE/MAE it is an absolute increase. Default 0.01.
#'
#' @param balance Class-imbalance handling for binary classification:
#'   `"auto"`, `"none"`, `"weight"`, `"up"`, `"down"`, `"smote"`,
#'   `"adasyn"`, `"rose"`, or `"compare"`. Balancing is applied only to the
#'   training portion of each resample. Test/validation observations are never
#'   resampled. `"auto"` uses class weighting when the minority class is below
#'   20 percent and otherwise leaves the sample unchanged. `"compare"` compares
#'   available approaches on development cross-validation.
#' @param balance_target Target minority proportion after sampling. Default 0.50.
#' @param neighbors Number of nearest neighbors for native SMOTE/ADASYN. Default 5.
#'
#' @param split Development/test split. A single number such as 0.80 means 80%
#'   training and 20% untouched test data. `FALSE` uses all development data for
#'   model development (appropriate when a separate `validation` data set is
#'   supplied). A length-three vector such as `c(.70,.15,.15)` creates training,
#'   internal validation, and test partitions; the internal validation portion
#'   is combined with training for final refitting after model decisions are
#'   completed, while the test portion remains untouched.
#' @param folds Number of cross-validation folds in the development training
#'   sample. Default 10. Classification folds are stratified when possible.
#' @param repeats Number of repeated cross-validation repetitions. Default 1.
#' @param nested Logical; if TRUE, hyperparameter tuning is repeated within each
#'   outer cross-validation fold. This is computationally expensive but gives a
#'   less optimistic development estimate. Regardless of this option, the final
#'   held-out test evaluation remains untouched by tuning.
#' @param seed Optional random seed used for splitting, resampling, tuning, and bootstrap. The default `NULL` does not set a seed.
#'
#' @param method Algorithms to fit. `"auto"` deliberately uses a compact
#'   low-dependency R4VN core set so that a routine analysis does not require
#'   installation of an ML ecosystem: logistic regression plus a decision tree
#'   for binary outcomes, linear regression plus a decision tree for continuous
#'   outcomes, and multinomial regression plus a decision tree and R4VN-native
#'   KNN for multiclass outcomes. The native KNN fallback keeps multiclass auto
#'   usable even if recommended modelling packages are unavailable. `"all"` tries every
#'   supported engine that is already installed. Optional engines are never
#'   installed automatically. Methods may also be supplied explicitly:
#'   `"logistic"`, `"linear"`, `"multinom"`, `"lasso"`, `"ridge"`, `"elastic"`,
#'   `"tree"`, `"rf"`, `"xgb"`, `"svm"`, `"knn"`, `"naive"`.
#' @param tune Hyperparameter tuning. `TRUE`/`"grid"` uses compact clinically
#'   practical grids; `"random"` samples candidate combinations; FALSE uses
#'   defaults. Penalized `glmnet` engines select lambda internally by CV.
#' @param tune_n Maximum random-search combinations when `tune = "random"`.
#' @param metric Primary model-selection metric. `"auto"` uses ROC-AUC for most
#'   binary problems, PR-AUC when the event is uncommon, macro F1 for multiclass,
#'   and RMSE for regression. Other supported names include `auc`, `pr_auc`,
#'   `accuracy`, `balanced_accuracy`, `sensitivity`, `specificity`, `f1`,
#'   `brier`, `logloss`, `macro_f1`, `weighted_f1`, `macro_auc`,
#'   `macro_pr_auc`, `rmse`, `mae`, `r2`, and `mape`.
#'
#' @param threshold Binary classification threshold. A numeric value between 0 and 1 fixes the
#'   cutoff. `"auto"`/`"youden"` maximizes sensitivity + specificity - 1 using
#'   out-of-fold predictions from the development training set. `"f1"` maximizes
#'   F1. `"sens"` chooses the most specific threshold achieving at least
#'   `target_sens`; `"spec"` chooses the most sensitive threshold achieving at
#'   least `target_spec`. The test set is never used to choose the threshold.
#' @param target_sens Target sensitivity used when `threshold = "sens"`.
#' @param target_spec Target specificity used when `threshold = "spec"`.
#'
#' @param ci Logical; calculate 95% confidence intervals (or the level supplied
#'   by `ci_level`) for performance measures whenever an implemented interval is
#'   statistically meaningful. Default TRUE.
#' @param ci_level Confidence level, default 0.95.
#' @param boot Number of bootstrap replicates for performance measures whose
#'   interval has no preferred closed-form method. Default 1000. For final
#'   publication analyses, 2000 or more may be preferred when runtime permits.
#' @param calibration Logical; for binary classification, calculate calibration
#'   intercept, calibration slope, Brier score, and calibration-curve data.
#' @param decision Logical; for binary classification, calculate decision-curve
#'   net benefit for the final model over `decision_thresholds`.
#' @param decision_thresholds Probability thresholds for decision-curve analysis.
#'   Default `seq(.01, .99, .01)`.
#' @param learning Logical; calculate training-size learning-curve summaries for
#'   the final model. Default FALSE because it can be computationally expensive.
#'
#' @param importance Logical; calculate permutation importance for the final
#'   model. Default TRUE. Importance is grouped back to original predictor names
#'   when dummy variables were created.
#' @param importance_repeats Number of repeated permutations used to stabilize
#'   permutation importance. Default 20.
#' @param explain Logical; retain explanation data and display the principal
#'   importance/calibration information in the report. Default TRUE.
#' @param shap Logical; calculate SHAP-like contribution output when a supported
#'   engine is available. Native XGBoost `predcontrib` is used for XGBoost;
#'   other engines are left as `NULL` unless a future R4VN explainer is added.
#' @param pdp Logical or character vector. TRUE calculates partial-dependence
#'   data for up to the five most important original numeric predictors;
#'   a character vector requests specific predictors.
#'
#' @param validation Optional external validation data frame. It must contain the
#'   same outcome and required predictors. It is never used for preprocessing,
#'   selection, tuning, balancing, threshold selection, or final model fitting.
#' @param predict Optional new data frame for predictions after the final model
#'   is fitted. Predictions are returned in `$predictions`.
#' @param id Optional identifier variable to copy into prediction output.
#'
#' @param show Logical; open the publication-style HTML report in the Viewer.
#'   Default TRUE.
#' @param plot Logical or character vector controlling figures embedded in the
#'   HTML Viewer. TRUE embeds every principal figure available for the analysis.
#'   Character values may include `"roc"`, `"pr"`, `"calibration"`,
#'   `"threshold"`, `"confusion"`, `"importance"`, `"decision"`, `"learning"`,
#'   `"observed"`, `"residual"`, `"pdp"`, and `"shap"`; `"all"` requests every
#'   available figure.
#'   Any figure sent to the Plot pane through `plot_display` is also retained in
#'   the Viewer whenever plotting is enabled.
#' @param plot_display Figure type(s) also drawn in the interactive R/RStudio
#'   Plot pane and Plot history. `"auto"` (default) draws the most useful
#'   diagnostic figures for the task. Use `"all"` for every available figure,
#'   `"none"` to keep figures only in the Viewer, or a character vector of figure
#'   names. This option never changes model fitting or performance estimates.
#' @param plot_args Named list of base-graphics options used for Viewer and Plot
#'   pane figures. Options may be common (for example
#'   `list(font_family="Arial")`) or nested, for example
#'   `list(all=list(font_family="Arial"), roc=list(lwd=3))`.
#' @param strict Logical. If FALSE (default), a non-essential figure that cannot
#'   be drawn is skipped with a warning while the analysis result is retained.
#'   If TRUE, such a plotting error stops the call.
#' @param console Logical; also print a compact console summary. Default FALSE.
#' @param digit Number of digits displayed for estimates. Default 3.
#' @param title Optional report title.
#' @param ai FALSE, TRUE, or a named R4VN AI endpoint. When R4VN `aiask()` is
#'   available, interpretation is requested only after the quantitative analysis
#'   object has been created. `tabmachine()` never calls AI when `ai = FALSE`;
#'   when enabling AI, the privacy behavior of the configured `aiask()` endpoint
#'   should be reviewed because fitted model objects may contain development
#'   information needed for prediction.
#' @param ... Reserved for future model-engine options.
#'
#' @details
#' ## 1. What `tabmachine()` regards as a complete ML workflow
#'
#' A typical call performs the following sequence:
#'
#' 1. resolve active/explicit data and variable labels;
#' 2. validate the outcome and candidate predictors;
#' 3. create an untouched test partition;
#' 4. inside training resamples, learn imputation/transformation/encoding rules;
#' 5. remove structural problems such as constant predictors;
#' 6. optionally create polynomial/interaction features and/or training-fold PCA;
#' 7. apply feature selection inside the training portion of each resample;
#' 8. apply class balancing only inside the training portion of each resample;
#' 9. tune and compare candidate algorithms;
#' 10. choose the final algorithm using development data only;
#' 11. choose a binary classification threshold from training out-of-fold
#'     predictions only;
#' 12. refit the selected pipeline on the full development training sample;
#' 13. evaluate the untouched test set and calculate confidence intervals;
#' 14. optionally validate on a completely external data set;
#' 15. assess calibration, decision-curve utility, and variable importance;
#' 16. store a prediction blueprint for future `predict()` calls.
#'
#' ## 2. Confidence intervals
#'
#' `ci = TRUE` is the default because R4VN is intended for scientific reporting.
#' The implementation does not attach a made-up CI to a quantity merely because
#' a point estimate exists. Methods currently used are:
#'
#' * sensitivity, specificity, PPV, NPV, accuracy, and prevalence: Wilson
#'   binomial intervals;
#' * ROC-AUC: DeLong interval through `pROC` when available; otherwise a
#'   stratified nonparametric bootstrap interval;
#' * an optimized binary classification threshold: stratified bootstrap interval; a user-fixed threshold has no sampling CI because it is specified rather than estimated;
#' * balanced accuracy, F1, MCC, kappa, PR-AUC, Brier score, log loss and other
#'   derived binary metrics: paired-observation nonparametric bootstrap;
#' * calibration intercept and slope: model-based Wald intervals, with bootstrap
#'   fallback when the calibration model is unstable;
#' * RMSE, MAE, R-squared and MAPE: nonparametric bootstrap over test subjects;
#' * multiclass accuracy/balanced accuracy/macro-F1/weighted-F1, macro one-vs-rest
#'   AUC/PR-AUC, and log loss: nonparametric bootstrap over test subjects;
#' * class-specific multiclass one-vs-rest sensitivity, specificity, PPV, NPV,
#'   accuracy and prevalence: Wilson intervals where the denominator is fixed;
#'   class-specific AUC uses DeLong through `pROC` when available, while PR-AUC,
#'   F1 and other derived class measures use nonparametric bootstrap;
#' * decision-curve net benefit: pointwise nonparametric bootstrap when `ci=TRUE`.
#'
#' These intervals quantify uncertainty in performance on the evaluation sample
#' conditional on the fitted development procedure. They do not replace full
#' external validation or transportability assessment.
#'
#' ## 3. Leakage prevention
#'
#' The most important implementation rule is that no data-dependent
#' preprocessing action is estimated on the test set. Imputation values,
#' transformations, factor levels, scaling parameters, feature selection,
#' balancing, tuning, and threshold optimization are fitted using training data.
#' During CV, those operations are refitted inside each training fold before
#' predictions are made for the corresponding validation fold.
#'
#' ## 4. Class imbalance
#'
#' `balance = "weight"` is the preferred automatic strategy because it does not
#' fabricate observations. `up`, `down`, native numeric-space `smote`, native
#' `adasyn`, and optional `ROSE` are available for explicit experiments. Synthetic
#' sampling occurs after fold-specific numeric encoding and is therefore applied
#' only to the analysis portion of a resample. The original test prevalence is
#' preserved for evaluation, PPV/NPV, calibration, Brier score and decision curves.
#'
#' ## 5. Feature selection versus feature importance
#'
#' `select` controls which predictors are allowed into the fitted model.
#' `importance` explains which predictors contribute most to a fitted final
#' model. They are intentionally separate concepts. A variable may survive
#' selection yet have weak final permutation importance, and correlated
#' variables may share or exchange importance.
#'
#' ## 6. Model comparison
#'
#' Development CV estimates are useful for choosing an algorithm. The untouched
#' test estimate is the primary internal-validation result. When several models
#' are evaluated on the same test subjects, `$performance` includes a CI for each
#' available metric. `$model_difference` additionally stores paired bootstrap
#' differences in the primary metric between each model and the selected model.
#'
#' ## 7. External validation
#'
#' Supply `validation = external_data` to evaluate the finalized development
#' pipeline without refitting it. External performance and its CIs are stored
#' separately. If external validation is the principal evaluation, use
#' `split = FALSE` to use all development observations for model development.
#'
#' ## 8. Optional packages and the low-dependency default
#'
#' R4VN deliberately does not require every ML engine for every user. The
#' default `method = "auto"` is intentionally based on base/recommended R
#' engines and does not require `caret`, `tidymodels`, `recipes`, `yardstick`,
#' or a collection of boosting/forest packages. Optional modelling engines are
#' checked only when the user explicitly requests them or uses `method = "all"`.
#' `glmnet` supplies penalized models; `ranger` random forests; `xgboost`
#' gradient boosting; `e1071` SVM and naive Bayes; `Boruta` Boruta feature
#' selection; `ROSE` ROSE sampling; and `pROC` DeLong ROC intervals. When
#' `pROC` is absent, R4VN uses its native bootstrap AUC interval. The returned
#' `$engines` table records what was requested, installed, and actually used.
#'
#' @return
#' Invisibly returns an object of class `c("r4vn_machine", "r4vn_tab")` with
#' major components:
#' \describe{
#'   \item{overview}{Data/task overview.}
#'   \item{engines}{Supported algorithms, package requirements, availability, and engines actually used.}
#'   \item{preprocessing}{Auditable preprocessing decisions.}
#'   \item{balance}{Chosen imbalance strategy and comparison when requested.}
#'   \item{selection}{Chosen feature-selection strategy and selected predictors.}
#'   \item{tuning}{Selected hyperparameters for each candidate algorithm.}
#'   \item{cv_performance}{Development cross-validation performance.}
#'   \item{performance}{Held-out test performance in long format with CI columns.}
#'   \item{comparison}{Publication-ready model comparison table.}
#'   \item{model_difference}{Paired difference in the primary metric versus the
#'     selected model, with bootstrap CI where available.}
#'   \item{overfitting}{Development-versus-evaluation comparison for the primary metric.}
#'   \item{coefficients}{For an unpenalized logistic or linear final model, model coefficients with 95% CI; logistic coefficients are exponentiated to OR.}
#'   \item{best}{Name of the selected algorithm.}
#'   \item{final}{Final fitted pipeline/model object.}
#'   \item{threshold}{Training-derived classification threshold information.}
#'   \item{confusion}{Final binary or multiclass confusion matrix counts.}
#'   \item{class_performance}{For multiclass outcomes, one-vs-rest class-specific discrimination and classification measures with 95% CI.}
#'   \item{calibration}{Calibration statistics and curve data.}
#'   \item{decision}{Decision-curve data.}
#'   \item{importance}{Grouped permutation importance.}
#'   \item{explanation}{Convenience list collecting importance, SHAP and PDP outputs when `explain = TRUE`.}
#'   \item{shap}{Native XGBoost contribution matrix when requested/supported.}
#'   \item{pdp}{Partial-dependence data when requested.}
#'   \item{external_performance}{External-validation performance when supplied.}
#'   \item{external_confusion}{External binary or multiclass confusion matrix when applicable.}
#'   \item{external_class_performance}{External multiclass one-vs-rest performance with CI.}
#'   \item{external_calibration}{External binary calibration statistics/curve when requested.}
#'   \item{external_decision}{External binary decision-curve data when requested.}
#'   \item{predictions}{Predictions for `predict=` data when supplied.}
#'   \item{plots}{Data required to replay publication plots.}
#'   \item{plot_titles}{Stable publication titles for available figure types.}
#'   \item{tables}{Named publication-ready tables used by R4VN/Studio/export.}
#'   \item{html}{Finished HTML report.}
#'   \item{settings}{Complete analysis settings.}
#'   \item{notes}{Methodological notes and any optional-engine skips.}
#' }
#'
#' @section Examples - binary classification:
#' ```
#' set.seed(123)
#' n <- 500
#' d <- data.frame(
#'   patient_id = seq_len(n),
#'   age = rnorm(n, 45, 12),
#'   sex = factor(sample(c("Female", "Male"), n, TRUE)),
#'   bmi = rnorm(n, 23, 3.5),
#'   smoke = factor(sample(c("No", "Yes"), n, TRUE, c(.75, .25)))
#' )
#' lp <- -6 + .055*d$age + .10*d$bmi + .65*(d$smoke == "Yes")
#' d$hypertension <- factor(rbinom(n, 1, plogis(lp)), 0:1, c("No", "Yes"))
#'
#' m <- tabmachine(
#'   hypertension,
#'   x = vars(age, sex, bmi, smoke),
#'   data = d,
#'   event = "Yes",
#'   method = "logistic",
#'   tune = FALSE,
#'   boot = 200,
#'   show = FALSE,
#'   plot = FALSE
#' )
#' m$comparison
#' m$performance
#' ```
#'
#' @section Examples - automatic ML comparison:
#' ```
#' \donttest{
#' # The automatic comparison uses the low-dependency R4VN core set.
#' # Optional engines join only when explicitly requested or method = "all".
#' m <- tabmachine(
#'   hypertension,
#'   x = vars(age, sex, bmi, smoke),
#'   data = d,
#'   event = "Yes",
#'   method = "auto",
#'   select = "auto",
#'   balance = "auto",
#'   ci = TRUE,
#'   boot = 2000
#' )
#' plot(m, "roc")
#' plot(m, "calibration")
#' plot(m, "importance")
#' }
#' ```
#'
#' @section Examples - class imbalance and feature selection:
#' ```
#' \donttest{
#' m2 <- tabmachine(
#'   hypertension,
#'   x = .,
#'   data = d,
#'   exclude = vars(patient_id),
#'   event = "Yes",
#'   balance = "compare",
#'   select = "compare",
#'   metric = "pr_auc",
#'   nested = TRUE
#' )
#' m2$balance
#' m2$selection
#' }
#' ```
#'
#' @section Examples - continuous outcome:
#' ```
#' set.seed(321)
#' r <- data.frame(
#'   age = rnorm(400, 50, 14),
#'   bmi = rnorm(400, 24, 4),
#'   sex = factor(sample(c("Female", "Male"), 400, TRUE))
#' )
#' r$sbp <- 75 + .65*r$age + 1.15*r$bmi + 4*(r$sex == "Male") + rnorm(400, 0, 9)
#'
#' mr <- tabmachine(
#'   sbp,
#'   x = vars(age, bmi, sex),
#'   data = r,
#'   method = "linear",
#'   tune = FALSE,
#'   boot = 200,
#'   show = FALSE,
#'   plot = FALSE
#' )
#' mr$performance
#' ```
#'
#' @section Examples - active data and new predictions:
#' ```
#' \donttest{
#' usedf(d)
#' fit <- tabmachine(
#'   hypertension,
#'   x = vars(age, sex, bmi, smoke),
#'   event = "Yes",
#'   method = "logistic"
#' )
#'
#' newpatients <- data.frame(
#'   age = c(35, 68),
#'   sex = factor(c("Female", "Male"), levels = levels(d$sex)),
#'   bmi = c(21, 31),
#'   smoke = factor(c("No", "Yes"), levels = levels(d$smoke))
#' )
#' predict(fit, newpatients)
#' }
#' ```
#'
#' @section Examples - external validation:
#' ```
#' \donttest{
#' dev <- d[1:350, ]
#' ext <- d[351:500, ]
#' me <- tabmachine(
#'   hypertension,
#'   x = vars(age, sex, bmi, smoke),
#'   data = dev,
#'   validation = ext,
#'   split = FALSE,
#'   event = "Yes",
#'   method = "auto"
#' )
#' me$external_performance
#' }
#' ```
#'
#' @section Examples - Viewer and Plot pane together:
#' ```
#' \donttest{
#' # Every requested figure remains in the Viewer. The selected figures are
#' # also added to Plot history so Previous/Next can be used in RStudio.
#' mv <- tabmachine(
#'   hypertension, vars(age, sex, bmi, smoke), data = d, event = "Yes",
#'   method = "logistic", tune = FALSE, boot = 200,
#'   plot = TRUE,
#'   plot_display = c("roc", "calibration", "confusion", "importance")
#' )
#' names(mv$plots)
#' mv$settings$viewer_plots
#' mv$settings$display_plots
#'
#' # Keep all figures in Viewer but draw none in the Plot pane.
#' mv2 <- tabmachine(
#'   hypertension, vars(age, sex, bmi, smoke), data = d, event = "Yes",
#'   method = "logistic", tune = FALSE, boot = 200,
#'   plot = TRUE, plot_display = "none"
#' )
#' }
#' ```
#'
#' @section Examples - classification thresholds:
#' ```
#' \donttest{
#' myouden <- tabmachine(hypertension, vars(age, sex, bmi, smoke), data=d,
#'   event="Yes", method="logistic", tune=FALSE, threshold="youden",
#'   boot=200, show=FALSE, plot=FALSE)
#' mf1 <- tabmachine(hypertension, vars(age, sex, bmi, smoke), data=d,
#'   event="Yes", method="logistic", tune=FALSE, threshold="f1",
#'   boot=200, show=FALSE, plot=FALSE)
#' msens <- tabmachine(hypertension, vars(age, sex, bmi, smoke), data=d,
#'   event="Yes", method="logistic", tune=FALSE, threshold="sens",
#'   target_sens=.90, boot=200, show=FALSE, plot=FALSE)
#' mspec <- tabmachine(hypertension, vars(age, sex, bmi, smoke), data=d,
#'   event="Yes", method="logistic", tune=FALSE, threshold="spec",
#'   target_spec=.90, boot=200, show=FALSE, plot=FALSE)
#' mfixed <- tabmachine(hypertension, vars(age, sex, bmi, smoke), data=d,
#'   event="Yes", method="logistic", tune=FALSE, threshold=.20,
#'   boot=200, show=FALSE, plot=FALSE)
#' myouden$threshold
#' plot(myouden, "threshold")
#' }
#' ```
#'
#' @section Examples - missing data, correlation, outliers, and transformations:
#' ```
#' \donttest{
#' mp <- tabmachine(
#'   hypertension, vars(age, sex, bmi, smoke), data=d, event="Yes",
#'   method="logistic", tune=FALSE,
#'   missing="auto", corr=.90, outlier="winsor", transform="none",
#'   boot=200, show=FALSE, plot=FALSE
#' )
#' mp$preprocessing
#' }
#' ```
#'
#' @section Examples - imbalance without extra packages:
#' ```
#' \donttest{
#' # weight, up, down, SMOTE, and ADASYN are implemented inside R4VN.
#' mw <- tabmachine(hypertension, vars(age, sex, bmi, smoke), data=d,
#'   event="Yes", method="logistic", tune=FALSE, balance="weight",
#'   metric="pr_auc", boot=200, show=FALSE, plot=FALSE)
#' msmote <- tabmachine(hypertension, vars(age, sex, bmi, smoke), data=d,
#'   event="Yes", method="logistic", tune=FALSE, balance="smote",
#'   metric="pr_auc", boot=200, show=FALSE, plot=FALSE)
#' mw$balance
#' msmote$balance
#' }
#' ```
#'
#' @section Examples - feature engineering and PCA:
#' ```
#' \donttest{
#' mfeat <- tabmachine(hypertension, vars(age, bmi), data=d, event="Yes",
#'   method="logistic", tune=FALSE, feature="all", degree=2,
#'   boot=200, show=FALSE, plot=FALSE)
#' mfeat$selection
#'
#' mpdp <- tabmachine(hypertension, vars(age, bmi, smoke), data=d, event="Yes",
#'   method="logistic", tune=FALSE, pdp=c("age","bmi"), boot=200,
#'   plot=TRUE, plot_display="pdp")
#' mpdp$pdp
#' plot(mpdp, "pdp")
#'
#' set.seed(11)
#' hd <- as.data.frame(matrix(rnorm(180*60),180,60))
#' names(hd) <- paste0("x",1:60)
#' hd$y <- factor(rbinom(180,1,plogis(hd$x1-.7*hd$x2+.5*hd$x3)),0:1,c("No","Yes"))
#' mpca <- tabmachine(y, x=., data=hd, event="Yes", method="logistic",
#'   tune=FALSE, reduce="pca", variance=.90, select="filter",
#'   simplify=FALSE, boot=100, show=FALSE, plot=FALSE)
#' mpca$selection
#' }
#' ```
#'
#' @section Examples - feature selection and parsimony:
#' ```
#' \donttest{
#' mfilter <- tabmachine(hypertension, vars(age, sex, bmi, smoke), data=d,
#'   event="Yes", method="logistic", tune=FALSE, select="filter",
#'   simplify=TRUE, boot=200, show=FALSE, plot=FALSE)
#' mfilter$selection
#' mfilter$simplify
#'
#' # Penalized selection is optional and used only when glmnet is installed.
#' if (requireNamespace("glmnet", quietly=TRUE)) {
#'   mlasso <- tabmachine(hypertension, vars(age, sex, bmi, smoke), data=d,
#'     event="Yes", method="lasso", select="lasso", boot=200,
#'     show=FALSE, plot=FALSE)
#'   mlasso$selection
#' }
#' }
#' ```
#'
#' @section Examples - repeated and nested validation:
#' ```
#' \donttest{
#' mrep <- tabmachine(hypertension, vars(age, sex, bmi, smoke), data=d,
#'   event="Yes", method=c("logistic","tree"), tune=TRUE,
#'   folds=5, repeats=3, boot=200, show=FALSE, plot=FALSE)
#' mrep$cv_performance
#'
#' mnested <- tabmachine(hypertension, vars(age, sex, bmi, smoke), data=d,
#'   event="Yes", method=c("logistic","tree"), tune=TRUE,
#'   folds=5, nested=TRUE, boot=200, show=FALSE, plot=FALSE)
#' mnested$cv_performance
#' }
#' ```
#'
#' @section Examples - three-way split:
#' ```
#' \donttest{
#' msplit <- tabmachine(hypertension, vars(age, sex, bmi, smoke), data=d,
#'   event="Yes", split=c(.70,.15,.15), method="auto", boot=200,
#'   show=FALSE, plot=FALSE)
#' msplit$overview
#' }
#' ```
#'
#' @section Examples - calibration, clinical utility, and learning curve:
#' ```
#' \donttest{
#' mclin <- tabmachine(hypertension, vars(age, sex, bmi, smoke), data=d,
#'   event="Yes", method="logistic", tune=FALSE,
#'   calibration=TRUE, decision=TRUE, learning=TRUE, boot=200,
#'   plot=TRUE, plot_display="all")
#' mclin$calibration$statistics
#' head(mclin$decision)
#' mclin$learning
#' }
#' ```
#'
#' @section Examples - regression diagnostics:
#' ```
#' \donttest{
#' mr2 <- tabmachine(sbp, vars(age,bmi,sex), data=r, method="auto",
#'   boot=200, plot=TRUE, plot_display=c("observed","residual","importance"))
#' mr2$comparison
#' plot(mr2,"observed")
#' plot(mr2,"residual")
#' }
#' ```
#'
#' @section Examples - multiclass classification:
#' ```
#' \donttest{
#' ir <- iris
#' mm <- tabmachine(Species, vars(Sepal.Length,Sepal.Width,Petal.Length,Petal.Width),
#'   data=ir, method="auto", folds=5, boot=200,
#'   plot=TRUE, plot_display=c("confusion","importance"))
#' mm$performance
#' mm$class_performance
#' mm$confusion
#' plot(mm,"confusion")
#' }
#' ```
#'
#' @section Examples - predictions for new patients:
#' ```
#' \donttest{
#' newpatients <- data.frame(
#'   patient_id=c("P001","P002"), age=c(35,68),
#'   sex=factor(c("Female","Male"),levels=levels(d$sex)),
#'   bmi=c(21,31), smoke=factor(c("No","Yes"),levels=levels(d$smoke))
#' )
#' mpred <- tabmachine(hypertension, vars(age,sex,bmi,smoke), data=d,
#'   event="Yes", method="logistic", tune=FALSE, boot=200,
#'   predict=newpatients, id=patient_id, show=FALSE, plot=FALSE)
#' mpred$predictions
#' predict(mpred,newpatients,type="prob")
#' predict(mpred,newpatients,type="class")
#' }
#' ```
#'
#' @section Examples - optional advanced engines:
#' ```
#' \donttest{
#' # None of these packages is needed for the ordinary R4VN auto workflow.
#' if (requireNamespace("ranger",quietly=TRUE)) {
#'   mrf <- tabmachine(hypertension, vars(age,sex,bmi,smoke), data=d,
#'     event="Yes", method="rf", boot=200, show=FALSE, plot=FALSE)
#' }
#' if (requireNamespace("xgboost",quietly=TRUE)) {
#'   mxgb <- tabmachine(hypertension, vars(age,sex,bmi,smoke), data=d,
#'     event="Yes", method="xgb", shap=TRUE, boot=200, show=FALSE, plot=FALSE)
#'   mxgb$shap
#'   plot(mxgb, "shap")
#' }
#' if (requireNamespace("e1071",quietly=TRUE)) {
#'   msvm <- tabmachine(hypertension, vars(age,sex,bmi,smoke), data=d,
#'     event="Yes", method="svm", boot=200, show=FALSE, plot=FALSE)
#' }
#' if (requireNamespace("Boruta",quietly=TRUE)) {
#'   mb <- tabmachine(hypertension, vars(age,sex,bmi,smoke), data=d,
#'     event="Yes", method="logistic", select="boruta", tune=FALSE,
#'     boot=200, show=FALSE, plot=FALSE)
#' }
#' if (requireNamespace("ROSE",quietly=TRUE)) {
#'   mrose <- tabmachine(hypertension, vars(age,sex,bmi,smoke), data=d,
#'     event="Yes", method="logistic", balance="rose", tune=FALSE,
#'     boot=200, show=FALSE, plot=FALSE)
#' }
#' }
#' ```
#'
#' @examples
#' set.seed(99)
#' n <- 110
#' dd <- data.frame(
#'   x1 = rnorm(n),
#'   x2 = rnorm(n),
#'   group = factor(sample(c("A", "B"), n, TRUE))
#' )
#' pp <- plogis(-0.4 + 0.9 * dd$x1 - 0.5 * dd$x2)
#' dd$y <- factor(rbinom(n, 1, pp), 0:1, c("No", "Yes"))
#'
#' z <- tabmachine(
#'   y,
#'   x = vars(x1, x2, group),
#'   data = dd,
#'   event = "Yes",
#'   method = "logistic",
#'   tune = FALSE,
#'   folds = 2,
#'   ci = FALSE,
#'   boot = 50,
#'   calibration = FALSE,
#'   decision = FALSE,
#'   importance = FALSE,
#'   explain = FALSE,
#'   show = FALSE,
#'   plot = FALSE
#' )
#' z$comparison
#' @export
tabmachine <- function(
  outcome,
  x = NULL,
  data = NULL,
  exclude = NULL,
  task = c("auto", "binary", "multiclass", "regression"),
  event = NULL,
  preprocess = c("auto", "none"),
  missing = c("auto", "median", "mean", "mode", "complete"),
  missing_max = 0.50,
  encode = c("auto", "dummy"),
  standardize = c("auto", "none", "z", "minmax", "robust"),
  transform = c("none", "auto", "log", "yeojohnson"),
  outlier = c("none", "detect", "winsor", "robust"),
  corr = "auto",
  feature = c("none", "auto", "polynomial", "interaction", "all"),
  degree = 2,
  reduce = c("none", "auto", "pca"),
  variance = 0.95,
  select = c("auto", "none", "filter", "lasso", "stepwise",
             "importance", "rfe", "boruta", "compare"),
  nfeatures = "auto",
  simplify = TRUE,
  simplify_tol = 0.01,
  balance = c("auto", "none", "weight", "up", "down", "smote",
              "adasyn", "rose", "compare"),
  balance_target = 0.50,
  neighbors = 5,
  split = 0.80,
  folds = 10,
  repeats = 1,
  nested = FALSE,
  seed = NULL,
  method = "auto",
  tune = TRUE,
  tune_n = 10,
  metric = "auto",
  threshold = "auto",
  target_sens = 0.90,
  target_spec = 0.90,
  ci = TRUE,
  ci_level = 0.95,
  boot = 1000,
  calibration = TRUE,
  decision = TRUE,
  decision_thresholds = seq(.01, .99, .01),
  learning = FALSE,
  importance = TRUE,
  importance_repeats = 20,
  explain = TRUE,
  shap = FALSE,
  pdp = FALSE,
  validation = NULL,
  predict = NULL,
  id = NULL,
  show = TRUE,
  plot = TRUE,
  plot_display = "auto",
  plot_args = list(),
  strict = FALSE,
  console = FALSE,
  digit = 3,
  title = NULL,
  ai = FALSE,
  ...
) {
  seed <- .r4vn_seed_normalize(seed)
  task <- match.arg(task)
  preprocess <- match.arg(preprocess)
  missing <- match.arg(missing)
  encode <- match.arg(encode)
  standardize <- match.arg(standardize)
  transform <- match.arg(transform)
  outlier <- match.arg(outlier)
  feature <- match.arg(feature)
  reduce <- match.arg(reduce)
  select <- match.arg(select)
  balance <- match.arg(balance)

  .r4vn_machine_validate_plot_options(plot, plot_display, plot_args, strict)

  .r4vn_machine_validate_scalar(missing_max, "missing_max", 0, 1, include_low = TRUE, include_high = FALSE)
  .r4vn_machine_validate_scalar(variance, "variance", 0, 1, include_low = FALSE, include_high = TRUE)
  if (!is.numeric(degree) || length(degree) != 1L || is.na(degree) || !degree %in% c(2, 3)) stop("`degree` must be 2 or 3.", call. = FALSE)
  .r4vn_machine_validate_scalar(balance_target, "balance_target", 0, 1, include_low = FALSE, include_high = FALSE)
  .r4vn_machine_validate_scalar(simplify_tol, "simplify_tol", 0, Inf, include_low = TRUE, include_high = TRUE)
  .r4vn_machine_validate_scalar(ci_level, "ci_level", 0, 1, include_low = FALSE, include_high = FALSE)
  .r4vn_machine_validate_scalar(target_sens, "target_sens", 0, 1, include_low = FALSE, include_high = TRUE)
  .r4vn_machine_validate_scalar(target_spec, "target_spec", 0, 1, include_low = FALSE, include_high = TRUE)
  if (!is.numeric(folds) || length(folds) != 1L || folds < 2) stop("`folds` must be at least 2.", call. = FALSE)
  if (!is.numeric(repeats) || length(repeats) != 1L || repeats < 1) stop("`repeats` must be at least 1.", call. = FALSE)
  if (!is.numeric(boot) || length(boot) != 1L || boot < 50) stop("`boot` must be at least 50.", call. = FALSE)
  if (!is.numeric(neighbors) || length(neighbors) != 1L || neighbors < 1) stop("`neighbors` must be at least 1.", call. = FALSE)
  if (!is.numeric(importance_repeats) || length(importance_repeats) != 1L || importance_repeats < 1) stop("`importance_repeats` must be at least 1.", call. = FALSE)
  if (!is.numeric(tune_n) || length(tune_n) != 1L || is.na(tune_n) || tune_n < 1) stop("`tune_n` must be at least 1.", call. = FALSE)
  if (!(is.character(nfeatures) && length(nfeatures) == 1L && identical(tolower(nfeatures), "auto"))) {
    if (!is.numeric(nfeatures) || length(nfeatures) != 1L || is.na(nfeatures) || nfeatures < 1) stop("`nfeatures` must be `\"auto\"` or a positive number.", call. = FALSE)
  }

  dat <- .r4vn_machine_resolve_data(data)
  outcome_expr <- substitute(outcome)
  outcome_name <- .r4vn_machine_one_name(outcome_expr, names(dat), "outcome")
  if (!outcome_name %in% names(dat)) stop("Outcome not found in `data`: ", outcome_name, call. = FALSE)

  x_expr <- if (missing(x)) NULL else substitute(x)
  ex_expr <- if (missing(exclude)) NULL else substitute(exclude)
  # Evaluate selector objects only when the expression is not itself a data
  # column. Bare column names have priority over same-named caller objects.
  x_value <- if (
    is.null(x_expr) ||
    (is.symbol(x_expr) && as.character(x_expr) %in% c(".", names(dat)))
  ) {
    NULL
  } else {
    tryCatch(eval(x_expr, parent.frame()), error = function(e) NULL)
  }

  ex_value <- if (
    is.null(ex_expr) ||
    (is.symbol(ex_expr) && as.character(ex_expr) %in% names(dat))
  ) {
    NULL
  } else {
    tryCatch(eval(ex_expr, parent.frame()), error = function(e) NULL)
  }
  predictor_names <- .r4vn_machine_predictor_names(x_expr, x_value, dat, outcome_name)
  exclude_names <- .r4vn_machine_names(ex_expr, ex_value, names(dat), allow_dot = FALSE)
  predictor_names <- setdiff(unique(predictor_names), c(outcome_name, exclude_names))
  if (!length(predictor_names)) stop("No candidate predictors remain after exclusions.", call. = FALSE)
  bad_names <- setdiff(predictor_names, names(dat))
  if (length(bad_names)) stop("Predictors not found in `data`: ", paste(bad_names, collapse = ", "), call. = FALSE)

  labels <- .r4vn_machine_labels(dat, c(outcome_name, predictor_names))
  outcome_info <- .r4vn_machine_outcome(dat[[outcome_name]], task = task, event = event, name = outcome_name)
  task <- outcome_info$task
  dat[[outcome_name]] <- outcome_info$y
  dat <- dat[!is.na(dat[[outcome_name]]), , drop = FALSE]
  if (nrow(dat) < 20L) stop("At least 20 observations with non-missing outcome are required.", call. = FALSE)

  notes <- outcome_info$notes
  warnings0 <- .r4vn_machine_data_warnings(dat, predictor_names, outcome_name)
  notes <- c(notes, warnings0)

  split_obj <- .r4vn_machine_split(dat, outcome_name, task, split = split, seed = seed)
  train_raw <- split_obj$train
  test_raw <- split_obj$test
  internal_raw <- split_obj$internal
  if(!nrow(test_raw)) {
    notes_split <- "No held-out internal test set was created; development performance is apparent/resubstitution performance. Prefer external validation or a held-out test for final assessment."
  } else notes_split <- character()

  # Resolve task-sensitive defaults only from the development training sample.
  notes <- c(notes, notes_split)
  event_rate <- if (task == "binary") mean(train_raw[[outcome_name]] == 1L) else NA_real_
  primary_metric <- .r4vn_machine_metric(metric, task, event_rate)
  select_resolved <- .r4vn_machine_resolve_select(select, train_raw, predictor_names, task)
  balance_resolved <- .r4vn_machine_resolve_balance(balance, task, event_rate)
  methods <- .r4vn_machine_methods(method, task)
  engine_table <- .r4vn_machine_engine_table(method, task, methods)
  if (!length(methods)) {
    unavailable <- engine_table[engine_table$Requested == "Yes" & engine_table$Available == "No", , drop = FALSE]
    detail <- if (nrow(unavailable)) paste0(
      " Missing optional package(s): ",
      paste(unique(unavailable$Package), collapse = ", "),
      ". Use `method = 'auto'` for the low-dependency R4VN core workflow."
    ) else ""
    stop("No requested machine-learning engine is available.", detail, call. = FALSE)
  }
  skipped_engines <- engine_table$Method[engine_table$Requested == "Yes" & engine_table$Available == "No"]
  if (length(skipped_engines)) {
    notes <- c(notes, paste0(
      "Requested optional engine(s) unavailable and therefore skipped: ",
      paste(skipped_engines, collapse = ", "), "."
    ))
  }

  # Strategy comparison is performed only on development training data.
  selection_compare <- data.frame()
  if (identical(select, "compare")) {
    sc <- .r4vn_machine_compare_selection(
      train_raw, outcome_name, predictor_names, task, outcome_info,
      methods = methods, preprocess = preprocess, missing = missing,
      missing_max = missing_max, encode = encode, standardize = standardize,
      transform = transform, outlier = outlier, corr = corr,
      feature = feature, degree = degree, reduce = reduce, variance = variance,
      nfeatures = nfeatures, balance = balance_resolved,
      balance_target = balance_target, neighbors = neighbors,
      folds = min(5L, as.integer(folds)), metric = primary_metric,
      seed = seed + 101L
    )
    select_resolved <- sc$best
    selection_compare <- sc$table
    notes <- c(notes, sc$notes)
  }

  balance_compare <- data.frame()
  if (identical(balance, "compare") && task == "binary") {
    bc <- .r4vn_machine_compare_balance(
      train_raw, outcome_name, predictor_names, task, outcome_info,
      methods = methods, preprocess = preprocess, missing = missing,
      missing_max = missing_max, encode = encode, standardize = standardize,
      transform = transform, outlier = outlier, corr = corr,
      feature = feature, degree = degree, reduce = reduce, variance = variance,
      select = select_resolved, nfeatures = nfeatures,
      balance_target = balance_target, neighbors = neighbors,
      folds = min(5L, as.integer(folds)), metric = primary_metric,
      seed = seed + 202L
    )
    balance_resolved <- bc$best
    balance_compare <- bc$table
    notes <- c(notes, bc$notes)
  }

  # Fit/tune each algorithm using development data only.
  model_results <- list()
  tuning_rows <- list()
  cv_rows <- list()
  engine_notes <- character()

  for (i in seq_along(methods)) {
    meth <- methods[i]
    grid <- .r4vn_machine_tune_grid(meth, task, tune = tune, tune_n = tune_n, seed = seed + i)

    tuned <- .r4vn_machine_tune(
      method = meth,
      grid = grid,
      raw = train_raw,
      outcome_name = outcome_name,
      predictor_names = predictor_names,
      task = task,
      outcome_info = outcome_info,
      preprocess = preprocess,
      missing = missing,
      missing_max = missing_max,
      encode = encode,
      standardize = standardize,
      transform = transform,
      outlier = outlier,
      corr = corr,
      feature = feature,
      degree = degree,
      reduce = reduce,
      variance = variance,
      select = select_resolved,
      nfeatures = nfeatures,
      balance = balance_resolved,
      balance_target = balance_target,
      neighbors = neighbors,
      folds = as.integer(folds),
      repeats = as.integer(repeats),
      nested = nested,
      metric = primary_metric,
      seed = seed + 1000L * i
    )

    if (!isTRUE(tuned$ok)) {
      engine_notes <- c(engine_notes, paste0(meth, ": ", tuned$error))
      next
    }
    model_results[[meth]] <- tuned
    if (nrow(tuned$tuning)) tuning_rows[[meth]] <- tuned$tuning
    cv_rows[[meth]] <- tuned$cv_summary
  }

  if (!length(model_results)) {
    stop("All requested models failed. Details: ", paste(engine_notes, collapse = " | "), call. = FALSE)
  }
  notes <- c(notes, engine_notes)
  tuning_table <- .r4vn_machine_bind(tuning_rows)
  cv_performance <- .r4vn_machine_bind(cv_rows)

  best_method <- .r4vn_machine_choose_best(cv_performance, primary_metric, task)
  if (!best_method %in% names(model_results)) best_method <- names(model_results)[1L]

  # OOF predictions for training-derived threshold use only the selected model.
  best_cv <- model_results[[best_method]]$cv
  threshold_info <- NULL
  threshold_value <- NA_real_
  if (task == "binary") {
    threshold_cv <- .r4vn_machine_oof_unique_binary(best_cv)
    threshold_info <- .r4vn_machine_threshold(
      y = threshold_cv$truth,
      p = threshold_cv$prob,
      threshold = threshold,
      target_sens = target_sens,
      target_spec = target_spec,
      ci = ci,
      level = ci_level,
      boot = as.integer(boot),
      seed = seed + 505L
    )
    threshold_value <- threshold_info$threshold
  }

  # Optional parsimony search is training-only and can replace the selected
  # encoded feature set used by the final fit. It never inspects test outcomes.
  simplify_info <- NULL
  if (isTRUE(simplify)) {
    simplify_info <- .r4vn_machine_simplify(
      train_raw, outcome_name, predictor_names, task, outcome_info,
      method = best_method,
      params = model_results[[best_method]]$params,
      preprocess = preprocess, missing = missing, missing_max = missing_max,
      encode = encode, standardize = standardize, transform = transform,
      outlier = outlier, corr = corr, feature = feature, degree = degree, reduce = reduce, variance = variance,
      select = select_resolved,
      nfeatures = nfeatures, balance = balance_resolved,
      balance_target = balance_target, neighbors = neighbors,
      metric = primary_metric, folds = min(5L, as.integer(folds)),
      tolerance = simplify_tol, seed = seed + 303L
    )
  }

  # Recombine the optional internal-validation partition only after all model
  # decisions. The held-out test remains untouched.
  final_train_raw <- if (nrow(internal_raw)) rbind(train_raw, internal_raw) else train_raw
  final_selected_override <- if (!is.null(simplify_info) && length(simplify_info$selected_encoded)) simplify_info$selected_encoded else NULL

  # Fit each candidate on full development training for fair held-out comparison.
  fitted_models <- list()
  perf_rows <- list()
  test_predictions <- list()
  for (meth in names(model_results)) {
    ff <- .r4vn_machine_fit_pipeline(
      raw = final_train_raw,
      outcome_name = outcome_name,
      predictor_names = predictor_names,
      task = task,
      outcome_info = outcome_info,
      method = meth,
      params = model_results[[meth]]$params,
      preprocess = preprocess,
      missing = missing,
      missing_max = missing_max,
      encode = encode,
      standardize = standardize,
      transform = transform,
      outlier = outlier,
      corr = corr,
      feature = feature,
      degree = degree,
      reduce = reduce,
      variance = variance,
      select = select_resolved,
      nfeatures = nfeatures,
      balance = balance_resolved,
      balance_target = balance_target,
      neighbors = neighbors,
      selected_override = if (identical(meth, best_method)) final_selected_override else NULL,
      seed = seed + match(meth, names(model_results)) * 7000L
    )
    fitted_models[[meth]] <- ff

    eval_raw <- if (nrow(test_raw)) test_raw else train_raw
    pr <- .r4vn_machine_predict_pipeline(ff, eval_raw, outcome_name = outcome_name)
    test_predictions[[meth]] <- pr
    ev <- .r4vn_machine_evaluate(
      truth = eval_raw[[outcome_name]],
      pred = pr,
      task = task,
      method = meth,
      threshold = if (task == "binary") threshold_value else NULL,
      ci = ci,
      ci_level = ci_level,
      boot = as.integer(boot),
      seed = seed + 9000L + match(meth, names(model_results)),
      dataset = if (nrow(test_raw)) "Test" else "Development (apparent)"
    )
    perf_rows[[meth]] <- ev$performance
  }

  performance <- .r4vn_machine_bind(perf_rows)
  comparison <- .r4vn_machine_comparison(performance, primary_metric, digit = digit)
  if(nrow(comparison)) comparison$Selected <- ifelse(comparison$Model==best_method,"Yes","")
  final_fit <- fitted_models[[best_method]]
  final_pred <- test_predictions[[best_method]]
  eval_raw <- if (nrow(test_raw)) test_raw else train_raw

  model_difference <- .r4vn_machine_model_differences(
    truth = eval_raw[[outcome_name]],
    predictions = test_predictions,
    best = best_method,
    task = task,
    metric = primary_metric,
    threshold = threshold_value,
    ci = ci,
    level = ci_level,
    boot = as.integer(boot),
    seed = seed + 12000L
  )

  overfitting <- .r4vn_machine_overfitting(cv_performance, performance, primary_metric)
  coefficients <- .r4vn_machine_coefficients(final_fit, labels = labels, ci = ci, ci_level = ci_level)

  # Final clinical evaluation uses exactly the rows returned by the finalized
  # prediction pipeline. This matters for policies such as missing = "complete",
  # where prediction rows can be a strict subset of the evaluation data.
  final_rows <- final_pred$rows %||% seq_len(nrow(eval_raw))
  final_truth <- eval_raw[[outcome_name]][final_rows]
  confusion <- data.frame()
  calibration_out <- NULL
  decision_out <- data.frame()
  if (task == "binary") {
    bm <- .r4vn_machine_binary_metrics(final_truth, final_pred$prob, threshold_value)
    confusion <- data.frame(
      Actual = c(outcome_info$event, outcome_info$negative),
      `Predicted event` = c(unname(bm["tp"]), unname(bm["fp"])),
      `Predicted non-event` = c(unname(bm["fn"]), unname(bm["tn"])),
      check.names = FALSE, stringsAsFactors = FALSE
    )
    if (isTRUE(calibration)) {
      calibration_out <- .r4vn_machine_calibration(
        final_truth, final_pred$prob,
        ci = ci, level = ci_level, boot = as.integer(boot), seed = seed + 13000L
      )
    }
    if (isTRUE(decision)) {
      decision_out <- .r4vn_machine_decision(
        final_truth, final_pred$prob,
        thresholds = decision_thresholds,
        ci = ci, level = ci_level, boot = as.integer(boot), seed = seed + 14000L
      )
    }
  }

  class_performance <- data.frame()
  if(task=="multiclass") {
    class_performance <- .r4vn_machine_multiclass_class_performance(
      final_truth, final_pred, ci=ci, level=ci_level, boot=as.integer(boot), seed=seed+14500L
    )
    confusion <- .r4vn_machine_multiclass_confusion(
      final_truth, final_pred$class, outcome_info$levels
    )
  }

  importance_out <- data.frame()
  if (isTRUE(importance)) {
    importance_out <- .r4vn_machine_importance(
      final_fit, eval_raw, outcome_name,
      metric = primary_metric,
      threshold = threshold_value,
      repeats = as.integer(importance_repeats),
      seed = seed + 15000L,
      labels = labels,
      ci = ci, level = ci_level
    )
  }

  shap_out <- NULL
  if (isTRUE(shap)) shap_out <- .r4vn_machine_shap(final_fit, eval_raw, outcome_name)
  pdp_out <- NULL
  if (!identical(pdp, FALSE)) {
    pdp_out <- .r4vn_machine_pdp(final_fit, eval_raw, outcome_name, pdp, importance_out, labels)
  }

  learning_out <- data.frame()
  if (isTRUE(learning)) {
    learning_out <- .r4vn_machine_learning_curve(
      final_train_raw, outcome_name, predictor_names, task, outcome_info,
      method = best_method, params = model_results[[best_method]]$params,
      preprocess = preprocess, missing = missing, missing_max = missing_max,
      encode = encode, standardize = standardize, transform = transform,
      outlier = outlier, corr = corr, feature = feature, degree = degree, reduce = reduce, variance = variance,
      select = select_resolved,
      nfeatures = nfeatures, balance = balance_resolved,
      balance_target = balance_target, neighbors = neighbors,
      metric = primary_metric, threshold = threshold_value,
      ci = ci, ci_level = ci_level, seed = seed + 16000L
    )
  }

  external_performance <- data.frame()
  external_predictions <- NULL
  external_confusion <- data.frame()
  external_class_performance <- data.frame()
  external_calibration <- NULL
  external_decision <- data.frame()
  if (!is.null(validation)) {
    if (!is.data.frame(validation)) stop("`validation` must be a data frame.", call. = FALSE)
    needed <- unique(c(outcome_name, final_fit$recipe$keep))
    missv <- setdiff(needed, names(validation))
    if (length(missv)) stop("External validation data are missing: ", paste(missv, collapse = ", "), call. = FALSE)
    val <- validation
    val[[outcome_name]] <- .r4vn_machine_apply_outcome(val[[outcome_name]], outcome_info)
    val <- val[!is.na(val[[outcome_name]]), , drop = FALSE]
    external_predictions <- .r4vn_machine_predict_pipeline(final_fit, val, outcome_name = outcome_name)
    external_rows <- external_predictions$rows %||% seq_len(nrow(val))
    external_truth <- val[[outcome_name]][external_rows]
    external_performance <- .r4vn_machine_evaluate(
      val[[outcome_name]], external_predictions, task, best_method,
      threshold = if (task == "binary") threshold_value else NULL,
      ci = ci, ci_level = ci_level, boot = as.integer(boot),
      seed = seed + 17000L, dataset = "External validation"
    )$performance
    if (task == "binary") {
      ebm <- .r4vn_machine_binary_metrics(external_truth, external_predictions$prob, threshold_value)
      external_confusion <- data.frame(
        Actual = c(outcome_info$event, outcome_info$negative),
        `Predicted event` = c(unname(ebm["tp"]), unname(ebm["fp"])),
        `Predicted non-event` = c(unname(ebm["fn"]), unname(ebm["tn"])),
        check.names = FALSE, stringsAsFactors = FALSE
      )
      if (isTRUE(calibration)) external_calibration <- .r4vn_machine_calibration(
        external_truth, external_predictions$prob,
        ci = ci, level = ci_level, boot = as.integer(boot), seed = seed + 17100L
      )
      if (isTRUE(decision)) external_decision <- .r4vn_machine_decision(
        external_truth, external_predictions$prob,
        thresholds = decision_thresholds, ci = ci, level = ci_level,
        boot = as.integer(boot), seed = seed + 17200L
      )
    } else if (task == "multiclass") {
      external_class_performance <- .r4vn_machine_multiclass_class_performance(
        external_truth, external_predictions, ci = ci, level = ci_level,
        boot = as.integer(boot), seed = seed + 17300L
      )
      external_confusion <- .r4vn_machine_multiclass_confusion(
        external_truth, external_predictions$class, outcome_info$levels
      )
    }
  }

  prediction_output <- NULL
  if (!is.null(predict)) {
    if (!is.data.frame(predict)) stop("`predict` must be a data frame.", call. = FALSE)
    id_expr <- if (missing(id)) NULL else substitute(id)
    id_name <- if (is.null(id_expr)) NULL else .r4vn_machine_one_name(id_expr, names(predict), "id")
    prediction_output <- .r4vn_machine_prediction_frame(
      final_fit, predict, threshold = threshold_value,
      id_name = id_name, outcome_info = outcome_info
    )
  }

  preprocessing_summary <- .r4vn_machine_preprocessing_table(final_fit)
  selection_table <- .r4vn_machine_selection_table(final_fit, selection_compare, labels)
  balance_table <- .r4vn_machine_balance_table(balance_resolved, final_train_raw[[outcome_name]], balance_compare, task)
  overview <- .r4vn_machine_overview(
    dat, outcome_name, predictor_names, task, outcome_info,
    train_raw, internal_raw, test_raw, validation, primary_metric
  )

  settings <- list(
    outcome = outcome_name,
    outcome_label = labels[[outcome_name]],
    predictors = predictor_names,
    excluded = exclude_names,
    task = task,
    event = outcome_info$event,
    preprocess = preprocess,
    missing = missing,
    missing_max = missing_max,
    encode = encode,
    standardize = standardize,
    transform = transform,
    outlier = outlier,
    corr = corr,
    feature = feature,
    degree = degree,
    reduce = reduce,
    variance = variance,
    select_requested = select,
    select = select_resolved,
    nfeatures = nfeatures,
    simplify = simplify,
    simplify_tol = simplify_tol,
    balance_requested = balance,
    balance = balance_resolved,
    balance_target = balance_target,
    neighbors = neighbors,
    split = split,
    folds = folds,
    repeats = repeats,
    nested = nested,
    seed = seed,
    methods = names(model_results),
    tune = tune,
    metric = primary_metric,
    threshold = threshold,
    threshold_value = threshold_value,
    ci = ci,
    ci_level = ci_level,
    boot = boot,
    calibration = calibration,
    decision = decision,
    learning = learning,
    importance = importance,
    explain = explain,
    shap = shap,
    pdp = pdp,
    plot = plot,
    plot_display = plot_display,
    plot_args = plot_args,
    strict = strict,
    digit = digit,
    title = title
  )

  out <- list(
    overview = overview,
    engines = engine_table,
    preprocessing = preprocessing_summary,
    balance = balance_table,
    selection = selection_table,
    simplify = simplify_info,
    tuning = tuning_table,
    cv_performance = cv_performance,
    performance = performance,
    comparison = comparison,
    model_difference = model_difference,
    overfitting = overfitting,
    coefficients = coefficients,
    best = best_method,
    final = final_fit,
    models = fitted_models,
    threshold = threshold_info,
    confusion = confusion,
    class_performance = class_performance,
    calibration = calibration_out,
    decision = decision_out,
    importance = importance_out,
    explanation = if (isTRUE(explain)) list(importance = importance_out, shap = shap_out, pdp = pdp_out) else NULL,
    shap = shap_out,
    pdp = pdp_out,
    learning = learning_out,
    external_performance = external_performance,
    external_predictions = external_predictions,
    external_confusion = external_confusion,
    external_class_performance = external_class_performance,
    external_calibration = external_calibration,
    external_decision = external_decision,
    predictions = prediction_output,
    plots = list(),
    plot_titles = .r4vn_machine_plot_title_map(),
    tables = list(),
    table = comparison,
    html = NULL,
    settings = settings,
    notes = unique(notes[nzchar(notes)]),
    call = match.call()
  )
  class(out) <- c("r4vn_machine", "r4vn_tab")

  out$plots <- .r4vn_machine_plot_data(out, eval_raw[[outcome_name]], test_predictions)
  out$tables <- .r4vn_machine_tables(out)
  plot_enabled <- !identical(plot, FALSE) && !is.null(plot) &&
    !(is.character(plot) && "none" %in% tolower(plot))
  viewer_types <- .r4vn_machine_plot_types(out, plot)
  display_types <- if (plot_enabled) .r4vn_machine_display_types(out, plot_display) else character()
  if (plot_enabled) {
    viewer_types <- unique(c(viewer_types, display_types))
  }
  out$settings$viewer_plots <- viewer_types
  out$settings$display_plots <- display_types
  out$html <- .r4vn_machine_html(out, plot = viewer_types, plot_args = plot_args, title = title, digit = digit)

  if (!identical(ai, FALSE) && exists("aiask", mode = "function", inherits = TRUE)) {
    ai_fun <- get("aiask", mode = "function", inherits = TRUE)
    out$ai <- tryCatch({
      fml <- names(formals(ai_fun))
      if (isTRUE(ai) || !"api" %in% fml) ai_fun(out) else ai_fun(out, api = ai)
    }, error = function(e) paste0("AI interpretation unavailable: ", conditionMessage(e)))
  }

  if (interactive() && plot_enabled && length(display_types)) {
    .r4vn_machine_draw_display(out, display_types, plot_args = plot_args, strict = strict)
  }
  if (isTRUE(show)) .r4vn_machine_show_html(out$html)
  if (isTRUE(console)) print(out)
  invisible(out)
}

# ============================================================================
# Internal utilities: validation, NSE, labels, outcome, splitting
# ============================================================================

.r4vn_machine_validate_scalar <- function(x, name, low, high,
                                          include_low = FALSE,
                                          include_high = FALSE) {
  if (!is.numeric(x) || length(x) != 1L || is.na(x)) {
    stop("`", name, "` must be a single numeric value.", call. = FALSE)
  }
  low_ok <- if (include_low) x >= low else x > low
  high_ok <- if (include_high) x <= high else x < high
  if (!low_ok || !high_ok) {
    stop("`", name, "` is outside its allowed range.", call. = FALSE)
  }
  invisible(TRUE)
}

.r4vn_machine_resolve_data <- function(data) {
  if (!is.null(data)) {
    if (!is.data.frame(data)) stop("`data` must be a data frame.", call. = FALSE)
    return(data)
  }
  if (exists(".r4vn_get_active", mode = "function", inherits = TRUE)) {
    z <- get(".r4vn_get_active", mode = "function", inherits = TRUE)()
    if (is.data.frame(z)) return(z)
  }
  stop("Supply `data` or select active data with `usedf()`/`opendata(..., active = TRUE)`.", call. = FALSE)
}

.r4vn_machine_one_name <- function(expr, data_names, what = "variable") {
  if (is.character(expr) && length(expr) == 1L) return(expr)
  if (is.symbol(expr)) return(as.character(expr))
  if (is.call(expr) && identical(expr[[1L]], as.name("get"))) {
    val <- tryCatch(eval(expr, parent.frame()), error = function(e) NULL)
    if (is.character(val) && length(val) == 1L) return(val)
  }
  val <- tryCatch(eval(expr, parent.frame()), error = function(e) NULL)
  if (is.character(val) && length(val) == 1L && val %in% data_names) return(val)
  stop("Could not resolve `", what, "` to one variable name.", call. = FALSE)
}

.r4vn_machine_names <- function(expr, value, data_names, allow_dot = TRUE) {
  if (is.null(expr)) return(character())

  # A bare symbol that matches a data column always means that column.
  # This prevents caller objects with the same name from hijacking NSE.
  if (is.symbol(expr)) {
    nm <- as.character(expr)
    if (allow_dot && identical(nm, ".")) return(data_names)
    if (nm %in% data_names) return(nm)

    # A symbol that is not a data column may name a character selector
    # object in the caller environment.
    if (is.character(value)) {
      return(unique(as.character(value)))
    }

    return(nm)
  }

  # Literal character selector.
  if (is.character(expr)) {
    return(unique(as.character(expr)))
  }

  # Evaluated character selectors, e.g. x = c("age", "bmi").
  if (is.character(value)) {
    return(unique(as.character(value)))
  }

  walk <- function(e) {
    if (is.null(e)) return(character())

    if (is.symbol(e)) {
      nm <- as.character(e)
      if (allow_dot && identical(nm, ".")) return(data_names)
      return(nm)
    }

    if (is.character(e)) {
      return(as.character(e))
    }

    if (is.call(e)) {
      head <- as.character(e[[1L]])

      if (head %in% c("vars", "c", "list")) {
        return(unique(unlist(
          lapply(as.list(e)[-1L], walk),
          use.names = FALSE
        )))
      }
    }

    character()
  }

  unique(walk(expr))
}

.r4vn_machine_predictor_names <- function(expr, value, dat, outcome_name) {
  if (is.null(expr)) return(setdiff(names(dat), outcome_name))
  z <- .r4vn_machine_names(expr, value, names(dat), allow_dot = TRUE)
  if (!length(z)) return(setdiff(names(dat), outcome_name))
  setdiff(z, outcome_name)
}

.r4vn_machine_label <- function(x, fallback) {
  lab <- attr(x, "label", exact = TRUE)
  if (is.null(lab) || !length(lab) || is.na(lab[1L]) || !nzchar(as.character(lab[1L]))) fallback else as.character(lab[1L])
}

.r4vn_machine_labels <- function(dat, vars) {
  out <- stats::setNames(vars, vars)
  for (v in vars) out[[v]] <- .r4vn_machine_label(dat[[v]], v)
  out
}

.r4vn_machine_get_label <- function(labels, v) {
  if (!is.null(labels) && v %in% names(labels)) {
    z <- labels[[v]]
    if (!is.null(z) && length(z) && !is.na(z[1L]) && nzchar(as.character(z[1L]))) return(as.character(z[1L]))
  }
  v
}

.r4vn_machine_positive_guess <- function(values) {
  vals <- unique(as.character(values[!is.na(values)]))
  if (!length(vals)) return(NULL)
  norm <- tolower(trimws(vals))
  positive <- c("1", "true", "yes", "y", "positive", "pos", "case", "event", "co", "c\u00f3")
  hit <- match(positive, norm, nomatch = 0L)
  hit <- hit[hit > 0L]
  if (length(hit)) return(vals[hit[1L]])
  NULL
}

.r4vn_machine_outcome <- function(y, task = "auto", event = NULL, name = "outcome") {
  raw <- y
  obs <- y[!is.na(y)]
  uniq <- unique(obs)
  if (!length(uniq)) stop("Outcome `", name, "` has no observed values.", call. = FALSE)

  if (task == "auto") {
    if (length(uniq) == 2L) {
      task <- "binary"
    } else if (is.factor(y) || is.character(y) || is.logical(y)) {
      task <- "multiclass"
    } else if ((is.integer(y) || all(abs(obs - round(obs)) < .Machine$double.eps^0.5)) && length(uniq) <= 10L) {
      task <- "multiclass"
    } else if (is.numeric(y)) {
      task <- "regression"
    } else {
      stop("Could not infer outcome type. Set `task=` explicitly.", call. = FALSE)
    }
  }

  notes <- character()
  if (task == "binary") {
    vals <- unique(as.character(obs))
    if (length(vals) != 2L) stop("Binary outcome must have exactly two observed levels.", call. = FALSE)
    ev <- if (!is.null(event)) as.character(event)[1L] else .r4vn_machine_positive_guess(obs)
    if (is.null(ev) || !ev %in% vals) ev <- vals[2L]
    neg <- setdiff(vals, ev)[1L]
    yy <- rep(NA_integer_, length(y))
    yy[!is.na(y) & as.character(y) == neg] <- 0L
    yy[!is.na(y) & as.character(y) == ev] <- 1L
    notes <- c(notes, paste0("Binary event level: ", ev, "; non-event level: ", neg, "."))
    return(list(task = task, y = yy, event = ev, negative = neg,
                levels = c(neg, ev), original_class = class(raw), notes = notes))
  }

  if (task == "multiclass") {
    yy <- factor(y)
    if (nlevels(yy) < 3L) stop("Multiclass outcome must have at least three observed levels.", call. = FALSE)
    return(list(task = task, y = yy, event = NULL, negative = NULL,
                levels = levels(yy), original_class = class(raw), notes = notes))
  }

  if (!is.numeric(y)) stop("Regression outcome must be numeric.", call. = FALSE)
  yy <- as.numeric(y)
  if (length(unique(yy[!is.na(yy)])) < 5L) {
    notes <- c(notes, "Regression outcome has fewer than five distinct values; verify that regression is intended.")
  }
  list(task = "regression", y = yy, event = NULL, negative = NULL,
       levels = NULL, original_class = class(raw), notes = notes)
}

.r4vn_machine_apply_outcome <- function(y, info) {
  if (info$task == "binary") {
    z <- rep(NA_integer_, length(y))
    z[!is.na(y) & as.character(y) == info$negative] <- 0L
    z[!is.na(y) & as.character(y) == info$event] <- 1L
    return(z)
  }
  if (info$task == "multiclass") {
    return(factor(as.character(y), levels = info$levels))
  }
  suppressWarnings(as.numeric(y))
}

.r4vn_machine_data_warnings <- function(dat, predictors, outcome) {
  out <- character()
  for (v in predictors) {
    z <- dat[[v]]
    miss <- mean(is.na(z))
    if (miss > .5) out <- c(out, paste0(v, ": ", round(100 * miss, 1), "% missing."))
    nu <- length(unique(z[!is.na(z)]))
    if (nu <= 1L) out <- c(out, paste0(v, ": no variation."))
    if (grepl("(^id$|_id$|^id_|identifier|record.?id|patient.?id)", tolower(v)) &&
        sum(!is.na(z)) > 0 && length(unique(z[!is.na(z)])) / sum(!is.na(z)) > .95) {
      out <- c(out, paste0(v, ": looks like an identifier; consider `exclude = vars(", v, ")`."))
    }
  }
  unique(out)
}

.r4vn_machine_strata_indices <- function(y, prop, seed) {
  .r4vn_set_seed_if(seed)
  n <- length(y)
  train <- integer()
  lev <- unique(y[!is.na(y)])
  for (lv in lev) {
    idx <- which(y == lv)
    ni <- max(1L, min(length(idx) - as.integer(length(idx) > 1L), round(length(idx) * prop)))
    if (length(idx) == 1L) ni <- 1L
    train <- c(train, sample(idx, ni))
  }
  sort(unique(train))
}

.r4vn_machine_split <- function(dat, outcome, task, split = .8, seed = NA_integer_) {
  if (identical(split, FALSE)) {
    return(list(train = dat, internal = dat[0, , drop = FALSE], test = dat[0, , drop = FALSE]))
  }
  if (!is.numeric(split) || !length(split)) stop("`split` must be FALSE or numeric.", call. = FALSE)
  if (length(split) == 1L) {
    if (split <= 0 || split >= 1) stop("Single `split` must be strictly between 0 and 1.", call. = FALSE)
    if (task %in% c("binary", "multiclass")) {
      tr <- .r4vn_machine_strata_indices(dat[[outcome]], split, seed)
    } else {
      .r4vn_set_seed_if(seed); tr <- sample(seq_len(nrow(dat)), max(1L, floor(nrow(dat) * split)))
    }
    te <- setdiff(seq_len(nrow(dat)), tr)
    return(list(train = dat[tr, , drop = FALSE], internal = dat[0, , drop = FALSE], test = dat[te, , drop = FALSE]))
  }
  if (length(split) != 3L || any(split < 0) || abs(sum(split) - 1) > 1e-8) {
    stop("Length-three `split` must contain nonnegative proportions summing to 1.", call. = FALSE)
  }
  # First select test, then internal validation from the remaining observations.
  train_internal_prop <- split[1L] + split[2L]
  if (task %in% c("binary", "multiclass")) {
    ti <- .r4vn_machine_strata_indices(dat[[outcome]], train_internal_prop, seed)
  } else {
    .r4vn_set_seed_if(seed); ti <- sample(seq_len(nrow(dat)), floor(nrow(dat) * train_internal_prop))
  }
  test <- setdiff(seq_len(nrow(dat)), ti)
  rem <- dat[ti, , drop = FALSE]
  ptrain <- if (train_internal_prop == 0) 0 else split[1L] / train_internal_prop
  if (ptrain <= 0 || ptrain >= 1) {
    tr_local <- if (ptrain >= 1) seq_len(nrow(rem)) else integer()
  } else if (task %in% c("binary", "multiclass")) {
    tr_local <- .r4vn_machine_strata_indices(rem[[outcome]], ptrain, seed + 1L)
  } else {
    .r4vn_set_seed_if(seed + 1L); tr_local <- sample(seq_len(nrow(rem)), floor(nrow(rem) * ptrain))
  }
  in_local <- setdiff(seq_len(nrow(rem)), tr_local)
  list(train = rem[tr_local, , drop = FALSE], internal = rem[in_local, , drop = FALSE], test = dat[test, , drop = FALSE])
}

.r4vn_machine_make_folds <- function(y, k = 10L, repeats = 1L, task, seed = NA_integer_) {
  n <- length(y)
  k <- min(as.integer(k), n)
  if (task %in% c("binary", "multiclass")) {
    tb <- table(y, useNA = "no")
    if (!length(tb) || min(tb) < 2L) stop("Each outcome class needs at least two training observations for cross-validation.", call. = FALSE)
    k <- min(k, as.integer(min(tb)))
  }
  if (k < 2L) stop("At least two cross-validation folds are required.", call. = FALSE)
  out <- list(); ii <- 0L
  for (r in seq_len(as.integer(repeats))) {
    .r4vn_set_seed_if(seed + r)
    fold_id <- integer(n)
    if (task %in% c("binary", "multiclass")) {
      lev <- unique(y)
      for (lv in lev) {
        idx <- sample(which(y == lv))
        fold_id[idx] <- rep(seq_len(k), length.out = length(idx))
      }
    } else {
      idx <- sample(seq_len(n))
      fold_id[idx] <- rep(seq_len(k), length.out = n)
    }
    for (f in seq_len(k)) {
      val <- which(fold_id == f)
      if (!length(val)) next
      ii <- ii + 1L
      out[[ii]] <- list(train = setdiff(seq_len(n), val), val = val, rep_id = r, fold = f)
    }
  }
  out
}

.r4vn_machine_metric <- function(metric, task, event_rate = NA_real_) {
  m <- tolower(as.character(metric)[1L])
  alias <- c(sens="sensitivity",spec="specificity",roc_auc="auc",prauc="pr_auc",f1_score="f1")
  if(m %in% names(alias)) m <- unname(alias[[m]])
  if (m == "auto") {
    m <- if (task == "binary") {
      if (is.finite(event_rate) && min(event_rate, 1 - event_rate) < .10) "pr_auc" else "auc"
    } else if (task == "multiclass") "macro_f1" else "rmse"
  }
  allowed <- if(task=="binary") c("auc","pr_auc","accuracy","balanced_accuracy","sensitivity","specificity","ppv","npv","f1","mcc","kappa","brier","logloss") else if(task=="multiclass") c("accuracy","balanced_accuracy","macro_f1","weighted_f1","macro_auc","macro_pr_auc","logloss") else c("rmse","mae","r2","mape")
  if(!m %in% allowed) stop("Unsupported `metric` for ",task,": ",m,". Allowed: ",paste(allowed,collapse=", "),".",call.=FALSE)
  m
}

.r4vn_machine_metric_direction <- function(metric) {
  if (tolower(metric) %in% c("rmse", "mae", "mape", "brier", "logloss")) "min" else "max"
}

.r4vn_machine_resolve_select <- function(select, train, predictors, task) {
  if (select != "auto") return(select)
  # Keep the default deterministic and dependency-light. More aggressive
  # selectors remain available through explicit `select=` or `select="compare"`.
  "filter"
}

.r4vn_machine_resolve_balance <- function(balance, task, event_rate) {
  if (task != "binary") return("none")
  if (balance != "auto") return(balance)
  if (is.finite(event_rate) && min(event_rate, 1 - event_rate) < .20) "weight" else "none"
}

.r4vn_machine_method_package <- function(method) {
  switch(method,
         lasso = "glmnet", ridge = "glmnet", elastic = "glmnet",
         tree = "rpart", rf = "ranger", xgb = "xgboost",
         svm = "e1071", naive = "e1071", multinom = "nnet",
         NULL)
}

.r4vn_machine_method_ok <- function(method, task) {
  allowed <- switch(task,
                    binary = c("logistic", "lasso", "ridge", "elastic", "tree", "rf", "xgb", "svm", "knn", "naive"),
                    regression = c("linear", "lasso", "ridge", "elastic", "tree", "rf", "xgb", "svm", "knn"),
                    multiclass = c("multinom", "lasso", "ridge", "elastic", "tree", "rf", "xgb", "svm", "knn", "naive"))
  method %in% allowed
}

.r4vn_machine_methods <- function(method, task) {
  req <- unique(tolower(as.character(method)))
  if (identical(req, "auto")) {
    # R4VN auto is intentionally low-dependency. Optional ML engines join only
    # when the user asks for them explicitly or requests method = "all".
    req <- switch(task,
                  binary = c("logistic", "tree"),
                  regression = c("linear", "tree"),
                  multiclass = c("multinom", "tree", "knn"))
  } else if (identical(req, "all")) {
    req <- switch(task,
                  binary = c("logistic", "lasso", "ridge", "elastic", "tree", "rf", "xgb", "svm", "knn", "naive"),
                  regression = c("linear", "lasso", "ridge", "elastic", "tree", "rf", "xgb", "svm", "knn"),
                  multiclass = c("multinom", "lasso", "ridge", "elastic", "tree", "rf", "xgb", "svm", "knn", "naive"))
  } else {
    bad <- req[!vapply(req, .r4vn_machine_method_ok, logical(1), task = task)]
    if (length(bad)) {
      allowed <- switch(task,
        binary = c("logistic", "lasso", "ridge", "elastic", "tree", "rf", "xgb", "svm", "knn", "naive"),
        regression = c("linear", "lasso", "ridge", "elastic", "tree", "rf", "xgb", "svm", "knn"),
        multiclass = c("multinom", "lasso", "ridge", "elastic", "tree", "rf", "xgb", "svm", "knn", "naive")
      )
      stop(
        "Unsupported `method` for ", task, ": ", paste(bad, collapse = ", "),
        ". Available methods are: ", paste(allowed, collapse = ", "), ".",
        call. = FALSE
      )
    }
  }
  req <- req[vapply(req, .r4vn_machine_method_ok, logical(1), task = task)]
  keep <- logical(length(req))
  for (i in seq_along(req)) {
    pkg <- .r4vn_machine_method_package(req[i])
    keep[i] <- is.null(pkg) || requireNamespace(pkg, quietly = TRUE)
  }
  out <- req[keep]
  # Keep a decision-tree fallback when the compact auto set collapses because
  # optional engines are unavailable. Explicit method requests are not altered.
  if (identical(tolower(as.character(method))[1L], "auto") && length(out) <= 1L &&
      requireNamespace("rpart", quietly=TRUE) && .r4vn_machine_method_ok("tree",task)) out <- unique(c(out,"tree"))
  out
}

.r4vn_machine_engine_table <- function(method, task, used) {
  supported <- switch(task,
    binary = c("logistic", "tree", "knn", "lasso", "ridge", "elastic", "rf", "xgb", "svm", "naive"),
    regression = c("linear", "tree", "knn", "lasso", "ridge", "elastic", "rf", "xgb", "svm"),
    multiclass = c("multinom", "tree", "knn", "lasso", "ridge", "elastic", "rf", "xgb", "svm", "naive")
  )
  req0 <- unique(tolower(as.character(method)))
  requested <- if (identical(req0, "auto")) {
    switch(task, binary = c("logistic", "tree"), regression = c("linear", "tree"), multiclass = c("multinom", "tree", "knn"))
  } else if (identical(req0, "all")) supported else intersect(req0, supported)
  pkg <- vapply(supported, function(m) {
    p <- .r4vn_machine_method_package(m)
    if (is.null(p)) "Base R / R4VN" else p
  }, character(1))
  available <- vapply(supported, function(m) {
    p <- .r4vn_machine_method_package(m)
    is.null(p) || requireNamespace(p, quietly = TRUE)
  }, logical(1))
  core <- supported %in% switch(task,
    binary = c("logistic", "tree"),
    regression = c("linear", "tree"),
    multiclass = c("multinom", "tree", "knn")
  )
  data.frame(
    Method = supported,
    Role = ifelse(core, "R4VN core", "Optional"),
    Package = pkg,
    Requested = ifelse(supported %in% requested, "Yes", ""),
    Available = ifelse(available, "Yes", "No"),
    Used = ifelse(supported %in% used, "Yes", ""),
    stringsAsFactors = FALSE, check.names = FALSE
  )
}

# ============================================================================
# Preprocessing blueprint
# ============================================================================

.r4vn_machine_mode <- function(x) {
  z <- x[!is.na(x)]
  if (!length(z)) return(NA)
  tab <- table(z, useNA = "no")
  names(tab)[which.max(tab)]
}

.r4vn_machine_nzv <- function(x) {
  z <- x[!is.na(x)]
  if (length(z) < 2L) return(TRUE)
  tab <- sort(table(z), decreasing = TRUE)
  unique_pct <- 100 * length(tab) / length(z)
  freq_ratio <- if (length(tab) >= 2L) as.numeric(tab[1L] / tab[2L]) else Inf
  unique_pct <= 10 && freq_ratio >= 19
}

.r4vn_machine_yj <- function(x, lambda) {
  out <- rep(NA_real_, length(x))
  pos <- !is.na(x) & x >= 0
  neg <- !is.na(x) & x < 0
  if (abs(lambda) < 1e-8) out[pos] <- log1p(x[pos])
  else out[pos] <- ((x[pos] + 1)^lambda - 1) / lambda
  if (abs(lambda - 2) < 1e-8) out[neg] <- -log1p(-x[neg])
  else out[neg] <- -((1 - x[neg])^(2 - lambda) - 1) / (2 - lambda)
  out
}

.r4vn_machine_yj_lambda <- function(x) {
  z <- x[is.finite(x)]
  if (length(z) < 10L || stats::sd(z) == 0) return(1)
  grid <- seq(-2, 2, by = .1)
  ll <- vapply(grid, function(lam) {
    t <- .r4vn_machine_yj(z, lam)
    vv <- stats::var(t)
    if (!is.finite(vv) || vv <= 0) return(-Inf)
    jac <- (lam - 1) * sum(sign(z) * log1p(abs(z)))
    -length(z) / 2 * log(vv) + jac
  }, numeric(1))
  grid[which.max(ll)]
}

.r4vn_machine_corr_remove <- function(dat, vars, threshold) {
  num <- vars[vapply(dat[vars], is.numeric, logical(1))]
  if (length(num) < 2L || is.null(threshold) || !is.finite(threshold)) return(character())
  cc <- suppressWarnings(stats::cor(dat[num], use = "pairwise.complete.obs"))
  if (!is.matrix(cc) || nrow(cc) < 2L) return(character())
  diag(cc) <- 0
  remove <- character()
  while (TRUE) {
    aa <- abs(cc)
    aa[!is.finite(aa)] <- 0
    mx <- max(aa)
    if (mx < threshold) break
    pos <- which(aa == mx, arr.ind = TRUE)[1L, ]
    v1 <- rownames(cc)[pos[1L]]; v2 <- colnames(cc)[pos[2L]]
    mean1 <- mean(abs(cc[v1, ]), na.rm = TRUE)
    mean2 <- mean(abs(cc[v2, ]), na.rm = TRUE)
    kill <- if (mean1 >= mean2) v1 else v2
    remove <- c(remove, kill)
    cc[kill, ] <- 0; cc[, kill] <- 0
  }
  unique(remove)
}

.r4vn_machine_corr_threshold <- function(corr) {
  if (identical(corr, FALSE) || identical(corr, "none") || is.null(corr)) return(NULL)
  if (is.character(corr) && identical(tolower(corr[1L]), "auto")) return(.95)
  if (is.numeric(corr) && length(corr) == 1L && corr > 0 && corr < 1) return(as.numeric(corr))
  stop("`corr` must be 'auto', 'none'/FALSE, or a number strictly between 0 and 1.", call. = FALSE)
}

.r4vn_machine_feature_mode <- function(feature) {
  f <- tolower(as.character(feature)[1L])
  if(f == "auto") "none" else f
}

.r4vn_machine_feature_expand_fit <- function(mm, raw_map, types, feature="none", degree=2L) {
  mode <- .r4vn_machine_feature_mode(feature)
  spec <- list(mode=mode, degree=as.integer(degree), powers=list(), interactions=list())
  out <- mm; map <- raw_map
  if(mode == "none") return(list(x=out,map=map,spec=spec))
  numeric_cols <- names(raw_map)[vapply(names(raw_map), function(nm) {
    rv <- raw_map[[nm]]
    !is.na(rv) && rv %in% names(types) && identical(types[[rv]], "numeric")
  }, logical(1))]
  numeric_cols <- intersect(numeric_cols,colnames(mm))
  if(mode %in% c("polynomial","all") && length(numeric_cols)) {
    for(nm in numeric_cols) for(pow in 2L:as.integer(degree)) {
      nn <- paste0(nm,"__pow",pow); out <- cbind(out, stats::setNames(data.frame(mm[,nm]^pow),nn));
      spec$powers[[nn]] <- list(source=nm,power=pow); map[[nn]] <- raw_map[[nm]]
    }
  }
  if(mode %in% c("interaction","all") && length(numeric_cols)>=2L) {
    cmb <- utils::combn(numeric_cols,2L,simplify=FALSE)
    if(length(cmb)>200L) cmb <- cmb[seq_len(200L)]
    for(cc in cmb) {
      nn <- paste0(cc[1L],"__X__",cc[2L]); out <- cbind(out, stats::setNames(data.frame(mm[,cc[1L]]*mm[,cc[2L]]),nn))
      spec$interactions[[nn]] <- cc
      rv <- unique(na.omit(c(raw_map[[cc[1L]]],raw_map[[cc[2L]]]))); map[[nn]] <- paste(rv,collapse=" x ")
    }
  }
  out <- as.matrix(out); storage.mode(out) <- "double"
  list(x=out,map=map[colnames(out)],spec=spec)
}

.r4vn_machine_feature_expand_apply <- function(mm, spec) {
  out <- mm
  if(is.null(spec) || identical(spec$mode,"none")) return(out)
  if(length(spec$powers)) for(nn in names(spec$powers)) {
    z <- spec$powers[[nn]]; out <- cbind(out, stats::setNames(data.frame(out[,z$source]^z$power),nn))
  }
  if(length(spec$interactions)) for(nn in names(spec$interactions)) {
    cc <- spec$interactions[[nn]]; out <- cbind(out, stats::setNames(data.frame(out[,cc[1L]]*out[,cc[2L]]),nn))
  }
  out <- as.matrix(out); storage.mode(out) <- "double"; out
}

.r4vn_machine_pca_fit <- function(x, reduce="none", variance=.95) {
  mode <- tolower(as.character(reduce)[1L])
  if(mode == "auto") mode <- if(ncol(x)>=50L && ncol(x)>nrow(x)) "pca" else "none"
  if(mode == "none") return(list(method="none",x=x,center=NULL,scale=NULL,rotation=NULL,variance=NULL,k=ncol(x)))
  cen <- colMeans(x); sc <- apply(x,2L,stats::sd); sc[!is.finite(sc)|sc==0] <- 1
  xs <- sweep(sweep(x,2L,cen,"-"),2L,sc,"/")
  pc <- stats::prcomp(xs,center=FALSE,scale.=FALSE)
  vv <- pc$sdev^2; prop <- vv/sum(vv); cum <- cumsum(prop)
  k <- which(cum>=variance)[1L]; if(is.na(k)) k <- length(cum); k <- max(1L,k)
  rot <- pc$rotation[,seq_len(k),drop=FALSE]
  xx <- xs %*% rot; colnames(xx) <- paste0("PC",seq_len(k))
  list(method="pca",x=xx,center=cen,scale=sc,rotation=rot,variance=prop[seq_len(k)],k=k,cumulative=cum[k])
}

.r4vn_machine_pca_apply <- function(x, pca) {
  if(is.null(pca) || identical(pca$method,"none")) return(x)
  xs <- sweep(sweep(x,2L,pca$center,"-"),2L,pca$scale,"/")
  z <- xs %*% pca$rotation; colnames(z) <- paste0("PC",seq_len(ncol(z))); z
}

.r4vn_machine_recipe_fit <- function(raw, predictors, preprocess, missing, missing_max,
                                     encode, transform, outlier, corr, feature="none", degree=2L,
                                     reduce="none", variance=.95) {
  d <- raw[, predictors, drop = FALSE]
  actions <- list()
  removed <- character()

  miss_prop <- vapply(d, function(z) mean(is.na(z)), numeric(1))
  if (preprocess == "auto") {
    rm_miss <- names(miss_prop)[miss_prop > missing_max]
    if (length(rm_miss)) {
      removed <- c(removed, rm_miss)
      actions[[length(actions) + 1L]] <- data.frame(step = "High missingness", variable = rm_miss,
                                                    action = paste0("Removed (> ", round(100 * missing_max), "% missing)"))
    }
  }
  keep <- setdiff(names(d), removed)
  if (!length(keep)) stop("All predictors were removed by missingness filtering.", call. = FALSE)

  if (preprocess == "auto") {
    constant <- keep[vapply(d[keep], function(z) length(unique(z[!is.na(z)])) <= 1L, logical(1))]
    nzv <- setdiff(keep[vapply(d[keep], .r4vn_machine_nzv, logical(1))], constant)
    if (length(constant)) {
      removed <- c(removed, constant)
      actions[[length(actions) + 1L]] <- data.frame(step = "Zero variance", variable = constant, action = "Removed")
    }
    if (length(nzv)) {
      removed <- c(removed, nzv)
      actions[[length(actions) + 1L]] <- data.frame(step = "Near-zero variance", variable = nzv, action = "Removed")
    }
  }
  keep <- setdiff(names(d), removed)

  cth <- .r4vn_machine_corr_threshold(corr)
  if (preprocess == "auto" && length(keep) && !is.null(cth)) {
    cr <- .r4vn_machine_corr_remove(d, keep, cth)
    if (length(cr)) {
      removed <- c(removed, cr)
      actions[[length(actions) + 1L]] <- data.frame(step = "High correlation", variable = cr,
                                                    action = paste0("Removed (|r| >= ", cth, ")"))
    }
  }
  keep <- setdiff(names(d), removed)
  if (!length(keep)) stop("No predictors remain after preprocessing filters.", call. = FALSE)

  # Factor blueprints and imputation learned on training only.
  types <- vapply(d[keep], function(z) if (is.numeric(z)) "numeric" else "categorical", character(1))
  levels_map <- list(); impute <- list(); transform_map <- list(); outlier_map <- list()
  proc <- d[keep]
  complete_rows <- rep(TRUE, nrow(proc))

  for (v in keep) {
    z <- proc[[v]]
    if (types[[v]] == "categorical") {
      zz <- as.character(z)
      lev <- unique(zz[!is.na(zz)])
      if (!length(lev)) lev <- "Missing"
      levels_map[[v]] <- lev
      if (missing == "complete") {
        complete_rows <- complete_rows & !is.na(zz)
      } else {
        mv <- .r4vn_machine_mode(zz)
        if (is.na(mv)) mv <- lev[1L]
        impute[[v]] <- mv
        zz[is.na(zz)] <- mv
      }
      proc[[v]] <- factor(zz, levels = levels_map[[v]])
    } else {
      zz <- suppressWarnings(as.numeric(z))
      if (missing == "complete") {
        complete_rows <- complete_rows & is.finite(zz)
      } else {
        meth <- if (missing == "auto") "median" else missing
        if (meth == "mode") {
          mv <- suppressWarnings(as.numeric(.r4vn_machine_mode(zz)))
        } else if (meth == "mean") {
          mv <- mean(zz, na.rm = TRUE)
        } else {
          mv <- stats::median(zz, na.rm = TRUE)
        }
        if (!is.finite(mv)) mv <- 0
        impute[[v]] <- mv
        zz[!is.finite(zz)] <- mv
      }

      if (outlier %in% c("detect", "winsor", "robust")) {
        q <- stats::quantile(zz, c(.25, .75), na.rm = TRUE, names = FALSE)
        iqr <- q[2L] - q[1L]
        detect_low <- q[1L] - 1.5 * iqr; detect_high <- q[2L] + 1.5 * iqr
        n_out <- sum(zz < detect_low | zz > detect_high, na.rm = TRUE)
        if (outlier == "winsor") {
          cuts <- stats::quantile(zz, c(.01, .99), na.rm = TRUE, names = FALSE)
          outlier_map[[v]] <- list(method = "winsor", low = cuts[1L], high = cuts[2L], detected = n_out)
          zz <- pmin(pmax(zz, cuts[1L]), cuts[2L])
        } else if (outlier == "robust") {
          md <- stats::median(zz, na.rm = TRUE); ma <- stats::mad(zz, center = md, constant = 1, na.rm = TRUE)
          if (!is.finite(ma) || ma == 0) ma <- stats::sd(zz, na.rm = TRUE)
          cuts <- c(md - 5 * ma, md + 5 * ma)
          outlier_map[[v]] <- list(method = "robust", low = cuts[1L], high = cuts[2L], detected = n_out)
          zz <- pmin(pmax(zz, cuts[1L]), cuts[2L])
        } else {
          outlier_map[[v]] <- list(method = "detect", low = detect_low, high = detect_high, detected = n_out)
        }
      }

      tr <- transform
      if (tr == "auto") tr <- "none"
      if (tr == "log") {
        minz <- min(zz, na.rm = TRUE)
        shift <- if (minz <= 0) 1 - minz else 0
        transform_map[[v]] <- list(method = "log", shift = shift)
        zz <- log(zz + shift)
      } else if (tr == "yeojohnson") {
        lam <- .r4vn_machine_yj_lambda(zz)
        transform_map[[v]] <- list(method = "yeojohnson", lambda = lam)
        zz <- .r4vn_machine_yj(zz, lam)
      }
      proc[[v]] <- zz
    }
  }

  if (missing == "complete") proc <- proc[complete_rows, , drop = FALSE]
  if (!nrow(proc)) stop("No complete training rows remain after `missing = 'complete'`.", call. = FALSE)

  # model.matrix blueprint learned from training data only.
  # Build model.frame first so its expanded terms object is available and
  # consistent with the model matrix. Do not rely on attr(model.matrix(), "terms").
  frm <- stats::as.formula("~ .")
  mf <- stats::model.frame(
    frm,
    data = proc,
    na.action = stats::na.pass,
    drop.unused.levels = FALSE
  )
  trm <- attr(mf, "terms")
  if (is.null(trm)) {
    stop("Could not create the training design terms.", call. = FALSE)
  }

  mm_full <- stats::model.matrix(
    trm,
    data = mf,
    contrasts.arg = NULL
  )

  assign <- attr(mm_full, "assign")
  term_labels <- attr(trm, "term.labels")

  if (is.null(assign) || length(assign) != ncol(mm_full)) {
    stop("The training model matrix does not contain a valid term assignment.", call. = FALSE)
  }
  if (!length(term_labels) && any(assign > 0L)) {
    stop("Could not recover predictor term labels from the training design.", call. = FALSE)
  }

  # Map each encoded column back to its original predictor.
  raw_map_full <- stats::setNames(
    rep(NA_character_, ncol(mm_full)),
    colnames(mm_full)
  )

  for (j in seq_along(assign)) {
    term_id <- assign[[j]]
    if (is.na(term_id) || term_id == 0L) next
    if (term_id < 1L || term_id > length(term_labels)) {
      stop(
        "Invalid model-matrix term assignment while building the preprocessing blueprint.",
        call. = FALSE
      )
    }
    raw_map_full[[j]] <- gsub(
      "`", "", term_labels[[term_id]], fixed = TRUE
    )
  }

  keep_mm <- colnames(mm_full) != "(Intercept)"
  mm <- mm_full[, keep_mm, drop = FALSE]
  if (!ncol(mm)) {
    stop("No encoded predictor columns remain.", call. = FALSE)
  }

  mm_cols <- colnames(mm)
  raw_map <- raw_map_full[mm_cols]

  # Defensive fallback: with the base ~ . blueprint every non-intercept encoded
  # column should map to one original predictor.
  if (anyNA(raw_map)) {
    bad <- names(raw_map)[is.na(raw_map)]
    stop(
      "Could not map encoded predictor column(s) back to original variables: ",
      paste(bad, collapse = ", "),
      call. = FALSE
    )
  }

  base_cols <- mm_cols

  fx <- .r4vn_machine_feature_expand_fit(mm, raw_map, types, feature=feature, degree=degree)
  mm_feature <- fx$x; raw_map_feature <- fx$map
  fsd <- apply(mm_feature,2L,stats::sd,na.rm=TRUE); good_fx <- is.finite(fsd) & fsd > 0
  if(any(!good_fx)) mm_feature <- mm_feature[,good_fx,drop=FALSE]
  if(!ncol(mm_feature)) stop("No non-constant encoded/engineered features remain.",call.=FALSE)
  feature_cols <- colnames(mm_feature); raw_map_feature <- raw_map_feature[feature_cols]
  pca <- .r4vn_machine_pca_fit(mm_feature, reduce=reduce, variance=variance)
  mm_final <- pca$x; mm_cols <- colnames(mm_final)
  raw_map <- if(identical(pca$method,"pca")) stats::setNames(rep(NA_character_,length(mm_cols)),mm_cols) else raw_map_feature[mm_cols]

  if (length(actions)) actions_df <- do.call(rbind, actions) else actions_df <- data.frame(step=character(), variable=character(), action=character())
  for (v in keep) {
    if (!is.null(impute[[v]])) {
      actions_df <- rbind(actions_df, data.frame(step = "Imputation", variable = v,
                                                 action = paste0(if (types[[v]] == "numeric") {
                                                   if (missing == "auto") "Median" else tools::toTitleCase(missing)
                                                 } else "Mode", ": ", format(impute[[v]], digits = 6))))
    }
    if (!is.null(transform_map[[v]])) actions_df <- rbind(actions_df, data.frame(step="Transformation", variable=v, action=transform_map[[v]]$method))
    if (!is.null(outlier_map[[v]])) actions_df <- rbind(actions_df, data.frame(step="Outlier", variable=v,
      action=paste0(outlier_map[[v]]$method, "; detected=", outlier_map[[v]]$detected)))
  }
  actions_df <- rbind(actions_df, data.frame(step="Encoding", variable="Categorical predictors", action="Training-derived dummy/model matrix"))
  if(!identical(fx$spec$mode,"none")) actions_df <- rbind(actions_df,data.frame(step="Feature engineering",variable="Numeric predictors",action=paste0(fx$spec$mode,"; degree=",degree,"; encoded features=",ncol(mm_feature))))
  if(identical(pca$method,"pca")) actions_df <- rbind(actions_df,data.frame(step="Dimensionality reduction",variable="Encoded predictors",action=paste0("PCA retained ",pca$k," components; cumulative variance=",formatC(pca$cumulative,format="f",digits=3))))

  list(
    predictors = predictors,
    keep = keep,
    removed = unique(removed),
    types = types,
    levels = levels_map,
    impute = impute,
    transform = transform_map,
    outlier = outlier_map,
    base_cols = base_cols,
    feature_spec = fx$spec,
    feature_cols = feature_cols,
    pca = pca,
    mm_cols = mm_cols,
    raw_map = raw_map,
    missing = missing,
    complete_rows = if (missing == "complete") which(complete_rows) else seq_len(nrow(raw)),
    actions = actions_df
  )
}

.r4vn_machine_recipe_apply <- function(recipe, raw) {
  d <- raw[, recipe$keep, drop = FALSE]
  keep_rows <- rep(TRUE, nrow(d))
  unseen <- character()
  for (v in recipe$keep) {
    if (recipe$types[[v]] == "categorical") {
      z <- as.character(d[[v]])
      if (recipe$missing == "complete") {
        keep_rows <- keep_rows & !is.na(z)
      } else {
        z[is.na(z)] <- as.character(recipe$impute[[v]])
      }
      bad <- !is.na(z) & !z %in% recipe$levels[[v]]
      if (any(bad)) {
        unseen <- c(unseen, paste0(v, ": ", paste(unique(z[bad]), collapse = ", ")))
        z[bad] <- recipe$levels[[v]][1L]
      }
      d[[v]] <- factor(z, levels = recipe$levels[[v]])
    } else {
      z <- suppressWarnings(as.numeric(d[[v]]))
      if (recipe$missing == "complete") keep_rows <- keep_rows & is.finite(z)
      else z[!is.finite(z)] <- as.numeric(recipe$impute[[v]])
      oo <- recipe$outlier[[v]]
      if (!is.null(oo) && oo$method %in% c("winsor", "robust")) z <- pmin(pmax(z, oo$low), oo$high)
      tt <- recipe$transform[[v]]
      if (!is.null(tt) && tt$method == "log") z <- log(z + tt$shift)
      if (!is.null(tt) && tt$method == "yeojohnson") z <- .r4vn_machine_yj(z, tt$lambda)
      d[[v]] <- z
    }
  }
  if (recipe$missing == "complete") d2 <- d[keep_rows, , drop = FALSE] else d2 <- d
  mm <- stats::model.matrix(~ ., data = d2)
  if ("(Intercept)" %in% colnames(mm)) mm <- mm[, colnames(mm) != "(Intercept)", drop = FALSE]
  # Align base encoded columns with the training blueprint.
  miss_cols <- setdiff(recipe$base_cols, colnames(mm))
  if (length(miss_cols)) {
    add <- matrix(0, nrow(mm), length(miss_cols), dimnames = list(NULL, miss_cols))
    mm <- cbind(mm, add)
  }
  extra <- setdiff(colnames(mm), recipe$base_cols)
  if (length(extra)) mm <- mm[, setdiff(colnames(mm), extra), drop = FALSE]
  mm <- mm[, recipe$base_cols, drop = FALSE]
  mm <- .r4vn_machine_feature_expand_apply(mm, recipe$feature_spec)
  miss_fx <- setdiff(recipe$feature_cols,colnames(mm))
  if(length(miss_fx)){add<-matrix(0,nrow(mm),length(miss_fx),dimnames=list(NULL,miss_fx));mm<-cbind(mm,add)}
  mm <- mm[,recipe$feature_cols,drop=FALSE]
  mm <- .r4vn_machine_pca_apply(mm, recipe$pca)
  mm <- mm[,recipe$mm_cols,drop=FALSE]
  list(x = mm, rows = which(keep_rows), unseen = unique(unseen))
}

.r4vn_machine_scale_fit <- function(x, method) {
  if (method == "none") return(list(method = "none", center = NULL, scale = NULL))
  if (method == "z") {
    cen <- colMeans(x, na.rm = TRUE); sc <- apply(x, 2L, stats::sd, na.rm = TRUE)
  } else if (method == "minmax") {
    lo <- apply(x, 2L, min, na.rm = TRUE); hi <- apply(x, 2L, max, na.rm = TRUE)
    cen <- lo; sc <- hi - lo
  } else {
    cen <- apply(x, 2L, stats::median, na.rm = TRUE); sc <- apply(x, 2L, stats::mad, na.rm = TRUE)
  }
  sc[!is.finite(sc) | sc == 0] <- 1
  list(method = method, center = cen, scale = sc)
}

.r4vn_machine_scale_apply <- function(x, scaler) {
  if (is.null(scaler) || scaler$method == "none") return(x)
  sweep(sweep(x, 2L, scaler$center, "-"), 2L, scaler$scale, "/")
}

.r4vn_machine_needs_scale <- function(method) method %in% c("lasso", "ridge", "elastic", "svm", "knn", "multinom")

.r4vn_machine_scale_method <- function(standardize, model_method) {
  if (standardize == "none") return("none")
  if (standardize == "auto") return(if (.r4vn_machine_needs_scale(model_method)) "z" else "none")
  standardize
}

# ============================================================================
# Feature selection and class balancing
# ============================================================================

.r4vn_machine_auto_nfeatures <- function(p, n) {
  max(3L, min(p, max(5L, floor(sqrt(max(1, p * n / max(p, 1)))))))
}

.r4vn_machine_univariate_score <- function(x, y, task) {
  if (task == "regression") {
    return(abs(suppressWarnings(stats::cor(x, y, use = "complete.obs"))))
  }
  if (task == "binary") {
    return(abs(mean(x[y == 1L], na.rm = TRUE) - mean(x[y == 0L], na.rm = TRUE)) /
             max(stats::sd(x, na.rm = TRUE), 1e-8))
  }
  # multiclass eta-squared-like between/total variation
  grand <- mean(x, na.rm = TRUE)
  lev <- levels(y)
  ssb <- sum(vapply(lev, function(g) {
    z <- x[y == g]; sum(!is.na(z)) * (mean(z, na.rm = TRUE) - grand)^2
  }, numeric(1)), na.rm = TRUE)
  sst <- sum((x - grand)^2, na.rm = TRUE)
  if (!is.finite(sst) || sst == 0) 0 else ssb / sst
}

.r4vn_machine_feature_select <- function(x, y, task, method, nfeatures = "auto", seed = NA_integer_) {
  p <- ncol(x)
  if (!p) stop("No encoded predictors available for feature selection.", call. = FALSE)
  if (method %in% c("none", "filter")) return(list(selected = colnames(x), detail = data.frame()))
  nf <- if (is.character(nfeatures) && nfeatures[1L] == "auto") .r4vn_machine_auto_nfeatures(p, nrow(x)) else as.integer(nfeatures[1L])
  nf <- max(1L, min(p, nf))

  if (method == "lasso") {
    if (!requireNamespace("glmnet", quietly = TRUE)) {
      return(list(selected = colnames(x), detail = data.frame(), note = "glmnet unavailable; feature selection fell back to filter."))
    }
    fam <- if (task == "binary") "binomial" else if (task == "multiclass") "multinomial" else "gaussian"
    yy <- if (task == "multiclass") y else as.numeric(y)
    .r4vn_set_seed_if(seed)
    cv <- tryCatch(glmnet::cv.glmnet(x, yy, family = fam, alpha = 1,
                                     nfolds = max(3L, min(5L, nrow(x))), standardize = FALSE),
                   error = function(e) NULL)
    if (is.null(cv)) return(list(selected = colnames(x), detail = data.frame(), note = "LASSO selection failed; filter selection used."))
    cf <- tryCatch(stats::coef(cv, s = "lambda.1se"), error = function(e) NULL)
    sel <- character()
    detail <- data.frame()
    if (task == "multiclass" && is.list(cf)) {
      rn <- unique(unlist(lapply(cf, rownames)))
      rn <- setdiff(rn, "(Intercept)")
      keep <- vapply(rn, function(v) any(vapply(cf, function(m) v %in% rownames(m) && as.numeric(m[v, 1L]) != 0, logical(1))), logical(1))
      sel <- rn[keep]
    } else if (!is.null(cf)) {
      vv <- as.matrix(cf)
      sel <- rownames(vv)[rownames(vv) != "(Intercept)" & as.numeric(vv[, 1L]) != 0]
      detail <- data.frame(feature = setdiff(rownames(vv), "(Intercept)"),
                           coefficient = as.numeric(vv[rownames(vv) != "(Intercept)", 1L]), stringsAsFactors = FALSE)
    }
    sel <- intersect(sel, colnames(x))
    if (!length(sel)) sel <- colnames(x)[seq_len(min(nf, p))]
    return(list(selected = sel, detail = detail))
  }

  if (method == "stepwise") {
    df <- data.frame(.y = y, x, check.names = FALSE)
    fit <- tryCatch({
      if (task == "binary") stats::glm(.y ~ ., data = df, family = stats::binomial())
      else if (task == "regression") stats::lm(.y ~ ., data = df)
      else NULL
    }, error = function(e) NULL)
    if (is.null(fit)) return(list(selected = colnames(x), detail = data.frame(), note = "Stepwise selection is not available for this task/model matrix; filter used."))
    st <- tryCatch(stats::step(fit, trace = 0), error = function(e) NULL)
    if (is.null(st)) return(list(selected = colnames(x), detail = data.frame(), note = "Stepwise selection failed; filter used."))
    sel <- attr(stats::terms(st), "term.labels")
    sel <- intersect(sel, colnames(x))
    if (!length(sel)) sel <- colnames(x)
    return(list(selected = sel, detail = data.frame(feature = sel, stringsAsFactors = FALSE)))
  }

  if (method == "boruta") {
    if (!requireNamespace("Boruta", quietly = TRUE)) {
      return(list(selected = colnames(x), detail = data.frame(), note = "Boruta unavailable; filter selection used."))
    }
    df <- data.frame(x, check.names = FALSE)
    yy <- if (task == "regression") as.numeric(y) else factor(y)
    .r4vn_set_seed_if(seed)
    b <- tryCatch(Boruta::Boruta(x = df, y = yy, doTrace = 0), error = function(e) NULL)
    if (is.null(b)) return(list(selected = colnames(x), detail = data.frame(), note = "Boruta selection failed; filter used."))
    decision <- b$finalDecision
    sel <- names(decision)[decision %in% c("Confirmed", "Tentative")]
    if (!length(sel)) sel <- names(decision)[decision == "Confirmed"]
    if (!length(sel)) sel <- colnames(x)
    detail <- data.frame(feature = names(decision), decision = as.character(decision), stringsAsFactors = FALSE)
    return(list(selected = intersect(sel, colnames(x)), detail = detail))
  }

  scores <- vapply(seq_len(ncol(x)), function(j) .r4vn_machine_univariate_score(x[, j], y, task), numeric(1))
  scores[!is.finite(scores)] <- -Inf
  ord <- order(scores, decreasing = TRUE)

  if (method == "importance") {
    sel <- colnames(x)[ord[seq_len(nf)]]
    return(list(selected = sel,
                detail = data.frame(feature = colnames(x), score = scores, selected = colnames(x) %in% sel, stringsAsFactors = FALSE)))
  }

  if (method == "rfe") {
    # Lightweight fold-local RFE: rank within the training fold, then choose the
    # smallest top-k set within 0.5% of the best simple holdout score.
    sizes <- unique(sort(pmax(1L, pmin(p, c(3L, 5L, 10L, 20L, nf, p)))))
    .r4vn_set_seed_if(seed)
    idx <- sample(seq_len(nrow(x)))
    cut <- max(2L, floor(.8 * nrow(x)))
    tr <- idx[seq_len(cut)]; va <- setdiff(idx, tr)
    if (!length(va)) va <- tr
    vals <- rep(NA_real_, length(sizes))
    for (i in seq_along(sizes)) {
      cols <- colnames(x)[ord[seq_len(sizes[i])]]
      if (task == "binary") {
        df <- data.frame(.y = y[tr], x[tr, cols, drop = FALSE], check.names = FALSE)
        ft <- tryCatch(stats::glm(.y ~ ., data = df, family = stats::binomial()), error = function(e) NULL)
        pp <- if (is.null(ft)) rep(mean(y[tr]), length(va)) else suppressWarnings(as.numeric(stats::predict(ft, newdata = data.frame(x[va, cols, drop=FALSE], check.names=FALSE), type="response")))
        vals[i] <- .r4vn_machine_auc(y[va], pp)
      } else if (task == "regression") {
        df <- data.frame(.y = y[tr], x[tr, cols, drop = FALSE], check.names = FALSE)
        ft <- tryCatch(stats::lm(.y ~ ., data = df), error = function(e) NULL)
        pp <- if (is.null(ft)) rep(mean(y[tr]), length(va)) else suppressWarnings(as.numeric(stats::predict(ft, newdata=data.frame(x[va, cols, drop=FALSE], check.names=FALSE))))
        vals[i] <- sqrt(mean((y[va] - pp)^2, na.rm = TRUE))
      } else {
        vals[i] <- mean(scores[ord[seq_len(sizes[i])]], na.rm = TRUE)
      }
    }
    if (task == "regression") bestv <- min(vals, na.rm = TRUE) else bestv <- max(vals, na.rm = TRUE)
    eligible <- if (task == "regression") which(vals <= bestv * 1.005) else which(vals >= bestv - .005)
    ii <- eligible[which.min(sizes[eligible])]
    sel <- colnames(x)[ord[seq_len(sizes[ii])]]
    return(list(selected = sel,
                detail = data.frame(nfeatures = sizes, score = vals, selected_size = sizes == sizes[ii])))
  }

  list(selected = colnames(x), detail = data.frame())
}

.r4vn_machine_class_weights <- function(y) {
  tab <- table(y)
  wmap <- sum(tab) / (length(tab) * tab)
  as.numeric(wmap[as.character(y)])
}

.r4vn_machine_updown <- function(x, y, method, target, seed) {
  .r4vn_set_seed_if(seed)
  tab <- table(y)
  if (length(tab) != 2L) return(list(x=x, y=y, weights=NULL))
  lev <- names(tab); minlev <- lev[which.min(tab)]; majlev <- lev[which.max(tab)]
  imin <- which(as.character(y) == minlev); imaj <- which(as.character(y) == majlev)
  if (method == "up") {
    n_target_min <- ceiling(length(imaj) * target / (1 - target))
    add <- sample(imin, max(0L, n_target_min - length(imin)), replace = TRUE)
    idx <- c(seq_along(y), add)
  } else {
    n_target_maj <- floor(length(imin) * (1 - target) / target)
    keepmaj <- sample(imaj, min(length(imaj), max(1L, n_target_maj)))
    idx <- c(imin, keepmaj)
  }
  idx <- sample(idx)
  list(x=x[idx,,drop=FALSE], y=y[idx], weights=NULL)
}

.r4vn_machine_knn_dist <- function(x, row) {
  sqrt(rowSums((x - matrix(row, nrow(x), ncol(x), byrow = TRUE))^2))
}

.r4vn_machine_smote <- function(x, y, target, k, seed, adasyn = FALSE) {
  .r4vn_set_seed_if(seed)
  tab <- table(y)
  if (length(tab) != 2L) return(list(x=x, y=y, weights=NULL))
  lev <- names(tab); minlev <- lev[which.min(tab)]; majlev <- lev[which.max(tab)]
  imin <- which(as.character(y) == minlev); imaj <- which(as.character(y) == majlev)
  n_target <- ceiling(length(imaj) * target / (1 - target))
  need <- max(0L, n_target - length(imin))
  if (!need || length(imin) < 2L) return(list(x=x, y=y, weights=NULL))
  kmin <- min(k, length(imin)-1L)
  minx <- x[imin,,drop=FALSE]

  source_prob <- rep(1/length(imin), length(imin))
  if (isTRUE(adasyn)) {
    kall <- min(k, nrow(x)-1L)
    difficulty <- numeric(length(imin))
    for (i in seq_along(imin)) {
      ds <- .r4vn_machine_knn_dist(x, x[imin[i],])
      nn <- order(ds)[seq_len(min(kall + 1L, length(ds)))]
      nn <- setdiff(nn, imin[i])[seq_len(min(kall, length(setdiff(nn, imin[i]))))]
      difficulty[i] <- if (length(nn)) mean(as.character(y[nn]) == majlev) else 0
    }
    if (sum(difficulty) > 0) source_prob <- difficulty / sum(difficulty)
  }

  synth <- matrix(NA_real_, need, ncol(x), dimnames = list(NULL, colnames(x)))
  for (b in seq_len(need)) {
    i <- sample(seq_along(imin), 1L, prob = source_prob)
    ds <- .r4vn_machine_knn_dist(minx, minx[i,])
    nn <- order(ds)[order(ds) != i]
    nn <- nn[seq_len(min(kmin, length(nn)))]
    j <- sample(nn, 1L)
    gap <- stats::runif(1L)
    synth[b,] <- minx[i,] + gap * (minx[j,] - minx[i,])
  }
  yy <- if (is.factor(y)) factor(c(as.character(y), rep(minlev, need)), levels=levels(y)) else c(y, rep(type.convert(minlev, as.is=TRUE), need))
  list(x=rbind(x, synth), y=yy, weights=NULL)
}

.r4vn_machine_balance_apply <- function(x, y, method, target=.5, k=5, seed=NA_integer_) {
  if (method == "none") return(list(x=x, y=y, weights=NULL, method="none"))
  if (method == "weight") return(list(x=x, y=y, weights=.r4vn_machine_class_weights(y), method="weight"))
  if (method %in% c("up","down")) {
    z <- .r4vn_machine_updown(x,y,method,target,seed); z$method <- method; return(z)
  }
  if (method == "smote") {
    z <- .r4vn_machine_smote(x,y,target,k,seed,FALSE); z$method <- method; return(z)
  }
  if (method == "adasyn") {
    z <- .r4vn_machine_smote(x,y,target,k,seed,TRUE); z$method <- method; return(z)
  }
  if (method == "rose") {
    if (requireNamespace("ROSE", quietly = TRUE)) {
      df <- data.frame(.y = factor(y), x, check.names=FALSE)
      .r4vn_set_seed_if(seed)
      rose_args <- list(formula = .y ~ ., data = df, N = nrow(df), p = target)
      if (is.finite(seed)) rose_args$seed <- as.integer(seed)
      rr <- tryCatch(do.call(ROSE::ROSE, rose_args)$data, error=function(e) NULL)
      if (!is.null(rr)) {
        yy <- rr$.y; xx <- as.matrix(rr[, setdiff(names(rr), ".y"), drop=FALSE]); storage.mode(xx) <- "double"
        return(list(x=xx,y=yy,weights=NULL,method="rose"))
      }
    }
    z <- .r4vn_machine_smote(x,y,target,k,seed,FALSE); z$method <- "smote (ROSE unavailable)"; return(z)
  }
  list(x=x,y=y,weights=NULL,method="none")
}

# ============================================================================
# Model engines, tuning, cross-validation, and complete pipelines
# ============================================================================

.r4vn_machine_default_params <- function(method, task, p = 10L) {
  switch(method,
    logistic = list(),
    linear = list(),
    multinom = list(decay = 0, maxit = 500),
    lasso = list(alpha = 1, lambda_rule = "lambda.1se"),
    ridge = list(alpha = 0, lambda_rule = "lambda.1se"),
    elastic = list(alpha = .5, lambda_rule = "lambda.1se"),
    tree = list(cp = .01, minsplit = 20L, maxdepth = 30L),
    rf = list(num.trees = 500L, mtry = max(1L, floor(sqrt(p))), min.node.size = if (task == "regression") 5L else 1L),
    xgb = list(nrounds = 250L, max_depth = 4L, eta = .05, subsample = .8, colsample_bytree = .8, min_child_weight = 1),
    svm = list(kernel = "radial", cost = 1, gamma = 1 / max(1,p)),
    knn = list(k = 7L),
    naive = list(laplace = 0),
    list()
  )
}

.r4vn_machine_tune_grid <- function(method, task, tune = TRUE, tune_n = 10, seed = NA_integer_) {
  if (identical(tune, FALSE)) return(list(.r4vn_machine_default_params(method, task)))
  mode <- if (isTRUE(tune)) "grid" else tolower(as.character(tune)[1L])
  grid <- switch(method,
    tree = list(
      list(cp=.001, minsplit=10L, maxdepth=10L),
      list(cp=.005, minsplit=20L, maxdepth=20L),
      list(cp=.01, minsplit=20L, maxdepth=30L)
    ),
    rf = list(
      list(num.trees=500L, mtry=NULL, mtry_factor=.5, min.node.size=1L),
      list(num.trees=500L, mtry=NULL, mtry_factor=1, min.node.size=if (task=="regression") 5L else 1L),
      list(num.trees=700L, mtry=NULL, mtry_factor=2, min.node.size=if (task=="regression") 10L else 5L)
    ),
    xgb = list(
      list(nrounds=200L,max_depth=3L,eta=.05,subsample=.8,colsample_bytree=.8,min_child_weight=1),
      list(nrounds=300L,max_depth=4L,eta=.05,subsample=.8,colsample_bytree=.8,min_child_weight=1),
      list(nrounds=250L,max_depth=6L,eta=.03,subsample=.9,colsample_bytree=.8,min_child_weight=2)
    ),
    svm = list(
      list(kernel="radial",cost=.5,gamma=NULL,gamma_factor=.5),
      list(kernel="radial",cost=1,gamma=NULL,gamma_factor=1),
      list(kernel="radial",cost=4,gamma=NULL,gamma_factor=2)
    ),
    knn = lapply(c(3L,5L,7L,11L,15L), function(k) list(k=k)),
    multinom = lapply(c(0,.001,.01,.1), function(d) list(decay=d,maxit=500)),
    elastic = lapply(c(.25,.5,.75), function(a) list(alpha=a,lambda_rule="lambda.1se")),
    list(.r4vn_machine_default_params(method, task))
  )
  if (mode == "random" && length(grid) > tune_n) {
    .r4vn_set_seed_if(seed); grid <- grid[sample(seq_along(grid), tune_n)]
  }
  grid
}

.r4vn_machine_resolve_params <- function(params, p, task) {
  z <- params
  if (!is.null(z$mtry_factor)) {
    z$mtry <- max(1L, min(p, round(sqrt(p) * z$mtry_factor))); z$mtry_factor <- NULL
  }
  if (!is.null(z$gamma_factor)) {
    z$gamma <- z$gamma_factor / max(1,p); z$gamma_factor <- NULL
  }
  z
}

.r4vn_machine_fit_engine <- function(x, y, task, method, params, weights = NULL, seed = NA_integer_, standardize = "auto") {
  p <- ncol(x)
  params <- .r4vn_machine_resolve_params(params, p, task)
  sm <- .r4vn_machine_scale_method(standardize, method)
  scaler <- .r4vn_machine_scale_fit(x, sm)
  xx <- .r4vn_machine_scale_apply(x, scaler)
  .r4vn_set_seed_if(seed)

  if (method == "logistic") {
    df <- data.frame(.y=as.numeric(y), xx, check.names=FALSE)
    fit <- stats::glm(.y ~ ., data=df, family=stats::binomial(), weights=weights)
  } else if (method == "linear") {
    df <- data.frame(.y=as.numeric(y), xx, check.names=FALSE)
    fit <- stats::lm(.y ~ ., data=df, weights=weights)
  } else if (method == "multinom") {
    if (!requireNamespace("nnet", quietly=TRUE)) stop("Package `nnet` is required for multinomial regression.")
    df <- data.frame(.y=factor(y), xx, check.names=FALSE)
    fit <- nnet::multinom(.y ~ ., data=df, weights=weights, trace=FALSE,
                         decay=params$decay %||% 0, maxit=params$maxit %||% 500)
  } else if (method %in% c("lasso","ridge","elastic")) {
    if (!requireNamespace("glmnet", quietly=TRUE)) stop("Package `glmnet` is required for penalized models.")
    fam <- if (task=="binary") "binomial" else if (task=="multiclass") "multinomial" else "gaussian"
    yy <- if (task=="multiclass") factor(y) else as.numeric(y)
    alpha <- params$alpha %||% switch(method,lasso=1,ridge=0,elastic=.5)
    fit <- glmnet::cv.glmnet(xx, yy, family=fam, alpha=alpha, weights=weights,
                             nfolds=max(3L,min(10L,nrow(xx))), standardize=FALSE)
  } else if (method == "tree") {
    if (!requireNamespace("rpart", quietly=TRUE)) stop("Package `rpart` is required for decision trees.")
    yy <- if (task=="regression") as.numeric(y) else factor(y)
    df <- data.frame(.y=yy, xx, check.names=FALSE)
    control <- rpart::rpart.control(cp=params$cp %||% .01,
                                   minsplit=params$minsplit %||% 20L,
                                   maxdepth=params$maxdepth %||% 30L,
                                   xval=0)
    fit <- rpart::rpart(.y ~ ., data=df, weights=weights,
                        method=if(task=="regression") "anova" else "class", control=control)
  } else if (method == "rf") {
    if (!requireNamespace("ranger", quietly=TRUE)) stop("Package `ranger` is required for random forest.")
    yy <- if (task=="regression") as.numeric(y) else factor(y)
    df <- data.frame(.y=yy, xx, check.names=FALSE)
    fit <- ranger::ranger(.y ~ ., data=df,
                          probability=task != "regression",
                          num.trees=params$num.trees %||% 500L,
                          mtry=min(p, params$mtry %||% max(1L,floor(sqrt(p)))),
                          min.node.size=params$min.node.size %||% if(task=="regression") 5L else 1L,
                          case.weights=weights,
                          importance="permutation",
                          seed=if (is.finite(seed)) seed else NULL)
  } else if (method == "xgb") {
    if (!requireNamespace("xgboost", quietly=TRUE)) stop("Package `xgboost` is required for XGBoost.")
    label <- if (task=="multiclass") as.integer(factor(y))-1L else as.numeric(y)
    dm <- if (is.null(weights)) xgboost::xgb.DMatrix(data=xx, label=label) else
      xgboost::xgb.DMatrix(data=xx, label=label, weight=weights)
    obj <- if(task=="binary") "binary:logistic" else if(task=="multiclass") "multi:softprob" else "reg:squarederror"
    xp <- list(objective=obj, max_depth=params$max_depth %||% 4L, eta=params$eta %||% .05,
               subsample=params$subsample %||% .8, colsample_bytree=params$colsample_bytree %||% .8,
               min_child_weight=params$min_child_weight %||% 1, verbosity=0)
    if(task=="multiclass") xp$num_class <- length(levels(factor(y)))
    fit <- xgboost::xgb.train(params=xp,data=dm,nrounds=params$nrounds %||% 250L,verbose=0)
  } else if (method == "svm") {
    if (!requireNamespace("e1071", quietly=TRUE)) stop("Package `e1071` is required for SVM.")
    yy <- if(task=="regression") as.numeric(y) else factor(y)
    df <- data.frame(.y=yy,xx,check.names=FALSE)
    class.weights <- NULL
    if (task != "regression" && !is.null(weights)) {
      tmp <- tapply(weights, yy, mean); class.weights <- as.numeric(tmp); names(class.weights) <- names(tmp)
    }
    fit <- e1071::svm(.y ~ .,data=df,type=if(task=="regression") "eps-regression" else "C-classification",
                      kernel=params$kernel %||% "radial",cost=params$cost %||% 1,
                      gamma=params$gamma %||% 1/max(1,p),probability=task!="regression",
                      class.weights=class.weights,scale=FALSE)
  } else if (method == "knn") {
    fit <- list(x=xx,y=y,k=as.integer(params$k %||% 7L),task=task)
    class(fit) <- "r4vn_machine_knn"
  } else if (method == "naive") {
    if (!requireNamespace("e1071", quietly=TRUE)) stop("Package `e1071` is required for naive Bayes.")
    df <- data.frame(.y=factor(y),xx,check.names=FALSE)
    fit <- e1071::naiveBayes(.y ~ .,data=df,laplace=params$laplace %||% 0)
  } else stop("Unsupported method: ",method,call.=FALSE)

  list(model=fit, method=method, task=task, scaler=scaler, params=params,
       x_cols=colnames(x), y_levels=if(task=="multiclass") levels(factor(y)) else if(task=="binary") c("0","1") else NULL)
}

`%||%` <- function(x,y) if(is.null(x)) y else x

.r4vn_machine_knn_predict <- function(fit, newx) {
  tr <- fit$x; y <- fit$y; k <- min(fit$k,nrow(tr))
  if(fit$task=="regression") {
    return(vapply(seq_len(nrow(newx)),function(i){d<-rowSums((tr-matrix(newx[i,],nrow(tr),ncol(tr),byrow=TRUE))^2);mean(as.numeric(y)[order(d)[seq_len(k)]])},numeric(1)))
  }
  lev <- levels(factor(y)); probs <- matrix(0,nrow(newx),length(lev),dimnames=list(NULL,lev))
  for(i in seq_len(nrow(newx))){
    d<-rowSums((tr-matrix(newx[i,],nrow(tr),ncol(tr),byrow=TRUE))^2); nn<-order(d)[seq_len(k)]; tb<-table(factor(y[nn],levels=lev)); probs[i,]<-as.numeric(tb)/k
  }
  probs
}

.r4vn_machine_predict_engine <- function(fit, x) {
  xx <- .r4vn_machine_scale_apply(x, fit$scaler)
  method <- fit$method; task <- fit$task; model <- fit$model
  if(method=="logistic") return(list(prob=as.numeric(stats::predict(model,newdata=data.frame(xx,check.names=FALSE),type="response"))))
  if(method=="linear") return(list(value=as.numeric(stats::predict(model,newdata=data.frame(xx,check.names=FALSE)))))
  if(method=="multinom") {
    p<-stats::predict(model,newdata=data.frame(xx,check.names=FALSE),type="probs"); if(is.vector(p)) p<-matrix(p,nrow=1)
    if(is.null(colnames(p)) && ncol(p)==length(fit$y_levels)) colnames(p)<-fit$y_levels
    return(list(prob=p,class=colnames(p)[max.col(p)]))
  }
  if(method %in% c("lasso","ridge","elastic")) {
    s <- fit$params$lambda_rule %||% "lambda.1se"
    pp <- stats::predict(model,newx=xx,s=s,type="response")
    if(task=="regression") return(list(value=as.numeric(pp)))
    if(task=="binary") return(list(prob=as.numeric(pp)))
    ar <- as.array(pp); if(length(dim(ar))==3L) ar<-ar[,,1L,drop=FALSE][,,1L]; if(is.vector(ar)) ar<-matrix(ar,nrow=nrow(xx))
    return(list(prob=ar,class=colnames(ar)[max.col(ar)]))
  }
  if(method=="tree") {
    if(task=="regression") return(list(value=as.numeric(stats::predict(model,newdata=data.frame(xx,check.names=FALSE)))))
    p<-stats::predict(model,newdata=data.frame(xx,check.names=FALSE),type="prob"); if(task=="binary") return(list(prob=as.numeric(p[,ncol(p)]))); return(list(prob=p,class=colnames(p)[max.col(p)]))
  }
  if(method=="rf") {
    p<-predict(model,data=data.frame(xx,check.names=FALSE))$predictions
    if(task=="regression") return(list(value=as.numeric(p)))
    if(task=="binary") return(list(prob=as.numeric(p[,ncol(p)])))
    return(list(prob=p,class=colnames(p)[max.col(p)]))
  }
  if(method=="xgb") {
    p<-predict(model,xgboost::xgb.DMatrix(xx))
    if(task=="regression") return(list(value=as.numeric(p)))
    if(task=="binary") return(list(prob=as.numeric(p)))
    nc<-length(fit$y_levels); if(is.null(dim(p))) p<-matrix(p,ncol=nc,byrow=TRUE); colnames(p)<-fit$y_levels
    return(list(prob=p,class=colnames(p)[max.col(p)]))
  }
  if(method=="svm") {
    if(task=="regression") return(list(value=as.numeric(stats::predict(model,newdata=data.frame(xx,check.names=FALSE)))))
    pr<-stats::predict(model,newdata=data.frame(xx,check.names=FALSE),probability=TRUE); p<-attr(pr,"probabilities")
    if(is.null(p)) {lev<-levels(pr);p<-matrix(0,length(pr),length(lev),dimnames=list(NULL,lev));p[cbind(seq_along(pr),match(pr,lev))]<-1}
    if(task=="binary") return(list(prob=as.numeric(p[,ncol(p)])))
    return(list(prob=p,class=colnames(p)[max.col(p)]))
  }
  if(method=="knn") {
    p<-.r4vn_machine_knn_predict(model,xx); if(task=="regression") return(list(value=p)); if(task=="binary") return(list(prob=as.numeric(p[,ncol(p)]))); return(list(prob=p,class=colnames(p)[max.col(p)]))
  }
  if(method=="naive") {
    p<-stats::predict(model,newdata=data.frame(xx,check.names=FALSE),type="raw"); if(task=="binary") return(list(prob=as.numeric(p[,ncol(p)]))); return(list(prob=p,class=colnames(p)[max.col(p)]))
  }
  stop("Prediction is not implemented for method: ",method,call.=FALSE)
}

.r4vn_machine_fit_pipeline <- function(raw,outcome_name,predictor_names,task,outcome_info,
                                       method,params,preprocess,missing,missing_max,encode,
                                       standardize,transform,outlier,corr,feature,degree,reduce,variance,select,nfeatures,
                                       balance,balance_target,neighbors,selected_override=NULL,seed=NA_integer_) {
  recipe <- .r4vn_machine_recipe_fit(raw,predictor_names,preprocess,missing,missing_max,encode,transform,outlier,corr,feature,degree,reduce,variance)
  ap <- .r4vn_machine_recipe_apply(recipe,raw)
  yy <- raw[[outcome_name]][ap$rows]
  fs <- .r4vn_machine_feature_select(ap$x,yy,task,select,nfeatures,seed+1L)
  selected <- fs$selected
  if(!is.null(selected_override)) {
    ov <- intersect(selected_override,colnames(ap$x)); if(length(ov)) selected<-ov
  }
  xx <- ap$x[,selected,drop=FALSE]
  balance_engine <- balance
  if (task == "binary" && identical(balance, "weight") && method %in% c("knn", "naive")) {
    balance_engine <- "up"
  }
  bal <- if(task=="binary") .r4vn_machine_balance_apply(xx,yy,balance_engine,balance_target,neighbors,seed+2L) else list(x=xx,y=yy,weights=NULL,method="none")
  if (!identical(balance_engine, balance)) bal$method <- paste0(bal$method, " (fallback for ", method, ")")
  eng <- .r4vn_machine_fit_engine(bal$x,bal$y,task,method,params,bal$weights,seed+3L,standardize)
  list(recipe=recipe,selected=selected,selection_detail=fs$detail,selection_note=fs$note %||% NULL,
       balance=bal$method,engine=eng,method=method,task=task,outcome=outcome_name,
       predictors=predictor_names,outcome_info=outcome_info,standardize=standardize)
}

.r4vn_machine_predict_pipeline <- function(fit,raw,outcome_name=NULL) {
  if(!is.data.frame(raw)) stop("Prediction data must be a data frame.",call.=FALSE)
  miss <- setdiff(fit$recipe$keep,names(raw))
  if(length(miss)) stop("Prediction data are missing required predictors: ",paste(miss,collapse=", "),call.=FALSE)
  ap <- .r4vn_machine_recipe_apply(fit$recipe,raw)
  x <- ap$x[,fit$selected,drop=FALSE]
  pr <- .r4vn_machine_predict_engine(fit$engine,x)
  pr$rows <- ap$rows; pr$unseen <- ap$unseen
  pr
}

.r4vn_machine_cv_one <- function(raw,outcome_name,predictor_names,task,outcome_info,method,params,
                                 preprocess,missing,missing_max,encode,standardize,transform,outlier,corr,
                                 feature,degree,reduce,variance,select,nfeatures,balance,balance_target,neighbors,folds,repeats,seed) {
  sp <- .r4vn_machine_make_folds(raw[[outcome_name]], folds, repeats, task, seed)
  rows <- list()
  errors <- character()
  ii <- 0L

  remember_error <- function(stage, s, e) {
    msg <- paste0(
      stage, " failed in repeat ", s$rep_id, ", fold ", s$fold,
      ": ", conditionMessage(e)
    )
    errors <<- c(errors, msg)
    invisible(NULL)
  }

  for (s in sp) {
    tr <- raw[s$train, , drop = FALSE]
    va <- raw[s$val, , drop = FALSE]

    fit <- tryCatch(
      .r4vn_machine_fit_pipeline(
        tr, outcome_name, predictor_names, task, outcome_info, method, params,
        preprocess, missing, missing_max, encode, standardize, transform, outlier, corr,
        feature, degree, reduce, variance, select, nfeatures, balance, balance_target, neighbors,
        seed = seed + s$rep_id * 100L + s$fold
      ),
      error = function(e) {
        remember_error("Model fitting", s, e)
        NULL
      }
    )

    if (is.null(fit)) next

    pr <- tryCatch(
      .r4vn_machine_predict_pipeline(fit, va, outcome_name),
      error = function(e) {
        remember_error("Prediction", s, e)
        NULL
      }
    )

    if (is.null(pr)) next
    if (!length(pr$rows)) {
      errors <- c(
        errors,
        paste0("Prediction returned zero usable rows in repeat ", s$rep_id, ", fold ", s$fold, ".")
      )
      next
    }

    truth <- va[[outcome_name]][pr$rows]
    ii <- ii + 1L

    if (task == "binary") {
      if (is.null(pr$prob) || length(pr$prob) != length(truth)) {
        errors <- c(errors, paste0("Binary prediction shape mismatch in repeat ", s$rep_id, ", fold ", s$fold, "."))
        ii <- ii - 1L
        next
      }
      rows[[ii]] <- data.frame(
        truth = as.numeric(truth),
        prob = as.numeric(pr$prob),
        row_id = s$val[pr$rows],
        rep_id = s$rep_id,
        fold = s$fold
      )
    } else if (task == "regression") {
      if (is.null(pr$value) || length(pr$value) != length(truth)) {
        errors <- c(errors, paste0("Regression prediction shape mismatch in repeat ", s$rep_id, ", fold ", s$fold, "."))
        ii <- ii - 1L
        next
      }
      rows[[ii]] <- data.frame(
        truth = as.numeric(truth),
        value = as.numeric(pr$value),
        row_id = s$val[pr$rows],
        rep_id = s$rep_id,
        fold = s$fold
      )
    } else {
      if (is.null(pr$prob) || !is.matrix(pr$prob) || nrow(pr$prob) != length(truth)) {
        errors <- c(errors, paste0("Multiclass prediction shape mismatch in repeat ", s$rep_id, ", fold ", s$fold, "."))
        ii <- ii - 1L
        next
      }
      p <- as.data.frame(pr$prob, check.names = FALSE)
      names(p) <- paste0("prob__", names(p))
      rows[[ii]] <- cbind(
        data.frame(
          truth = as.character(truth),
          row_id = s$val[pr$rows],
          rep_id = s$rep_id,
          fold = s$fold
        ),
        p
      )
    }
  }

  if (!length(rows)) {
    detail <- if (length(errors)) {
      paste0(" First fold error: ", errors[[1L]])
    } else {
      ""
    }
    stop(
      paste0("Cross-validation produced no valid predictions.", detail),
      call. = FALSE
    )
  }

  out <- do.call(rbind, rows)
  row.names(out) <- NULL
  attr(out, "cv_errors") <- unique(errors)
  out
}

.r4vn_machine_oof_unique_binary <- function(cv) {
  if (!all(c("truth", "prob") %in% names(cv))) stop("Binary OOF predictions are incomplete.", call. = FALSE)
  if (!"row_id" %in% names(cv) || anyNA(cv$row_id)) return(cv[, c("truth", "prob"), drop = FALSE])
  ids <- unique(cv$row_id)
  data.frame(
    truth = vapply(ids, function(i) cv$truth[match(i, cv$row_id)], numeric(1)),
    prob = vapply(ids, function(i) mean(cv$prob[cv$row_id == i], na.rm = TRUE), numeric(1)),
    row_id = ids,
    stringsAsFactors = FALSE
  )
}

.r4vn_machine_metric_from_cv <- function(cv,task,metric) {
  if(task=="binary") return(.r4vn_machine_binary_metric_value(cv$truth,cv$prob,metric,.5))
  if(task=="regression") return(.r4vn_machine_reg_metric_value(cv$truth,cv$value,metric))
  pcols<-grep("^prob__",names(cv),value=TRUE); pp<-as.matrix(cv[pcols]); colnames(pp)<-sub("^prob__","",pcols)
  .r4vn_machine_multi_metric_value(factor(cv$truth,levels=colnames(pp)),pp,metric)
}

.r4vn_machine_tune <- function(method,grid,raw,outcome_name,predictor_names,task,outcome_info,
                               preprocess,missing,missing_max,encode,standardize,transform,outlier,corr,
                               feature,degree,reduce,variance,select,nfeatures,balance,balance_target,neighbors,folds,repeats,nested,metric,seed) {
  vals<-rep(NA_real_,length(grid)); cvs<-vector("list",length(grid)); errs<-character()
  for(i in seq_along(grid)){
    cv<-tryCatch(.r4vn_machine_cv_one(raw,outcome_name,predictor_names,task,outcome_info,method,grid[[i]],preprocess,missing,missing_max,encode,standardize,transform,outlier,corr,feature,degree,reduce,variance,select,nfeatures,balance,balance_target,neighbors,folds,repeats,seed+i*50L),error=function(e)e)
    if(inherits(cv,"error")){errs<-c(errs,conditionMessage(cv));next}
    cvs[[i]]<-cv; vals[i]<-.r4vn_machine_metric_from_cv(cv,task,metric)
  }
  good<-which(is.finite(vals)); if(!length(good)) return(list(ok=FALSE,error=paste(unique(errs),collapse=" | ")))
  dir<-.r4vn_machine_metric_direction(metric); besti<-if(dir=="max") good[which.max(vals[good])] else good[which.min(vals[good])]
  cv<-cvs[[besti]]
  if (isTRUE(nested) && length(grid) > 1L) {
    nested_cv <- .r4vn_machine_nested_cv(raw,outcome_name,predictor_names,task,outcome_info,method,grid,
      preprocess,missing,missing_max,encode,standardize,transform,outlier,corr,feature,degree,reduce,variance,select,nfeatures,balance,balance_target,neighbors,folds,metric,seed+777L)
    if (!is.null(nested_cv)) cv <- nested_cv
  }
  # summary for the primary metric plus common secondary metrics on OOF predictions
  metrics<-if(task=="binary") c("auc","pr_auc","accuracy","balanced_accuracy","sensitivity","specificity","f1","brier","logloss") else if(task=="regression") c("rmse","mae","r2","mape") else c("accuracy","balanced_accuracy","macro_f1","weighted_f1","macro_auc","macro_pr_auc","logloss")
  summ<-data.frame(method=method,metric=metrics,estimate=vapply(metrics,function(m).r4vn_machine_metric_from_cv(cv,task,m),numeric(1)),stringsAsFactors=FALSE)
  tuning<-data.frame(method=method,candidate=seq_along(grid),score=vals,selected=seq_along(grid)==besti,stringsAsFactors=FALSE)
  tuning$params<-vapply(grid,function(z)paste(paste(names(z),unlist(z),sep="="),collapse="; "),character(1))
  list(ok=TRUE,params=grid[[besti]],cv=cv,cv_summary=summ,tuning=tuning,nested=nested)
}

# ============================================================================
# Performance metrics and confidence intervals
# ============================================================================

.r4vn_machine_div <- function(a,b) ifelse(b==0,NA_real_,a/b)

.r4vn_machine_wilson <- function(x,n,level=.95) {
  if(!is.finite(x)||!is.finite(n)||n<=0) return(c(NA_real_,NA_real_))
  z<-stats::qnorm(1-(1-level)/2); ph<-x/n; den<-1+z^2/n
  ctr<-(ph+z^2/(2*n))/den; half<-z*sqrt(ph*(1-ph)/n+z^2/(4*n^2))/den
  pmax(0,pmin(1,c(ctr-half,ctr+half)))
}

.r4vn_machine_auc <- function(y,p) {
  ok<-is.finite(p)&!is.na(y); y<-as.numeric(y[ok]);p<-p[ok]
  n1<-sum(y==1);n0<-sum(y==0);if(n1==0||n0==0)return(NA_real_)
  r<-rank(p,ties.method="average"); (sum(r[y==1])-n1*(n1+1)/2)/(n1*n0)
}

.r4vn_machine_pr_auc <- function(y,p) {
  ok<-is.finite(p)&!is.na(y);y<-as.numeric(y[ok]);p<-p[ok];npos<-sum(y==1);if(npos==0)return(NA_real_)
  o<-order(p,decreasing=TRUE); yy<-y[o];tp<-cumsum(yy==1);fp<-cumsum(yy==0)
  rec<-tp/npos;prec<-tp/pmax(tp+fp,1)
  rec<-c(0,rec);prec<-c(1,prec)
  sum(diff(rec)*(head(prec,-1)+tail(prec,-1))/2,na.rm=TRUE)
}

.r4vn_machine_binary_metrics <- function(y,p,threshold=.5) {
  ok<-!is.na(y)&is.finite(p);y<-as.numeric(y[ok]);p<-p[ok];cl<-as.integer(p>=threshold)
  tp<-sum(cl==1&y==1);fp<-sum(cl==1&y==0);fn<-sum(cl==0&y==1);tn<-sum(cl==0&y==0);n<-length(y)
  sens<-.r4vn_machine_div(tp,tp+fn);spec<-.r4vn_machine_div(tn,tn+fp);ppv<-.r4vn_machine_div(tp,tp+fp);npv<-.r4vn_machine_div(tn,tn+fn)
  acc<-.r4vn_machine_div(tp+tn,n);f1<-.r4vn_machine_div(2*tp,2*tp+fp+fn);bal<-mean(c(sens,spec),na.rm=TRUE)
  mccden<-sqrt((tp+fp)*(tp+fn)*(tn+fp)*(tn+fn));mcc<-.r4vn_machine_div(tp*tn-fp*fn,mccden)
  po<-acc; pe<-.r4vn_machine_div((tp+fp)*(tp+fn)+(fn+tn)*(fp+tn),n^2);kappa<-.r4vn_machine_div(po-pe,1-pe)
  prev<-mean(y==1); brier<-mean((p-y)^2); eps<-1e-15;pp<-pmin(pmax(p,eps),1-eps);logloss<- -mean(y*log(pp)+(1-y)*log(1-pp))
  c(n=n,tp=tp,fp=fp,fn=fn,tn=tn,prevalence=prev,auc=.r4vn_machine_auc(y,p),pr_auc=.r4vn_machine_pr_auc(y,p),
    accuracy=acc,balanced_accuracy=bal,sensitivity=sens,specificity=spec,ppv=ppv,npv=npv,f1=f1,mcc=mcc,kappa=kappa,brier=brier,logloss=logloss)
}

.r4vn_machine_binary_metric_value <- function(y,p,metric,threshold=.5) {
  z<-.r4vn_machine_binary_metrics(y,p,threshold); nm<-tolower(metric)
  if(nm=="sens")nm<-"sensitivity";if(nm=="spec")nm<-"specificity";if(nm=="roc_auc")nm<-"auc"
  unname(z[nm] %||% NA_real_)
}

.r4vn_machine_reg_metrics <- function(y,p) {
  ok<-is.finite(y)&is.finite(p);y<-as.numeric(y[ok]);p<-as.numeric(p[ok]);n<-length(y);if(!n)return(c(n=0,rmse=NA,mae=NA,r2=NA,mape=NA))
  rmse<-sqrt(mean((y-p)^2));mae<-mean(abs(y-p));sst<-sum((y-mean(y))^2);r2<-if(sst>0)1-sum((y-p)^2)/sst else NA_real_
  nz<-abs(y)>.Machine$double.eps^.5;mape<-if(any(nz))mean(abs((y[nz]-p[nz])/y[nz])) else NA_real_
  c(n=n,rmse=rmse,mae=mae,r2=r2,mape=mape)
}
.r4vn_machine_reg_metric_value <- function(y,p,metric) unname(.r4vn_machine_reg_metrics(y,p)[tolower(metric)] %||% NA_real_)

.r4vn_machine_multi_metrics <- function(y,p) {
  if(!is.matrix(p))p<-as.matrix(p);lev<-colnames(p);if(is.null(lev))lev<-levels(factor(y));y<-factor(y,levels=lev)
  pred<-factor(lev[max.col(p,ties.method="first")],levels=lev);cm<-table(y,pred);n<-sum(cm);acc<-sum(diag(cm))/n
  sens<-diag(cm)/pmax(rowSums(cm),1);prec<-diag(cm)/pmax(colSums(cm),1);f1<-2*prec*sens/pmax(prec+sens,.Machine$double.eps)
  macro_f1<-mean(f1,na.rm=TRUE);weighted_f1<-sum(f1*rowSums(cm))/n;bal<-mean(sens,na.rm=TRUE)
  aucs<-praucs<-rep(NA_real_,length(lev))
  for(j in seq_along(lev)){
    yy<-as.integer(y==lev[j]);aucs[j]<-.r4vn_machine_auc(yy,p[,j]);praucs[j]<-.r4vn_machine_pr_auc(yy,p[,j])
  }
  macro_auc<-mean(aucs,na.rm=TRUE);macro_pr_auc<-mean(praucs,na.rm=TRUE)
  eps<-1e-15;p<-pmin(pmax(p,eps),1-eps);idx<-cbind(seq_len(nrow(p)),match(y,lev));ll<- -mean(log(p[idx]))
  c(n=n,accuracy=acc,balanced_accuracy=bal,macro_f1=macro_f1,weighted_f1=weighted_f1,macro_auc=macro_auc,macro_pr_auc=macro_pr_auc,logloss=ll)
}
.r4vn_machine_multi_metric_value <- function(y,p,metric) unname(.r4vn_machine_multi_metrics(y,p)[tolower(metric)] %||% NA_real_)

.r4vn_machine_boot_binary <- function(y,p,threshold,metrics,boot,level,seed,stratified=TRUE) {
  .r4vn_set_seed_if(seed);n<-length(y);mat<-matrix(NA_real_,boot,length(metrics),dimnames=list(NULL,metrics));idx0<-which(y==0);idx1<-which(y==1)
  for(b in seq_len(boot)){
    if(stratified&&length(idx0)&&length(idx1)) idx<-c(sample(idx0,length(idx0),TRUE),sample(idx1,length(idx1),TRUE)) else idx<-sample(seq_len(n),n,TRUE)
    zz<-.r4vn_machine_binary_metrics(y[idx],p[idx],threshold);mat[b,]<-zz[metrics]
  }
  alpha<-(1-level)/2;t(apply(mat,2,stats::quantile,probs=c(alpha,1-alpha),na.rm=TRUE,names=FALSE))
}

.r4vn_machine_boot_reg <- function(y,p,metrics,boot,level,seed) {
  .r4vn_set_seed_if(seed);n<-length(y);mat<-matrix(NA_real_,boot,length(metrics),dimnames=list(NULL,metrics))
  for(b in seq_len(boot)){idx<-sample(seq_len(n),n,TRUE);mat[b,]<-.r4vn_machine_reg_metrics(y[idx],p[idx])[metrics]}
  alpha<-(1-level)/2;t(apply(mat,2,stats::quantile,probs=c(alpha,1-alpha),na.rm=TRUE,names=FALSE))
}

.r4vn_machine_boot_multi <- function(y,p,metrics,boot,level,seed) {
  .r4vn_set_seed_if(seed);n<-length(y);mat<-matrix(NA_real_,boot,length(metrics),dimnames=list(NULL,metrics))
  for(b in seq_len(boot)){idx<-sample(seq_len(n),n,TRUE);mat[b,]<-.r4vn_machine_multi_metrics(y[idx],p[idx,,drop=FALSE])[metrics]}
  alpha<-(1-level)/2;t(apply(mat,2,stats::quantile,probs=c(alpha,1-alpha),na.rm=TRUE,names=FALSE))
}

.r4vn_machine_multiclass_class_performance <- function(truth,pred,ci=TRUE,level=.95,boot=1000L,seed=NA_integer_) {
  p <- as.matrix(pred$prob); lev <- colnames(p)
  if(is.null(lev) || !length(lev)) return(data.frame())
  rows <- pred$rows %||% seq_along(truth)
  # `truth` may already have been aligned to the prediction rows by the caller.
  if (length(truth) == nrow(p)) rows <- seq_len(nrow(p))
  y <- factor(truth[rows],levels=lev); pcl <- factor(lev[max.col(p,ties.method="first")],levels=lev)
  out <- list(); ii <- 0L; .r4vn_set_seed_if(seed)
  for(j in seq_along(lev)) {
    cl <- lev[j]; yb <- as.integer(y==cl); cb <- as.integer(pcl==cl); ps <- p[,j]
    tp<-sum(cb==1&yb==1);fp<-sum(cb==1&yb==0);fn<-sum(cb==0&yb==1);tn<-sum(cb==0&yb==0);n<-length(yb)
    est <- c(
      auc=.r4vn_machine_auc(yb,ps), pr_auc=.r4vn_machine_pr_auc(yb,ps),
      sensitivity=.r4vn_machine_div(tp,tp+fn), specificity=.r4vn_machine_div(tn,tn+fp),
      ppv=.r4vn_machine_div(tp,tp+fp), npv=.r4vn_machine_div(tn,tn+fn),
      accuracy=.r4vn_machine_div(tp+tn,n),
      f1=.r4vn_machine_div(2*tp,2*tp+fp+fn), prevalence=mean(yb==1)
    )
    tab <- data.frame(class=cl,metric=names(est),estimate=as.numeric(est),lower=NA_real_,upper=NA_real_,ci_method=NA_character_,stringsAsFactors=FALSE)
    if(isTRUE(ci)) {
      den <- list(sensitivity=c(tp,tp+fn),specificity=c(tn,tn+fp),ppv=c(tp,tp+fp),npv=c(tn,tn+fn),accuracy=c(tp+tn,n),prevalence=c(tp+fn,n))
      for(nm in names(den)){cc<-.r4vn_machine_wilson(den[[nm]][1L],den[[nm]][2L],level);k<-which(tab$metric==nm);tab[k,c("lower","upper")]<-cc;tab$ci_method[k]<-"Wilson"}
      ac <- .r4vn_machine_auc_ci(yb,ps,level,boot,seed+j*101L);k<-which(tab$metric=="auc");tab[k,c("lower","upper")]<-ac$ci;tab$ci_method[k]<-ac$method
      # PR-AUC and F1 are derived measures; bootstrap complete subject rows.
      bpr<-bf1<-rep(NA_real_,boot)
      idx0<-which(yb==0);idx1<-which(yb==1)
      for(b in seq_len(boot)){
        ix<-if(length(idx0)&&length(idx1))c(sample(idx0,length(idx0),TRUE),sample(idx1,length(idx1),TRUE))else sample(seq_len(n),n,TRUE)
        bpr[b]<-.r4vn_machine_pr_auc(yb[ix],ps[ix])
        tpb<-sum(cb[ix]==1&yb[ix]==1);fpb<-sum(cb[ix]==1&yb[ix]==0);fnb<-sum(cb[ix]==0&yb[ix]==1)
        bf1[b]<-.r4vn_machine_div(2*tpb,2*tpb+fpb+fnb)
      }
      aq<-c((1-level)/2,1-(1-level)/2)
      qfun<-function(v){v<-v[is.finite(v)];if(length(v)<30L)c(NA_real_,NA_real_)else stats::quantile(v,aq,na.rm=TRUE,names=FALSE)}
      qq<-qfun(bpr);k<-which(tab$metric=="pr_auc");tab[k,c("lower","upper")]<-qq;tab$ci_method[k]<-"Stratified bootstrap"
      qq<-qfun(bf1);k<-which(tab$metric=="f1");tab[k,c("lower","upper")]<-qq;tab$ci_method[k]<-"Stratified bootstrap"
    }
    ii<-ii+1L;out[[ii]]<-tab
  }
  .r4vn_machine_bind(out)
}

.r4vn_machine_auc_ci <- function(y,p,level,boot,seed) {
  ok<-!is.na(y)&is.finite(p);y<-as.numeric(y[ok]);p<-p[ok]
  if(length(unique(y))<2L)return(list(ci=c(NA_real_,NA_real_),method=NA_character_))
  if(requireNamespace("pROC",quietly=TRUE)){
    rr<-tryCatch(pROC::roc(y,p,levels=c(0,1),direction="<",quiet=TRUE),error=function(e)NULL)
    if(!is.null(rr)){
      cc<-tryCatch(pROC::ci.auc(rr,conf.level=level,method="delong"),error=function(e)NULL)
      if(!is.null(cc))return(list(ci=c(as.numeric(cc[1L]),as.numeric(cc[3L])),method="DeLong"))
    }
  }
  bb<-.r4vn_machine_boot_binary(y,p,.5,"auc",boot,level,seed,TRUE)
  list(ci=as.numeric(bb[1,]),method="Stratified bootstrap")
}

.r4vn_machine_evaluate <- function(truth,pred,task,method,threshold=NULL,ci=TRUE,ci_level=.95,boot=1000,seed=NA_integer_,dataset="Test") {
  if(task=="binary"){
    y<-as.numeric(truth[pred$rows %||% seq_along(truth)]);p<-pred$prob;th<-threshold %||% .5;z<-.r4vn_machine_binary_metrics(y,p,th)
    metrics<-c("auc","pr_auc","accuracy","balanced_accuracy","sensitivity","specificity","ppv","npv","f1","mcc","kappa","brier","logloss","prevalence")
    out<-data.frame(dataset=dataset,method=method,metric=metrics,estimate=as.numeric(z[metrics]),lower=NA_real_,upper=NA_real_,ci_method=NA_character_,stringsAsFactors=FALSE)
    if(isTRUE(ci)){
      # Wilson intervals where the metric is a direct binomial proportion.
      denoms<-list(accuracy=c(z["tp"]+z["tn"],z["n"]),sensitivity=c(z["tp"],z["tp"]+z["fn"]),specificity=c(z["tn"],z["tn"]+z["fp"]),
                   ppv=c(z["tp"],z["tp"]+z["fp"]),npv=c(z["tn"],z["tn"]+z["fn"]),prevalence=c(z["tp"]+z["fn"],z["n"]))
      for(nm in names(denoms)){cc<-.r4vn_machine_wilson(denoms[[nm]][1],denoms[[nm]][2],ci_level);ii<-which(out$metric==nm);out[ii,c("lower","upper")]<-cc;out$ci_method[ii]<-"Wilson"}
      ac<-.r4vn_machine_auc_ci(y,p,ci_level,boot,seed+1L);ii<-which(out$metric=="auc");out[ii,c("lower","upper")]<-ac$ci;out$ci_method[ii]<-ac$method
      derived<-setdiff(metrics,c(names(denoms),"auc"));bb<-.r4vn_machine_boot_binary(y,p,th,derived,boot,ci_level,seed+2L,TRUE)
      for(nm in derived){ii<-which(out$metric==nm);out[ii,c("lower","upper")]<-bb[nm,];out$ci_method[ii]<-"Bootstrap"}
    }
    return(list(performance=out,metrics=z))
  }
  if(task=="regression"){
    y<-as.numeric(truth[pred$rows %||% seq_along(truth)]);p<-pred$value;metrics<-c("rmse","mae","r2","mape");z<-.r4vn_machine_reg_metrics(y,p)
    out<-data.frame(dataset=dataset,method=method,metric=metrics,estimate=as.numeric(z[metrics]),lower=NA_real_,upper=NA_real_,ci_method=NA_character_)
    if(isTRUE(ci)){bb<-.r4vn_machine_boot_reg(y,p,metrics,boot,ci_level,seed);out[,c("lower","upper")]<-bb[metrics,,drop=FALSE];out$ci_method<-"Bootstrap"}
    return(list(performance=out,metrics=z))
  }
  rows<-pred$rows %||% seq_along(truth);y<-factor(truth[rows],levels=colnames(pred$prob));p<-pred$prob;metrics<-c("accuracy","balanced_accuracy","macro_f1","weighted_f1","macro_auc","macro_pr_auc","logloss");z<-.r4vn_machine_multi_metrics(y,p)
  out<-data.frame(dataset=dataset,method=method,metric=metrics,estimate=as.numeric(z[metrics]),lower=NA_real_,upper=NA_real_,ci_method=NA_character_)
  if(isTRUE(ci)){
    pred_class <- factor(colnames(p)[max.col(p,ties.method="first")],levels=levels(y))
    acc_ci <- .r4vn_machine_wilson(sum(pred_class==y,na.rm=TRUE),sum(!is.na(y)),ci_level)
    ia <- which(out$metric=="accuracy"); out[ia,c("lower","upper")] <- acc_ci; out$ci_method[ia] <- "Wilson"
    derived <- setdiff(metrics,"accuracy")
    bb <- .r4vn_machine_boot_multi(y,p,derived,boot,ci_level,seed)
    for(nm in derived){ii<-which(out$metric==nm);out[ii,c("lower","upper")]<-bb[nm,];out$ci_method[ii]<-"Bootstrap"}
  }
  list(performance=out,metrics=z)
}

.r4vn_machine_threshold_point <- function(y,p,method="youden",target_sens=.9,target_spec=.9) {
  cuts <- sort(unique(c(0, p[is.finite(p)], 1)))
  if (length(cuts) > 2000L) cuts <- unique(stats::quantile(cuts, probs=seq(0,1,length.out=2000), na.rm=TRUE))
  tab <- do.call(rbind, lapply(cuts, function(th) {
    m <- .r4vn_machine_binary_metrics(y,p,th)
    data.frame(threshold=th,sensitivity=m["sensitivity"],specificity=m["specificity"],
               f1=m["f1"],youden=m["sensitivity"]+m["specificity"]-1)
  }))
  if (method == "youden") i <- which.max(tab$youden)
  else if (method == "f1") i <- which.max(tab$f1)
  else if (method == "sens") {
    ok <- which(tab$sensitivity >= target_sens)
    i <- if(length(ok)) ok[which.max(tab$specificity[ok])] else which.max(tab$sensitivity)
  } else if (method == "spec") {
    ok <- which(tab$specificity >= target_spec)
    i <- if(length(ok)) ok[which.max(tab$sensitivity[ok])] else which.max(tab$specificity)
  } else stop("Unknown `threshold` method: ",method,call.=FALSE)
  list(threshold=tab$threshold[i], curve=tab)
}

.r4vn_machine_threshold <- function(y,p,threshold="auto",target_sens=.9,target_spec=.9,
                                    ci=TRUE,level=.95,boot=1000L,seed=NA_integer_) {
  ok <- !is.na(y) & is.finite(p)
  y <- as.numeric(y[ok]); p <- as.numeric(p[ok])
  if (is.numeric(threshold)) {
    th <- max(0,min(1,as.numeric(threshold[1])))
    # Retain the full sensitivity/specificity curve for plotting even when the
    # cutoff itself was fixed a priori by the user.
    curve <- .r4vn_machine_threshold_point(y,p,"youden",target_sens,target_spec)$curve
    return(list(method="fixed",threshold=th,lower=NA_real_,upper=NA_real_,ci_method="Not applicable (fixed by user)",
                training_metrics=.r4vn_machine_binary_metrics(y,p,th),curve=curve))
  }
  method <- tolower(as.character(threshold)[1]); if(method=="auto") method <- "youden"
  pt <- .r4vn_machine_threshold_point(y,p,method,target_sens,target_spec)
  th <- pt$threshold; lo <- hi <- NA_real_; ci_method <- NA_character_
  if (isTRUE(ci) && length(unique(y)) == 2L && boot >= 50L) {
    .r4vn_set_seed_if(seed); idx0 <- which(y==0); idx1 <- which(y==1); vals <- rep(NA_real_,boot)
    for (b in seq_len(boot)) {
      ix <- c(sample(idx0,length(idx0),replace=TRUE), sample(idx1,length(idx1),replace=TRUE))
      vals[b] <- tryCatch(.r4vn_machine_threshold_point(y[ix],p[ix],method,target_sens,target_spec)$threshold,
                          error=function(e) NA_real_)
    }
    vals <- vals[is.finite(vals)]
    if (length(vals) >= max(30L, ceiling(.5*boot))) {
      a <- (1-level)/2; qq <- stats::quantile(vals,c(a,1-a),na.rm=TRUE,names=FALSE)
      lo <- qq[1L]; hi <- qq[2L]; ci_method <- "Stratified bootstrap"
    }
  }
  list(method=method,threshold=th,lower=lo,upper=hi,ci_method=ci_method,
       training_metrics=.r4vn_machine_binary_metrics(y,p,th),curve=pt$curve)
}

.r4vn_machine_comparison <- function(performance,primary,digit=3) {
  if(!nrow(performance))return(data.frame())
  methods<-unique(performance$method);metrics<-unique(performance$metric);out<-data.frame(Model=methods,stringsAsFactors=FALSE)
  fmt<-function(est,lo,hi){if(is.finite(lo)&&is.finite(hi))paste0(formatC(est,format="f",digits=digit)," (",formatC(lo,format="f",digits=digit),"-",formatC(hi,format="f",digits=digit),")") else formatC(est,format="f",digits=digit)}
  for(m in metrics){nm<-.r4vn_machine_pretty_metric(m);out[[nm]]<-vapply(methods,function(md){r<-performance[performance$method==md&performance$metric==m,,drop=FALSE];if(!nrow(r))return(NA_character_);fmt(r$estimate[1],r$lower[1],r$upper[1])},character(1))}
  out
}

.r4vn_machine_choose_best <- function(cv,metric,task) {
  z<-cv[cv$metric==metric&is.finite(cv$estimate),,drop=FALSE];if(!nrow(z))return(unique(cv$method)[1L]);if(.r4vn_machine_metric_direction(metric)=="max")z$method[which.max(z$estimate)] else z$method[which.min(z$estimate)]
}

.r4vn_machine_prediction_map <- function(pr, task) {
  rows <- pr$rows
  if (is.null(rows)) {
    n <- if(task=="binary") length(pr$prob) else if(task=="regression") length(pr$value) else nrow(pr$prob)
    rows <- seq_len(n)
  }
  list(rows=rows,
       prob=if(task=="binary") pr$prob else if(task=="multiclass") pr$prob else NULL,
       value=if(task=="regression") pr$value else NULL)
}

.r4vn_machine_model_differences <- function(truth,predictions,best,task,metric,threshold,ci,level,boot,seed) {
  if(length(predictions)<2L)return(data.frame())
  bmap <- .r4vn_machine_prediction_map(predictions[[best]],task)
  methods <- setdiff(names(predictions),best); out <- list(); .r4vn_set_seed_if(seed)
  for(md in methods) {
    cmap <- .r4vn_machine_prediction_map(predictions[[md]],task)
    rows <- intersect(bmap$rows,cmap$rows); if(!length(rows))next
    ib <- match(rows,bmap$rows); ic <- match(rows,cmap$rows); yy <- truth[rows]; n <- length(rows)
    if(task=="binary") {
      p1 <- bmap$prob[ib]; p2 <- cmap$prob[ic]
      v1 <- .r4vn_machine_binary_metric_value(yy,p1,metric,threshold)
      v2 <- .r4vn_machine_binary_metric_value(yy,p2,metric,threshold)
    } else if(task=="regression") {
      p1 <- bmap$value[ib]; p2 <- cmap$value[ic]
      v1 <- .r4vn_machine_reg_metric_value(yy,p1,metric); v2 <- .r4vn_machine_reg_metric_value(yy,p2,metric)
    } else {
      p1 <- bmap$prob[ib,,drop=FALSE]; p2 <- cmap$prob[ic,,drop=FALSE]
      lev <- intersect(colnames(p1),colnames(p2)); p1 <- p1[,lev,drop=FALSE]; p2 <- p2[,lev,drop=FALSE]
      yy <- factor(yy,levels=lev)
      v1 <- .r4vn_machine_multi_metric_value(yy,p1,metric); v2 <- .r4vn_machine_multi_metric_value(yy,p2,metric)
    }
    d0 <- v1-v2; lo <- hi <- NA_real_
    if(isTRUE(ci)) {
      bb <- rep(NA_real_,boot)
      for(b in seq_len(boot)) {
        ix <- sample(seq_len(n),n,TRUE)
        if(task=="binary") bb[b] <- .r4vn_machine_binary_metric_value(yy[ix],p1[ix],metric,threshold)-.r4vn_machine_binary_metric_value(yy[ix],p2[ix],metric,threshold)
        else if(task=="regression") bb[b] <- .r4vn_machine_reg_metric_value(yy[ix],p1[ix],metric)-.r4vn_machine_reg_metric_value(yy[ix],p2[ix],metric)
        else bb[b] <- .r4vn_machine_multi_metric_value(factor(yy[ix],levels=levels(yy)),p1[ix,,drop=FALSE],metric)-.r4vn_machine_multi_metric_value(factor(yy[ix],levels=levels(yy)),p2[ix,,drop=FALSE],metric)
      }
      good <- bb[is.finite(bb)]
      if(length(good)>=30L) {qq<-stats::quantile(good,c((1-level)/2,1-(1-level)/2),na.rm=TRUE,names=FALSE);lo<-qq[1];hi<-qq[2]}
    }
    out[[md]]<-data.frame(best=best,comparison=md,metric=metric,difference=d0,lower=lo,upper=hi,
                          ci_method=if(isTRUE(ci))"Paired bootstrap" else NA_character_,stringsAsFactors=FALSE)
  }
  .r4vn_machine_bind(out)
}

.r4vn_machine_overfitting <- function(cv, performance, metric) {
  cvz <- cv[cv$metric == metric, c("method","estimate"), drop=FALSE]
  tez <- performance[performance$metric == metric, c("method","estimate"), drop=FALSE]
  if(!nrow(cvz) || !nrow(tez)) return(data.frame())
  names(cvz)[2L] <- "development"
  names(tez)[2L] <- "evaluation"
  z <- merge(cvz,tez,by="method",all=FALSE)
  if(!nrow(z)) return(z)
  dir <- .r4vn_machine_metric_direction(metric)
  z$optimism <- if(dir=="max") z$development-z$evaluation else z$evaluation-z$development
  z$flag <- ifelse(is.finite(z$optimism) & z$optimism > .05,"Possible overfitting","")
  z$metric <- metric
  z[,c("method","metric","development","evaluation","optimism","flag"),drop=FALSE]
}

.r4vn_machine_coefficients <- function(fit, labels=NULL, ci=TRUE, ci_level=.95) {
  method <- fit$method
  if(!method %in% c("logistic","linear")) return(data.frame())
  model <- fit$engine$model
  sm <- tryCatch(summary(model)$coefficients,error=function(e)NULL)
  if(is.null(sm) || !nrow(sm)) return(data.frame())
  est <- sm[,1L]; se <- sm[,2L]
  pcol <- grep("Pr\\(",colnames(sm)); pv <- if(length(pcol)) sm[,pcol[1L]] else rep(NA_real_,length(est))
  z <- stats::qnorm(1-(1-ci_level)/2); lo <- if(isTRUE(ci)) est-z*se else rep(NA_real_,length(est)); hi <- if(isTRUE(ci)) est+z*se else rep(NA_real_,length(est))
  term <- names(est); clean <- gsub("`","",term,fixed=TRUE)
  raw <- rep(NA_character_,length(clean))
  raw[clean=="(Intercept)"] <- "(Intercept)"
  mp <- fit$recipe$raw_map
  for(i in which(is.na(raw))) {
    if(clean[i] %in% names(mp)) raw[i] <- as.character(mp[[clean[i]]])
    if(is.na(raw[i]) || !nzchar(raw[i])) raw[i] <- clean[i]
  }
  lab <- vapply(raw,function(v)if(v=="(Intercept)")"Intercept" else .r4vn_machine_get_label(labels,v),character(1))
  if(method=="logistic") {
    data.frame(Predictor=lab,Term=term,OR=exp(est),Lower=exp(lo),Upper=exp(hi),p=pv,
               `CI method`=if(isTRUE(ci))"Wald" else NA_character_,check.names=FALSE,stringsAsFactors=FALSE)
  } else {
    data.frame(Predictor=lab,Term=term,Beta=est,Lower=lo,Upper=hi,p=pv,
               `CI method`=if(isTRUE(ci))"Wald" else NA_character_,check.names=FALSE,stringsAsFactors=FALSE)
  }
}

# ============================================================================
# Strategy comparison, parsimony, calibration, decision curves, explanations
# ============================================================================

.r4vn_machine_baseline_method <- function(task, methods) {
  pref <- if(task=="binary") c("logistic","tree","rf") else if(task=="regression") c("linear","tree","rf") else c("multinom","tree","rf")
  z<-intersect(pref,methods);if(length(z))z[1L] else methods[1L]
}

.r4vn_machine_compare_selection <- function(raw,outcome_name,predictor_names,task,outcome_info,methods,
                                            preprocess,missing,missing_max,encode,standardize,transform,outlier,corr,
                                            feature,degree,reduce,variance,nfeatures,balance,balance_target,neighbors,folds,metric,seed) {
  cand<-c("filter","importance","rfe");if(requireNamespace("glmnet",quietly=TRUE))cand<-c(cand,"lasso");if(requireNamespace("Boruta",quietly=TRUE))cand<-c(cand,"boruta")
  meth<-.r4vn_machine_baseline_method(task,methods);par<-.r4vn_machine_default_params(meth,task,length(predictor_names));rows<-list();notes<-character()
  for(i in seq_along(cand)){
    cv<-tryCatch(.r4vn_machine_cv_one(raw,outcome_name,predictor_names,task,outcome_info,meth,par,preprocess,missing,missing_max,encode,standardize,transform,outlier,corr,feature,degree,reduce,variance,cand[i],nfeatures,balance,balance_target,neighbors,folds,1L,seed+i*17L),error=function(e)e)
    if(inherits(cv,"error")){notes<-c(notes,paste0("Selection ",cand[i]," skipped: ",conditionMessage(cv)));next}
    rows[[cand[i]]]<-data.frame(strategy=cand[i],metric=metric,estimate=.r4vn_machine_metric_from_cv(cv,task,metric),stringsAsFactors=FALSE)
  }
  tab<-.r4vn_machine_bind(rows);if(!nrow(tab))return(list(best="filter",table=tab,notes=notes))
  best<-if(.r4vn_machine_metric_direction(metric)=="max")tab$strategy[which.max(tab$estimate)] else tab$strategy[which.min(tab$estimate)]
  list(best=best,table=tab,notes=c(notes,paste0("Feature-selection comparison selected `",best,"` using development CV only.")))
}

.r4vn_machine_compare_balance <- function(raw,outcome_name,predictor_names,task,outcome_info,methods,
                                          preprocess,missing,missing_max,encode,standardize,transform,outlier,corr,
                                          feature,degree,reduce,variance,select,nfeatures,balance_target,neighbors,folds,metric,seed) {
  cand<-c("none","weight","up","down","smote","adasyn");if(requireNamespace("ROSE",quietly=TRUE))cand<-c(cand,"rose")
  meth<-.r4vn_machine_baseline_method(task,methods);par<-.r4vn_machine_default_params(meth,task,length(predictor_names));rows<-list();notes<-character()
  for(i in seq_along(cand)){
    cv<-tryCatch(.r4vn_machine_cv_one(raw,outcome_name,predictor_names,task,outcome_info,meth,par,preprocess,missing,missing_max,encode,standardize,transform,outlier,corr,feature,degree,reduce,variance,select,nfeatures,cand[i],balance_target,neighbors,folds,1L,seed+i*19L),error=function(e)e)
    if(inherits(cv,"error")){notes<-c(notes,paste0("Balancing ",cand[i]," skipped: ",conditionMessage(cv)));next}
    rows[[cand[i]]]<-data.frame(strategy=cand[i],metric=metric,estimate=.r4vn_machine_metric_from_cv(cv,task,metric),stringsAsFactors=FALSE)
  }
  tab<-.r4vn_machine_bind(rows);if(!nrow(tab))return(list(best="none",table=tab,notes=notes))
  best<-if(.r4vn_machine_metric_direction(metric)=="max")tab$strategy[which.max(tab$estimate)] else tab$strategy[which.min(tab$estimate)]
  list(best=best,table=tab,notes=c(notes,paste0("Balancing comparison selected `",best,"` using development CV only.")))
}

.r4vn_machine_simplify <- function(raw,outcome_name,predictor_names,task,outcome_info,method,params,
                                   preprocess,missing,missing_max,encode,standardize,transform,outlier,corr,
                                   feature,degree,reduce,variance,select,nfeatures,balance,balance_target,neighbors,metric,folds,tolerance,seed) {
  # Fit preprocessing and selection once on development training data only to
  # obtain a candidate ranking. Final test data remain untouched.
  rec<-tryCatch(.r4vn_machine_recipe_fit(raw,predictor_names,preprocess,missing,missing_max,encode,transform,outlier,corr,feature,degree,reduce,variance),error=function(e)NULL)
  if(is.null(rec))return(NULL);ap<-.r4vn_machine_recipe_apply(rec,raw);y<-raw[[outcome_name]][ap$rows]
  fs<-.r4vn_machine_feature_select(ap$x,y,task,select,nfeatures,seed);sel<-fs$selected
  if(length(sel)<=5L)return(list(selected_encoded=sel,table=data.frame(nfeatures=length(sel),score=NA_real_,chosen=TRUE),note="Already compact."))
  scores<-vapply(sel,function(v).r4vn_machine_univariate_score(ap$x[,v],y,task),numeric(1));scores[!is.finite(scores)]<- -Inf;ranked<-sel[order(scores,decreasing=TRUE)]
  sizes<-unique(sort(pmax(3L,pmin(length(ranked),c(3L,5L,8L,10L,15L,length(ranked))))))
  # Use a fixed training-only holdout for a conservative parsimony screen.
  sp<-.r4vn_machine_split(raw,outcome_name,task,.8,seed+1L);tr<-sp$train;va<-sp$test;vals<-rep(NA_real_,length(sizes))
  for(i in seq_along(sizes)){
    ft<-tryCatch(.r4vn_machine_fit_pipeline(tr,outcome_name,predictor_names,task,outcome_info,method,params,preprocess,missing,missing_max,encode,standardize,transform,outlier,corr,feature,degree,reduce,variance,select,nfeatures,balance,balance_target,neighbors,selected_override=ranked[seq_len(sizes[i])],seed=seed+i),error=function(e)NULL)
    if(is.null(ft))next;pr<-tryCatch(.r4vn_machine_predict_pipeline(ft,va,outcome_name),error=function(e)NULL);if(is.null(pr))next;yy<-va[[outcome_name]][pr$rows]
    vals[i]<-if(task=="binary").r4vn_machine_binary_metric_value(yy,pr$prob,metric,.5) else if(task=="regression").r4vn_machine_reg_metric_value(yy,pr$value,metric) else .r4vn_machine_multi_metric_value(factor(yy,levels=colnames(pr$prob)),pr$prob,metric)
  }
  good<-which(is.finite(vals));if(!length(good))return(list(selected_encoded=sel,table=data.frame(),note="Parsimony screen failed; original selected set retained."))
  dir<-.r4vn_machine_metric_direction(metric);bestv<-if(dir=="max")max(vals[good]) else min(vals[good]);eligible<-if(dir=="max")which(vals>=bestv-tolerance) else which(vals<=bestv+tolerance)
  if(!length(eligible))eligible<-good;ii<-eligible[which.min(sizes[eligible])];chosen<-ranked[seq_len(sizes[ii])]
  list(selected_encoded=chosen,table=data.frame(nfeatures=sizes,score=vals,chosen=seq_along(sizes)==ii),
       note=paste0("Training-only parsimony screen retained ",length(chosen)," encoded predictors."))
}

.r4vn_machine_calibration <- function(y,p,ci=TRUE,level=.95,boot=1000,seed=NA_integer_,bins=10L) {
  ok<-!is.na(y)&is.finite(p);y<-as.numeric(y[ok]);p<-pmin(pmax(p[ok],1e-6),1-1e-6);lp<-stats::qlogis(p)
  fit_s<-tryCatch(stats::glm(y~lp,family=stats::binomial()),error=function(e)NULL)
  fit_i<-tryCatch(stats::glm(y~offset(lp),family=stats::binomial()),error=function(e)NULL)
  slope<-intercept<-NA_real_;slo<-shi<-ilo<-ihi<-NA_real_;smeth<-imeth<-NA_character_
  zcrit<-stats::qnorm(1-(1-level)/2)
  if(!is.null(fit_s)){
    cf<-tryCatch(summary(fit_s)$coefficients,error=function(e)NULL);slope<-unname(stats::coef(fit_s)["lp"])
    if(isTRUE(ci)&&!is.null(cf)&&"lp"%in%rownames(cf)&&is.finite(cf["lp","Std. Error"])){
      slo<-slope-zcrit*cf["lp","Std. Error"];shi<-slope+zcrit*cf["lp","Std. Error"];smeth<-"Wald"
    }
  }
  if(!is.null(fit_i)){
    cf<-tryCatch(summary(fit_i)$coefficients,error=function(e)NULL);intercept<-unname(stats::coef(fit_i)[1L])
    if(isTRUE(ci)&&!is.null(cf)&&nrow(cf)&&is.finite(cf[1L,"Std. Error"])){
      ilo<-intercept-zcrit*cf[1L,"Std. Error"];ihi<-intercept+zcrit*cf[1L,"Std. Error"];imeth<-"Wald"
    }
  }
  brier<-mean((p-y)^2);blo<-bhi<-NA_real_;bmeth<-NA_character_
  if(isTRUE(ci)){
    .r4vn_set_seed_if(seed);n<-length(y);bs<-bi<-bb<-rep(NA_real_,boot)
    for(b in seq_len(boot)){
      ix<-sample(seq_len(n),n,TRUE);yb<-y[ix];pb<-p[ix];lpb<-stats::qlogis(pb)
      bb[b]<-mean((pb-yb)^2)
      fs<-tryCatch(stats::glm(yb~lpb,family=stats::binomial()),error=function(e)NULL)
      fi<-tryCatch(stats::glm(yb~offset(lpb),family=stats::binomial()),error=function(e)NULL)
      if(!is.null(fs))bs[b]<-unname(stats::coef(fs)["lpb"])
      if(!is.null(fi))bi[b]<-unname(stats::coef(fi)[1L])
    }
    aq<-c((1-level)/2,1-(1-level)/2)
    qfun<-function(v){v<-v[is.finite(v)];if(length(v)<30L)c(NA_real_,NA_real_)else stats::quantile(v,aq,na.rm=TRUE,names=FALSE)}
    qb<-qfun(bb);blo<-qb[1];bhi<-qb[2];bmeth<-"Bootstrap"
    if(!is.finite(slo)||!is.finite(shi)){qs<-qfun(bs);slo<-qs[1];shi<-qs[2];if(all(is.finite(qs)))smeth<-"Bootstrap"}
    if(!is.finite(ilo)||!is.finite(ihi)){qi<-qfun(bi);ilo<-qi[1];ihi<-qi[2];if(all(is.finite(qi)))imeth<-"Bootstrap"}
  }
  qs<-unique(stats::quantile(p,probs=seq(0,1,length.out=bins+1L),na.rm=TRUE));if(length(qs)<3L)grp<-factor(rep(1,length(p))) else grp<-cut(p,breaks=qs,include.lowest=TRUE,labels=FALSE)
  curve<-aggregate(cbind(predicted=p,observed=y)~grp,FUN=mean);curve$n<-as.numeric(table(grp)[as.character(curve$grp)])
  events <- tapply(y, grp, sum); curve$events <- as.numeric(events[as.character(curve$grp)])
  curve$observed_lower <- curve$observed_upper <- NA_real_
  curve$ci_method <- NA_character_
  if (isTRUE(ci)) {
    for (i in seq_len(nrow(curve))) {
      cc <- .r4vn_machine_wilson(curve$events[i], curve$n[i], level)
      curve$observed_lower[i] <- cc[1L]; curve$observed_upper[i] <- cc[2L]; curve$ci_method[i] <- "Wilson"
    }
  }
  statsdf<-data.frame(metric=c("Calibration intercept","Calibration slope","Brier score"),estimate=c(intercept,slope,brier),
                      lower=c(ilo,slo,blo),upper=c(ihi,shi,bhi),ci_method=c(imeth,smeth,bmeth),stringsAsFactors=FALSE)
  list(statistics=statsdf,curve=curve)
}

.r4vn_machine_decision_nb <- function(y,p,t) {
  n<-length(y);pred<-p>=t;tp<-sum(pred&y==1);fp<-sum(pred&y==0);tp/n-(fp/n)*(t/(1-t))
}

.r4vn_machine_decision <- function(y,p,thresholds=seq(.01,.99,.01),ci=TRUE,level=.95,boot=1000,seed=NA_integer_) {
  ok<-!is.na(y)&is.finite(p);y<-as.numeric(y[ok]);p<-p[ok];th<-sort(unique(thresholds[thresholds>0&thresholds<1]));prev<-mean(y==1)
  out<-data.frame(threshold=th,net_benefit=vapply(th,function(t).r4vn_machine_decision_nb(y,p,t),numeric(1)),
                  treat_all=prev-(1-prev)*th/(1-th),treat_none=0,lower=NA_real_,upper=NA_real_,
                  ci_method=if(isTRUE(ci))"Pointwise bootstrap" else NA_character_)
  if(isTRUE(ci)&&length(y)>5L){.r4vn_set_seed_if(seed);mat<-matrix(NA_real_,boot,length(th));for(b in seq_len(boot)){ix<-sample(seq_along(y),length(y),TRUE);mat[b,]<-vapply(th,function(t).r4vn_machine_decision_nb(y[ix],p[ix],t),numeric(1))};alpha<-(1-level)/2;qq<-apply(mat,2,stats::quantile,probs=c(alpha,1-alpha),na.rm=TRUE);out$lower<-qq[1,];out$upper<-qq[2,]}
  out
}

.r4vn_machine_importance <- function(fit,raw,outcome_name,metric,threshold,repeats=20,seed=NA_integer_,labels=NULL,ci=TRUE,level=.95) {
  basepr<-.r4vn_machine_predict_pipeline(fit,raw,outcome_name);yy<-raw[[outcome_name]][basepr$rows]
  base<-if(fit$task=="binary").r4vn_machine_binary_metric_value(yy,basepr$prob,metric,threshold) else if(fit$task=="regression").r4vn_machine_reg_metric_value(yy,basepr$value,metric) else .r4vn_machine_multi_metric_value(factor(yy,levels=colnames(basepr$prob)),basepr$prob,metric)
  vars<-fit$recipe$keep;res<-list();.r4vn_set_seed_if(seed);dir<-.r4vn_machine_metric_direction(metric)
  for(v in vars){vals<-numeric(repeats);for(r in seq_len(repeats)){d<-raw;d[[v]]<-sample(d[[v]]);pr<-tryCatch(.r4vn_machine_predict_pipeline(fit,d,outcome_name),error=function(e)NULL);if(is.null(pr)){vals[r]<-NA;next};yt<-d[[outcome_name]][pr$rows];pm<-if(fit$task=="binary").r4vn_machine_binary_metric_value(yt,pr$prob,metric,threshold) else if(fit$task=="regression").r4vn_machine_reg_metric_value(yt,pr$value,metric) else .r4vn_machine_multi_metric_value(factor(yt,levels=colnames(pr$prob)),pr$prob,metric);vals[r]<-if(dir=="max")base-pm else pm-base}
    a<-(1-level)/2;lo<-if(isTRUE(ci))stats::quantile(vals,a,na.rm=TRUE,names=FALSE)else NA_real_;hi<-if(isTRUE(ci))stats::quantile(vals,1-a,na.rm=TRUE,names=FALSE)else NA_real_
    res[[v]]<-data.frame(variable=v,label=if(!is.null(labels)&&v%in%names(labels))labels[[v]] else v,importance=mean(vals,na.rm=TRUE),lower=lo,upper=hi,metric=metric,ci_method=if(isTRUE(ci))"Empirical permutation interval" else NA_character_,stringsAsFactors=FALSE)
  }
  out<-.r4vn_machine_bind(res);if(nrow(out))out<-out[order(out$importance,decreasing=TRUE),,drop=FALSE];row.names(out)<-NULL;out
}

.r4vn_machine_shap <- function(fit,raw,outcome_name) {
  if(fit$method!="xgb"||!requireNamespace("xgboost",quietly=TRUE))return(NULL)
  ap<-.r4vn_machine_recipe_apply(fit$recipe,raw);x<-ap$x[,fit$selected,drop=FALSE];xx<-.r4vn_machine_scale_apply(x,fit$engine$scaler)
  tryCatch(predict(fit$engine$model,xgboost::xgb.DMatrix(xx),predcontrib=TRUE),error=function(e)NULL)
}

.r4vn_machine_pdp <- function(fit,raw,outcome_name,pdp,importance,labels=NULL) {
  vars<-if(isTRUE(pdp)){if(nrow(importance))head(importance$variable[vapply(raw[importance$variable],is.numeric,logical(1))],5L) else head(fit$recipe$keep[vapply(raw[fit$recipe$keep],is.numeric,logical(1))],5L)} else as.character(pdp)
  vars<-intersect(vars,fit$recipe$keep);out<-list()
  for(v in vars){if(!is.numeric(raw[[v]]))next;grid<-unique(as.numeric(stats::quantile(raw[[v]],probs=seq(.05,.95,length.out=20),na.rm=TRUE)));rows<-list();for(i in seq_along(grid)){d<-raw;d[[v]]<-grid[i];pr<-.r4vn_machine_predict_pipeline(fit,d,outcome_name);if(fit$task=="binary")rows[[i]]<-data.frame(variable=v,label=.r4vn_machine_get_label(labels,v),value=grid[i],class=NA,prediction=mean(pr$prob,na.rm=TRUE)) else if(fit$task=="regression")rows[[i]]<-data.frame(variable=v,label=.r4vn_machine_get_label(labels,v),value=grid[i],class=NA,prediction=mean(pr$value,na.rm=TRUE)) else {mm<-colMeans(pr$prob,na.rm=TRUE);rows[[i]]<-data.frame(variable=v,label=.r4vn_machine_get_label(labels,v),value=grid[i],class=names(mm),prediction=as.numeric(mm))}};out[[v]]<-.r4vn_machine_bind(rows)}
  .r4vn_machine_bind(out)
}

.r4vn_machine_learning_curve <- function(raw,outcome_name,predictor_names,task,outcome_info,method,params,
                                         preprocess,missing,missing_max,encode,standardize,transform,outlier,corr,
                                         feature,degree,reduce,variance,select,nfeatures,balance,balance_target,neighbors,metric,threshold,
                                         ci=TRUE,ci_level=.95,seed=NA_integer_) {
  fracs<-c(.3,.5,.7,.9,1);rows<-list();ii<-0L;nr<-if(isTRUE(ci))5L else 3L
  for(f in fracs){
    vals<-rep(NA_real_,nr)
    for(r in seq_len(nr)){
      .r4vn_set_seed_if(seed+r+round(f*100));nsub<-max(20L,floor(nrow(raw)*f));idx<-sample(seq_len(nrow(raw)),min(nsub,nrow(raw)));sub<-raw[idx,,drop=FALSE]
      sp<-.r4vn_machine_split(sub,outcome_name,task,.8,seed+r)
      ft<-tryCatch(.r4vn_machine_fit_pipeline(sp$train,outcome_name,predictor_names,task,outcome_info,method,params,preprocess,missing,missing_max,encode,standardize,transform,outlier,corr,feature,degree,reduce,variance,select,nfeatures,balance,balance_target,neighbors,seed=seed+r),error=function(e)NULL)
      if(is.null(ft)||!nrow(sp$test))next
      pr<-.r4vn_machine_predict_pipeline(ft,sp$test,outcome_name);yy<-sp$test[[outcome_name]][pr$rows]
      vals[r]<-if(task=="binary").r4vn_machine_binary_metric_value(yy,pr$prob,metric,threshold) else if(task=="regression").r4vn_machine_reg_metric_value(yy,pr$value,metric) else .r4vn_machine_multi_metric_value(factor(yy,levels=colnames(pr$prob)),pr$prob,metric)
    }
    ii<-ii+1L;good<-vals[is.finite(vals)];est<-if(length(good))mean(good)else NA_real_;lo<-hi<-NA_real_;cm<-NA_character_
    if(isTRUE(ci)&&length(good)>=2L){se<-stats::sd(good)/sqrt(length(good));crit<-stats::qt(1-(1-ci_level)/2,df=length(good)-1L);lo<-est-crit*se;hi<-est+crit*se;cm<-"t interval across repeated holdouts"}
    rows[[ii]]<-data.frame(fraction=f,n=round(nrow(raw)*f),metric=metric,estimate=est,lower=lo,upper=hi,ci_method=cm,resamples=length(good),stringsAsFactors=FALSE)
  }
  .r4vn_machine_bind(rows)
}

# ============================================================================
# Reporting tables, prediction output, plots, HTML, and S3 methods
# ============================================================================

.r4vn_machine_bind <- function(x) {
  x<-x[vapply(x,function(z)!is.null(z)&&is.data.frame(z)&&nrow(z)>0,logical(1))]
  if(!length(x))return(data.frame())
  alln<-unique(unlist(lapply(x,names),use.names=FALSE))
  x<-lapply(x,function(z){miss<-setdiff(alln,names(z));for(m in miss)z[[m]]<-NA;z[,alln,drop=FALSE]})
  do.call(rbind,x)
}

.r4vn_machine_overview <- function(dat,outcome,predictors,task,outcome_info,train,internal,test,validation,metric) {
  n<-nrow(dat);missx<-sum(is.na(dat[predictors]));cells<-n*length(predictors)
  vals<-list(
    "Outcome"=outcome,
    "Task"=task,
    "Observations"=n,
    "Candidate predictors"=length(predictors),
    "Predictor missing cells"=paste0(missx," (",formatC(if(cells)100*missx/cells else 0,format="f",digits=1),"%)"),
    "Development training"=nrow(train),
    "Internal validation"=nrow(internal),
    "Untouched test"=nrow(test),
    "External validation"=if(is.null(validation))0 else nrow(validation),
    "Primary metric"=metric
  )
  if(task=="binary")vals<-append(vals,list("Event"=outcome_info$event,"Development event rate"=paste0(formatC(100*mean(train[[outcome]]==1),format="f",digits=1),"%")),after=2L)
  data.frame(Item=names(vals),Value=unlist(vals,use.names=FALSE),stringsAsFactors=FALSE)
}

.r4vn_machine_preprocessing_table <- function(fit) {
  a<-fit$recipe$actions
  if(!nrow(a))a<-data.frame(step="Preprocessing",variable="All predictors",action="No data-dependent filtering requested",stringsAsFactors=FALSE)
  a<-rbind(a,data.frame(step="Feature selection",variable="Encoded predictors",action=paste0("Method: ",if(length(fit$selection_detail))"selected" else "filter/none","; retained ",length(fit$selected)),stringsAsFactors=FALSE))
  sc<-fit$engine$scaler$method
  rbind(a,data.frame(step="Standardization",variable="Model matrix",action=if(sc=="none")"Not required for selected algorithm" else paste0(sc," scaling learned on training data"),stringsAsFactors=FALSE))
}

.r4vn_machine_selection_table <- function(fit,comparison,labels) {
  if (!is.null(fit$recipe$pca) && identical(fit$recipe$pca$method, "pca")) {
    selected <- data.frame(
      type = "Selected component", variable = fit$selected, label = fit$selected,
      value = "Retained after PCA", stringsAsFactors = FALSE
    )
    original <- data.frame(
      type = "Original predictor entering PCA", variable = fit$recipe$keep,
      label = vapply(fit$recipe$keep, function(v) .r4vn_machine_get_label(labels, v), character(1)),
      value = "Encoded before PCA", stringsAsFactors = FALSE
    )
  } else {
    map <- fit$recipe$raw_map[fit$selected]
    raw <- unique(na.omit(as.character(map)))
    if (!length(raw)) raw <- fit$recipe$keep
    selected <- data.frame(
      type = "Selected predictor", variable = raw,
      label = vapply(raw, function(v) .r4vn_machine_get_label(labels, v), character(1)),
      value = "Retained", stringsAsFactors = FALSE
    )
    original <- data.frame()
  }
  removed <- fit$recipe$removed
  rem <- if(length(removed)) data.frame(
    type="Preprocessing removal", variable=removed,
    label=vapply(removed,function(v).r4vn_machine_get_label(labels,v),character(1)),
    value="Removed", stringsAsFactors=FALSE
  ) else data.frame()
  cmp <- if(nrow(comparison)) data.frame(
    type="Strategy comparison", variable=comparison$strategy,label=comparison$strategy,
    value=paste0(comparison$metric,"=",formatC(comparison$estimate,format="f",digits=3)),
    stringsAsFactors=FALSE
  ) else data.frame()
  .r4vn_machine_bind(list(selected, original, rem, cmp))
}

.r4vn_machine_balance_table <- function(strategy,y,comparison,task) {
  if(task!="binary")return(data.frame(Item=c("Strategy","Reason"),Value=c("none","Class balancing is only applied to binary classification.")))
  tab<-table(y); n0<-if("0"%in%names(tab))as.numeric(tab[["0"]])else 0; n1<-if("1"%in%names(tab))as.numeric(tab[["1"]])else 0
  base<-data.frame(Item=c("Chosen strategy","Non-event n","Event n","Event proportion"),Value=c(strategy,n0,n1,formatC(mean(y==1),format="f",digits=3)),stringsAsFactors=FALSE)
  if(nrow(comparison)){
    cmp<-data.frame(Item=paste0("Compared: ",comparison$strategy),Value=paste0(comparison$metric," = ",formatC(comparison$estimate,format="f",digits=3)),stringsAsFactors=FALSE);base<-rbind(base,cmp)
  }
  base
}

.r4vn_machine_pretty_metric <- function(x) {
  map<-c(auc="ROC-AUC",pr_auc="PR-AUC",accuracy="Accuracy",balanced_accuracy="Balanced accuracy",sensitivity="Sensitivity",specificity="Specificity",ppv="PPV",npv="NPV",f1="F1",mcc="MCC",kappa="Kappa",brier="Brier score",logloss="Log loss",prevalence="Prevalence",rmse="RMSE",mae="MAE",r2="R-squared",mape="MAPE",macro_f1="Macro F1",weighted_f1="Weighted F1",macro_auc="Macro one-vs-rest AUC",macro_pr_auc="Macro one-vs-rest PR-AUC")
  z <- unname(map[x]); if(length(z)==1L && !is.na(z)) z else x
}

.r4vn_machine_final_table <- function(performance,best,digit=3) {
  z<-performance[performance$method==best,,drop=FALSE]
  if(!nrow(z))return(data.frame())
  fmt<-function(e,l,u)if(is.finite(l)&&is.finite(u))paste0(formatC(e,format="f",digits=digit)," (",formatC(l,format="f",digits=digit),"-",formatC(u,format="f",digits=digit),")") else formatC(e,format="f",digits=digit)
  data.frame(Measure=vapply(z$metric,.r4vn_machine_pretty_metric,character(1)),`Estimate (CI)`=mapply(fmt,z$estimate,z$lower,z$upper),`CI method`=z$ci_method,check.names=FALSE,stringsAsFactors=FALSE)
}

.r4vn_machine_tables <- function(out) {
  lst<-list(
    "Data and modelling strategy"=out$overview,
    "Algorithms and package requirements"=out$engines,
    "Preprocessing"=out$preprocessing,
    "Class balance"=out$balance,
    "Feature selection"=out$selection,
    "Development cross-validation"=out$cv_performance,
    "Model comparison"=out$comparison,
    "Final model performance"=.r4vn_machine_final_table(out$performance,out$best,out$settings$digit)
  )
  if(nrow(out$tuning))lst[["Hyperparameter tuning"]]<-out$tuning
  if(!is.null(out$simplify) && is.list(out$simplify) && is.data.frame(out$simplify$table) && nrow(out$simplify$table)) {
    st <- out$simplify$table
    if (!is.null(out$simplify$note) && length(out$simplify$note)) st$Note <- c(out$simplify$note, rep("", max(0L, nrow(st)-1L)))
    lst[["Parsimony screen"]] <- st
  }
  if(nrow(out$model_difference))lst[["Paired model differences"]]<-out$model_difference
  if(nrow(out$overfitting))lst[["Development versus evaluation"]]<-out$overfitting
  if(nrow(out$coefficients))lst[["Interpretable model coefficients"]]<-out$coefficients
  if(!is.null(out$threshold)&&length(out$threshold))lst[["Classification threshold"]]<-data.frame(Method=out$threshold$method,Threshold=out$threshold$threshold,Lower=out$threshold$lower,Upper=out$threshold$upper,`CI method`=out$threshold$ci_method,check.names=FALSE)
  if(nrow(out$confusion))lst[["Confusion matrix"]]<-out$confusion
  if(nrow(out$class_performance))lst[["Class-specific performance"]]<-out$class_performance
  if(!is.null(out$calibration))lst[["Calibration"]]<-out$calibration$statistics
  if(nrow(out$importance))lst[["Predictor importance"]]<-out$importance
  if(nrow(out$learning))lst[["Learning curve"]]<-out$learning
  if(!is.null(out$plots$shap) && is.data.frame(out$plots$shap) && nrow(out$plots$shap)) lst[["SHAP summary"]] <- out$plots$shap
  if(nrow(out$external_performance))lst[["External validation"]]<-.r4vn_machine_final_table(out$external_performance,out$best,out$settings$digit)
  if(nrow(out$external_confusion))lst[["External confusion matrix"]]<-out$external_confusion
  if(nrow(out$external_class_performance))lst[["External class-specific performance"]]<-out$external_class_performance
  if(!is.null(out$external_calibration))lst[["External calibration"]]<-out$external_calibration$statistics
  if(!is.null(out$predictions)&&is.data.frame(out$predictions))lst[["Predictions"]]<-out$predictions
  lst[vapply(lst,function(z)is.data.frame(z)&&nrow(z)>0,logical(1))]
}

.r4vn_machine_prediction_frame <- function(fit,newdata,threshold=.5,id_name=NULL,outcome_info=NULL) {
  pr<-.r4vn_machine_predict_pipeline(fit,newdata)
  idv<-if(!is.null(id_name)&&id_name%in%names(newdata))newdata[[id_name]][pr$rows] else pr$rows
  if(fit$task=="binary"){
    cl<-ifelse(pr$prob>=threshold,outcome_info$event,outcome_info$negative)
    return(data.frame(ID=idv,Probability=pr$prob,Predicted=cl,stringsAsFactors=FALSE))
  }
  if(fit$task=="regression")return(data.frame(ID=idv,Prediction=pr$value,stringsAsFactors=FALSE))
  pp<-as.data.frame(pr$prob,check.names=FALSE);names(pp)<-paste0("Probability: ",names(pp));data.frame(ID=idv,Predicted=pr$class,pp,check.names=FALSE,stringsAsFactors=FALSE)
}

.r4vn_machine_multiclass_confusion <- function(truth, predicted, levels) {
  truth <- factor(truth, levels = levels)
  predicted <- factor(predicted, levels = levels)
  tab <- table(truth, predicted, useNA = "no")
  out <- as.data.frame.matrix(tab, stringsAsFactors = FALSE)
  names(out) <- paste0("Predicted: ", names(out))
  data.frame(Actual = rownames(out), out, row.names = NULL, check.names = FALSE, stringsAsFactors = FALSE)
}

.r4vn_machine_confusion_matrix <- function(tab) {
  if (is.null(tab) || !is.data.frame(tab) || !nrow(tab) || ncol(tab) < 3L) return(NULL)
  z <- as.matrix(tab[, -1L, drop = FALSE])
  storage.mode(z) <- "numeric"
  rownames(z) <- as.character(tab[[1L]])
  colnames(z) <- sub("^Predicted: ?", "", sub("^Predicted ", "", colnames(z)))
  z
}

.r4vn_machine_roc_curve <- function(y,p) {
  cuts<-sort(unique(c(Inf,p,-Inf)),decreasing=TRUE);if(length(cuts)>1000L)cuts<-unique(stats::quantile(cuts[is.finite(cuts)],probs=seq(0,1,length.out=998),na.rm=TRUE));cuts<-c(Inf,sort(cuts,decreasing=TRUE),-Inf)
  rows<-lapply(cuts,function(t){m<-.r4vn_machine_binary_metrics(y,p,t);data.frame(threshold=t,fpr=1-m["specificity"],tpr=m["sensitivity"])})
  z<-do.call(rbind,rows);z[order(z$fpr,z$tpr),,drop=FALSE]
}

.r4vn_machine_pr_curve <- function(y,p) {
  o<-order(p,decreasing=TRUE);yy<-y[o];tp<-cumsum(yy==1);fp<-cumsum(yy==0);data.frame(recall=tp/max(1,sum(y==1)),precision=tp/pmax(tp+fp,1),threshold=p[o])
}

.r4vn_machine_plot_choices <- function() {
  c("roc", "pr", "calibration", "threshold", "confusion", "importance",
    "decision", "learning", "observed", "residual", "pdp", "shap")
}

.r4vn_machine_plot_title_map <- function() {
  c(
    roc = "ROC curve",
    pr = "Precision-recall curve",
    calibration = "Calibration plot",
    threshold = "Classification threshold",
    confusion = "Confusion matrix",
    importance = "Predictor importance",
    decision = "Decision curve analysis",
    learning = "Learning curve",
    observed = "Observed versus predicted",
    residual = "Residual plot",
    pdp = "Partial dependence",
    shap = "SHAP summary"
  )
}

.r4vn_machine_validate_plot_options <- function(plot, plot_display, plot_args, strict) {
  choices <- .r4vn_machine_plot_choices()
  validate_spec <- function(z, name, allow_auto = TRUE) {
    if (is.null(z) || identical(z, FALSE) || identical(z, TRUE)) return(invisible(TRUE))
    if (!is.character(z)) stop("`", name, "` must be TRUE, FALSE, or plot name(s).", call. = FALSE)
    q <- unique(tolower(z[!is.na(z) & nzchar(z)]))
    allowed <- c(choices, "all", "none", if (allow_auto) "auto")
    bad <- setdiff(q, allowed)
    if (length(bad)) stop("Unsupported `", name, "`: ", paste(bad, collapse = ", "), ".", call. = FALSE)
    invisible(TRUE)
  }
  validate_spec(plot, "plot")
  validate_spec(plot_display, "plot_display")
  if (!is.list(plot_args) || (length(plot_args) && (is.null(names(plot_args)) || any(!nzchar(names(plot_args)))))) {
    stop("`plot_args` must be a named list.", call. = FALSE)
  }
  nested <- intersect(names(plot_args), c("all", choices))
  if (length(nested) && any(!vapply(plot_args[nested], is.list, logical(1)))) {
    stop("Nested `plot_args` entries such as `all`, `roc`, or `calibration` must be lists.", call. = FALSE)
  }
  if (!is.logical(strict) || length(strict) != 1L || is.na(strict)) stop("`strict` must be TRUE or FALSE.", call. = FALSE)
  invisible(TRUE)
}

.r4vn_machine_plot_options <- function(plot_args, type) {
  if (!length(plot_args)) return(list())
  choices <- c("all", .r4vn_machine_plot_choices())
  nested <- intersect(names(plot_args), choices)
  common <- plot_args[setdiff(names(plot_args), nested)]
  if ("all" %in% names(plot_args)) common <- utils::modifyList(common, plot_args$all)
  if (type %in% names(plot_args)) common <- utils::modifyList(common, plot_args[[type]])
  common
}

.r4vn_machine_plot_available <- function(z) {
  if (is.null(z)) return(FALSE)
  if (is.data.frame(z) || is.matrix(z)) return(nrow(z) > 0L && ncol(z) > 0L)
  if (is.list(z)) return(length(z) > 0L && any(vapply(z, function(a) {
    if (is.data.frame(a) || is.matrix(a)) nrow(a) > 0L else length(a) > 0L
  }, logical(1))))
  length(z) > 0L
}

.r4vn_machine_plot_data <- function(out,truth,predictions) {
  z<-list()
  if(out$settings$task=="binary"){
    z$roc<-lapply(predictions,function(pr).r4vn_machine_roc_curve(as.numeric(truth[pr$rows%||%seq_along(truth)]),pr$prob))
    z$pr<-lapply(predictions,function(pr).r4vn_machine_pr_curve(as.numeric(truth[pr$rows%||%seq_along(truth)]),pr$prob))
    z$calibration<-out$calibration$curve%||%data.frame()
    z$threshold<-out$threshold$curve%||%data.frame()
    z$decision<-out$decision
    z$confusion<-.r4vn_machine_confusion_matrix(out$confusion)
  } else if(out$settings$task=="regression"){
    pr<-predictions[[out$best]]
    z$observed<-data.frame(observed=as.numeric(truth[pr$rows%||%seq_along(truth)]),predicted=pr$value)
    z$residual<-transform(z$observed,residual=observed-predicted)
  } else if(out$settings$task=="multiclass") {
    z$confusion<-.r4vn_machine_confusion_matrix(out$confusion)
  }
  z$importance<-out$importance
  z$learning<-out$learning
  z$pdp<-out$pdp
  if (!is.null(out$shap) && (is.matrix(out$shap) || is.data.frame(out$shap))) {
    sm <- as.matrix(out$shap)
    if (ncol(sm)) {
      keep <- !grepl("^(BIAS|Bias|bias)$", colnames(sm) %||% rep("", ncol(sm)))
      if (any(keep)) {
        sm <- sm[, keep, drop = FALSE]
        z$shap <- data.frame(
          feature = colnames(sm) %||% paste0("Feature ", seq_len(ncol(sm))),
          mean_abs_shap = colMeans(abs(sm), na.rm = TRUE),
          stringsAsFactors = FALSE, check.names = FALSE
        )
        z$shap <- z$shap[order(z$shap$mean_abs_shap, decreasing = TRUE), , drop = FALSE]
      }
    }
  }
  z
}

.r4vn_machine_available_plot_types <- function(out) {
  choices <- .r4vn_machine_plot_choices()
  intersect(choices, names(out$plots)[vapply(out$plots, .r4vn_machine_plot_available, logical(1))])
}

.r4vn_machine_plot_types <- function(out,plot) {
  avail <- .r4vn_machine_available_plot_types(out)
  if(identical(plot,FALSE)||is.null(plot))return(character())
  default <- if(out$settings$task=="binary") {
    c("roc","pr","calibration","threshold","confusion","importance","decision","learning","pdp","shap")
  } else if(out$settings$task=="regression") {
    c("observed","residual","importance","learning","pdp","shap")
  } else c("confusion","importance","learning","pdp","shap")
  if(isTRUE(plot))return(intersect(default,avail))
  q <- unique(tolower(as.character(plot)))
  if("none" %in% q)return(character())
  if("all" %in% q)return(avail)
  if("auto" %in% q)return(intersect(default,avail))
  intersect(q,avail)
}

.r4vn_machine_display_types <- function(out,plot_display) {
  avail <- .r4vn_machine_available_plot_types(out)
  if(is.null(plot_display)||identical(plot_display,FALSE))return(character())
  if(isTRUE(plot_display)) plot_display <- "auto"
  q <- unique(tolower(as.character(plot_display)))
  if("none" %in% q)return(character())
  if("all" %in% q)return(avail)
  if("auto" %in% q) {
    q <- if(out$settings$task=="binary") c("roc","calibration","importance") else
      if(out$settings$task=="regression") c("observed","residual","importance") else c("confusion","importance")
  }
  intersect(q,avail)
}

.r4vn_machine_metric_legend <- function(x, methods, metric, label) {
  vapply(methods, function(md) {
    r <- x$performance[x$performance$method == md & x$performance$metric == metric, , drop = FALSE]
    if (!nrow(r) || !is.finite(r$estimate[1L])) md else paste0(md, " (", label, " ", formatC(r$estimate[1L], format="f", digits=x$settings$digit), ")")
  }, character(1))
}

.r4vn_machine_plot_call <- function(fun, fixed, dots = list(), protected = names(fixed)) {
  protected <- unique(c(protected, names(fixed)))
  if (length(dots)) {
    dn <- names(dots)
    keep <- is.null(dn) || !(dn %in% protected)
    if (!is.null(dn)) dots <- dots[keep]
  }
  do.call(fun, c(fixed, dots))
}

#' Plot a tabmachine result
#' @param x A `r4vn_machine` object.
#' @param type Plot type: `roc`, `pr`, `calibration`, `threshold`, `confusion`,
#'   `importance`, `decision`, `learning`, `observed`, `residual`, `pdp`, or
#'   `shap`.
#' @param title Optional figure title. A publication-ready default is supplied.
#' @param font_family Base-graphics font family. Default `"sans"`, which is also
#'   the safest choice for HTML/SVG Viewer rendering.
#' @param ... Additional base-graphics arguments where applicable.
#' @return The input `r4vn_machine` object, invisibly. The requested figure is
#'   drawn on the current graphics device as a side effect; the fitted machine-
#'   learning result itself is not modified.
#' @export
plot.r4vn_machine <- function(x,type=NULL,title=NULL,font_family="sans",...) {
  choices<-.r4vn_machine_plot_choices()
  if(is.null(type)) type<-if(x$settings$task=="binary")"roc" else if(x$settings$task=="regression")"observed" else "confusion"
  type<-match.arg(type,choices);pd<-x$plots[[type]]
  if(!.r4vn_machine_plot_available(pd))stop("Plot data are not available for `",type,"`.",call.=FALSE)
  dots <- list(...)
  line_width <- if (!is.null(dots$lwd) && is.numeric(dots$lwd) && length(dots$lwd)) as.numeric(dots$lwd[1L]) else 2
  dots$lwd <- NULL
  if (is.null(title) && !is.null(dots$main)) { title <- dots$main; dots$main <- NULL }
  if (is.null(title)) title <- unname(.r4vn_machine_plot_title_map()[type])
  oldpar <- graphics::par(no.readonly = TRUE)
  on.exit(graphics::par(oldpar), add = TRUE)
  try(graphics::par(family = font_family), silent = TRUE)

  if(type=="roc"){
    mods<-names(pd); .r4vn_machine_plot_call(graphics::plot,
      list(x=c(0,1),y=c(0,1),type="n",xlab="1 - Specificity",ylab="Sensitivity",xlim=c(0,1),ylim=c(0,1),main=title),dots,
      protected=c("x","y","type","xlab","ylab","xlim","ylim","main"))
    graphics::abline(0,1,lty=2)
    lty<-seq_along(mods);for(i in seq_along(mods))graphics::lines(pd[[i]]$fpr,pd[[i]]$tpr,lty=lty[i],lwd=line_width)
    graphics::legend("bottomright",legend=.r4vn_machine_metric_legend(x,mods,"auc","AUC"),lty=lty,lwd=line_width,bty="n")
    return(invisible(x))
  }
  if(type=="pr"){
    mods<-names(pd); .r4vn_machine_plot_call(graphics::plot,
      list(x=c(0,1),y=c(0,1),type="n",xlab="Recall (Sensitivity)",ylab="Precision (PPV)",xlim=c(0,1),ylim=c(0,1),main=title),dots,
      protected=c("x","y","type","xlab","ylab","xlim","ylim","main"))
    lty<-seq_along(mods);for(i in seq_along(mods))graphics::lines(pd[[i]]$recall,pd[[i]]$precision,lty=lty[i],lwd=line_width)
    graphics::legend("bottomleft",legend=.r4vn_machine_metric_legend(x,mods,"pr_auc","PR-AUC"),lty=lty,lwd=line_width,bty="n")
    return(invisible(x))
  }
  if(type=="calibration"){
    .r4vn_machine_plot_call(graphics::plot,
      list(x=pd$predicted,y=pd$observed,xlim=c(0,1),ylim=c(0,1),xlab="Mean predicted probability",ylab="Observed event proportion",pch=19,main=title),dots,
      protected=c("x","y","xlim","ylim","xlab","ylab","main"))
    graphics::abline(0,1,lty=2)
    if(all(c("observed_lower","observed_upper")%in%names(pd))) graphics::segments(pd$predicted,pd$observed_lower,pd$predicted,pd$observed_upper)
    graphics::lines(pd$predicted,pd$observed,lwd=line_width);return(invisible(x))
  }
  if(type=="threshold"){
    .r4vn_machine_plot_call(graphics::plot,
      list(x=pd$threshold,y=pd$sensitivity,type="l",lwd=line_width,xlim=c(0,1),ylim=c(0,1),xlab="Probability threshold",ylab="Performance",main=title),dots,
      protected=c("x","y","type","xlim","ylim","xlab","ylab","main"))
    graphics::lines(pd$threshold,pd$specificity,lty=2,lwd=line_width)
    if(!is.null(x$threshold)&&is.finite(x$threshold$threshold))graphics::abline(v=x$threshold$threshold,lty=3)
    graphics::legend("bottomleft",c("Sensitivity","Specificity","Selected threshold"),lty=c(1,2,3),lwd=c(2,2,1),bty="n")
    return(invisible(x))
  }
  if(type=="confusion"){
    m<-pd; nr<-nrow(m);nc<-ncol(m); mx<-max(m,na.rm=TRUE); den<-if(is.finite(mx)&&mx>0)mx else 1
    graphics::plot(c(.5,nc+.5),c(.5,nr+.5),type="n",axes=FALSE,xlab="Predicted",ylab="Actual",main=title,xaxs="i",yaxs="i")
    for(i in seq_len(nr))for(j in seq_len(nc)){
      shade <- 1 - .35 * (m[i,j]/den)
      graphics::rect(j-.5,nr-i+.5,j+.5,nr-i+1.5,col=grDevices::gray(shade),border="white")
      graphics::text(j,nr-i+1,labels=format(m[i,j],trim=TRUE))
    }
    graphics::axis(1,at=seq_len(nc),labels=colnames(m),las=2)
    graphics::axis(2,at=seq_len(nr),labels=rev(rownames(m)),las=2)
    graphics::box();return(invisible(x))
  }
  if(type=="importance"){
    d<-head(pd[order(pd$importance,decreasing=TRUE),,drop=FALSE],20L)
    .r4vn_machine_plot_call(graphics::barplot,
      list(height=rev(d$importance),names.arg=rev(d$label),horiz=TRUE,las=1,xlab=paste("Permutation importance (",unique(d$metric)[1],")",sep=""),main=title),dots,
      protected=c("height","names.arg","horiz","las","xlab","main"));return(invisible(x))
  }
  if(type=="decision"){
    .r4vn_machine_plot_call(graphics::plot,
      list(x=pd$threshold,y=pd$net_benefit,type="l",lwd=line_width,xlab="Threshold probability",ylab="Net benefit",main=title),dots,
      protected=c("x","y","type","xlab","ylab","main"))
    if(all(c("lower","upper")%in%names(pd))&&any(is.finite(pd$lower))&&any(is.finite(pd$upper))){graphics::lines(pd$threshold,pd$lower,lty=3);graphics::lines(pd$threshold,pd$upper,lty=3)}
    graphics::lines(pd$threshold,pd$treat_all,lty=2);graphics::abline(h=0,lty=3);graphics::legend("topright",c("Model","Treat all","Treat none"),lty=c(1,2,3),bty="n");return(invisible(x))
  }
  if(type=="learning"){
    .r4vn_machine_plot_call(graphics::plot,
      list(x=pd$n,y=pd$estimate,type="b",pch=19,xlab="Training observations",ylab=.r4vn_machine_pretty_metric(unique(pd$metric)[1]),main=title),dots,
      protected=c("x","y","type","xlab","ylab","main"))
    if(all(is.finite(pd$lower))&&all(is.finite(pd$upper)))graphics::segments(pd$n,pd$lower,pd$n,pd$upper);return(invisible(x))
  }
  if(type=="observed"){
    .r4vn_machine_plot_call(graphics::plot,
      list(x=pd$observed,y=pd$predicted,xlab="Observed",ylab="Predicted",pch=19,main=title),dots,
      protected=c("x","y","xlab","ylab","main"));graphics::abline(0,1,lty=2);return(invisible(x))
  }
  if(type=="residual"){
    .r4vn_machine_plot_call(graphics::plot,
      list(x=pd$predicted,y=pd$residual,xlab="Predicted",ylab="Residual",pch=19,main=title),dots,
      protected=c("x","y","xlab","ylab","main"));graphics::abline(h=0,lty=2);return(invisible(x))
  }
  if(type=="pdp"){
    panels <- pd[vapply(pd, function(d) is.data.frame(d) && nrow(d) && all(c("value","prediction") %in% names(d)), logical(1))]
    if (!length(panels)) stop("Partial-dependence data are unavailable.", call. = FALSE)
    np <- length(panels); nc <- if (np == 1L) 1L else 2L; nr <- ceiling(np / nc)
    graphics::par(mfrow = c(nr, nc), mar = c(4.2, 4.2, 3.1, 1.1))
    for (nm in names(panels)) {
      d <- panels[[nm]]
      lab <- if ("label" %in% names(d) && any(nzchar(as.character(d$label)))) as.character(d$label[which(nzchar(as.character(d$label)))[1L]]) else nm
      cls <- if ("class" %in% names(d)) unique(as.character(d$class[!is.na(d$class)])) else character()
      panel_title <- if (np == 1L) paste0(title, ": ", lab) else lab
      if (!length(cls)) {
        graphics::plot(d$value, d$prediction, type="l", lwd=line_width, xlab=lab, ylab="Partial dependence", main=panel_title)
      } else {
        yr <- range(d$prediction, finite=TRUE); xr <- range(d$value, finite=TRUE)
        graphics::plot(xr, yr, type="n", xlab=lab, ylab="Partial dependence", main=panel_title)
        lt <- seq_along(cls)
        for (i in seq_along(cls)) {
          dd <- d[as.character(d$class) == cls[i], , drop=FALSE]
          graphics::lines(dd$value, dd$prediction, lty=lt[i], lwd=line_width)
        }
        graphics::legend("topright", legend=cls, lty=lt, lwd=line_width, bty="n")
      }
    }
    return(invisible(x))
  }
  if(type=="shap"){
    d <- head(pd[order(pd$mean_abs_shap, decreasing=TRUE), , drop=FALSE], 20L)
    .r4vn_machine_plot_call(graphics::barplot,
      list(height=rev(d$mean_abs_shap),names.arg=rev(d$feature),horiz=TRUE,las=1,xlab="Mean absolute SHAP value",main=title),dots,
      protected=c("height","names.arg","horiz","las","xlab","main"))
    return(invisible(x))
  }
  invisible(x)
}

.r4vn_machine_draw_display <- function(out, types, plot_args = list(), strict = FALSE) {
  for (tp in types) {
    opts <- .r4vn_machine_plot_options(plot_args, tp)
    ok <- tryCatch({ do.call(graphics::plot, c(list(x=out,type=tp),opts)); TRUE }, error=function(e) {
      if (isTRUE(strict)) stop(e)
      warning("Machine-learning Plot-pane figure `",tp,"` was unavailable: ",conditionMessage(e),call.=FALSE)
      FALSE
    })
  }
  invisible(out)
}

.r4vn_machine_html_escape <- function(x) {
  x<-as.character(x);x<-gsub("&","&amp;",x,fixed=TRUE);x<-gsub("<","&lt;",x,fixed=TRUE);x<-gsub(">","&gt;",x,fixed=TRUE);x<-gsub('"',"&quot;",x,fixed=TRUE);x
}

.r4vn_machine_html_table <- function(d,digit=3) {
  if(is.null(d)||!is.data.frame(d)||!nrow(d))return("")
  z<-d
  # Wide-and-shallow result tables are hard to read in the Viewer. When the
  # first column identifies a few models/groups, rotate only the HTML display;
  # the returned R data frame remains unchanged for reproducibility/export.
  if(ncol(z) > 8L && nrow(z) <= 6L && (is.character(z[[1L]]) || is.factor(z[[1L]]))) {
    ids <- make.unique(as.character(z[[1L]]))
    body <- z[, -1L, drop = FALSE]
    zz <- data.frame(Measure = names(body), stringsAsFactors = FALSE, check.names = FALSE)
    for(i in seq_len(nrow(z))) zz[[ids[i]]] <- as.character(unlist(body[i, , drop = FALSE], use.names = FALSE))
    z <- zz
  }
  for(j in seq_along(z))if(is.numeric(z[[j]]))z[[j]]<-ifelse(is.na(z[[j]]),"",formatC(z[[j]],format="fg",digits=digit,flag="#"))
  head<-paste0("<tr>",paste0("<th>",.r4vn_machine_html_escape(names(z)),"</th>",collapse=""),"</tr>")
  body<-vapply(seq_len(nrow(z)),function(i){vals<-unlist(z[i,,drop=FALSE],use.names=FALSE);paste0("<tr>",paste0("<td>",.r4vn_machine_html_escape(vals),"</td>",collapse=""),"</tr>")},character(1))
  paste0("<div class='table-wrap'><table><thead>",head,"</thead><tbody>",paste(body,collapse=""),"</tbody></table></div>")
}

.r4vn_machine_svg <- function(out,type,plot_args=list()) {
  f<-tempfile(fileext=".svg");ok<-FALSE
  opts <- .r4vn_machine_plot_options(plot_args, type)
  tryCatch({
    grDevices::svg(f,width=7.2,height=5.2,onefile=TRUE,family="sans")
    do.call(graphics::plot,c(list(x=out,type=type),opts))
    grDevices::dev.off();ok<-TRUE
  },error=function(e){try(grDevices::dev.off(),silent=TRUE)})
  if(!ok||!file.exists(f))return("")
  ln<-readLines(f,warn=FALSE,encoding="UTF-8");unlink(f)
  i<-grep("<svg",ln,fixed=TRUE)[1L];if(is.na(i))return("");paste(ln[i:length(ln)],collapse="\n")
}

.r4vn_machine_html <- function(out,plot=TRUE,plot_args=list(),title=NULL,digit=3) {
  ttl<-title%||%"Machine learning analysis"
  sections<-character()
  add<-function(name,content){if(nzchar(content))sections<<-c(sections,paste0("<section><h2>",.r4vn_machine_html_escape(name),"</h2>",content,"</section>"))}
  if(length(out$tables)) {
    for (nm in names(out$tables)) {
      heading <- if (identical(nm,"Final model performance")) paste0("Final model: ",out$best) else nm
      add(heading,.r4vn_machine_html_table(out$tables[[nm]],digit))
    }
  }
  pts<-.r4vn_machine_plot_types(out,plot)
  if(length(pts))for(tp in pts){
    sv<-.r4vn_machine_svg(out,tp,plot_args=plot_args)
    if(nzchar(sv)) {
      heading <- out$plot_titles[[tp]] %||% tools::toTitleCase(tp)
      add(paste0("Figure: ",heading),paste0("<div class='plot'>",sv,"</div>"))
    }
  }
  if(length(out$notes))add("Methodological notes",paste0("<ul>",paste0("<li>",.r4vn_machine_html_escape(out$notes),"</li>",collapse=""),"</ul>"))
  paste0("<!doctype html><html><head><meta charset='utf-8'><style>",
         "body{font-family:Arial,Helvetica,sans-serif;margin:24px;color:#222;line-height:1.45}h1{font-size:25px;margin-bottom:4px}h2{font-size:18px;margin-top:28px;border-bottom:1px solid #ddd;padding-bottom:6px}.subtitle{color:#666;margin-bottom:18px}.table-wrap{overflow-x:auto}table{border-collapse:collapse;width:100%;font-size:13px}th,td{border:1px solid #ddd;padding:7px 9px;text-align:left;vertical-align:top}th{background:#f5f5f5;font-weight:600}tr:nth-child(even)td{background:#fafafa}.plot{max-width:900px;margin:10px auto}.plot svg{display:block;max-width:100%;height:auto;margin:auto}ul{padding-left:22px}",
         "</style></head><body><h1>",.r4vn_machine_html_escape(ttl),"</h1><div class='subtitle'>R4VN comprehensive machine-learning report</div>",paste(sections,collapse="\n"),"</body></html>")
}

.r4vn_machine_show_html <- function(html) {
  if(!interactive())return(invisible(FALSE));f<-tempfile(fileext=".html");writeLines(html,f,useBytes=TRUE);viewer<-getOption("viewer");if(is.function(viewer))viewer(f) else utils::browseURL(f);invisible(TRUE)
}

#' Print a tabmachine result
#' @param x A `r4vn_machine` object.
#' @param ... Unused.
#' @return The input `r4vn_machine` object, invisibly, after printing a compact
#'   summary of the analysis and the final-model performance table.
#' @export
print.r4vn_machine <- function(x,...) {
  cat("Machine learning analysis\n")
  cat("Task: ",x$settings$task," | Best model: ",x$best," | Primary metric: ",x$settings$metric,"\n",sep="")
  if(x$settings$task=="binary")cat("Event: ",x$settings$event," | Threshold: ",formatC(x$settings$threshold_value,format="f",digits=3),"\n",sep="")
  print(.r4vn_machine_final_table(x$performance,x$best,x$settings$digit),row.names=FALSE)
  invisible(x)
}

#' Predict from a tabmachine model
#' @param object A fitted `r4vn_machine` object.
#' @param newdata New predictor data.
#' @param type `"response"`, `"prob"`, or `"class"`.
#' @param threshold Optional binary threshold overriding the training-derived
#'   final threshold.
#' @param ... Unused.
#' @return Predictions whose structure depends on the task and `type`. For a
#'   regression task, a numeric vector of predicted outcomes is returned. For a
#'   binary task, `type = "response"` or `"prob"` returns a numeric vector of
#'   probabilities for the modeled event, while `type = "class"` returns a
#'   factor of predicted classes. For a multiclass task, `type = "response"`
#'   or `"prob"` returns a numeric matrix of class probabilities and
#'   `type = "class"` returns a factor of predicted classes.
#' @method predict r4vn_machine
#' @export
predict.r4vn_machine <- function(object,newdata,type=c("response","prob","class"),threshold=NULL,...) {
  type<-match.arg(type);pr<-.r4vn_machine_predict_pipeline(object$final,newdata);task<-object$settings$task
  if(task=="regression")return(pr$value)
  if(task=="binary"){
    th<-threshold%||%object$settings$threshold_value
    if(type%in%c("response","prob"))return(pr$prob)
    return(factor(ifelse(pr$prob>=th,object$settings$event,object$final$outcome_info$negative),levels=object$final$outcome_info$levels))
  }
  if(type=="prob")return(pr$prob);if(type=="class")return(factor(pr$class,levels=object$final$outcome_info$levels));pr$prob
}

#' Convert a tabmachine result to a data frame
#' @param x A `r4vn_machine` object.
#' @param ... Unused.
#' @return A data frame containing the final evaluation performance stored in
#'   `x$performance`. Each row represents a model/metric combination and reports
#'   the evaluation data set, method, metric, estimate, confidence limits when
#'   available, and the confidence-interval method.
#' @export
as.data.frame.r4vn_machine <- function(x,...) x$performance

# Nested-CV helper used only when `nested = TRUE`.
.r4vn_machine_nested_cv <- function(raw,outcome_name,predictor_names,task,outcome_info,method,grid,
                                    preprocess,missing,missing_max,encode,standardize,transform,outlier,corr,
                                    feature,degree,reduce,variance,select,nfeatures,balance,balance_target,neighbors,folds,metric,seed) {
  outer<-.r4vn_machine_make_folds(raw[[outcome_name]],folds,1L,task,seed);rows<-list();ii<-0L
  for(s in outer){
    otr<-raw[s$train,,drop=FALSE];ova<-raw[s$val,,drop=FALSE];scores<-rep(NA_real_,length(grid))
    for(g in seq_along(grid)){
      inn<-tryCatch(.r4vn_machine_cv_one(otr,outcome_name,predictor_names,task,outcome_info,method,grid[[g]],preprocess,missing,missing_max,encode,standardize,transform,outlier,corr,feature,degree,reduce,variance,select,nfeatures,balance,balance_target,neighbors,min(5L,max(2L,folds-1L)),1L,seed+s$fold*1000L+g*20L),error=function(e)NULL)
      if(!is.null(inn))scores[g]<-.r4vn_machine_metric_from_cv(inn,task,metric)
    }
    good<-which(is.finite(scores));if(!length(good))next;gsel<-if(.r4vn_machine_metric_direction(metric)=="max")good[which.max(scores[good])] else good[which.min(scores[good])]
    fit<-tryCatch(.r4vn_machine_fit_pipeline(otr,outcome_name,predictor_names,task,outcome_info,method,grid[[gsel]],preprocess,missing,missing_max,encode,standardize,transform,outlier,corr,feature,degree,reduce,variance,select,nfeatures,balance,balance_target,neighbors,seed=seed+s$fold*4000L),error=function(e)NULL);if(is.null(fit))next
    pr<-tryCatch(.r4vn_machine_predict_pipeline(fit,ova,outcome_name),error=function(e)NULL);if(is.null(pr))next;truth<-ova[[outcome_name]][pr$rows];ii<-ii+1L
    if(task=="binary")rows[[ii]]<-data.frame(truth=as.numeric(truth),prob=pr$prob,row_id=s$val[pr$rows],rep_id=1L,fold=s$fold)
    else if(task=="regression")rows[[ii]]<-data.frame(truth=as.numeric(truth),value=pr$value,row_id=s$val[pr$rows],rep_id=1L,fold=s$fold)
    else {pp<-as.data.frame(pr$prob,check.names=FALSE);names(pp)<-paste0("prob__",names(pp));rows[[ii]]<-cbind(data.frame(truth=as.character(truth),row_id=s$val[pr$rows],rep_id=1L,fold=s$fold),pp)}
  }
  if(!length(rows))return(NULL);do.call(rbind,rows)
}