INLASPDE.jl

INLASPDE

SPDE–Matérn machinery for LatentGaussianModels. Provides FEM-assembled sparse precision matrices Q(τ, κ) on triangulated meshes and exposes SPDE2 <: AbstractLatentComponent for use in LatentGaussianModel stacks.

See plans/plan.md for the package roadmap and /plans/decisions.md ADR-001, ADR-007, ADR-013 for load-bearing design decisions.

Module layout (populated incrementally per milestone)

• assembly/ — per-element FEM matrices and Q(α, τ, κ) (M1).
• components/ — SPDE2 <: AbstractLatentComponent (M2).
• priors/ — PC-Matérn prior on (range, σ) (M2).
• mesh/ — inla_mesh_2d and refinement helpers (M3).
• projector.jl — MeshProjector for observation-to-vertex maps (M4).

FEMMatrices{T, SC, SG}

Container for the precomputed FEM matrices (C, G₁, C̃, G₂) on a fixed triangular mesh. These matrices depend only on mesh geometry, not on the SPDE hyperparameters (τ, κ), so they are assembled once and reused for every precision evaluation.

Fields

• C::SC — full P1 mass matrix.
• G1::SG — P1 stiffness matrix (∫ ∇φ_i · ∇φ_j).
• C_lumped::SC — diagonal mass matrix from lumped_mass(C).
• G2::SG — G₁ · C̃⁻¹ · G₁, used for α = 2.

Construct via FEMMatrices(points, triangles).

FEMMatrices(points, triangles)

Assemble C, G₁, C̃, and G₂ from a 2D triangular mesh given as raw arrays. See assemble_fem_matrices for the argument conventions.

FEMMatrices(points::AbstractVector, segments::AbstractMatrix)

1D variant: assemble C, G₁, C̃, and G₂ from a 1D segment mesh given as raw arrays. See assemble_fem_matrices_1d for the argument conventions.

For α = 2 in 1D the smoothness is ν = 3/2 (well-defined under PC-Matérn); the same G₂ = G₁ · C̃⁻¹ · G₁ is used as in 2D.

FEMMatrices(mesh::INLAMesh1D)

Convenience overload: assemble FEM matrices from a 1D mesh.

FEMMatrices(mesh::INLAMesh)

Assemble the FEM matrices (C, G₁, C̃, G₂) directly from an INLA mesh. Forwards to FEMMatrices(points, triangles).

GaussianBasisPrior(; mean, prec)
GaussianBasisPrior(p::Integer; mean = 0.0, prec = 1.0)

Independent Gaussian prior on the basis-coefficient vector θ = [θ_τ; θ_κ] of SPDE2NonStationary. Each coefficient has its own mean and precision:

log π(θ) = ∑_k -½ λ_k (θ_k - μ_k)^2 + ½ log(λ_k / 2π)

The normalising constant is included so log_basis_prior_density returns the proper log density — required for marginal-likelihood comparison against R-INLA's inla.spde2.matern, which uses the same parameterisation under theta.prior.mean / theta.prior.prec.

Fields

• mean::Vector{T} — prior means μ_k, one per basis coefficient.
• prec::Vector{T} — prior precisions λ_k > 0, one per basis coefficient. prec[k] = 0 is rejected; users wanting an improper flat prior on a coefficient should choose a small but positive value (e.g. 1.0e-3), mirroring R-INLA practice.

Length contract

length(prior.mean) == length(prior.prec) is enforced at construction and equals the total number of basis coefficients (p_τ + p_κ for SPDE2NonStationary).

Defaults

mean = 0, prec = 1 for every coefficient — R-INLA's documented inla.spde2.matern defaults. Override per-coefficient for canonical non-stationary fits (a wide intercept prior + tighter spline-basis priors is the standard pattern).

Example

# Two τ-coefficients (intercept + spline) and three κ-coefficients
# (intercept + two region-indicators):
prior = GaussianBasisPrior(
    mean = [0.0, 0.0, 0.0, 0.0, 0.0],
    prec = [1.0e-3, 1.0, 1.0e-3, 1.0, 1.0]   # wide on intercepts
)

See ADR-028 for the design rationale (R-INLA-parity choice over a unit-Gaussian default and over PC-on-basis-norm).

INLAMesh{T, TRI}

A 2D triangular mesh in INLA-SPDE form: raw points and triangles matrices suitable for FEMMatrices / SPDE2, plus the underlying DelaunayTriangulation.Triangulation for advanced use (further refinement, point-in-triangle queries, etc.).

Fields

• points::Matrix{T} — n × 2 vertex coordinates.
• triangles::Matrix{Int} — m × 3 one-based vertex indices.
• boundary::Vector{Int} — ordered boundary vertex indices (open loop in the CCW convention; the loop is closed implicitly).
• triangulation::TRI — the raw DT object.

Construct via inla_mesh_2d.

INLAMesh1D{T}

A 1D segment mesh in INLA-SPDE form: a sorted vector of vertex coordinates plus the implicit segment topology (each row of segments connects consecutive vertices). Suitable for the 1D FEMMatrices overload and SPDE1D.

Fields

• points::Vector{T} — vertex coordinates in ascending order.
• segments::Matrix{Int} — (n - 1) × 2 one-based vertex indices, segments[k, :] == [k, k + 1].
• boundary::NTuple{2, Int} — endpoint vertex indices (1, n).

Construct via inla_mesh_1d.

MeshProjector{T, M, S}

Sparse barycentric projection from mesh-vertex values to observation locations. For a field u defined on mesh vertices, the interpolated value at observation i is (A * u)[i], where A is the n_obs × n_vertices sparse matrix stored in the A field. Each row of A has at most three nonzeros: the barycentric weights of observation i in its enclosing mesh triangle, which sum to 1.

Fields

• mesh::INLAMesh — the source mesh.
• locations::Matrix{T} — n_obs × 2 observation coordinates.
• A::SparseMatrixCSC{T, Int} — the sparse projection matrix.

Construct via MeshProjector(mesh, locations).

Example

mesh = inla_mesh_2d(; boundary = [0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0],
                    max_edge = 0.25)
locs = rand(50, 2)
P = MeshProjector(mesh, locs)
u  = mesh.points[:, 1]              # field defined on mesh vertices
y  = P * u                          # interpolate at locations

MeshProjector(mesh::INLAMesh, locations; outside = :error, atol = 0.0)

Build a sparse barycentric projector from mesh to the points in locations (n_obs × 2 matrix).

Keyword arguments

• outside::Symbol = :error — policy for locations falling outside the mesh domain:
  • :error (default) — throw ArgumentError on the first outside point.
  • :zero — emit an empty row (zero interpolation weight); the location still occupies an output row.
• atol::Real = 0.0 — extra tolerance when classifying a location as outside. A point with signed barycentric coordinates ≥ -atol in its nearest triangle is accepted. Useful when observation coordinates lie slightly outside the numerical mesh boundary due to rounding.

MeshProjector1D{T, S}

Sparse linear-interpolation projector from a 1D segment mesh to observation locations on the line. For a field u defined on mesh vertices, the interpolated value at observation i is (A * u)[i], where each row of A has at most two nonzeros — the 1D barycentric weights of observation i in its enclosing segment, summing to 1.

Fields

• mesh::INLAMesh1D — the source mesh.
• locations::Vector{T} — observation coordinates.
• A::SparseMatrixCSC{T, Int} — n_obs × n_vertices sparse projector.

Construct via MeshProjector1D(mesh, locations).

MeshProjector1D(mesh, locations; outside = :error, atol = 0.0) -> MeshProjector1D

Build a sparse linear-interpolation projector from mesh to the points in locations.

Keyword arguments

• outside::Symbol = :error — policy for locations falling outside [first(mesh.points), last(mesh.points)]:
  • :error — throw ArgumentError on the first outside point.
  • :zero — emit an empty row.
• atol::Real = 0.0 — extra tolerance when classifying as outside; a location within atol of either endpoint is accepted and clamped.

PCMatern{D}(; range_U, range_α, sigma_U, sigma_α)
PCMatern(...)                    # alias for PCMatern{2}(...)

Joint penalised-complexity (PC) prior on the Matérn range ρ and marginal standard deviation σ, following Fuglstad, Simpson, Lindgren and Rue (2019). The prior is elicited by the two tail probabilities

P(ρ < range_U) = range_α,     P(σ > sigma_U) = sigma_α,

with range_U, sigma_U > 0 and both range_α, sigma_α ∈ (0, 1).

The dimension D is carried as a type parameter so the density formula specialises at compile time. D = 2 is the default and matches SPDE2; D = 1 is used by SPDE1D.

Scope

This prior assumes positive smoothness ν = α - D/2 > 0:

• D = 1 ⇒ valid for α ∈ {1, 2, …} (ν ∈ {0.5, 1.5, …}).
• D = 2 ⇒ valid for α ∈ {2, 3, …} (ν ∈ {1, 2, …}). The case α = 1, D = 2 (ν = 0) is degenerate and rejected by SPDE2.

Fractional-α support (Bolin–Kirchner) is deferred to v0.3.

Density

Marginally, the PC-Matern prior is a product of two independent pieces:

• Range: π(ρ) = (D/2) · λ_ρ · ρ^(-D/2-1) · exp(-λ_ρ · ρ^(-D/2)), with λ_ρ = -log(range_α) · range_U^(D/2).
• Sigma: π(σ) = λ_σ · exp(-λ_σ · σ), with λ_σ = -log(sigma_α) / sigma_U.

See pc_matern_log_density for the evaluator used by SPDE2 / SPDE1D.

Defaults

Defaults (range_U = 1.0, range_α = 0.05, sigma_U = 1.0, sigma_α = 0.01) follow R-INLA's inla.spde2.pcmatern. Always override these with problem-specific scales.

SPDE1D{α}(points, segments; pc = PCMatern{1}())
SPDE1D{α}(mesh::INLAMesh1D; pc = PCMatern{1}())

Integer-α SPDE–Matérn component on a 1D segment mesh, implementing the AbstractLatentComponent contract from LatentGaussianModels.

The latent field has length num_vertices(mesh) and follows a GMRF with sparse precision Q(τ, κ) assembled from the 1D FEM matrices (C, G₁, C̃, G₂). Internal hyperparameters are θ = [log τ, log κ] per ADR-013; the user-facing scale is the Matérn (ρ, σ) pair with mapping

ρ = √(8ν) / κ,
σ² = Γ(ν) / (Γ(ν + 0.5) · (4π)^(1/2) · κ^(2ν) · τ²)         (d = 1)

where ν = α - d/2. Both α = 1 (⇒ ν = 0.5, exponential covariance) and α = 2 (⇒ ν = 1.5) are supported; the 1D PC-Matern density is non-degenerate at both.

Fields

• fem::FEMMatrices — the 1D FEM matrices, assembled once at construction.
• graph::GMRFGraph — mesh adjacency, derived from the off-diagonal pattern of C. Used for diagnostics and component-machinery compatibility.
• pc::PCMatern{1} — joint PC prior on (ρ, σ) with D = 1.

Example

mesh = inla_mesh_1d([0.0, 1.0, 2.0]; max_edge = 0.1)
spde = SPDE1D(mesh; α = 2,
              pc = PCMatern{1}(range_U=0.5, range_α=0.05,
                               sigma_U=1.0, sigma_α=0.01))
Q = precision_matrix(spde, [0.0, 0.0])   # (log τ, log κ) = (0, 0)

SPDE2{α}(mesh::INLAMesh; pc = PCMatern())
SPDE2{α}(points, triangles; pc = PCMatern())

Integer-α SPDE–Matérn component on a 2D triangular mesh, implementing the AbstractLatentComponent contract from LatentGaussianModels.

The latent field has length n_vertices(mesh) and follows a GMRF with sparse precision Q(τ, κ) assembled from the FEM matrices (C, G₁, C̃, G₂) (milestone M1). Internal hyperparameters are θ = [log τ, log κ] per ADR-013; the user-facing scale is the Matérn (ρ, σ) pair with mapping

ρ = √(8ν) / κ,        σ² = 1 / (4π · ν · κ^(2ν) · τ²)         (d = 2)

where ν = α - d/2. In v0.1 only α = 2 (⇒ ν = 1) is supported alongside the PC-Matern prior; α = 1 is deferred (ν = 0 in 2D is degenerate for the PC-Matern density).

Fields

• fem::FEMMatrices — the M1 FEM matrices (C, G₁, C̃, G₂), assembled once at construction.
• graph::GMRFGraph — mesh adjacency, derived from the off-diagonal pattern of C. Used for diagnostics, visualisation and compatibility with the LGM component machinery.
• pc::PCMatern — joint PC prior on (ρ, σ).
• mesh::Union{INLAMesh, Nothing} — the source mesh when the component was constructed via SPDE2(::INLAMesh; …); nothing when constructed from raw (points, triangles). The retained mesh is what MeshProjector(spde.mesh, locations) and the LGMFormula f((east, north), spde) extension consume.

Example

mesh = inla_mesh_2d(loc; max_edge = 0.1)
spde = SPDE2(mesh;
             pc = PCMatern(range_U=0.5, range_α=0.05,
                           sigma_U=1.0, sigma_α=0.01))
Q = precision_matrix(spde, [0.0, 0.0])   # (log τ, log κ) = (0, 0)

# Raw constructor — mesh field is `nothing`, projector calls then need
# an explicit mesh argument.
spde_raw = SPDE2(mesh.points, mesh.triangles)

SPDE2(mesh::INLAMesh; α = 2, pc = PCMatern())

Assemble an SPDE2 component directly on an INLA mesh and retain mesh in the component's mesh field. Downstream projector-aware code (MeshProjector(spde.mesh, locations), the LGMFormula f((east, north), spde) extension) requires the mesh be captured this way; the raw SPDE2(points, triangles; …) constructor sets mesh = nothing for back-compat.

SPDE2NonStationary(points, triangles; α=2, B_τ, B_κ, prior=GaussianBasisPrior(p_τ + p_κ))
SPDE2NonStationary(mesh::INLAMesh; α=2, B_τ, B_κ, prior=...)

Non-stationary integer-α SPDE-Matérn component on a 2D triangular mesh. Per-vertex log τ_v = (B_τ θ_τ)_v and log κ_v = (B_κ θ_κ)_v let the Matérn precision and range vary across the domain — the parameterisation behind R-INLA's inla.spde2.matern(..., B.tau, B.kappa, theta.prior.mean, theta.prior.prec).

Currently restricted to α = 2 (Matérn smoothness ν = 1 in 2D); α = 1 is deferred to a follow-up.

Hyperparameter layout

The internal hyperparameter vector concatenates τ-coefficients and κ-coefficients in that order:

θ = [θ_τ_1, …, θ_τ_{p_τ}, θ_κ_1, …, θ_κ_{p_κ}]

with length(θ) = p_τ + p_κ = size(B_τ, 2) + size(B_κ, 2). The GaussianBasisPrior on θ (ADR-028) factorises independently across coefficients with per-coefficient mean and prec, matching R-INLA's theta.prior.mean / theta.prior.prec.

Fields

• fem::FEMMatrices — precomputed M1 FEM matrices.
• graph::GMRFGraph — mesh adjacency.
• B_τ::Matrix{T} — (n_v, p_τ) τ-basis matrix.
• B_κ::Matrix{T} — (n_v, p_κ) κ-basis matrix.
• prior::GaussianBasisPrior — independent Gaussian on basis coefficients.
• mesh::Union{INLAMesh, Nothing} — the source mesh when the component was constructed via SPDE2NonStationary(::INLAMesh; …); nothing when constructed from raw (points, triangles). The retained mesh is what MeshProjector(spde.mesh, locations) and raster-side predict_raster(model, res, template) consume.

Stationary limit

If B_τ = ones(n_v, 1) and B_κ = ones(n_v, 1) and θ = [log τ, log κ], the precision reduces to SPDE2's stationary form τ²·(κ⁴ C̃ + 2κ² G₁ + G₂). This makes stationary-recovery a clean sub-test of the LRL §3.2 oracle fixture.

Example

mesh = inla_mesh_2d(loc; max_edge = 0.1)
n_v = num_vertices(mesh)

# Two-region piecewise-constant κ, constant τ:
region = [x < 0.5 ? 1 : 2 for (x, _) in eachrow(mesh.points)]
B_κ = hcat([region .== 1, region .== 2]...) .|> Float64
B_τ = ones(n_v, 1)

prior = GaussianBasisPrior(
    mean = [0.0, 0.0, 0.0],
    prec = [1.0e-3, 1.0, 1.0]            # wide on intercept, tight on regions
)
spde = SPDE2NonStationary(mesh.points, mesh.triangles;
                          B_τ = B_τ, B_κ = B_κ, prior = prior)

GMRFs.constraints(::SPDE1D) -> NoConstraint

The SPDE precision is strictly positive-definite for all (τ, κ) > 0 on a non-degenerate 1D mesh, so no hard linear constraint is required.

GMRFs.constraints(::SPDE2NonStationary) -> NoConstraint

The non-stationary SPDE precision is SPD as long as τ_v, κ_v > 0, which the construction enforces.

GMRFs.constraints(::SPDE2) -> NoConstraint

The SPDE precision is strictly positive-definite for all (τ, κ) > 0 (κ² C + G₁ and κ⁴ C̃ + 2κ² G₁ + G₂ are SPD), so no hard linear constraint is required.

_as_boundary_matrix(b) -> Matrix{Float64} | Nothing

Coerce the public boundary argument of inla_mesh_2d to the internal k × 2 matrix form.

• nothing passes through (boundary defaults to the convex hull of loc).
• AbstractMatrix{<:Real} is copied to Matrix{Float64}.
• GeoInterface PolygonTrait / LineStringTrait / LinearRingTrait geometries are handled by INLASPDEGeoInterfaceExt once GeoInterface is loaded.

_as_location_matrix(l) -> Matrix{Float64} | Nothing

Coerce the public loc / locations arguments of inla_mesh_2d and MeshProjector to the internal n × 2 matrix form.

• nothing passes through.
• AbstractMatrix{<:Real} is copied to Matrix{Float64}.
• GeoInterface MultiPointTrait / PointTrait geometries and vectors of PointTrait geometries are handled by INLASPDEGeoInterfaceExt once GeoInterface is loaded.

_barycentric(x, y, x1, y1, x2, y2, x3, y3) -> (λ1, λ2, λ3)

Barycentric coordinates of (x, y) with respect to the triangle (p1, p2, p3). Sums to 1 exactly by construction.

_interior_only(pts, polygon) -> Matrix{Float64}

Rows of pts strictly inside the CCW polygon (not on the boundary).

_mesh_graph_from_C(C) -> GMRFGraph

Build the mesh adjacency graph from the sparsity pattern of the FEM mass matrix. Two vertices are adjacent iff they share a triangle — equivalently, iff the corresponding off-diagonal entry of C is nonzero.

_strip_hull_points(pts, hull) -> Matrix{Float64}

Rows of pts strictly inside hull; rows on the hull are dropped.

assemble_fem_matrices(points, triangles) -> (C, G1)

Assemble the P1 finite-element mass matrix C and stiffness matrix G₁ on a 2D triangular mesh.

The mass matrix has entries C[i,j] = ∫_Ω φ_i φ_j dx and the stiffness matrix G₁[i,j] = ∫_Ω ∇φ_i · ∇φ_j dx, where φ_i is the piecewise-linear hat basis function at vertex i on the mesh defined by points and triangles.

Arguments

• points::AbstractMatrix{<:Real} — n × 2 vertex coordinates.
• triangles::AbstractMatrix{<:Integer} — m × 3 one-based vertex indices, one row per triangle. Orientation (CW / CCW) is irrelevant; the unsigned triangle area is used.

Returns

• C::SparseMatrixCSC — n × n symmetric positive-definite mass matrix.
• G1::SparseMatrixCSC — n × n symmetric positive-semidefinite stiffness matrix. G₁ · 1 = 0 (constant-preserving).

Notes

A Meshes.jl overload accepting Meshes.SimpleMesh is introduced in M3 alongside inla_mesh_2d. This low-level method is what that overload ultimately calls.

assemble_fem_matrices_1d(points, segments) -> (C, G1)

Assemble the P1 finite-element mass matrix C and stiffness matrix G₁ on a 1D segment mesh.

The mass matrix has entries C[i,j] = ∫_Ω φ_i φ_j dx and the stiffness matrix G₁[i,j] = ∫_Ω (dφ_i/dx) · (dφ_j/dx) dx, where φ_i is the piecewise-linear hat basis function at vertex i.

Arguments

• points::AbstractVector{<:Real} — vertex coordinates (length n).
• segments::AbstractMatrix{<:Integer} — m × 2 one-based vertex indices, one row per segment. Orientation is irrelevant; segment length |x_{i_2} - x_{i_1}| is used.

Returns

• C::SparseMatrixCSC — n × n symmetric positive-definite mass matrix.
• G1::SparseMatrixCSC — n × n symmetric positive-semidefinite stiffness matrix. G₁ · 1 = 0 (constant-preserving).

Element-level formulas

For a segment of length h with endpoint indices (i, j):

• mass: C_loc = (h/6) · [2 1; 1 2]
• stiffness: G_loc = (1/h) · [1 -1; -1 1]

convex_hull_polygon(points) -> Matrix{Float64}

Return the convex-hull polygon of points (n × 2 matrix) as a k × 2 matrix in counter-clockwise order, without repeating the first vertex at the end. Uses DelaunayTriangulation.convex_hull, which emits CCW-ordered vertices.

Collinear points are excluded from the hull vertex set. Throws ArgumentError if the hull is degenerate (fewer than three distinct vertices).

cutoff_dedup(points, cutoff) -> Matrix{Float64}

Remove near-duplicate rows of points so that the returned set has pairwise distance ≥ cutoff. First occurrence is kept (greedy). Runs in O(n²) — adequate for mesh vertex counts; for 10⁶-point clouds a KD-tree version would replace this.

cutoff ≤ 0 returns Matrix{Float64}(points) unchanged.

expand_polygon(polygon, offset) -> Matrix{Float64}

Expand a simple convex polygon outward by a perpendicular distance offset. polygon is a k × 2 matrix of CCW-ordered vertices (no closing duplicate). Each vertex is shifted along the outward bisector such that each edge of the output polygon lies offset away from the corresponding input edge.

Used to construct the outer mesh domain from the convex hull of the observation locations (R-INLA's offset argument to inla.mesh.2d).

inla_mesh_1d(loc; max_edge, cutoff = 0.0, boundary = nothing) -> INLAMesh1D

Build a 1D segment mesh covering the interval implied by loc and boundary, refined so every segment is no longer than max_edge. The 1D analogue of inla_mesh_2d.

Arguments

• loc::AbstractVector{<:Real} — observation locations on the line. Need not be sorted; duplicates and points within cutoff of each other are collapsed (first-wins).
• max_edge::Real — upper bound on segment length. Required.
• cutoff::Real = 0 — minimum allowed spacing; closer points are merged.
• boundary::Union{Nothing, NTuple{2, <:Real}} = nothing — explicit domain endpoints (a, b) with a < b. Must contain all loc values. If omitted, the domain defaults to (minimum(loc), maximum(loc)).

Returns

An INLAMesh1D with sorted vertex coordinates and consecutive-segment topology.

Example

loc  = randn(50)
mesh = inla_mesh_1d(loc; max_edge = 0.1, cutoff = 0.01,
                    boundary = (-3.0, 3.0))
fem  = FEMMatrices(mesh)

inla_mesh_2d([loc]; boundary, max_edge, offset, cutoff, min_angle,
             subdivide_boundary) -> INLAMesh

Julia-native analogue of R-INLA's inla.mesh.2d: builds a constrained Delaunay triangulation refined to satisfy a maximum edge length and a minimum triangle angle.

Arguments

• loc::AbstractMatrix{<:Real} (optional, n × 2) — observation locations. If boundary is not supplied, the mesh domain defaults to the convex hull of loc, optionally expanded by offset.
• boundary::AbstractMatrix{<:Real} — k × 2 CCW polygon defining the mesh domain (no repeating closing vertex). If omitted, the domain is derived from loc.
• max_edge::Real or NTuple{2, <:Real} — upper bound on triangle edge length. A scalar applies one bound globally; a tuple (inner, outer) matches R-INLA's max.edge=c(inner, outer): inner triangles (centroid inside the inner hull) refine to inner, triangles in the buffer band refine to outer. The tuple form requires loc (and an NTuple{2} offset); use the scalar form with explicit boundary. Bound is implemented as a triangle-area constraint max_area = √3 / 4 · max_edge² (equilateral reference).
• offset::Real = 0 or NTuple{2, <:Real} — extension buffer(s) applied to the loc convex hull when boundary is not given. Scalar = single buffer (R-INLA's offset scalar). Tuple (o₁, o₂) expands the hull by o₁ to form the inner domain and by o₁ + o₂ to form the outer domain (R-INLA's offset=c(o₁, o₂)).
• cutoff::Real = 0 — minimum pairwise distance between interior points; closer points are collapsed (first-wins) before triangulation. Has no effect on boundary vertices.
• min_angle::Real = 21.0 — minimum interior triangle angle in degrees. Must be below ≈33.8° for Ruppert's algorithm to terminate (Shewchuk 2002); values much below 20° give fewer refinement steps but coarser meshes near acute corners.
• subdivide_boundary::Bool = false — when true, subdivide each boundary segment so that no input edge exceeds the prevailing max_edge. Strict per-edge enforcement on the boundary; the Ruppert area bound continues to govern the interior. Without this, long input edges of an explicit boundary survive at their input length.

Returns

An INLAMesh with points, triangles, boundary, and the underlying DelaunayTriangulation object.

Example

loc = randn(50, 2)
mesh = inla_mesh_2d(loc; max_edge = 0.3, offset = 0.5, min_angle = 25.0)

# Two-region: tight near data, coarser in the buffer band
mesh2 = inla_mesh_2d(loc; max_edge = (0.2, 0.6), offset = (0.3, 1.0))

log_basis_prior_density(prior::GaussianBasisPrior, θ) -> Real

Log density of the independent-Gaussian basis prior at coefficient vector θ, including the normalising constant ½ ∑_k log(λ_k / 2π).

Throws ArgumentError on length mismatch.

lumped_mass(C) -> C_lumped

Diagonal mass lumping: replace the full mass matrix C with the diagonal matrix whose i-th entry is the sum of row i of C. This preserves the total mass (sum(C_lumped) == sum(C), which equals the domain area for P1 FEM) while making C_lumped trivially invertible.

Used in SPDE-Matérn precision assembly for α = 2 to keep G₂ = G₁ C̃⁻¹ G₁ sparse, following Lindgren, Rue and Lindström (2011, Appendix C).

Arguments

• C::AbstractSparseMatrix — the full mass matrix from assemble_fem_matrices.

Returns

• C_lumped::SparseMatrixCSC — diagonal sparse matrix with the same size and element type as C.

nonconvex_hull_polygon(points; α) -> Matrix{Float64}

Return an α-shape (Edelsbrunner-Mücke) of points (n × 2) as a CCW polygon, without repeating the first vertex. Larger α yields the convex hull in the limit; smaller α follows the point cloud more tightly.

Arguments

• points::AbstractMatrix{<:Real} — n × 2 point cloud.
• α::Real — circumradius bound. A Delaunay triangle is kept iff its circumradius is ≤ α. If α is omitted it defaults to twice the median nearest-neighbour distance, a heuristic that gives a "reasonable concavity" on roughly uniform point clouds.

Errors

• ArgumentError if points has fewer than 3 rows or is not 2D.
• ArgumentError if the alpha-shape is empty (α too small) or multi-component (the simply-connected boundary tracer fails); raise α and try again.

See also: convex_hull_polygon for the convex-hull case.

pc_matern_log_density(pc::PCMatern{D}, log_ρ, log_σ) -> Real

Log density of a D-dimensional PC-Matern prior evaluated on the (log ρ, log σ) scale. Includes the log ρ + log σ Jacobian that converts from density on (ρ, σ) to density on (log ρ, log σ).

Used internally by SPDE2 (D = 2) and SPDE1D (D = 1); the change of variables from (log ρ, log σ) to the internal (log τ, log κ) has unit Jacobian and does not contribute an extra term (ADR-013).

scimloperator(P::MeshProjector)

Wrap the sparse projection matrix as a SciMLOperators.MatrixOperator for use with the SciMLOperators algebra (lazy adjoint / composition with other operators).

spde_internal_scale(c::SPDE1D{α}, ρ, σ) -> (log τ, log κ)

Inverse of spde_user_scale for d = 1, α ∈ {1, 2}.

spde_internal_scale(c::SPDE2{α}, ρ, σ) -> (log τ, log κ)

Inverse of spde_user_scale. Returns the internal (log τ, log κ) pair corresponding to user-scale Matérn (ρ, σ).

spde_precision(α, τ, κ, C, G1[, C_lumped, G2]) -> Q

Stateless form: assemble Q(α, τ, κ) directly from the raw FEM matrices. Missing C_lumped and G2 are derived on the fly. Prefer spde_precision(::FEMMatrices, ...) in hot loops — the FEMMatrices constructor precomputes C̃ and G₂ once.

spde_precision(fem::FEMMatrices, α, τ, κ) -> Q

Assemble the SPDE-Matérn precision matrix on user-scale parameters (τ, κ). Supported orders are α ∈ {1, 2}.

• α = 1: Q = τ² · (κ² C + G₁) (Matérn smoothness ν = 0).
• α = 2: Q = τ² · (κ⁴ C̃ + 2κ² G₁ + G₂) (ν = 1), using the lumped mass matrix C̃ — this matches R-INLA's implementation per Lindgren-Rue-Lindström (2011, Appendix C).

Fractional α is deferred to v0.3 (Bolin–Kirchner 2020 rational approximation).

spde_precision_nonstationary(fem::FEMMatrices, α, τ_v, κ_v) -> Q

Sparse SPDE-Matérn precision with per-vertex τ_v, κ_v (vectors of length n_vertices). Implements R-INLA's inla.spde2.matern formula:

• α = 2: Q = D_τ · (D_κ²·C̃·D_κ² + D_κ²·G₁ + G₁·D_κ² + G₂) · D_τ, where D_τ = diag(τ_v), D_κ² = diag(κ²_v), and C̃ is the lumped (diagonal) mass matrix. D_κ²·C̃·D_κ² is just Diagonal(κ⁴_v ⊙ c̃_v); the cross-terms D_κ²·G₁ + G₁·D_κ² symmetrise to G₁_{ij}·(κ²_i + κ²_j).

In the stationary limit (τ_v ≡ τ, κ_v ≡ κ) this reduces to τ² · (κ⁴ C̃ + 2κ² G₁ + G₂), matching the SPDE2 stationary precision. Only α = 2 is supported for the non-stationary form; α = 1 and fractional α are deferred (no R-INLA-parity oracle for either yet).

spde_user_scale(c::SPDE1D{α}, θ) -> (ρ, σ)

Convert internal θ = [log τ, log κ] to the user-facing Matérn pair (ρ, σ) for d = 1, α ∈ {1, 2}.

spde_user_scale(c::SPDE2{α}, θ) -> (ρ, σ)

Convert the internal θ = [log τ, log κ] to the user-facing Matérn pair (ρ, σ) using the fixed mapping for α = 2 in 2D.

stiffness_squared(G1, C_lumped) -> G2

Construct G₂ = G₁ · C̃⁻¹ · G₁, where C̃ is a lumped (diagonal) mass matrix. This is the sparse approximation of G₁ · C⁻¹ · G₁ used for α = 2 SPDE precision.

Throws ArgumentError if any diagonal entry of C_lumped is zero — this indicates a vertex with zero associated area, i.e. a degenerate mesh.

subdivide_polygon(polygon, max_edge) -> Matrix{Float64}

Insert Steiner points along each edge of polygon (a k × 2 CCW matrix) so that every consecutive pair of returned vertices is at most max_edge apart. Returns the augmented polygon (no closing duplicate).

Used by inla_mesh_2d(...; subdivide_boundary = true) to enforce a strict per-edge max_edge bound on the input boundary, before Ruppert refinement runs over the interior. Without this, the boundary edges survive at their input length and the area-based Ruppert bound only softly limits interior edges.

initial_hyperparameters(c::SPDE2NonStationary) -> Vector

Default initial point for the BFGS θ-mode finder: the prior mean. This is the natural starting point for a Gaussian prior and matches R-INLA's theta.initial = theta.prior.mean default.

log_hyperprior(c::SPDE1D{α}, θ) -> Real

PC-Matern log-prior density evaluated at θ = [log τ, log κ]. Maps (log τ, log κ) → (log ρ, log σ) per the d=1 formulas and calls pc_matern_log_density. The (log τ, log κ) ↔ (log ρ, log σ) Jacobian has absolute determinant 1 (ADR-013).

log_hyperprior(c::SPDE2NonStationary, θ) -> Real

Independent-Gaussian log density on θ = [θ_τ; θ_κ] via log_basis_prior_density.

log_hyperprior(c::SPDE2, θ) -> Real

PC-Matern log-prior density evaluated at θ = [log τ, log κ] on the internal scale. Maps (log τ, log κ) → (log ρ, log σ) and calls pc_matern_log_density. The (log τ, log κ) ↔ (log ρ, log σ) Jacobian has absolute determinant 1 and contributes no extra term (ADR-013).

log_normalizing_constant(c::SPDE1D, θ) -> Real

R-INLA-style log normalizing constant for the proper Gaussian SPDE prior, evaluated at internal θ = [log τ, log κ]:

log NC = -½ d log(2π) + ½ log|Q(θ)|

where d = length(c) is the number of mesh vertices and Q(θ) is the SPDE precision.

log_normalizing_constant(c::SPDE2, θ) -> Real

R-INLA-style log normalizing constant for the proper Gaussian SPDE prior, evaluated at internal θ = [log τ, log κ]:

log NC = -½ d log(2π) + ½ log|Q(θ)|

where d = length(c) is the number of mesh vertices and Q(θ) is the SPDE precision (κ⁴ C̃ + 2κ² G₁ + G₂ scaled by τ², for α = 2). Required by the marginal-likelihood formula in LatentGaussianModels.jl's Laplace inference (commit 635cbb9 / ADR-style change).

log_normalizing_constant(c::SPDE2NonStationary, θ) -> Real

Proper-Gaussian normalizing constant -½ d log(2π) + ½ log|Q(θ)|, matching SPDE2. Required for marginal-likelihood comparison with R-INLA.

precision_matrix(c::SPDE1D{α}, θ) -> SparseMatrixCSC

Sparse SPDE precision at θ = [log τ, log κ], delegating to spde_precision(fem, α, τ, κ).

precision_matrix(c::SPDE2NonStationary, θ) -> SparseMatrixCSC

Sparse non-stationary precision at θ = [θ_τ; θ_κ] of length p_τ + p_κ. Computes per-vertex τ_v = exp(B_τ θ_τ) and κ_v = exp(B_κ θ_κ), then dispatches to spde_precision_nonstationary.

precision_matrix(c::SPDE2{α}, θ) -> SparseMatrixCSC

Sparse SPDE precision at θ = [log τ, log κ], delegating to the M1 assembly spde_precision(fem, α, τ, κ).