---
title: "Getting Started with scanr"
format:
  html:
    toc: true
    toc-depth: 2
    number-sections: true
vignette: >
  %\VignetteIndexEntry{Getting Started with scanr}
  %\VignetteEngine{quarto::html}
  %\VignetteEncoding{UTF-8}
---

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

# Overview

The `scanr` package detects multiple change points in long univariate time series
using window-based SCAN statistics. The method compares local windows along a
series, calibrates evidence for changes using resampling, and combines
information across several window sizes.

This vignette demonstrates a typical workflow:

1. simulate a univariate time series with known change points;
2. run `scan_cpd()` with multiple window sizes;
3. compare estimated change points with the known truth;
4. inspect the result using diagnostic visualizations;
5. repeat the workflow for distributional changes;
6. apply the workflow to hourly Bitcoin closing prices;
7. use the SWAL statistic for single change-point localization.


# Load the Package

Users can install the package from CRAN and load it as follows:

```{r}
#| eval: false
install.packages("scanr")
```

```{r}
library(scanr)
```

Alternatively, the development version of `scanr` can be installed directly
from GitHub:

```{r}
#| eval: false
if (!requireNamespace("pak", quietly = TRUE)) { install.packages("pak") }

pak::pak("Prabashoka/scanr")
```

# Example 1 - Changes in Mean

We first simulate a white-noise time series of length 20,000 with 20 change
points and 21 piecewise-constant mean regimes. The noise variance is fixed, so
this example focuses on changes in the mean.

```{r}
set.seed(1234)

n <- 20000

change_points <- c(
  952, 1905, 2858, 3810, 4763, 5715, 6668, 7620, 8573, 9525,
  10478, 11430, 12383, 13335, 14288, 15240, 16193, 17145, 18098,
  19050
)

means <- c(
  0, 2, -1, 3, 0.5, -2, 2, 5, -0.5, 2.5, 0, -2.5, -1.5, 1.5,
  3, 1, 0, 1.25, -2, 3.5, -1.5
)

segment_starts <- c(1L, change_points + 1L)
segment_ends <- c(change_points, n)
x_mean <- numeric(n)

for (j in seq_along(means)) {
  segment_index <- segment_starts[j]:segment_ends[j]
  x_mean[segment_index] <- rnorm(length(segment_index), mean = means[j], sd = 1)
}

change_points
```

The plot below shows the simulated series.

```{r}
#| echo: false
#| warning: false
#| fig.width: 10
#| fig.height: 4
vis_time_series(
  x_mean,
  x_label = "Time",
  y_label = "Value",
  title = ""
)
```

## Run SCAN Change-Point Detection

The `window_sizes` argument controls the local scales used by the scan. The
window sizes should be smaller than the spacing between nearby change points,
while still being large enough to estimate the local distributions on each side
of a candidate split. Smaller values improve localization and sensitivity to
closely spaced changes, while larger values provide more stable local estimates.

`default_window_sizes()` randomly samples the requested number of distinct
scales. By default, the largest is `floor(n^(2 / 3))`; it is always capped at
`floor(n / 2)` so that both sides of a candidate split can contain a full
window. Its optional `seed` argument makes the sampled grid reproducible
without changing the caller's random-number generator state. Supplying
problem-informed bounds is useful when the approximate minimum segment length
is known. Omitting `window_sizes` makes `scan_cpd()` use the same helper with
its `min_window`, `max_window`, and `n_windows` arguments.

```{r}
mean_window_sizes <- default_window_sizes(
  n = length(x_mean),
  min_window = 100,
  max_window = floor(length(x_mean)^(2 / 3)),
  n_windows = 7,
  seed = 52
)
mean_window_sizes

fit_mean <- scan_cpd(
  x_mean,
  window_sizes = mean_window_sizes,
  n_boot = 400,
  random_state = 1234,
  change_type = "mean",
  n_jobs = -1
)

fit_mean
```

::: {.callout-note}
## Remark on bootstrap calibration

The adaptive threshold for declaring a change point is derived through a
**tapered block bootstrap** across the combined windows. Increasing the number
of bootstrap samples (`n_boot`) produces a more accurate null distribution at
the cost of longer computation time; a value of 400–1,000 is typically
sufficient. The method supports parallel processing via
[Rust's](https://rust-lang.org/) **Rayon** library for computational
efficiency.
:::

The estimated change points can be extracted directly from the fitted object.

```{r}
fit_mean$change_points
```

## Evaluate Detection Accuracy

For simulated data, estimated change points can be compared with the known
truth. The `tolerance` argument controls how close an estimated change point
must be to a true change point to count as a match.

```{r}
cpd_metrics(
  true_cps = change_points,
  estimated_cps = fit_mean$change_points,
  n = length(x_mean),
  tolerance = 20
)
```

## Visualize the Result

The package includes visual helpers for inspecting detected change points,
window-level votes, and the overall vote scree.

```{r}
#| fig.width: 10
#| fig.height: 4.5
vis_change_points(
  x_mean,
  fit_mean,
  true_change_points = change_points,
  x_label = "Time",
  y_label = "Value"
)
```

## Tune the Voting Threshold

The ensemble vote threshold, `vote_threshold`, controls how much agreement is
required across window sizes before a candidate is retained. The default value
is 0.5, meaning that a candidate must be supported by at least half of the
available window-level votes.

```{r}
vis_vote_scree(fit_mean)
```

```{r}
#| fig.width: 10
#| fig.height: 4.5
vis_window_votes(fit_mean)
```

# Example 2 - Changes in Distribution

The next example uses segments from different distribution families, including
normal, exponential, Poisson, gamma, uniform, Student's t, lognormal, beta,
Weibull, and chi-square distributions. Each segment is standardized before
being rescaled, so the signal is not only a simple shift in the mean. This
creates a more challenging distributional change-point problem.

```{r}
set.seed(1234)

n <- 20000

change_points <- c(952, 1905, 2858, 3810, 4763, 5715, 6668, 7620, 8573, 9525, 10478, 11430, 12383, 13335, 14288, 15240, 16193, 17145, 18098, 19050)

segment_starts <- c(1, change_points + 1)
segment_ends <- c(change_points, n)

families <- c("normal", "exponential", "poisson", "t", "gamma", "uniform", "lognormal", "weibull", "chisq", "beta", "normal", "poisson", "exponential", "uniform", "gamma", "t", "lognormal", "beta", "weibull", "chisq", "normal")

scale_factors <- c(0.7, 1.4, 0.6, 1.8, 0.75, 2.5, 0.65, 3.0, 0.7, 2.6, 0.6, 1.7, 0.65, 2.9, 0.7, 2.8, 0.6, 2.6, 0.65, 3.0, 0.7)

simulate_base <- function(m, family) {
  z <- switch(family, normal = rnorm(m, mean = 0, sd = 1), exponential = rexp(m, rate = 1), poisson = rpois(m, lambda = 2), gamma = rgamma(m, shape = 2, rate = 1), uniform = runif(m, min = -sqrt(3), max = sqrt(3)), t = rt(m, df = 3), lognormal = rlnorm(m, meanlog = 0, sdlog = 0.7), beta = rbeta(m, shape1 = 2, shape2 = 5), weibull = rweibull(m, shape = 1.5, scale = 1), chisq = rchisq(m, df = 5))
  as.numeric((z - mean(z)) / sd(z))
}

x_dist <- numeric(n)

for (j in seq_along(families)) {
  segment_index <- segment_starts[j]:segment_ends[j]
  x_dist[segment_index] <-scale_factors[j] * simulate_base(length(segment_index), families[j])
}

change_points
```

The simulated series contains changes in distributional shape and scale.

```{r}
#| echo: false
#| warning: false
#| fig.width: 10
#| fig.height: 4
vis_time_series(
  x_dist,
  x_label = "Time",
  y_label = "Value",
  title = ""
)
```

## Run SCAN Change-Point Detection

For distributional changes, set `change_type = "distribution"`. This uses the
distribution-sensitive local statistic rather than a mean-only statistic.

```{r}
distribution_window_sizes <- default_window_sizes(
  n = length(x_dist),
  min_window = 100,
  max_window = floor(length(x_mean)^(2 / 3)),
  n_windows = 15,
  seed = 52
)

fit_dist <- scan_cpd(
  x_dist,
  window_sizes = distribution_window_sizes,
  n_boot = 1000,
  random_state = 1234,
  change_type = "distribution",
  vote_threshold = 0.5,
  n_jobs = -1
)
```

The set of estimated change points is:

```{r}
fit_dist$change_points
```

## Evaluate Detection Accuracy

```{r}
cpd_metrics(
  true_cps = change_points,
  estimated_cps = fit_dist$change_points,
  n = length(x_dist),
  tolerance = 20
)
```

## Visualize the Result

```{r}
#| fig.width: 10
#| fig.height: 4.5
vis_change_points(
  x_dist,
  fit_dist,
  true_change_points = change_points,
  x_label = "Time",
  y_label = "Value"
)
```

## Real Data Example

This example applies change-point detection to the PIT-502 pressure sensor in
the [Secure Water Treatment (SWaT) dataset](https://www.sutd.edu.sg/itrust/itrust-labs/datasets/dataset-characteristics/swat/).

## Load the Data

```{r}
#| echo: false
swat_path <- system.file(
  "extdata",
  "PIT-502.csv",
  package = "scanr"
)
swat <- read.csv(
  swat_path,
  check.names = FALSE,
  strip.white = TRUE
)
```
Load the dataset and resample it at one-minute intervals for change-point detection.
```{r}
head(swat)
```

## Choose Window Sizes

The largest candidate window is $\lfloor n^{2/3}\rfloor$. The package helper
selects 15 default window sizes between 20 minutes and that upper bound.

```{r}
n <- nrow(swat)
window_sizes <- default_window_sizes(
  n,
  min_window = 100,
  max_window = floor(n^(2 / 3)),
  n_windows = 10,
  seed = 400
  
)
window_sizes
```

The PIT-502 measurements are standardized before detection. Standardization
changes their units but preserves the locations of mean shifts.

```{r}
x <- as.numeric(scale(swat[["PIT502"]]))
```

## Detect change-points

The R interface uses `n_boot` and `random_state` for bootstrap calibration and
reproducibility.

```{r}
fit_swat <- scan_cpd(
  x,
  window_sizes = window_sizes,
  n_boot = 400,
  vote_threshold = 0.15,
  random_state = 100,
  change_type = "mean",
  n_jobs = -1
)

fit_swat
```

The fitted object stores change points as positions in the one-minute series.
We can map those positions back to timestamps, observed pressures, and the
dataset's operating-state labels.

```{r}
detected_changes <- swat[
  fit_swat$change_points,
  c("Timestamp", "PIT502", "Normal/Attack"),
  drop = FALSE
]

cat("Number of change-points:", nrow(detected_changes), "\n")
```

## Visualize the Result

```{r}
#| fig-width: 11
#| fig-height: 5

vis_change_points(
  x,
  fit_swat,
  index = seq_along(x),
  x_label = "Time",
  y_label = "Standardized PIT-502 pressure",
  title = "Mean changes in the SWaT PIT-502 pressure sensor"
)
```

## Inspect the Vote Scree

The vote scree shows how much support each candidate location receives across
the selected window sizes. The horizontal threshold corresponds to the
`vote_threshold = 0.15` used above.

```{r}
#| fig-width: 11
#| fig-height: 5

vis_vote_scree(fit_swat)
```

# SWAL Statistic for Single Change-Point Detection

The SCAN procedure above combines local evidence across many windows. The
package also exposes lower-level localization tools for a single region that is
believed to contain one change point.

For a mean change, the mean version of the SWAL statistic aligns with the CUSUM
localizer. Both methods search for the split that best separates the left and
right sides of the local region.

```{r}
set.seed(1234)

true_single_cp <- 150
mean_region <- c(
  rnorm(true_single_cp, mean = 0, sd = 1),
  rnorm(true_single_cp, mean = 2, sd = 1)
)

c(
  truth = true_single_cp,
  cusum = ts_cusum(mean_region),
  swal_distribution = swal_statistic(mean_region, change_type = "distribution")
)
```

The distributional version is more general. It can localize changes where the
mean is approximately unchanged but other aspects of the distribution, such as
variance, change sharply. In the example below, both segments have mean zero,
but the second segment has a much larger standard deviation.

```{r}
set.seed(1234)

var_region <- c(
  rnorm(true_single_cp, mean = 0, sd = 0.5),
  rnorm(true_single_cp, mean = 0, sd = 2)
)

c(
  truth = true_single_cp,
  cusum = ts_cusum(var_region),
  swal_distribution = swal_statistic(var_region, change_type = "distribution")
)
```

The `ts_wasserstein()` function returns both the estimated split and the full
sequence of split statistics. The plotting helper `vis_swal_curve()` displays
that localization curve, making it easier to inspect where the statistic is
maximized.

```{r}
vis_swal_curve(
  var_region,
  start = 1,
  end = length(var_region),
  x_label = "Candidate split",
  y_label = "Scaled Wasserstein statistic",
  title = ""
)
```

# Session Information

Recording session information helps make the vignette reproducible across R
versions and operating systems.

```{r}
sessionInfo()
```