topolow

• sanity_check_subsample(): Validates subsample suitability for cross-validation
• prune_sparse_matrix(): Prunes sparse dissimilarity matrices to a well-connected subset
• plot_mcmc_diagnostics(): Replacement for create_diagnostic_plots() (see Deprecations)
• plot_ll_improvement(), plot_performance_trace(): Diagnostics for the parameter search
• plot_euclidify_diagnostics(), create_diagnostic_report(): Diagnostics for Euclidify() output
• **Enhanced Functions**:
• initial_parameter_optimization(): Now accepts opt_subsample parameter
• run_adaptive_sampling(): Now accepts opt_subsample parameter
• adaptive_MC_sampling(): Now accepts opt_subsample parameter (internal)
• Euclidify(): Now accepts opt_subsample parameter

How Subsampling Works

When opt_subsample is specified:

1. Each parameter evaluation uses a random subsample of the specified size

2. Connectivity is automatically validated; disconnected subsamples are rejected

3. If connectivity fails, sample size needs to be increased

4. Different parameter evaluations use different subsamples for robustness

5. **Final embedding always uses the full dataset**

The opt_subsample parameter is optional (default: NULL = use full data).

Performance Benefits

• Substantially reduces the cost of parameter optimization on large datasets (>500 points), since each evaluation embeds a subsample rather than the full matrix
• Reduces memory usage proportional to subsample size
• Parameters found on subsamples generalize well to full data

Recommendations

• Datasets < 500 points: Use full data (opt_subsample = NULL)
• Datasets > 500 points: Recommended opt_subsample = 200-500
• Always ensure opt_subsample >= folds for reliable cross-validation

Epoch-Based Parameter Search

• initial_parameter_optimization() gained an epochs argument (default: 1). With more than one epoch the search runs as a simple evolutionary strategy: 1. **Initialization**: start from the user-provided parameter ranges. 2. **Epoch loop**, for each epoch: a. Generate num_samples points by Latin hypercube sampling within the current ranges. b. If opt_subsample is set, each evaluation uses its own random subsample. c. Evaluate the parameter sets by cross-validation, in parallel batches. d. **Range update** (after every epoch but the last): sort by NLL, keep the top 50%, and set the next epoch's range to 0.75 * min(survivors) through 1.25 * max(survivors). This lets the search drift toward and zoom in on promising regions. 3. **Finalization**: log-transform the results of the final epoch for direct use with adaptive sampling.

Other changes

• Package gridExtra is a required import now.

Bug Fixes

• initial_parameter_optimization() (and therefore Euclidify()) could fail on Windows with invalid connection whenever the parameter search ran in more than one batch. A cluster shutdown handler was registered per batch, and all of them referred to the last cluster, so it was stopped repeatedly while the earlier clusters leaked. Each batch's cluster is now stopped as soon as that batch finishes. The failure was reliably triggered by max_cores = 2, which is the configuration used on CRAN check machines.
• Conversion of matrices to numeric in R/adaptive_sampling.R is now handled by the package's extract_numeric_values() function.
• run_adaptive_sampling() no longer sleeps for a fixed 1.5 seconds on every call. Two Sys.sleep() calls guarded file visibility on networked filesystems (NFS, Lustre); they now poll for the expected state instead, so the wait is effectively free on a local filesystem while networked ones still get up to 5 seconds. This removed about a second from every Euclidify() run and has no effect on results.

Documentation

• Clarified that the *_range arguments of Euclidify() are a **starting range for the search, not a hard bound**. Initial sampling draws inside the range, but adaptive refinement samples from a kernel density estimate of the best parameter sets found so far, which can propose values outside it. So ndim_range = c(3, 6) does not guarantee a 3- to 6-dimensional result; only the sanity limits (for example 1 <= ndim <= 50) are enforced. This behaviour is unchanged from earlier versions and is now documented, in a new section of ?Euclidify. Use euclidean_embedding() with an explicit ndim when a hard cap is required.

Improvements

• Enhanced connectivity checking using igraph
• Better error messages for disconnected data
• Adaptive strategies for handling sparse data
• Comprehensive logging of subsampling operations
• New diagnostic plots including MCMC exploration and parameter fit traces

C++ Backend for Core Optimization

• **This is the first release of topolow containing compiled code.** The core optimization loop of euclidean_embedding() has been rewritten in C++ using Rcpp, which speeds up the embedding substantially, most visibly on larger matrices. The remaining R-level loops in that function were replaced with vectorized operations.
• **New parameter for euclidean_embedding()**:
• convergence_check_freq: How often, in iterations, to test for convergence (default: 3). Larger values reduce the bookkeeping overhead; smaller values stop more promptly.
• **Implementation details**:
• **COO format**: edge data is held as a coordinate list, which avoids the zero-dropping behaviour of sparse matrix classes
• **Shuffled sweep order**: the sweep order is permuted each iteration with a C++ Mersenne Twister (std::mt19937), which is what lets the configuration escape local optima
• **Immediate updates**: Gauss-Seidel style in-place position updates, as in the original R implementation
• **Single-pass error calculation**: the MAE over active constraints is accumulated in one pass over the edge list at each convergence check
• **Cache-friendly layout**: edge data is stored in contiguous arrays
• **Pre-computed factors**: degree-based normalization factors are computed once, before the iteration begins
• **Direct memory access**: the inner loop writes to the raw column-major position buffer
• **Return value enhancement**: the convergence field of the returned topolow object now includes:
• achieved: whether convergence was reached
• error: final MAE on active constraints
• final_k: final spring constant after cooling
• **Dependencies**: added Rcpp to LinkingTo (compile-time only; no new runtime dependency). The backend is Rcpp-only and links no BLAS or LAPACK routines, so no src/Makevars is required on any platform.

Reproducibility and numerical equivalence

• The algorithm, its force model, cooling schedule and convergence criterion are unchanged from 2.0.1. Results from this version are **statistically equivalent** to 2.0.1 but **not bit-identical**, for two reasons:
• The sweep order is randomized by a C++ generator seeded from std::random_device, independent of R's RNG. set.seed() fixes the starting configuration (drawn with stats::runif()) but not the sweep order, so repeated runs differ slightly. The R implementation in 2.0.1 shuffled the sweep order in the same way.
• The convergence MAE (result$convergence$error) is now accumulated in package code rather than by a BLAS routine, so it can differ in the last representable digit from earlier development builds. This is far below the run-to-run variation produced by the shuffle, and it makes the value consistent across platforms, which it previously was not. Embedded positions are unaffected.
• See the new Reproducibility section of ?euclidean_embedding for what is and is not seeded, and for how to obtain a single canonical map.

Backward Compatibility

• All existing code using euclidean_embedding() will work without modification
• Default parameter values preserve the original algorithm's behavior
• Output structure remains compatible with downstream functions (Euclidify(), parameter optimization, etc.)