---
title: "soilVAE Workflow"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{soilVAE Workflow}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4)
```

This vignette predicts soil carbon from NIR spectra with **soilVAE** and
compares it with a **PLS** baseline.

The data are `NIRsoil` from the **prospectr** package. They are loaded from
**prospectr** when the vignette runs; **soilVAE** itself contains code only and
bundles no data.

> The model runs in Python. Install the backend once with `vae_install()`;
> soilVAE then finds it in every session. The model-fitting chunks below run
> only when `vae_available()` is `TRUE` (it never installs anything), and the
> whole vignette needs the suggested packages **prospectr** and **pls**.

```{r packages, message=FALSE}
library(soilVAE)
has_pkgs <- requireNamespace("prospectr", quietly = TRUE) &&
  requireNamespace("pls", quietly = TRUE)
has_tf <- has_pkgs && vae_available()
knitr::opts_chunk$set(eval = has_pkgs)
# The tuning example fits several models; it runs only outside CRAN checks.
run_tuning <- has_tf && identical(Sys.getenv("NOT_CRAN"), "true")
```

```{r attach, message=FALSE}
library(prospectr)
library(pls)
```

## Data

`NIRsoil` holds 825 soil samples with NIR absorbance spectra (1100–2498 nm)
and laboratory values. The target here is total carbon (`Ciso`, %). The data
come with a calibration/test split (`train`).

```{r data}
data("NIRsoil", package = "prospectr")
d <- NIRsoil[is.finite(NIRsoil$Ciso), ]
table(set = ifelse(d$train == 1, "calibration", "test"))
```

## Preprocessing

Standard normal variate (SNV) removes scatter. No further scaling is needed:
`vae_fit()` centres each band and divides all bands by one overall SD (the
default `scale_x = "global"`, which keeps the shape of the spectrum), and
stores these statistics in the model, so `vae_predict()` scales new spectra the
same way. (PLS centres the data itself.)

```{r preprocess}
X <- standardNormalVariate(d$spc)
cal <- d$train == 1

y <- d$Ciso
X_cal <- X[cal, ];  y_cal <- y[cal]
X_te  <- X[!cal, ]; y_te  <- y[!cal]

wl <- as.numeric(colnames(d$spc))
matplot(wl, t(standardNormalVariate(d$spc[1:30, ])), type = "l", lty = 1,
        col = "grey40", xlab = "Wavelength (nm)", ylab = "SNV absorbance")
```

## Metrics

```{r metrics}
metrics <- function(obs, pred) {
  rmse <- sqrt(mean((obs - pred)^2))
  c(R2   = 1 - sum((obs - pred)^2) / sum((obs - mean(obs))^2),
    RMSE = rmse,
    RPIQ = IQR(obs) / rmse)
}
```

## Baseline: PLS

The number of components is chosen by 10-fold cross-validation on the
calibration set.

```{r pls}
set.seed(1)
pls_fit <- plsr(y_cal ~ X_cal, ncomp = 30, validation = "CV")
ncomp <- which.min(RMSEP(pls_fit)$val[1, 1, -1])
pls_pred <- predict(pls_fit, X_te, ncomp = ncomp)[, 1, 1]
round(metrics(y_te, pls_pred), 3)
```

## soilVAE

A fifth of the calibration samples is held out for early stopping: training
stops when the **prediction** error on these samples stops improving, and the
best weights are kept. Early stopping starts after a warm-up of `patience`
epochs. (Passing `validation_split = 0.2` instead of
`X_val`/`y_val` does the same split inside `vae_fit()`.)

```{r vae, eval = has_tf}
set.seed(1)
va <- sample(nrow(X_cal), round(0.2 * nrow(X_cal)))

model <- vae_build(input_dim = ncol(X_cal), seed = 1L)
vae_fit(model, X_cal[-va, ], y_cal[-va], X_cal[va, ], y_cal[va],
        epochs = 300L, patience = 20L)

vae_pred <- vae_predict(model, X_te)
round(metrics(y_te, vae_pred), 3)
```

The training history shows the validation prediction loss that early stopping
watched:

```{r history, eval = has_tf}
h <- vae_history(model)
plot(h$val_y_loss, type = "l", xlab = "Epoch", ylab = "Validation y loss (standardised)")
```

### Tuning

`tune_vae_train_val()` fits a grid of configurations on a train/validation
split; `select_best_from_grid()` picks one. Each configuration gets its own
seed, so a tuning run is reproducible. (This chunk is skipped during CRAN
checks because it fits several models.)

```{r tune, eval = run_tuning}
grid <- data.frame(latent_dim = c(16L, 32L), dropout = 0.1, lr = 5e-4,
                   beta_kl = 1e-3, alpha_y = 1, epochs = 300L,
                   batch_size = 64L, patience = 20L)
grid$hidden_enc <- list(c(512L, 256L), c(512L, 256L))
grid$hidden_dec <- list(c(256L, 512L), c(256L, 512L))

tuned <- tune_vae_train_val(X_cal[-va, ], y_cal[-va], X_cal[va, ], y_cal[va],
                            seed = 1, grid_vae = grid)
tuned$tuning_df[, c("latent_dim", "R2_val", "RMSE_val")]
select_best_from_grid(tuned$tuning_df, "rmse")$best$latent_dim
```

## Comparison on the test set

```{r compare, eval = has_tf}
round(rbind(PLS = metrics(y_te, pls_pred), soilVAE = metrics(y_te, vae_pred)), 3)
plot(y_te, vae_pred, xlab = "Observed Ciso (%)", ylab = "Predicted Ciso (%)",
     pch = 16, col = "grey30")
abline(0, 1, lty = 2)
```

## Saving a model

A fitted model, including the predictor and target scaling, is saved to one
`.keras` file:

```{r save, eval = has_tf}
f <- tempfile(fileext = ".keras")
vae_save(model, f)
model2 <- vae_load(f)
all.equal(vae_predict(model2, X_te), vae_pred)
```