---
title: "Introduction to densemlp"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to densemlp}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

```{r setup}
library(densemlp)
```

## What densemlp is

`densemlp` is a compact dense feedforward neural network (multilayer perceptron) for regression,
classification and survival analysis on tabular data. Forward propagation, backpropagation, and
Adam optimization are implemented natively in C++ via `RcppArmadillo` - there is no
`torch`/`libtorch` dependency, which is what makes it fast to install and fast to train.

Hidden layers are `Linear -> BatchNorm -> ReLU -> [gate] -> [dropout] -> [residual]` (set
`batch_norm = FALSE` to drop the normalization step), with a plain `Linear -> task activation`
output layer (linear for regression, sigmoid for binary classification, softmax for multiclass,
a linear risk score for survival). Numeric predictors are centered/scaled and categorical
predictors are one-hot encoded internally, so `x` can be a mixed-type data frame.

## Fitting a single model

```{r}
set.seed(1)
n <- 300
x <- data.frame(a = rnorm(n), b = rnorm(n), c = rnorm(n))
y <- 2 * x$a - x$b + 0.5 * x$a * x$c + rnorm(n, sd = 0.2)

train <- sample.int(n, floor(0.8 * n))
test <- setdiff(seq_len(n), train)

fit <- densemlp(x[train, ], y[train], hidden_units = c(32, 16), epochs = 60)
pred <- predict(fit, x[test, ])
sqrt(mean((pred - y[test])^2))
```

`task = "auto"` (the default) infers `"regression"` from a numeric outcome, or `"binary"` /
`"multiclass"` from a factor/character outcome with 2 or more than 2 levels respectively.

```{r}
y_class <- factor(ifelse(y > median(y), "high", "low"))
fit_class <- densemlp(x[train, ], y_class[train], hidden_units = c(32, 16), epochs = 60)
predict(fit_class, x[test, ], type = "prob")[1:5, ]
predict(fit_class, x[test, ], type = "class")[1:5]
```

## Architecture options

All accuracy-oriented options are opt-in and default to off, so a plain `densemlp(x, y)` call keeps
the smallest, fastest architecture.

```{r}
fit_opts <- densemlp(
  x[train, ], y[train], hidden_units = c(32, 16),
  residual = TRUE,          # learned skip connection per hidden block
  gated = TRUE,              # learned sigmoid gate per hidden block
  dropout = 0.05,            # inverted dropout per hidden layer
  batch_norm = TRUE,         # batch normalization inside each hidden block
  input_projection = 8,      # linear map of the encoded inputs before layer 1
  ema_decay = 0.99,          # moving-average weights for eval/final predictions
  lr_schedule = "cosine",
  epochs = 60
)
```

`input_projection = k` inserts a bare linear layer (no activation, no batch normalization) that
maps the encoded predictors to `k` dimensions before the first hidden block; it cannot be
combined with `interaction`.

`interaction = TRUE` (an efficient learned cross-feature layer) is also available, but is
currently numerically unstable on small-n/high-p regression data and is not recommended - see
`NEWS.md` for details. It is left out of `tune_densemlp()`'s search grid for the same reason.

## Internal ensembles

Setting `ensemble > 1` fits several bootstrap-resampled members internally and averages their
predictions transparently - `predict()` still returns a single vector/matrix. Use `ncores` to
fit members in parallel (via `parallel::mclapply` on Unix-alikes; falls back to serial on
Windows).

```{r}
fit_ens <- densemlp(
  x[train, ], y[train], hidden_units = c(32, 16), epochs = 60,
  ensemble = 5, ncores = 2
)
sqrt(mean((predict(fit_ens, x[test, ]) - y[test])^2))
```

## Cross-validation

`cv_densemlp()` fits `densemlp()` on each of `folds` splits and scores the held-out fold with
`densemlp_metrics()` (RMSE/NRMSE/R² for regression; accuracy/balanced accuracy/macro AUC/log loss for
classification). Extra arguments are forwarded to every fold's `densemlp()` call.

```{r}
cv <- cv_densemlp(x, y, folds = 5, hidden_units = c(32, 16), epochs = 40)
cv
```

## Hyperparameter tuning

`tune_densemlp()` grid-searches over `hidden_units`, `dropout`, `residual`, `gated`, `ema_decay`,
`lr_schedule`, `epochs`, `batch_size`, and `lr`. Each candidate is repeated (`repeats`) with
successive seeds and ranked by mean best-epoch validation loss; the best configuration is refit on
the full data by default.

```{r}
tuned <- tune_densemlp(
  x, y, repeats = 2,
  grid = list(
    hidden_units = list(c(16), c(32, 16)),
    lr = c(1e-3, 3e-3)
  )
)
tuned$best_config
```

## Survival outcomes

With `task = "survival"` the response is a `survival::Surv(time, event)` object (or a two-column
`(time, event)` matrix). Two losses are available: `loss = "cox"` (the default) trains a single
linear risk score with a batch-wise Breslow-tie Cox partial likelihood, and `loss = "brier"`
trains a discrete-time hazard head against the IPCW integrated Brier score.

```{r, eval = requireNamespace("survival", quietly = TRUE)}
library(survival)
data(lung, package = "survival")
lung <- na.omit(lung[, c("time", "status", "age", "sex", "ph.ecog", "ph.karno", "wt.loss")])
sy <- Surv(lung$time, lung$status == 2)
sx <- lung[, c("age", "sex", "ph.ecog", "ph.karno", "wt.loss")]

sfit <- densemlp(sx, sy, task = "survival", hidden_units = c(16, 8), epochs = 60)
risk <- predict(sfit, sx, type = "response")
densemlp_metrics(sy, risk, task = "survival")
```

## Variable importance

`perm_importance()` gives model-agnostic permutation importance: each predictor is shuffled in
turn and the drop in `densemlp_metrics()` performance is recorded.

```{r}
imp <- perm_importance(fit, x[test, ], y[test])
imp$data
```

## Visualizing training

```{r, fig.width = 6, fig.height = 4}
plot(fit)          # or, equivalently, plot_history(fit)
```

## Learn more

See `?densemlp`, `?cv_densemlp`, `?tune_densemlp`, `?densemlp_metrics` and `?perm_importance` for
full argument documentation, and `NEWS.md` for the changelog.