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

## ----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")

## ----attach, message=FALSE----------------------------------------------------
library(prospectr)
library(pls)

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

## ----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------------------------------------------------------------------
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)
}

## ----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)

## ----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)

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

## ----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

## ----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)

## ----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)