---
title: "tidymodels workflow with psvr"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{tidymodels workflow with psvr}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

This vignette demonstrates the full tidymodels pipeline with **psvr**: data
splitting, preprocessing, hyperparameter tuning by cross-validation, and final
model evaluation.  We use `psvr_rmspe_rbf()` (LS-SVR with RMSPE loss, RBF
kernel) and tune the regularisation parameter `cost` ($\Gamma$) against MAPE.

```{r libs}
library(psvr)
library(parsnip)
library(rsample)
library(recipes)
library(workflows)
library(tune)
library(dials)
library(yardstick)
library(dplyr)
```

## Data

The synthetic even-function dataset from the package README:
$y = 2 + x_1^2 + 0.5\,x_2^2 + \varepsilon$, $\varepsilon \sim \mathcal{N}(0,\,0.1^2)$.
Targets are strictly positive by construction ($y > 0$).

```{r 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)
```

## 1 — Split

```{r split}
set.seed(1)
split <- initial_split(dat, prop = 0.75)
train <- training(split)
test  <- testing(split)
```

## 2 — Preprocessing recipe

Centre and scale all predictors so the RBF kernel operates on a
standardised feature space.

```{r recipe}
rec <- recipe(y ~ x1 + x2, data = train) |>
  step_normalize(all_predictors())
```

## 3 — Model spec with `tune()`

Both `cost` (maps to $\Gamma$) and `rbf_sigma` (the RBF bandwidth $\sigma$)
are `tune()` placeholders; the grid search will explore all combinations.

```{r spec}
spec <- psvr_rmspe_rbf(cost = tune(), rbf_sigma = tune()) |>
  set_engine("psvr")
```

## 4 — Workflow

```{r workflow}
wf <- workflow() |>
  add_recipe(rec) |>
  add_model(spec)
```

## 5 — Tune with 5-fold CV

We search over a 15-point Latin hypercube of `cost` and `rbf_sigma` values and
evaluate each fold by MAPE.  Both search ranges are set from the data, and
**neither happens automatically** — `dials` cannot finalize either one, so both
have to be passed explicitly through `param_info`:

- `rbf_sigma_psvr_data()` centres the bandwidth range on the median pairwise
  distance, so it is computed on the **baked** predictors: the heuristic only
  means anything on the scale the model is actually fitted on.
- `cost_psvr_ls_data()` widens the `cost` range.  Here `cost` is $\Gamma$, and
  its registered default of $[-2, 10]$ on the log2 scale ($\Gamma \le 1024$) is
  the $\epsilon$-SVR range — far too low for LS-SVR, where the optimum scales
  with `var(y) * n`.  Left at the default the grid tops out at $\Gamma = 1024$,
  which on this dataset is an order of magnitude below the value selected once
  the range is widened — compare the `cost` column printed below.  The search
  is boundary-trapped: it cannot reach the optimum at all, and nothing warns
  you, because every candidate it did evaluate was legal.
  This one cannot be automated even in principle,
  because `tune` finalizes parameters from the predictors alone and never passes
  the outcome to `dials::finalize()`.

The RMSPE LS-SVR only solves an `(N+1) × (N+1)`
linear system — no iterative solver is involved — so 75
fits complete in seconds.

```{r 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)
)
```

Cross-validated MAPE for each candidate (lower is better):

```{r tune-results}
collect_metrics(tune_res)[, c("cost", "rbf_sigma", "mean", "std_err")]
```

## 6 — Select best

```{r best}
best_params <- select_best(tune_res, metric = "mape")
best_params
```

## 7 — Final fit and test-set evaluation

`last_fit()` refits on the full training set with the chosen `cost` and
evaluates once on the held-out test data.

```{r 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 on the test set:

```{r predictions}
preds <- collect_predictions(final_fit)
head(preds[, c(".row", "y", ".pred")])
```

The fitted workflow can also be used directly for new data:

```{r 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)
```

## 8 — Inspecting the fitted psvr model

The tidymodels layer wraps a `psvr_rmspe` object (returned by the engine fit wrapper `psvr_rmspe_rbf_fit()`). Extract it to use `print()` and `coef()` directly.

```{r 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)
```

```{r 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)))
```