Estimate distances on a declared quadratic graph surface

Internal common interface for paired endpoint queries on \(F(u)=(u,u^T A u)\). A method must be selected explicitly.

Usage

quadform_geodesics(
  A,
  from,
  to,
  domain,
  method,
  control = list(),
  return.paths = FALSE
)

Arguments

A

Finite exactly symmetric 2 by 2 matrix. The Delaunay and radius graph references also accept 3 by 3 and 4 by 4 matrices, representing a single quadratic height over the corresponding domain dimension.

from, to

Two-element numeric vectors or matching two-column matrices of domain coordinates. Rows define pairs; no recycling or Cartesian product.

domain

A list with kind = "ball", center and radius, or kind = "box", lower and upper. The domain is closed and convex.

method

One of "single_point", "local_network", "best_of_six", "grid_dijkstra", "paraboloid_clairaut", "geodesic_shooting" or "geodesic_collocation", "boundary_optimization", "polyhedral_mesh", "delaunay_graph" or "radius_graph". No partial matching.

control

Named list of method-specific controls; unknown names are errors. Refinements accept the corresponding lower-level function's controls, with cache_edges = 0 and trace = FALSE by default. Grid controls are grid_size = c(33,33) (3–257 per axis), direction_radius = 2 (1–8), endpoint_neighbors = 8 (1–32), include_direct = TRUE, max_edges = 1000000 (0–2000000), keep_graph = FALSE and max_seconds = Inf. Clairaut controls are iterations = 80 (1–256), path_samples = 257 (3–4097), domain_check_depth = 16 (0–20), angle_tolerance = 1e-12 (1e-15–1e-8) and max_seconds = Inf. The two continuous methods accept ode_tolerance = 1e-6 (1e-10–1e-2), endpoint_tolerance = 1e-8 (1e-12–1e-3), path_tolerance = 1e-7 (1e-10–1e-2), iterations = 30 (1–100), continuation_steps = 8 (1–128), continuation_attempts = 128 (1–1024), initial_nodes = 17 (3–257), max_nodes = 1025 (initial_nodes–4097), max_path_vertices = 4097 (3–65537), max_evaluations = 500000 (0–10000000), domain_check_depth = 20 (0–24), initial_bends = c(0,-1,1) (1–9 distinct values between -4 and 4), and max_seconds = Inf. Shooting additionally accepts integration_tolerance = 1e-10 (1e-13–1e-5). Boundary optimization accepts initial_edges = 8 (2–64), levels = 2 (1–3), evaluations_per_start = 2000 (1–100000), max_evaluations = 20000 (0–1000000), position_tolerance = 1e-6 (1e-10–1e-2), initial_step = 0.05 (1e-6–0.5), initial_bends = c(0,-0.5,0.5) (1–9 distinct values in -4–4), initial_path = NULL (otherwise a feasible two-column path with 2–65 rows and exactly matching endpoints), and max_seconds = Inf. Reference graphs and meshes accept resolution (2–33; defaults 17, 5, 3 in domain dimensions 2, 3, 4), vertices = NULL (otherwise a finite domain-coordinate matrix), max_vertices, max_edges = 200000 (0–2000000), max_pair_checks = 1000000 (0–10000000), max_seconds = Inf and keep_graph = FALSE. Default vertex caps are 4096 for radius graphs, 2048/512/128 for 2D/3D/4D Delaunay graphs, and 1024 for meshes; the Delaunay and mesh caps cannot be raised. Radius graphs accept radius, defaulting to twice the largest domain bounding-box width divided by resolution minus one, in domain units. Mesh controls additionally include max_faces = 2048 (1–4096), max_propagations = 100000 (0–1000000), max_intervals = 200000 (1–1000000), and keep_mesh = FALSE. Time and calculation limits apply to each pair, not the whole batch.

return.paths

Whether to retain paths, including those in nested refinement results. Default FALSE. No files or checkpoints are written.

Value

A list with method, geometry, domain, resolved control, row-preserving summary, and detailed pair results. Status is "candidate", "partial", "unsupported" or "failed". Summary length is the numerical distance estimate; length_error_estimate concerns numerical length evaluation, not discretization or optimization error. approximation_error is missing and global_optimality_certified is FALSE. Detailed results retain backend_result, including all six outcomes and total costs for the best-of-six strategy. A completed budget is not an accuracy certificate. User interrupts and input errors propagate.

Details

Grid edges use the native analytic lifted-segment lengths. Primitive lattice offsets up to direction_radius supply directions; radius one is the eight-neighbor stencil. Finer spacing alone does not remove its directional bias. Off-grid endpoints attach to nearest retained vertices. The optional direct connector is an explicit graph edge, not a fallback. Search uses strict binary64 comparisons; the selected path total is recomputed with the native exact accumulator.

Clairaut requires A = a I exactly. It solves radial-monotone or turning branches of the continuous surface metric, with separate flat, apex, radial and antipodal cases. Origin-centered disks use a whole-path radial envelope; other domains require conservative radial/angular section enclosures. Leaving the domain or unresolved containment returns unsupported, not a boundary-constrained replacement. Curved polar cases outside the initial numeric scope (minimum endpoint radius below 1e-100, maximum above 1e100, or twice absolute curvature times maximum radius above 1e6) are explicitly unsupported. These restrictions do not apply to flat, coincident or apex cases.

Boundary optimization uses NLopt COBYLA to vary interior vertices of a lifted domain polyline. Bounds constrain a rectangle; a quadratic inequality per vertex constrains a disk. All accepted vertices must also pass the strict floating-point domain check, independently of NLopt's constraint status. Convexity then keeps every lifted domain segment in D. The direct connector is an explicit initial candidate. Each accepted replacement must be shorter by more than the sum of the two numerical length-error estimates. A feasible supplied initial path competes with the direct connector, rather than silently replacing it.

Optimization coordinates are centered and divided by the domain radius or largest box width. The position tolerance and initial step use these normalized units. Starts add sine-shaped normal displacements controlled by initial_bends, projected into D only to construct feasible initial guesses. Subsequent levels bisect the retained path, up to 257 vertices. Objective-evaluation or time limits, and incomplete individual searches, retain the best feasible path as partial. A position-tolerance stop is not a stationarity or global-optimality certificate. No states are saved.

Delaunay and radius references insert endpoints into either the supplied vertex set or a regular lattice clipped to D, remove exact duplicates, and sort vertices lexicographically. Delaunay uses geometry::delaunayn with Qhull options Qt Qbb Qc Qz, without randomized perturbation. The package geometry is an optional dependency for Delaunay and mesh methods. Its simplices supply all one-skeleton edges. Radius graphs instead join all pairs at Euclidean domain distance at most radius. There is no direct-edge override or silent connectivity repair. Edge weights are analytic lifted straight-segment lengths, not ambient chords. Boost Dijkstra selects a graph path, whose length is recomputed with the exact accumulator. In 3D and 4D, the same analytic integral is used with higher-dimensional endpoint slopes and propagated numerical error estimates. The entire triangulation and pair preparation count toward total time, but Qhull is an indivisible call: its time allowance is checked before and after it, not enforced as an operating-system deadline.

The polyhedral method uses Kirsanov's implementation of the Mitchell–Mount–Papadimitriou triangle-interior propagation algorithm on a two-dimensional Delaunay triangulation lifted at its vertices. Its domain is the triangulated convex hull of these points, a subset of D, not the whole disk when D is circular. The upstream small-interval cutoff is 1e-6 times edge length. Live intervals and propagation counts are bounded; user interrupts release owned intervals before propagating to R. Degenerate triangles or unresolved numerical invariants return no path. In particular, inserting very close endpoints can create angles below the mesh engine's 1e-5-radian limit. Rejections retain the offending angle and the endpoint separation relative to mesh scale; endpoints are never merged or moved to make the triangulation pass.

For polyhedral_surface_polyline, surface_path is a path across flat mesh triangles and path is its domain projection. length measures that mesh path, not a smooth-surface path. The separate diagnostic smooth_lifted_length measures the lifted domain polyline obtained from its projection. These two numbers must not be interchanged in comparisons. Mesh length_error_estimate is missing, not zero; the algorithm's name does not imply a rigorous certificate. Kept mesh vertices and triangles are available under backend_result$mesh when keep_mesh = TRUE.

Shooting adjusts the initial velocity of the geodesic equation using its variational equations and damped Newton iterations. Integration uses Boost's controlled Dormand–Prince 5(4) stepper. Collocation independently solves the four-state boundary-value equations with fourth-order Lobatto collocation, a sparse analytic Jacobian and damped Newton steps. Both use the equation $$u''=-\frac{2Au}{1+\|2Au\|^2}(u'^T(2A)u').$$ Coordinates are translated to the canonical first endpoint and divided by the domain-coordinate distance between endpoints. Thus endpoint tolerance is relative to that distance, not to the domain's diameter.

A zero initial bend starts with a flat surface and continues through increasing multiples of A to the requested surface. The initial step is 1/continuation_steps; successful steps grow by 1.5 and unsuccessful steps halve. Collocation refines and checks the equation residual at each continuation stage before advancing. Nonzero bends are additional direct starts on the full surface. If delta joins the normalized endpoints and J rotates by 90 degrees, shooting starts with velocity delta + bend * J delta; collocation starts with the curve a + t * delta + bend * sin(pi*t) * J delta, with its derivative. initial_nodes sets the initial collocation mesh; shooting uses adaptive integration nodes. The shortest accepted start is retained, with overlapping numerical length-error intervals resolved by start order.

ode_tolerance bounds a sampled, componentwise scaled residual of cubic Hermite state interpolation, not a rigorous solution error. Intervals are tested at fractions 0.2113248654, 0.5 and 0.7886751346. Only intervals that fail this check are bisected. Collocation resolves the enlarged system; shooting integrates from each affected interval's left state to its midpoint, retaining the existing nodes. Shooting's relative integration tolerance is integration_tolerance, with absolute tolerance one hundredth of it. iterations limits each Newton solve; max_evaluations counts equation evaluations across all starts, continuation and residual checks. Collocation's Newton stopping test divides each differential-equation row by one plus the Simpson-weighted absolute right-hand side; endpoint rows retain scale one. The Newton step and backtracking merit function remain unscaled. This avoids imposing an unattainable absolute equation tolerance on steep surfaces without relaxing endpoint accuracy.

Continuous candidates are converted to lifted domain polylines. Recursive cubic Bezier control-hull checks establish containment of the interpolated domain curve using floating-point arithmetic. Midpoint and quarter-point subdivision compares lifted-connector lengths using a total allowance of path_tolerance times the direct connector length. Each accepted section emits its coarse connector; midpoint and quarter points are used for the check, not automatically added to the output. This is a refinement diagnostic, not a certified discretization bound. length_error_estimate concerns only the returned polyline's length. Endpoint residuals are measured before fixing the returned endpoints exactly; reported residuals have the original domain-coordinate units.

A stationary path need not be the shortest path. The diagnostics retain each start's outcome, the selected start, equation and endpoint residuals, and the direct connector length for comparison. There is no silent fallback to a different method and no boundary-following solver. A retained candidate with unresolved starts or an exhausted overall budget is partial; no usable path gives failed, or unsupported when all completed attempts leave the domain and none remains unresolved. A mixture of outside curves and unresolved starts is failed, not evidence that the geometry is unsupported. The first observed outside point is retained for each rejected start. These methods do not handle more than one quadratic form or domain dimension other than two. Nontrivial pairs require endpoint separation between 1e-100 and 1e100, maximum absolute entry of 2 * separation * A at most 2048 and maximum absolute component of 2 * A * canonical_from at most 1e6. Out-of-scope cases return unsupported; flat and coincident pairs bypass these restrictions.

For analytic_clairaut_curve, path is a sampled display of the curve, not a polyline whose length equals length. curve_parameters describe the authoritative curve; the numerical endpoint residual is retained. Parameterize it by t = t_from + s * t_span, for s between zero and one; use the retained t_span, not subtraction of rounded endpoints. With \(r = \sqrt{c^2 + t^2}\), its polar angle is theta_turn plus angular_sign * sign(t) * H(c, abs(t)), where \(H(c,t)=kc\,\mathrm{asinh}(kt/\sqrt{1+k^2c^2})+ \mathrm{atan2}(t,c\sqrt{1+k^2c^2+k^2t^2})\). caller_reversed indicates whether to traverse this parameterization backward to follow the caller's endpoint order. For lifted_domain_polyline, consecutive domain points specify lifted straight segments. All numerical error estimates are nonrigorous. The historical development adapter registry is independent of this package interface and is not activated by these methods. Grid, Clairaut and continuous results are symmetric under endpoint reversal when calculation limits, rather than wall-clock timing, determine completion. A single randomized refinement run is not required to be: its internal orientation is randomized relative to the supplied endpoints. The wrapper preserves the lower-level solver's behavior and seed semantics.