---
title: "Synthetic datasets for geometric smoothing"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Synthetic datasets for geometric smoothing}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

Simulated data with a known true response are useful for checking a
smoother. This vignette shows how to describe and generate such data
reproducibly. For an overview of the whole package, see
[Finding your way around geosmooth](function-guide.html).

`geosmooth` describes a synthetic dataset by four components: the
*geometry* (the shape on which the points lie and how it sits in the
observed coordinates), the *sampling* (how points are drawn on that shape),
the *truth* (the true response at each point), and the *response* (the noise
added to the truth). [`synthetic.spec()`](https://pgajer.github.io/geosmooth/reference/synthetic.spec.html) combines the four components into a
specification without drawing any data. [`materialize.synthetic()`](https://pgajer.github.io/geosmooth/reference/materialize.synthetic.html) then draws
the data, using the seed and the random-number policy you give it, and
returns a validated `synthetic_dataset`.

```{r circle}
library(geosmooth)

circle.spec <- synthetic.spec(
  geometry = dgraphs::synthetic.circle(radius = 1, ambient.dim = 2),
  sampling = dgraphs::synthetic.sampling.grid.interval(
    0, 2 * pi, endpoints = "exclude.lower"),
  truth = synthetic.truth.polynomial(c(b0 = 0)),
  response = synthetic.response.gaussian(sd = 0.05)
)
circle.data <- materialize.synthetic(circle.spec, n = 80, seed = 2026)
circle.data
```

The circle has one latent coordinate, the angle, and two observed
coordinates, the position in the plane. The truth is deliberately constant,
so the response in this example is pure noise around zero; the example is
not meant to test whether a smoother recovers a varying signal.

```{r circle-views, fig.width=8, fig.height=4, out.width="100%", fig.cap="Left: the circle in its two observed coordinates. Right: responses versus latent angle, with the constant zero truth shown by a dashed line. These views describe the same 80 observations.", fig.alt="A circle of 80 points beside noisy responses scattered around a dashed zero line over angles from zero to two pi."}
old.par <- par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
plot(circle.data, view = "geometry", col = "#0072B2", asp = 1,
     main = "Observed geometry")
plot(circle.data, view = "response", col = "grey40", xlab = "Angle (radians)",
     main = "Response and truth")
par(old.par)
```

The second kind of geometry, `dgraphs::synthetic.quadform()`, is a surface
of a chosen intrinsic dimension inside a space of a chosen, larger ambient
dimension. Its shape is given by quadratic forms: a flat surface has none,
and a paraboloid or a saddle has one or more.

```{r quadform}
surface.spec <- synthetic.spec(
  geometry = dgraphs::synthetic.quadform(
    intrinsic.dim = 2,
    ambient.dim = 5,
    forms = list(diag(c(1, -1))),
    frame = "random.orthonormal"
  ),
  sampling = dgraphs::synthetic.sampling.uniform.box(-1, 1),
  truth = synthetic.truth.polynomial(c(b0 = 0, b1 = 1, b22 = -0.5)),
  response = synthetic.response.gaussian(sd = 0.1)
)
surface.data <- materialize.synthetic(surface.spec, n = 100, seed = 17)
dim(surface.data$latent)
dim(surface.data$predictors)
head(as.data.frame(surface.data))
```

This surface has two latent coordinates and five observed coordinates; a
random orthonormal frame places it in five dimensions. A plot of two of the
observed coordinates is a projection, so points that look close in the plot
may be far apart in the coordinates that are not shown.

```{r surface-view, fig.width=6, fig.height=4.5, out.width="80%", fig.cap="Two-coordinate projection of a surface embedded in five dimensions. Point colors indicate the known, nonconstant conditional mean; the key uses its actual numerical range.", fig.alt="A projected surface point cloud colored by its known conditional mean, with a numeric color key."}
plot(surface.data, view = "geometry", color = "truth", main = "Projected surface")
```

In this example the truth is a polynomial that varies over the surface. A
smoother is given only the observed coordinates and the noisy response, not
the latent coordinates or the truth:

```{r surface-fit}
surface.fit <- fit.lps(surface.data$predictors, surface.data$response,
  foldid = rep(1:5, length.out = surface.data$n),
  support.grid = c(20L, 30L), degree.grid = 1L,
  kernel.grid = "gaussian", coordinate.method = "local.pca", chart.dim = 2L)
surface.fit$selected
head(data.frame(response = surface.data$response, fitted = fitted(surface.fit),
                truth = surface.data$truth))
```

The fitted values are at the observed points. The truth is known only
because the data are simulated. The cross-validation score that selected the
settings is not an independent estimate of prediction accuracy. The
[first-fit guide](function-guide.html#a-first-regression-fit) explains that distinction.

The package also contains a registry of predefined specifications, listed
by [`synthetic.registry.ids()`](https://pgajer.github.io/geosmooth/reference/synthetic.registry.ids.html). The identifiers S01--S16 name signal shapes, and
the suffixes V1--V3 name combinations of sampling and noise; together they
form the test suite used for Hessian-energy regression. For example, S01.V1
is a signal with two separated, narrow peaks of equal height, sampled
uniformly, with Gaussian noise. The identifiers G1--G7 are kept from earlier
versions: G1 is a polynomial in the observed coordinates, and G2 is a
polynomial on a flat surface. The `dataset.id` of a generated dataset
includes a hash of its complete specification.

```{r registry}
ids <- c("G1", "G2", "S01.V1")
do.call(rbind, lapply(ids, function(id) {
  spec <- synthetic.registry.spec(id)
  data.frame(id, geometry = spec$geometry.spec$family,
             truth = spec$truth.spec$parameters$estimand)
}))
registered <- synthetic.registry.spec("S01.V1")
registered.data <- materialize.synthetic(
  registered, n = 160,
  seed = synthetic.registry.seed("S01.V1", replicate = 1),
  rng.policy = "legacy"
)
synthetic.dataset.checksum(registered.data)
```

For new work, use `rng.policy = "named.stream.v1"`. The `"legacy"` policy
reproduces the random numbers of earlier versions exactly and exists only so
that earlier results can be reproduced. Noise in the predictors, noise in
the response, and the normalization of the truth over the sampled points are
specified separately, and changing one does not change the others.