--- title: "Fermat distances: sampling density, curvature, and geodesic recovery" output: rmarkdown::html_vignette: toc: true self_contained: true math_method: mathml css: compact-guides.css vignette: > %\VignetteIndexEntry{Fermat distances: sampling density, curvature, and geodesic recovery} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup,include=FALSE} knitr::opts_chunk$set(collapse=TRUE,comment='#>',message=FALSE,warning=FALSE, out.width='100%',fig.align='center') # Prefer this checkout's snapshot when rendering the maintained source. bundle <- '../inst/extdata/fermat-density-curvature/benchmark.rds' if(!file.exists(bundle)) bundle <- system.file('extdata','fermat-density-curvature', 'benchmark.rds',package='dgraphs',mustWork=TRUE) b <- readRDS(bundle); m <- b$metrics; s <- b$summary; power_data <- b$power_tuning surface_names <- c(paraboloid_a1='Paraboloid, $a=1$',paraboloid_a2='Paraboloid, $a=2$', paraboloid_a4='Paraboloid, $a=4$',saddle_a1='Saddle, $a=1$',saddle_a2='Saddle, $a=2$', saddle_a4='Saddle, $a=4$',swiss_roll='Swiss roll') helper <- '../inst/doc-tools/compact-figures.R' if(!file.exists(helper)) helper <- system.file('doc-tools','compact-figures.R',package='dgraphs',mustWork=TRUE) source(helper,local=TRUE) figure <- function(name) knitr::include_graphics(dgraphs_document_figure( paste0('figures/fermat-density-curvature/',name,'.png')), rel_path=FALSE) ``` Fermat distances can improve on ambient straight-line distances, but the useful power depends on geometry and sampling density. This article examines that interaction before any embedding into two or three dimensions. It compares complete-graph Fermat distances, ambient chords, and ordinary Euclidean nearest-neighbor graph distances on area-uniform quadratic surfaces and a Swiss roll. The main findings are: - Intermediate powers matter: a grid containing only integer powers can miss substantial improvements already at small sample sizes. - Denser samples improve many higher-power estimates, but no tested power is best across all surfaces. - Ordinary Euclidean kNN shortest paths are more accurate at the largest sample size on these uniformly sampled examples. Density-sensitive geometry and recovery of ordinary surface geodesics are different objectives. - Calibration-only tuning at 300 points yields useful powers across sample sizes, but does not identify a power that remains optimal as sampling grows denser. - Numerical reference error matters on the steepest paraboloid. This complements [Fermat distances and their supporting graphs on surfaces](fermat-surfaces.html) and the [Fermat distance introduction](fermat-distances.html). The benchmark results were completed on `r b$provenance$completion$completed`; power-selection results were completed on `r power_data$completion$completed`; this reading copy was built on `r format(Sys.time(),tz='America/New_York',usetz=TRUE)`. Rendering reads the frozen benchmark, saved figures, and saved surface coordinates. The teaching examples below are displayed without execution; rendering does not rerun the numerical analyses or reference searches. ## What distance is being estimated? For observations $X_1,\ldots,X_n\in\mathbb R^3$, the unrooted Fermat distance is $$ D_p(i,j)=\min_{\pi:i\leadsto j} \sum_{(u,v)\in\pi}\lVert X_u-X_v\rVert_2^p,\qquad p\geq1. $$ The minimum ranges over paths through any observed points, including the direct edge. At $p=1$, the triangle inequality gives $D_1(i,j)=\lVert X_i-X_j\rVert_2$: adding intermediate observations cannot shorten the ambient chord. For $p>1$, several short steps can cost less than a long step. More intermediate points can therefore change both the preferred route and its cost. The units are length to power $p$; this article takes no $p$-th root. A surface geodesic measures length along the underlying surface. We compare Fermat distances to that intrinsic reference after fitting one multiplicative scale. These are finite-sample experiments, not a proof of asymptotic convergence or an estimate of its rate. ## A compact example using the existing area sampler Consider the saddle $(u,v,a(u^2-v^2))$ over $u^2+v^2\leq1$. For a quadratic form matrix $A$, the surface-area element is $\sqrt{1+\lVert2A(u,v)^T\rVert_2^2}$. Uniform draws in the parameter disk are therefore not uniform on the surface. Use the existing dgraphs sampler: ```{r small-sample,eval=FALSE} library(dgraphs) A <- diag(c(2,-2)) geometry <- synthetic.quadform(2L,3L,forms=list(A)) sampling <- synthetic.sampling.quadform.lab('disk',extent=1,mode='area') sample <- sample.synthetic.geometry(geometry,sampling,n=300L,seed=62117L) X <- sample$predictors endpoints <- 1:12 D150 <- fermat.distances(points=X[1:150,],p=1.25,backend='implicit', sources=endpoints,targets=endpoints,rooted=FALSE) D300 <- fermat.distances(points=X,p=1.25,backend='implicit', sources=endpoints,targets=endpoints,rooted=FALSE) stopifnot(all(D300 <= D150 + 1e-10)) ``` Only the requested endpoint distances are returned, but paths may use every observation in the supplied matrix. This is an exact restriction of the queries, not a landmark restriction of the paths. The check illustrates that unrooted Fermat distances cannot increase when the available path vertices expand. These small teaching observations are separate from the frozen benchmark below. For the ordinary neighborhood baseline, construct a union-symmetric graph: an edge exists if either endpoint selects the other among its $k$ nearest neighbors. Its weights remain Euclidean lengths, without a Fermat power. ```{r small-graph,eval=FALSE} graph <- create.sknn.graph(X,k=10,neighbor.method='exact', connect.components=FALSE) edges <- graph.edges(graph) ig <- igraph::graph_from_data_frame(edges[,1:2],directed=FALSE, vertices=data.frame(name=seq_len(nrow(X)))) Dknn <- igraph::distances(ig,v=endpoints,to=endpoints,weights=edges$length) igraph::components(ig)$no ``` If the graph is disconnected, some entries can be infinite. This study records that limitation rather than adding repair edges. This is an **ordinary Euclidean kNN graph**, distinct from `create.fermat.sknn.graph()`, which selects neighbors using Fermat distances. ## The seven-surface benchmark The quadratic surfaces are $$ S_a^+(u,v)=(u,v,a(u^2+v^2)),\qquad S_a^-(u,v)=(u,v,a(u^2-v^2)),\qquad u^2+v^2\leq1, $$ with $a\in\{1,2,4\}$. Principal curvature magnitudes at the origin are $2a$. Increasing $a$ also changes the surface area and curvature away from the origin; it does not uniformly multiply curvature everywhere. Both patches have area $$ A(a)=\frac{\pi}{6a^2}\left[(1+4a^2)^{3/2}-1\right]. $$ ```{r surfaces,echo=FALSE,fig.alt='Paraboloid and saddle patches at three coefficients, with a common vertical range within each family.'} figure('surfaces') ``` **Figure 1.** The six quadratic patches. A common vertical range within each family makes the increased bending and height visible. The seventh surface is $(t\cos t,h,t\sin t)$, where $1.5\pi\leq t\leq4.5\pi$ and $0\leq h\leq20$. Its arc coordinate is $$ s(t)=\frac12\left[t\sqrt{1+t^2}+\operatorname{asinh}(t)\right]. $$ Uniform draws in the rectangle with coordinates $(s,h)$ give area-uniform samples. Its exact intrinsic distance is $g_{ij}=\sqrt{(s_i-s_j)^2+(h_i-h_j)^2}$. ### Interactive surface views An optional ivue view complements Figure 1 and also shows the Swiss roll. In interactive builds, choose a surface, drag to rotate, and scroll to zoom; **Reset view** restores its initial camera. Orange points are the 64 fixed endpoints, and blue points are 236 additional observations from the original first sampling repetition at $n=300$. These are saved area-uniform samples, not newly generated points. The translucent mesh is a drawing approximation of the target surface. These are target coordinates, not metric-MDS embeddings. Equal data-unit scales preserve their geometry, and common viewing bounds within each quadform family allow comparison of the three coefficients. The mesh resolution affects only the display; it is not used to estimate distances or reference geodesics. ```{r interactive-surfaces,echo=FALSE} have.ivue <- !identical(Sys.getenv('DGRAPHS_STATIC_VIGNETTES','true'),'true') && requireNamespace('ivue',quietly=TRUE) && requireNamespace('rgl',quietly=TRUE) && requireNamespace('htmlwidgets',quietly=TRUE) && requireNamespace('htmltools',quietly=TRUE) && all(c('plot3D.plain','layer3D.callback','layer3D.axes','camera.zup') %in% getNamespaceExports('ivue')) && 'limits' %in% names(formals(ivue::plot3D.plain)) if(have.ivue) { options(rgl.useNULL=TRUE) helper <- '../inst/doc-tools/fermat-surface-views.R' if(!file.exists(helper)) helper <- system.file('doc-tools','fermat-surface-views.R',package='dgraphs',mustWork=TRUE) source(helper,local=TRUE) surface_views <- readRDS(file.path(dirname(bundle),'surface-views.rds')) fermat_surface_gallery(surface_views) } else { cat('The installed guide uses the static surface views in Figure 1. Interactive views are available when rebuilding with DGRAPHS_STATIC_VIGNETTES=false and the optional visualization packages.') } ``` Interactive HTML needs a WebGL-capable browser, with no running R or Shiny session. To rebuild it interactively, install [ivue](https://github.com/pgajer/ivue) and its rgl backend and set `DGRAPHS_STATIC_VIGNETTES=false`. Installed guides use compact static illustrations by default; the website uses interactive views. The original figure files, benchmark tables and numerical results are unchanged. ### Sampling and distance comparisons Each surface has 64 fixed sampled endpoints. Of the 2,016 unordered endpoint pairs, 500 are assigned to scale calibration and 1,000 different pairs to assessment. Pairs may share endpoints. Five independent pools of additional observations give nested total sample sizes $n\in\{150,300,600,1200,2400\}$, always including the same 64 endpoints. All observations may support paths. Repetitions measure variation from the added path vertices, conditional on one endpoint set per surface. The power grid is $p\in\{1,1.25,1.5,1.75,2,3,4\}$. Ordinary kNN graphs use $k\in\{5,10,20,40\}$ and the prespecified rule $k=\lceil2\log n\rceil$. The logarithm in that rule is natural. Quadform references use `quadform_geodesics(method="best_of_six")` with 64 search rounds, 16 initial edges, 32 candidate draws and a 257-vertex limit. Paraboloid pairs also use `method="paraboloid_clairaut"`. This is currently an internal dgraphs research interface, accessed by the reproduction scripts; it is not presented as a stable exported API. Best-of-six returns finite-search candidate paths, not certified globally shortest geodesics. All **2,100 distance conditions** and **10,500 main reference pairs** completed. There were no failed main reference pairs. Four ordinary kNN conditions were disconnected, all at $k=5$; 238 evaluation-pair occurrences were unreachable across those conditions. These counts are retained in the saved evidence. ## Compare distance structure on a separate pair set For calibration pairs $C$, evaluation pairs $T$, estimates $D_q$, and reference lengths $g_q$, fit the scale and evaluate the normalized error: $$ \widehat b=\frac{\sum_{q\in C}D_qg_q}{\sum_{q\in C}D_q^2},\qquad E_T=\sqrt{\frac{\sum_{q\in T}(\widehat bD_q-g_q)^2} {\sum_{q\in T}g_q^2}}. $$ Zero means proportional agreement with the reference. Plots show $100E_T$ as a percentage. Raw errors with $\widehat b=1$ are saved too; they answer a different question because powered lengths have different units. Unreachable graph pairs are excluded from their finite-distance score, and coverage is reported. Scale fitting likewise uses only finite calibration pairs. Paired differences use ambient scores on the same reachable evaluation pairs. A small error with incomplete coverage does not demonstrate recovery of the entire surface. ```{r coverage,echo=FALSE} q <- subset(m,family=='euclidean_knn') coverage <- do.call(rbind,lapply(split(q,q$policy),function(x) data.frame(Rule=x$policy[1],Minimum.coverage=100*min(x$test_pairs/x$requested_test), Disconnected=sum(x$components>1)))) knitr::kable(coverage,row.names=FALSE,digits=1, col.names=c('$k$ rule','Minimum evaluation coverage (%)','Disconnected conditions')) ``` **Table 1.** Coverage over all surfaces, repetitions and sample sizes. “Adaptive” denotes the fixed logarithmic rule above. A disconnected component may contain no evaluated endpoints, so connectivity and pair coverage are reported separately. ## Intermediate powers and denser samples ```{r paraboloid,echo=FALSE,fig.alt='Fermat distance errors across sample sizes on three paraboloids; intermediate powers outperform ambient chords.'} figure('paraboloid_fermat') ``` **Figure 2.** Paraboloid errors against best-of-six references. Thin lines are individual repetitions; bold lines are means. Shading gives pointwise 95% Student-$t$ intervals with four degrees of freedom. The horizontal axis is logarithmic. All intervals here and below are conditional on the fixed endpoints and references, not simultaneous over the tested grid, and omit solver uncertainty. At $a=1$, $p=1.25$ already improves on ambient chords at $n=150$: mean error is 3.73% versus 8.51%. For $p=2$, error falls from 11.72% to 5.67% as the sample increases to 2,400. Thus both the choice of power and the availability of intermediate observations matter. On these three paraboloids, $p=1.25$ is the best tested complete power throughout the sample-size grid. ```{r saddle,echo=FALSE,fig.alt='Fermat distance errors on three saddles, showing the benefit of a larger intermediate power on the strongest saddle.'} figure('saddle_fermat') ``` **Figure 3.** The same comparison on saddles. On the strongest saddle at $n=2400$, $p=1.5$ has 2.55% error versus 28.32% for ambient and 11.99% for $p=1.25$. A power useful on a mild patch need not be sufficient on a more strongly bent patch. ```{r roll,echo=FALSE,fig.alt='Swiss-roll errors across sample sizes: larger Fermat powers suppress shortcuts more effectively than powers close to one.'} figure('swiss_roll_fermat') ``` **Figure 4.** Swiss-roll errors against exact intrinsic lengths. Ambient error is 55.84%; at $n=2400$, $p=3$ achieves 7.70%, whereas $p=1.25$ still has 51.99%. The best tested power changes from $p=4$ at $n=150$ to $p=3$ at $n=2400$. ```{r winners,echo=FALSE} z <- subset(s,n==2400 & family=='complete' & p>1) w <- do.call(rbind,lapply(split(z,z$id),function(x)x[which.min(x$mean),])) knitr::kable(data.frame(Surface=unname(surface_names[w$id]),Power=w$p, Fermat=100*w$mean,Ambient=100*w$ambient,Wins=w$wins),row.names=FALSE,digits=2, col.names=c('Surface','Best tested $p$','Fermat error (%)','Ambient error (%)','Wins out of 5')) ``` **Table 2.** Descriptive minima over the tested powers at $n=2400$. Errors are means over five repetitions; wins count paired repetitions with lower error than ambient. These powers were identified using the displayed evaluation scores, so this table is not an independent assessment of a tuned selection procedure. A separate calibration-only power search at $n=300$ and transfer to the other sample sizes is reported below, after the original fixed-grid comparisons. ## How much do ambient chords underestimate surface lengths? ```{r ratios,echo=FALSE,fig.alt='Geodesic-to-chord ratios on seven surfaces, with a logarithmic ratio axis and a long Swiss-roll tail.'} figure('chord_distortion') ``` **Figure 5.** Distributions of $R_q=g_q/c_q$ on the 1,000 evaluation pairs per surface, where $c_q$ is the ambient chord. One denotes equality. The axis has $\log_2$ spacing but is labeled in ratio units; two means twice the chord length. Violin widths show density on the log-ratio scale, with medians and quartiles inside. Shared endpoints make these descriptive pair distributions. The rolled surface exposes very large shortcuts. Quadform coefficients change the distribution too, but need not increase every pair's distortion: curvature varies spatially, and area-uniform sampling puts more mass in steeper regions. ## Ordinary kNN distances are a strong baseline ```{r knn,echo=FALSE,fig.alt='Euclidean symmetric kNN errors on paraboloids, saddles and the Swiss roll, with ambient baselines and conditional intervals.'} figure(c('paraboloid_knn','saddle_knn','swiss_roll_knn')) ``` **Figure 6.** Ordinary Euclidean kNN errors on paraboloids (top), saddles (middle) and the Swiss roll (bottom). Thin lines are repetitions; bold means and error bars use conditional 95% $t$ intervals. Intervals can extend below zero, although errors cannot. The dashed ambient baseline uses all evaluation pairs; consult Table 1 before comparing it with a restricted finite-pair score. At $n=2400$, $k=20$ is more accurate than every tested complete Fermat power on each of the seven examples. Its errors on the mild paraboloid, mild saddle and Swiss roll are 0.62%, 0.49% and 0.35%. Larger neighborhoods are not automatically better: on the roll, $k=40$ has 13.03% error, consistent with shortcut edges. Connectivity alone is not a test of geodesic accuracy. These findings concern ordinary geodesic recovery under uniform sampling. They do not settle whether density-sensitive distances are useful on nonuniform biological data, where the desired geometry itself may depend on density. ## Reference accuracy can change the apparent error floor ```{r clairaut,echo=FALSE,fig.alt='Paraboloid errors using Clairaut references, showing reduced apparent error on the strongest patch.'} figure('paraboloid_clairaut') ``` **Figure 7.** The paraboloid comparison using Clairaut lengths. Observations, pairs and estimated distances are identical to Figure 2; only the reference and its fitted calibration scale change. For $a=4$, the $p=1.25$ error at $n=2400$ changes from 4.82% against best-of-six to 2.05% against Clairaut. The $k=20$ error changes from 4.24% to 0.60%. Some apparent estimator error is therefore reference-method disagreement. The original best-of-six results remain part of the evidence. A separate saddle diagnostic increased search rounds from 64 to 256 on 50 regularly spaced pair indices per saddle, selected without inspecting errors. The largest relative length reduction over those 150 pairs was 0.137%. This supports stability for those pairs and budgets, but does not certify global minimality or remove reference uncertainty elsewhere. ## Paired improvements over ambient distances ```{r paired,echo=FALSE,fig.alt='Paired Fermat-minus-ambient errors at 2,400 observations, showing replicate differences and conditional intervals.'} figure('paired_difference') ``` **Figure 8.** Scaled error minus ambient error at $n=2400$, in percentage points. Negative values favor Fermat. Black points are individual repetitions; blue means and intervals show conditional 95% $t$ intervals. Paired comparisons hold the endpoint set and added sample pool fixed within each repetition. ## Selecting a power at 300 observations and transferring it This separate follow-up tunes each of the six quadforms at $n=300$, then keeps that power fixed at $n\in\{150,300,600,1200,2400\}$. It reuses the original area-uniform samples, fixed endpoints, five repetitions, and numerical reference lengths. No surface samples or reference geodesics were recomputed. The Swiss roll is outside this tuning experiment; its original results above are unchanged. For each surface, the 500 calibration pairs were randomly assigned to five fixed folds of 100 pairs. The same folds were used for every power, repetition, and reference method. For repetition $r$, power $p$, and validation fold $F_j$, the scale was fitted on the other 400 pairs: $$ \widehat c_{r,p,-j}= \frac{\sum_{q\in\mathcal C\setminus F_j}D_{r,p,q}g_q} {\sum_{q\in\mathcal C\setminus F_j}D_{r,p,q}^{2}}, \qquad E_{r,p}^{\mathrm{CV}}= \sqrt{\frac{\sum_{j=1}^{5}\sum_{q\in F_j} (\widehat c_{r,p,-j}D_{r,p,q}-g_q)^2} {\sum_{q\in\mathcal C}g_q^2}}. $$ Here $\mathcal C$ is the calibration-pair set, $D_{r,p,q}$ the unrooted complete-graph Fermat distance, and $g_q$ the reference geodesic length. Selection minimizes $\overline E_p^{\mathrm{CV}}=\frac15\sum_{r=1}^{5} E_{r,p}^{\mathrm{CV}}$. Pooling the out-of-fold squared residuals before taking the square root gives each pair its contribution to the same normalized error; we do not average five differently normalized fold scores. The coarse grid has spacing $0.025$ on $[1,1.5]$ for all paraboloids and the coefficient-one saddle, and on $[1,2]$ for the two stronger saddles. The best coarse power is refined over its $\pm0.025$ neighborhood in steps of $0.005$, clipped to the search interval. The reported choice is the **best power found at this resolution**, not a continuous optimum. The ambient baseline $p=1$ participated in selection. Exact ties would favor the smaller power; none of the nine selected powers lies on a search boundary. All nine choices were saved before evaluation: six using best-of-six references and three independent paraboloid selections using Clairaut references. The follow-up completed 301 surface/reference/power combinations, each with five repetitions and five folds (7,525 fold scores), and 1,275 transfer comparison records. Every transfer record uses all 500 calibration and 1,000 evaluation pairs. No distance condition failed or was excluded. The number of records includes baseline comparisons and scoring each paraboloid selection against both reference methods; it is not a count of distinct distance searches. ```{r power-full,echo=FALSE,fig.alt='Full cross-validation curves for all six quadforms, including ambient power one and the selected interior minima.'} figure('power_tuning_full') ``` **Figure 9.** Full tuning curves at $n=300$ using best-of-six references. Thin lines show the five repetitions; blue curves show their mean and pale bands pointwise conditional 95% $t$ intervals. Red diamonds and dotted lines mark the selected power. Panels use different horizontal ranges because the two stronger saddles were searched through $p=2$; vertical ranges also differ. All coarse and refined grid points are included. Smaller error is better. ```{r power-near,echo=FALSE,fig.alt='Enlarged neighborhoods of each tuning minimum, displaying relative excess error and a one-percent threshold.'} figure('power_tuning_near') ``` **Figure 10.** Neighborhoods within $0.030$ of the selected power. The horizontal coordinate is $p-p_{\mathrm{chosen}}$; the vertical coordinate is $100(\overline E_p^{\mathrm{CV}}/\overline E_{p_{\mathrm{chosen}}}^{\mathrm{CV}}-1)$, the percentage increase relative to the minimum mean error. The dashed line marks a 1% relative increase. Bands use paired differences between each candidate and the selected power within each repetition, with pointwise conditional 95% $t$ intervals divided by the observed minimum mean error. They are not confidence intervals for the optimum's location. Dots are evaluated powers and lines only connect them. ```{r power-choices,echo=FALSE} x<-subset(power_data$selection,reference=='best_of_six') z<-subset(power_data$transfer_summary,reference=='best_of_six' & method=='selected_best_of_six') at_n<-function(n)z$mean[match(paste(x$id,n),paste(z$id,z$n))] knitr::kable(data.frame(Surface=surface_names[x$id],`Selected $p$`=sprintf('%.3f',x$p), `Within 1% of minimum`=sprintf('%.3f–%.3f',x$near_min,x$near_max), `Test error, $n=300$ (%)`=round(100*at_n(300),2), `Test error, $n=2400$ (%)`=round(100*at_n(2400),2),check.names=FALSE),row.names=FALSE) ``` **Table 3.** Best powers using best-of-six references, a descriptive near-minimum range, and mean held-out errors after refitting scale on all calibration pairs. The range contains the smallest and largest **tested** powers with mean cross-validation error no more than 1% above its minimum; it is not a confidence interval and does not classify untested powers. The 1% threshold is relative to the minimum error, not one percentage point of distance error. The minima span two to eight qualifying grid points. Their breadth varies with the surface, so the third decimal place records search resolution rather than certainty about a uniquely preferred power. ```{r power-transfer,echo=FALSE,fig.alt='Frozen selected-power error across five sample sizes, with ambient, fixed powers, and ordinary Euclidean neighborhood-graph comparisons.'} figure('power_transfer') ``` **Figure 11.** Transfer of the power selected at $n=300$ to the five sample sizes, scored against best-of-six references. Scale is refitted on the 500 calibration pairs at each size and repetition; the 1,000 evaluation pairs are unchanged. Ambient chords, fixed powers $1.25$ and $2$, and ordinary Euclidean symmetric-kNN shortest paths with $k=20$ provide comparisons. Thin lines show repetitions; thicker lines and points show means, with pointwise conditional 95% $t$ intervals. The sample-size axis is logarithmic and panel vertical ranges differ. All displayed neighborhood graphs are connected. The selected powers improve on ambient chords on all six surfaces at every sample size, in all five repetitions. Their mean errors decrease as the sample size increases. At $n=2400$, errors range from 1.00% on the mild saddle to 4.78% on the steepest paraboloid against best-of-six references. This supports useful transfer, but does not mean that selection at $n=300$ remains optimal at every size. ```{r power-difference,echo=FALSE,fig.alt='Paired selected-power minus fixed power 1.25 errors, showing a reversal on the intermediate saddle as sample size increases.'} figure('power_transfer_difference') ``` **Figure 12.** Selected-power error minus fixed-$p=1.25$ error on the same evaluation pairs and sample pools, in percentage points. Negative values favor the selected power. Thin lines show the five paired differences; blue means and pointwise conditional 95% $t$ intervals summarize them. For example, on the intermediate saddle the selected $p=1.335$ has 3.95% error at $n=300$, versus 5.91% for $p=1.25$. At $n=2400$ the ordering reverses: 2.00% versus 1.58%. On the strongest saddle, transferred $p=1.680$ reaches 3.25%, whereas the original grid's $p=1.5$ reaches 2.55%. Selection at one sample size therefore identifies a useful finite-sample choice, not a universal power. Ordinary Euclidean $k=20$ graph distances remain more accurate at $n=2400$ than these selected complete-graph Fermat distances on all six examples. ### Clairaut sensitivity on the paraboloids ```{r power-clairaut-tuning,echo=FALSE,fig.alt='Paraboloid cross-validation curves against best-of-six and Clairaut references, with separate selected powers.'} figure('power_clairaut_tuning') ``` **Figure 13.** The same paraboloid distance estimates, calibration folds and repetitions scored against the two reference methods. Each method fits its own scales and selects its own power. Thin lines show repetitions, thicker curves their means, and symbols their selected powers. The strongest paraboloid shows the largest shift in selected power and apparent error floor. ```{r power-clairaut-table,echo=FALSE} x<-subset(power_data$selection,reference=='clairaut') a<-subset(power_data$selection,reference=='best_of_six') z<-subset(power_data$transfer_summary,reference=='clairaut' & n==2400) value<-function(method) {q<-z[z$method==method,];round(100*q$mean[match(x$id,q$id)],2)} knitr::kable(data.frame(Surface=surface_names[x$id], `Best-of-six $p$`=sprintf('%.3f',a$p[match(x$id,a$id)]), `Clairaut $p$`=sprintf('%.3f',x$p), `Best-of-six choice error (%)`=value('selected_best_of_six'), `Clairaut choice error (%)`=value('selected_clairaut'),check.names=FALSE),row.names=FALSE) ``` **Table 4.** Separate paraboloid choices. Both error columns use Clairaut reference lengths at $n=2400$, with scales refitted against Clairaut; only the selected power differs. This separates the effect of power selection from the effect of changing the scoring reference. ```{r power-clairaut-transfer,echo=FALSE,fig.alt='Sample-size transfer curves scored against Clairaut, comparing powers selected with each reference and the unchanged baselines.'} figure('power_clairaut_transfer') ``` **Figure 14.** Paraboloid transfer scored against Clairaut references. Red curves use Clairaut-selected powers, blue curves the powers selected with best-of-six. All methods refit scale against Clairaut at every size. Means, repetitions and pointwise conditional intervals have the same meaning as in Figure 11. Clairaut selects $p=1.130,1.190,1.210$ for the three paraboloids, compared with $1.135,1.200,1.240$ under best-of-six. On the steepest paraboloid, the Clairaut-selected power has 1.77% error at $n=2400$ against Clairaut; the best-of-six-selected power scores 1.98% against that same reference and 4.78% against best-of-six. Reference disagreement has a larger effect on the apparent error floor than this modest change of power. The split protects the 1,000 evaluation pairs from this power-selection rule, but those pairs have already been examined in the earlier exploratory study. They are not a newly collected confirmatory test set. Calibration and evaluation pairs also share endpoints, and all five repetitions share the same 64 endpoints on each surface. Intervals are conditional on these endpoints and the selected power; they do not include uncertainty from repeating selection, choosing new endpoints, or numerical reference error. The experiment transfers the power, not the scale, and does not establish a population-wide optimum. ## Reproduce the evidence without slowing article builds The installed snapshot can be inspected directly: ```{r read-snapshot,eval=FALSE} path <- system.file('extdata','fermat-density-curvature','benchmark.rds', package='dgraphs',mustWork=TRUE) benchmark <- readRDS(path) nrow(benchmark$metrics) # 2,100 condition-level records nrow(benchmark$references) # 10,500 main reference pairs benchmark$power_tuning$selection nrow(benchmark$power_tuning$transfer_metrics) # 1,275 follow-up records benchmark$provenance ``` The companion `surface-views.rds` contains only the seven saved 300-point coordinate sets used in the interactive display; `surface-views-provenance.R` records their source checksums and sampling repetition. The benchmark snapshot contains condition-level errors and coverage, summaries, reference comparisons, the saddle sensitivity diagnostic, the power-selection fold scores and transfer results, and provenance. All data are synthetic. The companion `provenance.R` records source-file checksums and the research revision from which the snapshot and figures were imported. The full research experiment is named `fermat_density_curvature_01_oct_2026`. To reproduce the expensive calculations, copy the installed scripts to a writable directory: ```{r copy-workflow,eval=FALSE} src <- system.file('examples','fermat-density-curvature',package='dgraphs', mustWork=TRUE) work <- file.path(tempdir(),'fermat-recovery-reproduction') dir.create(work,showWarnings=FALSE) file.copy(list.files(src,full.names=TRUE),work,overwrite=FALSE) ``` From that directory run the following commands. This is intentionally separate from rendering the article. It requires dgraphs, igraph and jsonlite; use one worker on platforms without forked multiprocessing. ```sh Rscript --vanilla prepare.R Rscript --vanilla references.R 1 Rscript --vanilla distances.R 1 Rscript --vanilla reference_sensitivity.R Rscript --vanilla evaluate.R Rscript --vanilla validate.R Rscript --vanilla power_tuning.R 1 Rscript --vanilla power_transfer.R 1 Rscript --vanilla power_validate.R ``` Reference searches can be expensive. Completed chunks are reused; use a new working directory when changing the design. The scripts use an internal numerical reference interface, so retain the recorded package versions and source revision when interpreting differences between reruns. The packaged snapshot remains the fixed evidence for this article. The follow-up tuning script refuses to overwrite a frozen selection. Its recorded checksum is checked again during evaluation. The additional `power_summarize.R` script regenerates follow-up figures using ggplot2; article builds read saved figures and do not repeat power selection. The practical conclusion is to compare useful intermediate powers, sample sizes, and an ordinary neighborhood-graph baseline before attributing a result to Fermat geometry. Keep distance-estimation error separate from embedding distortion, and show reference uncertainty when it can change the comparison.