## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  message  = FALSE,
  warning  = FALSE
)

## ----libs---------------------------------------------------------------------
library(psvr)
library(parsnip)
library(rsample)
library(recipes)
library(workflows)
library(tune)
library(dials)
library(yardstick)
library(dplyr)

## ----data---------------------------------------------------------------------
set.seed(42)
n   <- 200
x1  <- runif(n, -3, 3)
x2  <- runif(n, -3, 3)
y   <- 2 + x1^2 + 0.5 * x2^2 + rnorm(n, sd = 0.1)
dat <- data.frame(y = y, x1 = x1, x2 = x2)

## ----split--------------------------------------------------------------------
set.seed(1)
split <- initial_split(dat, prop = 0.75)
train <- training(split)
test  <- testing(split)

## ----recipe-------------------------------------------------------------------
rec <- recipe(y ~ x1 + x2, data = train) |>
  step_normalize(all_predictors())

## ----spec---------------------------------------------------------------------
spec <- psvr_rmspe_rbf(cost = tune(), rbf_sigma = tune()) |>
  set_engine("psvr")

## ----workflow-----------------------------------------------------------------
wf <- workflow() |>
  add_recipe(rec) |>
  add_model(spec)

## ----tune---------------------------------------------------------------------
set.seed(2)
folds <- vfold_cv(train, v = 5)

# Data-driven rbf_sigma range centred on median pairwise distance
train_baked      <- rec |> prep() |> bake(new_data = train)
rbf_sigma_custom <- rbf_sigma_psvr_data(train_baked |> select(-y))

wf_params <- extract_parameter_set_dials(wf) |>
  update(
    cost      = cost_psvr_ls_data(train$y),
    rbf_sigma = rbf_sigma_custom
  )

tune_res <- tune_grid(
  wf,
  resamples  = folds,
  grid       = 15,
  param_info = wf_params,
  metrics    = metric_set(yardstick::mape)
)

## ----tune-results-------------------------------------------------------------
collect_metrics(tune_res)[, c("cost", "rbf_sigma", "mean", "std_err")]

## ----best---------------------------------------------------------------------
best_params <- select_best(tune_res, metric = "mape")
best_params

## ----final--------------------------------------------------------------------
final_wf  <- finalize_workflow(wf, best_params)
final_fit <- last_fit(final_wf, split, metrics = metric_set(yardstick::mape))

collect_metrics(final_fit)

## ----predictions--------------------------------------------------------------
preds <- collect_predictions(final_fit)
head(preds[, c(".row", "y", ".pred")])

## ----new-data-----------------------------------------------------------------
new_obs <- data.frame(x1 = c(0, 1, -2), x2 = c(0, 1, 2))
predict(extract_workflow(final_fit), new_data = new_obs)

## ----engine-fit---------------------------------------------------------------
# extract_fit_engine() unwraps the parsnip/workflow layer to the raw psvr
# object -- the same class psvr_rmspe() returns when called directly
engine_fit <- extract_fit_engine(extract_workflow(final_fit))
print(engine_fit)

## ----engine-coef--------------------------------------------------------------
cf <- coef(engine_fit)
# alpha:        N dual variables; weight each training point in
#               f(x) = sum_k alpha_k K(x_k, x) + b
# b:            bias / intercept term
# support_data: all N training inputs (LS-SVR has no sparsity — every training
#               point contributes, so despite the name this is not a subset)
cat(sprintf("b = %.4f  |  alpha range: [%.4f, %.4f]\n",
            cf$b, min(cf$alpha), max(cf$alpha)))