---
title: "Simulating Trajectories"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Simulating Trajectories}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  fig.width  = 6,
  fig.height = 5
)
```

```{r load}
if (requireNamespace("pkgload", quietly = TRUE)) {
  pkgload::load_all("..", export_all = FALSE, helpers = FALSE, quiet = TRUE)
} else if (requireNamespace("radiatR", quietly = TRUE)) {
  library(radiatR)
} else {
  stop("Package 'radiatR' not installed and 'pkgload' not available.")
}
library(ggplot2)
```

## Overview

`simulate_tracks()` generates synthetic trajectories in circular space without
any tracking files. Its main uses are:

- Verifying that a new analysis pipeline works before real data arrives.
- Reproducing plots and examples in documentation.
- Teaching: demonstrating how concentration, tortuosity, and directional bias
  each affect the appearance of paths.
- Informal power analysis: generate data under a plausible effect size and check
  whether the expected statistical pattern is visible.

---

## Quick Start

Called with no arguments, `simulate_tracks()` returns a tibble of three
conditions (60 trials total, 200 frames each):

```{r quick-start}
sim <- simulate_tracks(seed = 1)
dim(sim)
names(sim)
```

Each row is one frame. The key columns are:

| Column | Description |
|---|---|
| `condition` | Condition label |
| `trial_id` | Unique trial identifier (`<condition>_<index>`) |
| `frame` | Frame number within the trial |
| `predictor` | Per-trial continuous covariate value |
| `concentration` | Effective von Mises κ for the trial |
| `tortuosity` | Effective angular noise σ for the trial |
| `final_heading` | Drawn heading (unit-circle radians) |
| `abs_x`, `abs_y` | Absolute position on the unit disc |
| `rel_x`, `rel_y` | Position re-centred on the heading direction |
| `abs_theta`, `rel_theta` | Corresponding angular coordinates |

---

## Output Formats

The `output` argument controls what is returned.

```{r output-formats, eval=FALSE}
# Default: long-form tibble
sim_tbl  <- simulate_tracks(seed = 1, output = "tibble")

# Tracks directly (ready for radiate(), derive_headings(), etc.)
sim_ts   <- simulate_tracks(seed = 1, output = "trajset")

# Both representations in a named list
sim_both <- simulate_tracks(seed = 1, output = "both")
names(sim_both)  # "tibble" "trajset"
```

The `output = "trajset"` path wraps the tibble in a `Tracks` (absolute
coordinates, no normalisation). Use it whenever you want to plug straight into
the plotting or analysis functions.

---

## Plotting the Default Simulation

```{r plot-default}
ts_default <- simulate_tracks(seed = 1, output = "trajset")
radiate(ts_default,
        group_col    = "trial_id",
        colour_cycle = 10,
        panel_by     = "condition",
        ncol         = 3,
        show_labels  = FALSE,
        show_arrow   = TRUE)
```

The three default conditions differ in directional bias, concentration, and
tortuosity (see below).

---

## The Conditions Table

All simulation behaviour is controlled by a data frame passed to `conditions`.
When `conditions = NULL` the function fills in a three-condition template.
Supply your own data frame to override any subset of columns; missing columns
receive sensible defaults.

| Column | Default | Description |
|---|---|---|
| `condition` | `"condition_N"` | Label for the condition |
| `n_trials` | `10` | Number of trajectories |
| `ref_mean` | `0` | Baseline reference direction (unit-circle radians) |
| `concentration_base` | `5` | Baseline von Mises κ |
| `concentration_slope` | `0` | Slope of κ on the per-trial predictor |
| `tortuosity_base` | `0.06` | Baseline angular noise σ |
| `tortuosity_slope` | `0` | Slope of σ on the predictor |
| `tortuosity_sd` | `0.01` | Trial-to-trial variability in σ |
| `predictor_mean` | `0` | Mean of the per-trial predictor distribution |
| `predictor_sd` | `0.2` | SD of the predictor distribution |
| `predictor_values` | (none) | Optional list-column of explicit predictor values |

---

## Concentration (κ)

`concentration_base` is the von Mises κ passed to `rvonmises()`. Higher values
produce tightly clustered headings; values near zero produce near-uniform
heading distributions.

```{r concentration, fig.width=9, fig.height=4}
conds_kappa <- data.frame(
  condition          = c("low (κ=0.1)", "medium (κ=1)", "high (κ=12)"),
  n_trials           = 20L,
  concentration_base = c(0.1, 1, 12),
  tortuosity_base    = 0.05
)
ts_kappa <- simulate_tracks(conditions = conds_kappa, seed = 42,
                             output = "trajset")
radiate(ts_kappa,
        group_col    = "trial_id",
        colour_cycle = 10,
        panel_by     = "condition",
        ncol         = 3,
        show_labels  = FALSE,
        show_arrow   = TRUE)
```

At κ = 1 the mean arrow is short (low resultant length); at κ = 12 nearly all
trials point in the same direction.

---

## Tortuosity and Path Smoothness

Two parameters control within-trial path shape:

- `tortuosity_base` — the SD of the per-step angular noise. Larger values
  produce more sinuous paths.
- `phi` — autocorrelation of successive angular deviations (0 = white noise,
  values near 1 produce smooth, sweeping curves). Passed directly to
  `simulate_tracks()`, not via the conditions table.

```{r tortuosity, fig.width=9, fig.height=4}
conds_tort <- data.frame(
  condition          = c("smooth (φ=0.95)", "moderate (φ=0.7)", "noisy (φ=0)"),
  n_trials           = 15L,
  concentration_base = 6,
  tortuosity_base    = c(0.12, 0.12, 0.12)
)
phi_vals <- c(0.95, 0.7, 0)

# Simulate each phi separately and bind
parts <- lapply(seq_along(phi_vals), function(i) {
  d <- conds_tort[i, , drop = FALSE]
  simulate_tracks(conditions = d, n_points = 200,
                  phi = phi_vals[i], seed = 7 + i)
})
ts_tort <- tracks(
  do.call(rbind, parts),
  id = "trial_id", time = "frame",
  x = "abs_x", y = "abs_y", angle = "abs_theta",
  angle_unit = "radians", normalize_xy = FALSE
)
radiate(ts_tort,
        group_col    = "trial_id",
        colour_cycle = 10,
        panel_by     = "condition",
        ncol         = 3,
        show_labels  = FALSE,
        show_arrow   = FALSE)
```

---

## Directional Bias (ref_mean)

`ref_mean` shifts the reference direction away from East (0 rad). Use it to
simulate experiments where the reference is not along the positive x-axis.

```{r stim-mean, fig.width=9, fig.height=4}
conds_bias <- data.frame(
  condition          = c("East (0°)", "North (90°)", "Northwest (135°)"),
  n_trials           = 20L,
  concentration_base = 5,
  tortuosity_base    = 0.06,
  ref_mean          = c(0, pi / 2, 3 * pi / 4)
)
ts_bias <- simulate_tracks(conditions = conds_bias, seed = 99,
                           output = "trajset")
radiate(ts_bias,
        group_col    = "trial_id",
        colour_cycle = 10,
        panel_by     = "condition",
        ncol         = 3,
        show_labels  = FALSE,
        show_arrow   = TRUE)
```

---

## Multi-Condition Gradient

A common design is a monotone gradient in cue strength. The example below
simulates four levels with increasing concentration and decreasing tortuosity —
mimicking an orientation experiment where stronger cues produce tighter,
smoother paths.

```{r gradient, fig.width=9, fig.height=4}
conds_grad <- data.frame(
  condition          = paste0(c(15, 30, 60, 120), "° arc"),
  n_trials           = 25L,
  concentration_base = c(1.5, 3, 6, 10),
  tortuosity_base    = c(0.14, 0.10, 0.06, 0.03),
  tortuosity_sd      = 0.015
)
ts_grad <- simulate_tracks(conditions = conds_grad, n_points = 200,
                            seed = 55, output = "trajset")
radiate(ts_grad,
        group_col    = "trial_id",
        colour_col   = "condition",
        panel_by     = "condition",
        ncol         = 4,
        show_labels  = FALSE,
        show_arrow   = TRUE)
```

---

## Continuous Predictor

When `concentration_slope` or `tortuosity_slope` is non-zero, the final kappa
and tortuosity of each trial depend on its sampled predictor value. This lets
you model experiments where individual differences (e.g. body size, prior
exposure) modulate the response.

```{r predictor}
conds_pred <- data.frame(
  condition            = "gradient",
  n_trials             = 40L,
  concentration_base   = 3,
  concentration_slope  = 4,    # kappa rises with predictor
  tortuosity_base      = 0.10,
  tortuosity_slope     = -0.04, # smoother at high predictor
  predictor_mean       = 0,
  predictor_sd         = 1
)
sim_pred <- simulate_tracks(conditions = conds_pred, seed = 7, output = "both")

# Each trial gets its own predictor sample
range(sim_pred$tibble$predictor)
range(sim_pred$tibble$concentration)
```

Plotting concentration against the predictor confirms the linear relationship:

```{r predictor-plot, fig.width=5, fig.height=3}
trial_summary <- unique(sim_pred$tibble[, c("trial_id", "predictor",
                                             "concentration")])
ggplot(trial_summary, aes(predictor, concentration)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = FALSE, colour = "steelblue") +
  labs(x = "Predictor", y = "Effective κ") +
  theme_minimal()
```

---

## Connecting to Circular Analysis

Simulated `Tracks` objects feed directly into `derive_headings()` and
`circ_summarise()`, making it straightforward to run the same analysis pipeline
on synthetic data.

```{r analysis}
conds_analysis <- data.frame(
  condition          = c("diffuse", "concentrated"),
  n_trials           = 30L,
  concentration_base = c(1.5, 8),
  tortuosity_base    = c(0.12, 0.04)
)
sim_analysis <- simulate_tracks(conditions = conds_analysis, seed = 21,
                                 output = "both")

hd <- derive_headings(sim_analysis$trajset, rule = "crossing",
                       circ0 = 0.2, circ1 = 0.4,
                       coords = "absolute",
                       angle_convention = "unit_circle",
                       return_coords = TRUE)
names(hd)[names(hd) == "id"] <- "trial_id"

# derive_headings returns only id/heading/coords; join condition back
cond_map <- unique(sim_analysis$tibble[, c("trial_id", "condition")])
hd <- merge(hd, cond_map, by = "trial_id")

compute_circ_mean(hd, colour_col = "condition")[,
  c("condition", "mean_dir", "resultant_R")]
```

Visualise the headings alongside the trajectories:

```{r analysis-plot}
p <- radiate(sim_analysis$trajset,
             group_col    = "trial_id",
             colour_col   = "condition",
             panel_by     = "condition",
             ncol         = 2,
             show_labels  = FALSE,
             show_arrow   = FALSE) +
  add_heading_points(hd, colour_col = "condition", size = 2, alpha = 0.7) +
  add_heading_arrow(hd, colour_col = "condition", colour = "black")
p
```

---

## Known Modality and Shape

`simulate_tracks()` can also attach a known directional ground truth, which is
handy for sanity-checking the circular statistics: you know what the analysis
*should* recover. The `modality` of a condition sets how its per-trial principal
headings are distributed — `uniform` (no preferred direction), `unimodal` (one
von Mises mode), `axial` (two antipodal modes), or `multimodal` (`n_modes`
evenly spaced modes). Here one condition per modality, fed through a directional
heading rule (`net`) and `circ_model_select()`, recovers the generating model:

```{r modality-recovery}
conds <- tibble::tibble(
  condition          = c("uniform", "unimodal", "axial"),
  n_trials           = 60L,
  ref_mean           = 0.4,
  concentration_base = 6,
  modality           = c("uniform", "unimodal", "axial")
)
sim <- simulate_tracks(conditions = conds, n_points = 80,
                       output = "both", seed = 7)

# net heading per trial, with the condition label joined back on
hd <- derive_headings(sim$trajset, rule = "net")
names(hd)[names(hd) == "id"] <- "trial_id"
hd <- merge(hd, unique(sim$tibble[, c("trial_id", "condition")]),
            by = "trial_id")

# best-supported model per condition (AICc; lowest dAICc = best)
ms <- circ_model_select(hd, group_col = "condition")
do.call(rbind, lapply(split(ms, ms$condition), function(d) {
  d[which.min(d$AICc), c("condition", "model", "weight")]
}))
```

The `unimodal` and `axial` conditions are picked out cleanly (Akaike weight near
1); the `uniform` condition leaves *uniform* and *unimodal* in a near-tie, as it
should when there is no preferred direction to fit.

The `track_shape` of a condition sets the *within-track* geometry instead.
A `track_shape = "oscillatory"` condition produces back-and-forth tracks along a
principal axis — the bidirectional pattern the per-track axial methods are built
for. The oscillatory tracks form a genuinely axial position cloud, so the
position-based estimator `pca_axis` recovers the axis at the default sampling
density, while a directional rule (`net`) sees the to-and-fro directions cancel:

```{r oscillatory-axis}
osc <- tibble::tibble(condition = "osc", n_trials = 40L, ref_mean = 0.6,
                      concentration_base = 50, track_shape = "oscillatory",
                      n_reversals = 4L)
ts_osc <- simulate_tracks(conditions = osc, output = "trajset", seed = 11)

# pca_axis recovers the ~0.6 rad axis (read it as an axial mean) at defaults
pa <- derive_headings(ts_osc, rule = "pca_axis")
circ_summarise(pa, "heading", axial = TRUE, units = "radians",
               stats = "mean_dir")

# net cannot: the directions cancel to a tiny resultant length
net <- derive_headings(ts_osc, rule = "net")
circ_summarise(net, "heading", units = "radians", stats = "resultant_R")
```

The oscillatory line-width (`line_width`) is a fixed fraction of the amplitude,
independent of `tortuosity_base`, so the axis is recoverable at the default
settings used here without any signal-to-noise tuning. Note the position-based
methods (`pca_axis`, `ransac_straight`) are the robust choice: the step-based
`velocity_axis` recovers the same axis only at coarse sampling, because its
per-step axial signal shrinks with `n_points` while the perpendicular jitter
step does not, so dense sampling lets the jitter flip the estimate toward the
perpendicular.

---

## Reproducibility and Saving

Pass `seed` for exact reproducibility. Pass `write_path` to write the tibble
to a CSV alongside (or instead of) returning it — useful when you want to
hand the data to a colleague or another tool.

```{r save, eval=FALSE}
simulate_tracks(
  conditions = conds_grad,
  n_points   = 200,
  seed       = 55,
  write_path = "sim_gradient.csv"
)
```