---
title: "Finding your way around geosmooth"
output:
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{Finding your way around geosmooth}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

<style>
td, th { vertical-align: top; }
td code { white-space: normal; overflow-wrap: anywhere; }
</style>

```{r standalone-width, echo=FALSE, results='asis'}
# Let pkgdown own the site layout; widen only the standalone installed vignette.
if (!identical(Sys.getenv("IN_PKGDOWN"), "true")) {
  cat("<style>body { max-width: 1120px; }</style>")
}
```

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
read.recipe <- getExportedValue("knitr", "read_chunk")
read.recipe(system.file("doc-tools/first-fit.R", package = "geosmooth"))
read.recipe(system.file("doc-tools/graph-workflow.R", package = "geosmooth"))
```

`geosmooth` estimates a smooth response, or a probability mass, from data
given as a matrix of coordinates or as a weighted graph. This guide is
organized by task. The table below names the function to start with for each
task, the next section walks through a first fit, and the catalog that
follows lists every exported function with a link to its help page.

Many of the exported functions build operators or generate synthetic data.
You do not need any of them to fit your first model.

## Choose a starting point

| Your task | Start with | What to know before choosing |
|---|---|---|
| Regress a response on coordinates | [`fit.lps()`](https://pgajer.github.io/geosmooth/reference/fit.lps.html) | Local polynomial smoothing; selects neighborhood size, degree, and kernel by cross-validation. |
| Average predictions from overlapping local models | [`fit.malps()`](https://pgajer.github.io/geosmooth/reference/fit.malps.html) | Model-averaged local polynomial smoothing; supports coordinate and graph-geodesic neighborhoods. |
| Couple local polynomial predictions across overlaps | [`fit.ps.lps()`](https://pgajer.github.io/geosmooth/reference/fit.ps.lps.html) | Prediction-synchronized local polynomial smoothing; experimental R reference implementation. |
| Fit a locally adaptive trend with a sparse difference penalty | [`fit.lpl.tf()`](https://pgajer.github.io/geosmooth/reference/fit.lpl.tf.html) | Local polynomial lifting trend filtering; requires optional package `genlasso`. |
| Fit a one-dimensional continuous response | `fit.gaussian.trend.filtering()` | Order-2 Gaussian trend filtering with squared-error loss, explicit penalties and convergence diagnostics. |
| Fit a one-dimensional binary or binomial response | `fit.logistic.trend.filtering()` | Order-2 logistic trend filtering with a native C++ solver, explicit penalties, optional covariates/offsets, and convergence diagnostics. |
| Fit adaptive curves along fixed or learned projections | `fit.projection.trend.filtering()` | Gaussian or binomial additive projection models; mean-loss penalties and separate joint-fit and search diagnostics. |
| Add synchronization to a lifting trend filter | [`fit.slpl.tf()`](https://pgajer.github.io/geosmooth/reference/fit.slpl.tf.html) | Synchronized local polynomial lifting trend filtering; adds a second, quadratic penalty. |
| Smooth a response already attached to graph vertices | [`fit.metric.graph.lowpass()`](https://pgajer.github.io/geosmooth/reference/fit.metric.graph.lowpass.html) or [`fit.graph.trend.filtering()`](https://pgajer.github.io/geosmooth/reference/fit.graph.trend.filtering.html) | Choose spectral low-pass smoothing or a graph-difference penalty. Edge-weight meanings differ; see below. |
| Penalize local second derivatives on a point cloud | [`fit.ssrhe.hessian.regression()`](https://pgajer.github.io/geosmooth/reference/fit.ssrhe.hessian.regression.html) | Hessian-energy regression, following semi-supervised regression using Hessian energy (SSRHE); use `lambda.selection = "cv"` or `"gcv"` to select penalties. |
| Use an absolute-value Hessian penalty | [`fit.ssrhe.hessian.l1.regression()`](https://pgajer.github.io/geosmooth/reference/fit.ssrhe.hessian.l1.regression.html) | Different regularization from quadratic Hessian energy; has its own solver and tuning controls. |
| Interpolate graph values while fixing a region boundary | [`perform.harmonic.smoothing()`](https://pgajer.github.io/geosmooth/reference/perform.harmonic.smoothing.html) | Iterative harmonic relaxation; [`harmonic.smoother()`](https://pgajer.github.io/geosmooth/reference/harmonic.smoother.html) also records extrema and topology diagnostics. |
| Estimate probability mass over a fixed support | [`fit.density()`](https://pgajer.github.io/geosmooth/reference/fit.density.html) | Takes nonnegative counts or masses; returns mass over the rows of `X`. |
| Estimate where one subject visits that support | [`fit.subject.od()`](https://pgajer.github.io/geosmooth/reference/fit.subject.od.html) | Takes visit indices, including repeats; handles occupation-density fitting and optional visit-level validation. |
| Build reproducible synthetic examples | [`synthetic.spec()`](https://pgajer.github.io/geosmooth/reference/synthetic.spec.html) then [`materialize.synthetic()`](https://pgajer.github.io/geosmooth/reference/materialize.synthetic.html) | Combine geometry, sampling, truth, and response noise. |

The rows of this table are different estimators, listed in no particular
order of predictive performance. Two experimental methods,
parallel-transport trend filtering and transported graph Hessians, are
described in their own [section](#experimental-transport-methods).

## A first regression fit

`fit.lps()` takes a matrix `X` with one observation per row and a response
`y` with one value per row. By default it works in the coordinates you
supply. It can also fit each local polynomial in a *local chart*, a
coordinate system computed near the point of interest, usually by principal
component analysis (PCA). For a first fit, keep the supplied coordinates and
the default numerical settings.

```{r first-fit-model}
```

```{r first-fit-plot, echo=FALSE, fig.width=7.2, fig.height=4, out.width="100%", fig.cap="The points are noisy observations. The solid blue line is the selected fit; the dashed orange line is known only because these data are synthetic.", fig.alt="Noisy observations along a sine curve, with a fitted solid curve close to the dashed synthetic truth."}
```

Here CV selected `r fit$selected$support.size`-point neighborhoods, degree
`r fit$selected$degree`, and the Gaussian kernel, with a tuning RMSE of
`r format(fit$selected$cv.rmse.observed, digits = 3)` response units. A larger neighborhood
averages over more observations, and the degree determines the shape that
is fitted within a neighborhood. Because this RMSE was used to choose the
settings, it is not an independent estimate of prediction error.

```{r first-predictions}
head(fitted(fit))
new.X <- matrix(seq(0.1, 0.9, length.out = 5), ncol = 1)
predict(fit, newdata = new.X)
```

Some terms used throughout the package: cross-validation (CV) chooses
tuning parameters by predicting responses that were held out of the fit. A
*support* is the neighborhood used for one local fit, and a *grid* is a list
of candidate values. The support size and the polynomial degree together
determine whether a local fit has enough observations. For grouped
observations, use
[`lps.grouped.foldid()`](https://pgajer.github.io/geosmooth/reference/lps.grouped.foldid.html) to keep each group in one fold. Use [`lps.nested.cv()`](https://pgajer.github.io/geosmooth/reference/lps.nested.cv.html) when
you need an outer assessment of the whole tuning procedure.

### Prediction and reusable fits

| Fitted object | How to get predictions | Scope |
|---|---|---|
| Local polynomial smoother | `predict(fit)` or `predict(fit, newdata = Xnew)` | Existing evaluation points or new coordinate rows. |
| Model-averaged local polynomial smoother | `predict(fit)` or `predict(fit, newdata = Xnew)` | New coordinates are supported for coordinate-distance fits; graph-geodesic fits currently reject new targets. |
| Either lifting trend filter | `predict(fit)` | Training points only; `newdata` must be `NULL`. |
| Prediction-synchronized, chart-kernel, local-likelihood, graph low-pass, graph trend-filtering, Hessian, or parallel-transport fit | `fitted(fit)` | No registered `predict()` method for these classes in this version. Chart-kernel and local-likelihood fits can receive `X.eval` when fitting. |
| Density or occupation-density result | `fit$rho` | Probability mass on the evaluation support. |
| Tracked harmonic result (`harmonic.smoother()`) | `fitted(fit)` | Values on the original graph vertices. |
| Basic harmonic result (`perform.harmonic.smoothing()`) | `fit$harmonic_predictions` | An unclassed list; retains the basic function's original contract. |

`summary(fit)` shows the selected settings and the diagnostics stored in the
fit, and says so when no convergence information is available.
`fitted(fit)` returns the fitted values as a vector or a matrix, as stored
by the fit. `residuals(fit)` returns the response minus the fitted values,
with `NA` where the response is missing. Residuals exist only at the
training points, so a fit that was evaluated at other coordinates has none.
They are on the scale of the response, also for binomial fits; they are not
Pearson or deviance residuals. Probability mass is a different quantity and
is read from `fit$rho`. See `help("model-values", package = "geosmooth")` and
`help("model-summary", package = "geosmooth")` for details.

Three further terms: a *refit* fits a new response while reusing the
neighborhoods or the operator of an existing fit. It does not repeat the
model selection of the original fit; the help page of each refit method says
whether the penalties stay fixed or can be selected again. An *operator* is
the matrix of differences or derivatives that a method penalizes; it is not
a fitted response. A *spectral basis* is a stored set of graph eigenvectors,
which allows many responses to be smoothed without recomputing them.

### Refit and extract smoother matrices

`refit(fit, y = new.y, ...)` works for the six kinds of model in the table
below. R chooses the method from the class of the fitted object. The
arguments named in the last column belong to that method only.

| Fitted model | Accepted response | What is reused and what can be changed |
|---|---|---|
| MALPS | Vector; omitted or `NULL` reuses the stored response | Supports and averaging weights stay fixed. Optional `weights` affect local coefficient estimation; robust residual weights are recomputed. |
| Lifting trend filter | Vector required | Reuses the residual operator and stored penalty, or accepts a fixed `lambda`. |
| Synchronized lifting trend filter | Vector required | Reuses both operators and stored penalties, or accepts fixed `lambda1` and `lambda2`. |
| Metric graph low-pass | Vector or matrix required | Reuses the spectral basis and smoothing parameter; `per.column.gcv = TRUE` selects a parameter for each response column. |
| Quadratic Hessian-energy regression | Vector or matrix; omitted or `NULL` reuses the stored response | Reuses the operator and fixed penalties; accepts `lambda1`, `lambda2`, `weights`, and `ridge`. CV/GCV fits dispatch here and retain their selected penalties by default. |
| Hessian L1 regression | One response; omitted or `NULL` reuses the stored response | Reuses the derivative operator; accepts `lambda.grid`, fixed or CV `lambda.selection`, and solver controls. |

A refit of a metric graph low-pass fit or of a quadratic Hessian fit
returns a summary that cannot itself be refitted, so for further responses
keep calling `refit()` on the original fit. The other four kinds return
fitted objects that can be refitted again.

`smoother.matrix(fit, ...)` returns the matrix that maps responses to
predictions, for LPS and MALPS fits. For LPS, the fit must have a single
candidate configuration; the matrix maps the training responses to the
predictions at the evaluation points, and `check.tol` sets the tolerance of
an internal check. For MALPS, the matrix maps training responses to
predictions at the training points, given the stored supports and weights;
`max.n` limits the size of the matrix. For a robust MALPS fit,
`allow.robust = TRUE` returns the matrix for the final robust weights held
fixed. It does not reflect how those weights depend on the response.

```{r shared-generics}
malps.fit <- fit.malps(
  X, y, degree = 1L, support.type = "knn", support.size = 12L
)
new.y <- cos(2 * pi * X[, 1])
malps.refit <- refit(malps.fit, y = new.y)
S <- smoother.matrix(malps.fit)
max(abs(as.vector(S %*% new.y) - predict(malps.refit)))
```

If you are updating code written for version 0.1.0: the six `refit.*()`
functions and the two `*.smoother.matrix()` functions are no longer
exported. Call `refit()` and `smoother.matrix()` instead, and in refit calls
replace `fitted.model =` with `object =` and `y.new =` with `y =`. The other
arguments and the classes of the results are unchanged.
`help("refit", package = "geosmooth")` and
`help("smoother.matrix", package = "geosmooth")` link to the arguments of
each method, and the old help names still work.

## Nearest-neighbor probabilities for a binary outcome

`fit.neighbor.probability()` estimates the probability of a binary outcome
at a point from the outcomes of its nearest neighbors. It offers three
estimators, selected by `method`. `"knn"` is the ordinary proportion among
the `k` nearest people. `"case_raw"` and `"case_corrected"` instead extend
the neighborhood until it contains `k` cases; the second corrects the
proportion for this stopping rule. So `k` counts people for `"knn"` and
cases for the other two. Give `k` directly, or give fold labels in `foldid`
to select `k` by the Brier score.

The function is designed for data in which a participant can have several
records:

* `group` identifies the participant of each record. A participant is
  represented by the closest of their records that has an outcome, and all
  records of a participant are left out together when predicting for that
  participant.
* `id` identifies the records and is used to break ties in distance, so that
  results do not depend on row order.
* A record with a missing outcome is never counted as a control.
* `predict(fit)` leaves each participant out of their own prediction. For
  new points that belong to a participant in the training data, pass
  `query.group`.
* If too few neighbors can be reached, the prediction is `NA`; no overall
  rate is substituted.

Three settings control how neighbors are found. With Euclidean coordinates,
`search = "auto"` uses an exact tree search (ANN with `eps = 0`), and
`search = "brute"` uses an exhaustive search that serves as a reference.
`search = "ann"` with, for example, `eps = 0.1` allows approximate
neighbors. The tolerance `eps` bounds the error in distance, not the error
in the estimated probability, and an approximate search can miss
participants or tied records, so the diagnostics for ties are complete only
with an exact search. Chebyshev, intrinsic simplex, and precomputed
distances always use the exhaustive search; asking for ANN or a positive
`eps` with them is an error.

For compositional data normalized so that the largest coordinate is one,
`metric = "linf.simplex"` with `normalization = "max"` measures distance
along the surface of the simplex, using the solver in `dgraphs`. Paths may
cross intermediate faces. Distances are computed for one query point at a
time, so neither a matrix of all pairwise distances nor a graph of the
samples is needed. With `normalization = "none"`, the rows must already be
max-normalized to within `sqrt(.Machine$double.eps)`, and
`normalization = "sum"` is not allowed for this metric. Note that max
normalization combined with `metric = "euclidean"` measures straight-line
distance instead, which can give different neighbors. The search structures
are reused within one prediction call or cross-validation split and rebuilt
on the next call, so a fitted object can be saved and reloaded.

```{r neighbor-probability, eval=FALSE}
fit <- fit.neighbor.probability(X, y, method = "case_corrected", k = 3,
                                id = record_id, group = participant_id,
                                normalization = "max", metric = "linf.simplex")
predict(fit, newdata = query_X, query.group = query_participant_id,
        type = "diagnostics")
```

When the neighborhood stops at a number of cases, it depends on the
outcomes, so the estimator is not a linear smoother and `smoother.matrix()`
is not available for it. To assess a model whose `k` was tuned, use a
separate set of outer folds, grouped by participant.

## Function catalog

The tables in this section list every function exported by `geosmooth`,
each with a one-line description. To open the help page of a function in R,
use, for example, `help("fit.lps", package = "geosmooth")`. The last two
tables list related functions that belong to the `dgraphs` package. They are
not exported by `geosmooth`, and their help pages are in `dgraphs`.

### Shared model operations

| Function | Purpose |
|---|---|
| [`refit()`](https://pgajer.github.io/geosmooth/reference/refit.html) | Refit a supported model with a new response using stored geometry and method-specific tuning controls. |
| [`smoother.matrix()`](https://pgajer.github.io/geosmooth/reference/smoother.matrix.html) | Extract the conditional linear map for a supported LPS or MALPS fit. |

### Local regression and prediction synchronization

| Function | Purpose |
|---|---|
| `fit.neighbor.probability()` | Fit C++ ordinary kNN, raw nearest-case or corrected nearest-case binary estimators, with deterministic ties and grouped tuning. |
| [`fit.lps()`](https://pgajer.github.io/geosmooth/reference/fit.lps.html) | Fit and tune local polynomial regression, including supported Gaussian and binary response modes. |
| [`fit.malps()`](https://pgajer.github.io/geosmooth/reference/fit.malps.html) | Fit and optionally tune a model-averaged collection of local polynomial models. |
| [`fit.ps.lps()`](https://pgajer.github.io/geosmooth/reference/fit.ps.lps.html) | Fit experimental local polynomial models whose predictions are synchronized on overlaps. |
| [`fit.chart.kernel()`](https://pgajer.github.io/geosmooth/reference/fit.chart.kernel.html) | Fit a chart-kernel field, with optional reference-measure quadrature weights. |
| [`fit.local.likelihood()`](https://pgajer.github.io/geosmooth/reference/fit.local.likelihood.html) | Fit local density or Bernoulli likelihood models; returns a fitted field. |
| `unfold.simplex.chart()` | Certify and unfold a local one/two-face simplex neighborhood before constructing charts. |
| `fit.masked.graph.field()` | Fit a graph field from partially observed responses. |
| `support.masked.graph.field()` | Identify graph components supported by observed responses. |
| `cv.masked.graph.field()` | Assess masked fields using participant-grouped nested cross-validation. |
| `compare.masked.graph.fields()` | Compare candidate graph fields with independent workers. |
| `select.masked.graph.field()` | Assemble nested graph-selection predictions from inner-fold choices. |

### Model assessment and uncertainty

| Function | Purpose |
|---|---|
| `assess.logistic.trend.filtering()` | Test constant binary probability by retuned outcome permutations and measure nested-CV improvement over a training-only constant. Summarize absolute probability deviation and total variation; optionally test them using the same fits, with Holm adjustment. Requires independent unit-weight binary observations without covariates or offsets. |
| `assess.gaussian.trend.filtering()` | Assess continuous outcomes with retuned squared-error permutation tests, optional absolute-deviation/total-variation tests and nested-CV MSE/MAE improvement. Requires independent unit-weight observations without covariates or offsets, exchangeable under a common-distribution null. |
| [`lps.grouped.foldid()`](https://pgajer.github.io/geosmooth/reference/lps.grouped.foldid.html) | Assign whole groups of observations to validation folds. |
| [`lps.nested.cv()`](https://pgajer.github.io/geosmooth/reference/lps.nested.cv.html) | Assess a tuned local polynomial smoother using outer and inner validation folds. |
| [`lps.backend.diagnostics()`](https://pgajer.github.io/geosmooth/reference/lps.backend.diagnostics.html) | Inspect the computational backend, chart dimensions, and selected controls. |
| [`lps.pointwise.band()`](https://pgajer.github.io/geosmooth/reference/lps.pointwise.band.html) | Calculate pointwise uncertainty for a supported fixed-configuration Gaussian fit. |
| [`malps.gcv()`](https://pgajer.github.io/geosmooth/reference/malps.gcv.html) | Calculate generalized cross-validation and optional leave-one-out diagnostics for a stored model-averaged fit. |
| [`bootstrap.malps()`](https://pgajer.github.io/geosmooth/reference/bootstrap.malps.html) | Measure refit variability by bootstrap case weights while holding the selected geometry fixed. |

Generalized cross-validation (GCV) estimates prediction error from the
effective degrees of freedom of a smoother, without fitting a separate model
for each left-out observation. It is a different criterion from
cross-validation with folds, and the two need not select the same settings.

The uncertainty functions for LPS and for MALPS rest on different
assumptions. `lps.pointwise.band()` requires a Gaussian fit made with the R
backend and the orthogonal polynomial basis, with a single candidate
configuration and a fixed chart dimension, and evaluated at the training
points. `malps.gcv()` and `bootstrap.malps()` treat the selected supports
and weights as fixed and do not repeat the model selection. Neither gives
simultaneous coverage, and neither accounts for the uncertainty of choosing
the settings from the data.

### Lifting trend filters

| Function | Purpose |
|---|---|
| [`lpl.tf.operator()`](https://pgajer.github.io/geosmooth/reference/lpl.tf.operator.html) | Build a local polynomial residual operator. |
| [`fit.lpl.tf()`](https://pgajer.github.io/geosmooth/reference/fit.lpl.tf.html) | Fit a trend filter with an absolute-value penalty on those residuals. |
| [`slpl.tf.operator()`](https://pgajer.github.io/geosmooth/reference/slpl.tf.operator.html) | Build the residual operator together with an overlap synchronization operator. |
| [`fit.slpl.tf()`](https://pgajer.github.io/geosmooth/reference/fit.slpl.tf.html) | Fit the residual penalty plus quadratic synchronization. |

Both fitting functions accept an operator that was built beforehand, or
build one themselves. For a single fit, call the fitting function directly.
Build the operator yourself when you want to fit several responses, examine
the neighborhoods, or avoid recomputing an expensive operator. Both fitting
functions require the optional package `genlasso`, and both predict only at
the training points.

### Graph spectral smoothing and graph trend filtering

| Function | Purpose |
|---|---|
| [`fit.metric.graph.lowpass()`](https://pgajer.github.io/geosmooth/reference/fit.metric.graph.lowpass.html) | Select a graph-spectral smoothing parameter by GCV and fit a response. |
| [`metric.graph.lowpass.operator()`](https://pgajer.github.io/geosmooth/reference/metric.graph.lowpass.operator.html) | Construct a graph Laplacian from edge lengths or explicit conductances. |
| [`metric.graph.lowpass.basis()`](https://pgajer.github.io/geosmooth/reference/metric.graph.lowpass.basis.html) | Construct and retain a partial or complete graph spectrum for repeated smoothing. |
| [`apply.metric.graph.lowpass.path()`](https://pgajer.github.io/geosmooth/reference/apply.metric.graph.lowpass.path.html) | Apply a supplied basis over a parameter grid; retain the candidate fits. |
| [`metric.graph.heat.eta.grid()`](https://pgajer.github.io/geosmooth/reference/metric.graph.heat.eta.grid.html) | Choose heat-time candidates using spectral information. |
| [`metric.graph.heat.extend.lower()`](https://pgajer.github.io/geosmooth/reference/metric.graph.heat.extend.lower.html) | Propose a smaller heat time for an external search, subject to resolution safeguards. |
| [`graph.trend.filtering.operator()`](https://pgajer.github.io/geosmooth/reference/graph.trend.filtering.operator.html) | Construct graph differences or experimental path divided differences. |
| [`fit.graph.trend.filtering()`](https://pgajer.github.io/geosmooth/reference/fit.graph.trend.filtering.html) | Fit an absolute-value graph-difference penalty; requires `genlasso`. |

`geosmooth` does not construct graphs or compute shortest paths; the
`dgraphs` package does. A graph is passed as an adjacency list, in which
element `i` holds the indices (starting from 1) of the neighbors of vertex
`i`, and a weight list of the same shape that holds a value for each edge.

**The meaning of the edge values depends on the method.** A *length* is a
distance between two vertices. A *conductance* is the strength of their
connection, and larger values mean a stronger connection. Metric low-pass
smoothing expects lengths by default and converts them to conductances; pass
`weight.type = "conductance"` if your values are already conductances. Graph
trend filtering with the default recursive operator treats the values as
conductances, whereas its path divided-difference operator treats them as
lengths.

The heat time `eta` controls the amount of smoothing in graph low-pass
filtering. A basis that contains only part of the spectrum may be unable to
represent the filter for small positive heat times. If you build or extend a
grid of heat times yourself, consult the spectral diagnostics returned with
the basis.

With a complete spectrum, heat smoothing also accepts missing responses, as
long as every connected component of the graph has at least one observed
response. The [observed-data heat guide](verified-heat-fallback.html)
describes the criterion that is minimized in this case and the optional
numerical cross-checks requested with
`numerical.policy = "verified.fallback"`. Under that policy, a heat time is
not selected automatically when the GCV scores cannot distinguish the
candidates, either because the scores are constant by construction or
because their differences are smaller than the numerical accuracy. The
resolution bands reported there describe numerical accuracy; they are not
confidence intervals.

#### From a constructed graph to a fitted response

This example constructs a graph with `dgraphs` and smooths a response on it.
The 20 locations lie on a line, and the radius used to construct the graph
joins only adjacent locations, so the graph is a chain. Vertex *i*, row *i*
of the coordinates, and entry *i* of the response must refer to the same
observation. The smoother cannot detect a response that is in the wrong
order, so keep your identifiers in a table and match the response to the
vertices by identifier before fitting.

```{r graph-workflow}
```

With `conductance.alpha = 1`, an edge of length 2 receives the conductance
`1 / (2 + 1e-8)`, about 0.5, so longer edges connect their endpoints less
strongly. The conversion from lengths to conductances affects the result of
the smoothing. The example computes all eigenvectors of the graph Laplacian,
so nothing is lost by truncating the spectrum. For a large graph you would
keep only some of them; `help("metric.graph.lowpass.basis")` describes the
diagnostics to consult in that case. The second response is smoothed with
the same eigenvectors and the heat time already selected. The GCV score that
is printed was used for tuning and is not an estimate of held-out prediction
error.

| Where the values are used | What they mean | Consequence |
|---|---|---|
| Coordinate neighborhoods | Distances in the supplied coordinate units | Rescaling coordinates changes supports. |
| Graph-geodesic neighborhoods | Sums of edge lengths along paths | Unreachable vertices have infinite distance; support rules must handle connectivity. |
| Metric low-pass, `weight.type = "length"` (default) | Positive lengths transformed to conductances | Reciprocal edges and matching lengths are required; zero lengths are rejected. |
| Metric low-pass, `weight.type = "conductance"` | Supplied nonnegative connection strengths | Reciprocal weights must agree; zero-weight edges are removed and positive weights are used unchanged. |
| Recursive graph differences | Conductance-like weights | Greater weight penalizes differences more strongly. |
| Path divided differences | Edge lengths | The operator uses metric spacing, not conductance. |

For example, `metric.graph.lowpass.operator(list(2L, integer()), list(1, numeric()))`
is rejected, because the edge from vertex 1 to vertex 2 is not also listed
from vertex 2 to vertex 1. If your input has directed edges or edges of
length zero, decide what they should mean before you construct an undirected
graph. A conductance of zero means that there is no connection. A length of
zero is not allowed, because it cannot be converted to a conductance.

### One-dimensional Gaussian and logistic trend filtering

| Function | Purpose |
|---|---|
| `fit.logistic.trend.filtering()` | Fit a univariate piecewise-quadratic logit curve using the native C++ backend. |
| `fit.gaussian.trend.filtering()` | Fit a univariate piecewise-quadratic mean curve with squared-error loss using the same order-2 geometry and native C++ solvers. |
| `fit.projection.trend.filtering()` | Fit centered order-2 curves on fixed directions or learn directions by projection pursuit; supports Gaussian and binomial losses. |
| `cv.projection.trend.filtering()` | Select component counts and mean-loss penalties using supplied, optionally grouped folds; repeat all direction learning within each training split. |

For a continuous outcome, fit the model over a list of penalties that you
supply, and then assess both the association with the predictor and the
held-out prediction:

```{r gaussian-trend-filtering-assessment}
x.continuous <- seq(0, 1, length.out = 40)
y.continuous <- sin(2*pi*x.continuous) + 0.2*cos(seq_along(x.continuous))
gaussian.fit <- fit.gaussian.trend.filtering(
  x.continuous, y.continuous, lambda = c(Inf, 0.01), backend = "banded")
# Nine permutations keep this API example short; use more for inference.
gaussian.assessment <- assess.gaussian.trend.filtering(
  gaussian.fit, permutations = 9, inner.folds = 3, outer.folds = 3, seed = 42,
  additional.tests = c("absolute.deviation", "total.variation"))
gaussian.assessment$tests
gaussian.assessment$prediction$scores
```

The test statistic is the reduction in the residual sum of squares relative
to the sample mean. The tuning is repeated for every permutation. The test
assumes that the outcomes are exchangeable under the null hypothesis; equal
means with variances that depend on the predictor are not enough. The
nested cross-validation is a separate assessment. It compares the smoother
with the mean of the outer training data, using MSE and MAE on the same
held-out rows. The curve summaries are in the units of the response. Note
that `lambda = Inf` gives a quadratic curve, whereas the baseline of the
assessment is a constant. The example shows how the functions are used; with
nine permutations it says nothing about the calibration or the power of the
test.

The [logistic trend-filtering vignette](logistic-trend-filtering.html) walks
through synthetic probabilities, trial weights and held-out Pima predictions.

`fit.logistic.trend.filtering()` accepts binary observations, or
proportions together with the numbers of trials as weights. You must supply
the penalties in `lambda`. They refer to the predictor after it has been
rescaled to the range of the training data, not to its original units.
Observations with exactly equal predictor values share one value of the
curve but keep their own covariates and offsets. To select a penalty by
cross-validation, supply `foldid`; give repeated observations of one subject
the same fold label.

Always inspect `fit$diagnostics`. A fit that reached the iteration limit,
detected separation, or failed numerically is not a successful fit.
`predict()` applies the transformations learned from the training data and
holds the curve constant outside the training range. No ridge penalty is
added unless you request one.

Two solvers are available. `backend = "dense"` is the reference solver for
small problems. The banded solver does not form dense matrices over all
knots, yet its stopping checks still examine every knot coefficient.
`backend = "auto"` chooses the banded solver when there are more than 256
distinct predictor values or when the dense solver would need too much
memory. Both solvers check an estimate of their memory use before
allocating.

Cross-validation for Gaussian, logistic, local polynomial lifting, and graph
trend filtering uses `cv.penalty.scale = "training_mass"` by default. A fold
is fitted to less data than the full fit, so its penalty is multiplied by
the ratio of the training weight to the full weight. This keeps the penalty
in the same proportion to the mean loss. Set `cv.penalty.scale = "none"` to
use the same numerical penalty in every fold, as earlier versions did. The
one-dimensional fits scale explicit ridge penalties in the same way and
record the settings used in each fold. The assessment functions apply the
same rule to their outer and inner training sets. Fits saved by earlier
versions keep the earlier convention, and fits with a fixed penalty are not
affected.

### Hessian-energy regression

| Function | Purpose |
|---|---|
| [`ssrhe.hessian.operator()`](https://pgajer.github.io/geosmooth/reference/ssrhe.hessian.operator.html) | Estimate local derivative operators and Hessian-energy matrices from point-cloud neighborhoods. |
| [`fit.ssrhe.hessian.regression()`](https://pgajer.github.io/geosmooth/reference/fit.ssrhe.hessian.regression.html) | Fit quadratic Hessian-energy regression using fixed penalties, label CV, or GCV; optionally select supports. |
| [`fit.ssrhe.hessian.l1.regression()`](https://pgajer.github.io/geosmooth/reference/fit.ssrhe.hessian.l1.regression.html) | Fit an absolute-value derivative penalty with fixed or CV-selected regularization. |
| [`ssrhe.support.grid()`](https://pgajer.github.io/geosmooth/reference/ssrhe.support.grid.html) | Build candidate adaptive-radius support profiles for tuning. |

The Hessian of a function is its matrix of second derivatives. These
methods estimate it locally from the points near each observation and
penalize its size. A quadratic penalty reduces curvature everywhere by a
moderate amount, whereas an absolute-value penalty can set some second
derivatives exactly to zero. The two penalties therefore give different
estimates.

Cross-validation works when some points have no observed response
(semi-supervised regression). GCV requires a single response vector with no
missing values and strictly positive observation weights. The trace needed
for GCV can be computed exactly or estimated with random probes; the exact
computation costs more, and the estimate depends on a random seed. For the
absolute-value penalty, `solver = "admm"` uses the built-in alternating
direction method of multipliers, and the `genlasso` solver requires that
optional package.

To select the penalties of the quadratic method, call the same fitting
function with `lambda.selection = "cv"` or `"gcv"`, and put the settings of
the selection in `cv.control` or `gcv.control`. A single penalty value and a
grid of candidates cannot be given together.

```{r hessian-selection}
hessian.fit <- fit.ssrhe.hessian.regression(
  X, y, k = 12L, tangent.dim = 1L,
  lambda.selection = "cv", lambda1.grid = c(0.01, 0.1),
  cv.control = list(foldid = rep(1:5, length.out = nrow(X)))
)
hessian.fit$lambda
```

The separate `.cv()` and `.gcv()` fitting functions of version 0.1.0 have
been removed. Their old help names now open the help page of
`fit.ssrhe.hessian.regression()`, which explains how to update a call. Calls
with fixed `lambda1` and optional `lambda2` work as before.

### Harmonic interpolation

| Function | Purpose |
|---|---|
| [`perform.harmonic.smoothing()`](https://pgajer.github.io/geosmooth/reference/perform.harmonic.smoothing.html) | Relax interior values while holding the region boundary fixed; return convergence diagnostics. |
| [`harmonic.smoother()`](https://pgajer.github.io/geosmooth/reference/harmonic.smoother.html) | Run relaxation while recording extrema, prediction history, and topology stability. |
| [`get.region.boundary()`](https://pgajer.github.io/geosmooth/reference/get.region.boundary.html) | Find region vertices with a neighbor outside the region. |

Harmonic smoothing replaces the values inside a region of the graph by
values that are, at each vertex, the weighted average of the neighbors'
values, while the values on the boundary of the region stay fixed. A
boundary vertex is a vertex of the region that has a neighbor outside it. A
vertex of degree one is therefore not a boundary vertex when the region is
the whole graph.

The two functions behave differently when the region has no boundary, so
they are not interchangeable. `perform.harmonic.smoothing()` returns its
input unchanged with `status = "no_boundary"`, whereas `harmonic.smoother()`
still runs the iteration. Both report numerical convergence only when the
change between iterations and the harmonic residual are both small.

`harmonic.smoother()` returns a `harmonic_smoother` object with `print()`,
`summary()`, and `plot()` methods. It also reports whether the local extrema
of the values stopped changing: `stability_detected` says whether a stable
stretch of iterations was found, and `stable_iteration` identifies it, or is
`NA` if none was found. This is a separate diagnostic from convergence. For
example, on a path whose values alternate, recording every second iteration
can look stable while the values are still oscillating. Check `converged`
and `max_residual`, and do not take stable extrema as a sign that the amount
of smoothing is right.

### Probability mass and subject occupation

| Function | Purpose |
|---|---|
| [`fit.density()`](https://pgajer.github.io/geosmooth/reference/fit.density.html) | Fit empirical or graph-random-walk mass estimates from nonnegative weights on a support. |
| [`normalize.density()`](https://pgajer.github.io/geosmooth/reference/normalize.density.html) | Convert a numeric field or supported fitted object to probability mass, with normalization diagnostics. |
| [`fit.subject.od()`](https://pgajer.github.io/geosmooth/reference/fit.subject.od.html) | Convert repeated visit indices into counts and fit a subject-occupation distribution. |

```{r occupation}
support <- matrix(seq(0, 1, length.out = 6), ncol = 1)
occupation <- fit.subject.od(
  support, subject.index = c(1L, 3L, 3L, 6L), method = "empirical"
)
occupation$rho
sum(occupation$rho)
```

In this example the third support point is visited twice and the first and
sixth once each. `rho` is a vector of probability **masses**, one for each
support point, that sums to one. It is not a density per unit volume.
[`normalize.density()`](https://pgajer.github.io/geosmooth/reference/normalize.density.html) converts a fitted field to such a vector; it does not
fit anything. If the field has negative values, they are set to zero before
normalizing, and the result reports how much was changed.

[`fit.subject.od()`](https://pgajer.github.io/geosmooth/reference/fit.subject.od.html) offers more methods than [`fit.density()`](https://pgajer.github.io/geosmooth/reference/fit.density.html), because it can
also fit one of the regression smoothers and normalize the result. Two kinds
of cross-validation occur here and should not be confused. Visit-level
cross-validation holds out visits and takes its folds from `visit.foldid`.
The cross-validation of a regression smoother holds out rows of the support
and takes its folds from `foldid`, or from the argument of that name in the
chosen method.

### Experimental transport methods

| Function | Purpose |
|---|---|
| [`pttf.geometry()`](https://pgajer.github.io/geosmooth/reference/pttf.geometry.html) | Build local tangent frames and their transport across a graph. |
| [`pttf.operator()`](https://pgajer.github.io/geosmooth/reference/pttf.operator.html) | Construct parallel-transport trend-filtering differences from that geometry. |
| [`pttf.operator.filter.rows()`](https://pgajer.github.io/geosmooth/reference/pttf.operator.filter.rows.html) | Select operator rows while maintaining provenance and related operator data. |
| [`fit.pttf.trend.filtering()`](https://pgajer.github.io/geosmooth/reference/fit.pttf.trend.filtering.html) | Fit a parallel-transport trend filter with an absolute-value or quadratic penalty. |
| [`transported.graph.hessian.operator()`](https://pgajer.github.io/geosmooth/reference/transported.graph.hessian.operator.html) | Build a separate experimental graph Hessian diagnostic with several direction-matching rules. |

Derivatives estimated at neighboring points are expressed in different
local coordinate frames. Parallel transport aligns the frames so that the
derivatives can be compared. The functions in this table are experimental,
and they return detailed diagnostics about the frames, the operator, and the
solver. `transported.graph.hessian.operator()` serves a different purpose
from [`pttf.operator()`](https://pgajer.github.io/geosmooth/reference/pttf.operator.html): it builds a diagnostic operator and does not lead to
a fitted response.

### Synthetic datasets: specifications, responses, and reproducibility

| Function | Purpose |
|---|---|
| [`synthetic.spec()`](https://pgajer.github.io/geosmooth/reference/synthetic.spec.html) | Combine geometry, sampling, truth, and response descriptors without drawing data. |
| [`materialize.synthetic()`](https://pgajer.github.io/geosmooth/reference/materialize.synthetic.html) | Draw a dataset from a specification with an explicit seed and random-number policy. |
| [`materialize.synthetic.instance()`](https://pgajer.github.io/geosmooth/reference/materialize.synthetic.instance.html) | Reproduce a registered frozen instance and optionally validate its checksum. |
| [`synthetic.registry.ids()`](https://pgajer.github.io/geosmooth/reference/synthetic.registry.ids.html) | List maintained recipe identifiers. |
| [`synthetic.registry.spec()`](https://pgajer.github.io/geosmooth/reference/synthetic.registry.spec.html) | Look up a recipe specification; some recipe families allow parameter overrides. |
| [`synthetic.registry.seed()`](https://pgajer.github.io/geosmooth/reference/synthetic.registry.seed.html) | Reproduce legacy SSRHE recipe seeds; this is not a seed generator for every registry family. |
| [`validate.synthetic.dataset()`](https://pgajer.github.io/geosmooth/reference/validate.synthetic.dataset.html) | Check the structure and consistency of a materialized dataset. |
| [`compare.synthetic.dataset()`](https://pgajer.github.io/geosmooth/reference/compare.synthetic.dataset.html) | Compare two datasets with explicit numerical tolerances. |
| [`synthetic.dataset.checksum()`](https://pgajer.github.io/geosmooth/reference/synthetic.dataset.checksum.html) | Compute a content checksum for reproducibility checks. |
| [`synthetic.sampling.stratified()`](https://pgajer.github.io/geosmooth/reference/synthetic.sampling.stratified.html) | Describe allocation between the supported synthetic strata. |
| [`synthetic.truth.named()`](https://pgajer.github.io/geosmooth/reference/synthetic.truth.named.html) | Select a named truth specification with optional parameters. |
| [`synthetic.truth.polynomial()`](https://pgajer.github.io/geosmooth/reference/synthetic.truth.polynomial.html) | Specify polynomial truth in the chosen evaluation coordinates. |
| [`synthetic.truth.gaussian.mixture()`](https://pgajer.github.io/geosmooth/reference/synthetic.truth.gaussian.mixture.html) | Specify Gaussian-mixture truth using centers, scales, and weights. |
| [`synthetic.truth.occupation.mixture()`](https://pgajer.github.io/geosmooth/reference/synthetic.truth.occupation.mixture.html) | Specify occupation-mixture truth with covariance and probability controls. |
| [`synthetic.truth.logit()`](https://pgajer.github.io/geosmooth/reference/synthetic.truth.logit.html) | Specify a logit-based truth with prevalence and clipping controls. |
| [`synthetic.response.gaussian()`](https://pgajer.github.io/geosmooth/reference/synthetic.response.gaussian.html) | Specify Gaussian response noise with a supplied standard deviation. |
| [`synthetic.response.heteroskedastic.gaussian()`](https://pgajer.github.io/geosmooth/reference/synthetic.response.heteroskedastic.gaussian.html) | Specify Gaussian noise whose scale depends on truth. |
| [`synthetic.response.clustered.gaussian()`](https://pgajer.github.io/geosmooth/reference/synthetic.response.clustered.gaussian.html) | Specify within-cluster correlated Gaussian noise. |
| [`synthetic.response.laplace.outlier()`](https://pgajer.github.io/geosmooth/reference/synthetic.response.laplace.outlier.html) | Specify Laplace noise and an outlier component. |
| [`synthetic.response.bernoulli()`](https://pgajer.github.io/geosmooth/reference/synthetic.response.bernoulli.html) | Specify binary response draws and optional positive-count requirements. |

The functions whose names begin with `synthetic.` only describe a data set;
no random numbers are drawn until you call [`materialize.synthetic()`](https://pgajer.github.io/geosmooth/reference/materialize.synthetic.html). The
result is a `synthetic_dataset` object with `print()`, `plot()`, and
`as.data.frame()` methods. `vignette("synthetic-datasets", package =
"geosmooth")` gives complete examples. For new work, use the random-number
policy `"named.stream.v1"`. The other seed and algorithm settings exist so
that data sets generated by earlier versions can be reproduced exactly;
leave them at their defaults unless that is your aim.

### Geometry constructors and helpers supplied by dgraphs

The following functions were exported by `geosmooth` 0.1.0 and are now
provided by `dgraphs`. In existing code, replace `geosmooth::` with
`dgraphs::`. Their help pages are in `dgraphs`; see, for example,
`help("synthetic.circle", package = "dgraphs")`.

| Function | Purpose |
|---|---|
| `dgraphs::synthetic.circle()` | Describe a circle or circular arc and its embedding. |
| `dgraphs::synthetic.helix()` | Describe a helical curve. |
| `dgraphs::synthetic.trefoil()` | Describe a trefoil curve. |
| `dgraphs::synthetic.torus.patch()` | Describe a patch of a torus. |
| `dgraphs::synthetic.sphere.cap()` | Describe a spherical cap. |
| `dgraphs::synthetic.quadform()` | Describe a flat or quadratic geometry with intrinsic and ambient dimensions. |
| `dgraphs::synthetic.simplex()` | Describe a simplex with a given number of parts. |
| `dgraphs::synthetic.stratified()` | Combine strata and optional junctions. |
| `dgraphs::synthetic.stratum.point()` | Describe a point stratum. |
| `dgraphs::synthetic.stratum.segment()` | Describe a line-segment stratum. |
| `dgraphs::synthetic.stratum.rectangle()` | Describe a rectangular stratum. |
| `dgraphs::synthetic.point.line.junction()` | Describe a point-and-line junction geometry. |
| `dgraphs::embed.synthetic.geometry()` | Map latent coordinates to the ambient space. |
| `dgraphs::edge.lengths.synthetic.geometry()` | Evaluate geometry-aware lengths between latent endpoints. |
| `dgraphs::quadform.gradient()` | Evaluate the gradient of a quadratic geometry. |
| `dgraphs::quadform.metric()` | Evaluate its induced metric. |

### Sampling constructors supplied by dgraphs

| Function | Purpose |
|---|---|
| `dgraphs::synthetic.sampling.uniform.interval()` | Describe uniform draws on an interval. |
| `dgraphs::synthetic.sampling.uniform.box()` | Describe uniform draws in a box. |
| `dgraphs::synthetic.sampling.uniform.rectangle()` | Describe uniform draws in a rectangle with its specified draw algorithm. |
| `dgraphs::synthetic.sampling.uniform.disk()` | Describe uniform draws in a disk. |
| `dgraphs::synthetic.sampling.grid.interval()` | Describe an interval grid with an endpoint convention. |
| `dgraphs::synthetic.sampling.gapped.uniform()` | Describe sampling across separated intervals. |
| `dgraphs::synthetic.sampling.truncated.normal()` | Describe bounded Gaussian draws. |
| `dgraphs::synthetic.sampling.clustered()` | Describe clustered observations around random centers. |
| `dgraphs::synthetic.sampling.dirichlet.zeros()` | Describe Dirichlet sampling with structural zeros. |

Two sampling descriptions of similar regions do not necessarily draw the
same random numbers, because each description also fixes the order of the
draws and the algorithm.

## Finding detailed help

The function links in this guide open the online reference. The same help
pages are installed with the package and need no internet connection: use
`help("fit.lps", package = "geosmooth")`, with the name of the function you
want. `help("geosmooth-package", package = "geosmooth")` lists the main
functions by task, and `browseVignettes("geosmooth")` lists the installed
guides.

The [synthetic datasets guide](synthetic-datasets.html) develops reproducible
geometry, sampling, truth, and response specifications in more detail.


```{r help-navigation, eval=FALSE}
help(package = "geosmooth")
help("geosmooth-package", package = "geosmooth")
help("fit.lps", package = "geosmooth")
help("predict.lps", package = "geosmooth")
vignette("function-guide", package = "geosmooth")
vignette("synthetic-datasets", package = "geosmooth")
```

The help pages give the exact arguments, defaults, return values, and
examples. Use this guide to find the right function, and read its help page
before changing its more specialized settings.

## Tuning and numerical settings

The help page of each main fitting function groups its arguments into those
needed for a first fit and those that control model selection, the
neighborhoods, the numerical solver, and the diagnostics that are returned.
Start by tuning over a small grid of the statistical settings. Tightening
numerical tolerances is not a substitute for that.

| Method | Arguments for folds and selection | What is selected, and the main cost |
|---|---|---|
| LPS / MALPS | `foldid`, `cv.folds`, support/degree/kernel grids | Held-out response selection; candidate combinations multiply fitting work. |
| Quadratic Hessian | `lambda.selection`, `cv.control`, `gcv.control` | Label CV or GCV; exact trace calculation is more expensive than sampled traces. |
| Hessian L1 | `fold.id`, `nfolds`, `lambda.grid`, `admm.*` | CV excludes unconverged solves; each candidate and fold needs a solve. |
| Graph low-pass | `eta.grid`, `n.eigenpairs` | GCV over heat times; a full spectral basis uses quadratic storage. |

The argument names differ among the methods because the selection criteria
differ; giving them the same names would not make the criteria equivalent.

`smoother.matrix()` returns a dense matrix with one row per evaluation point
and one column per training point. Stored candidate fits and solver
histories also take memory, so request them only when you intend to examine
them.

The catalog above covers every function exported by `geosmooth` 0.2.0. The
geometry and sampling functions belong to `dgraphs`. Methods for generics
such as `print()`, `predict()`, and `summary()` are not listed individually;
the package index, `help(package = "geosmooth")`, lists them all.

For intrinsic simplex charts, weighted PCA/PLS and their limits, see the
[local simplex charts guide](local-simplex-charts.html).

For worked synthetic comparisons, held-out prediction figures and reproducible
PCA/PLS examples, see [PCA and PLS charts](pca-pls-chart-comparison.html).