← All skills

gaudi

analysis

Unsupervised multi-omics integration with GAUDI — two-stage UMAP, HDBSCAN clustering, and SHAP metagenes, in R or Python. Use when two or more omics layers share samples and you need sample clusters plus the features driving them.

GAUDI: multi-omics integration via UMAP + HDBSCAN

GAUDI (Group Aggregation via UMAP Data Integration) is an unsupervised, non-linear multi-omics integration method. It embeds each omics layer with UMAP independently, concatenates those embeddings, runs a second UMAP over the concatenation, clusters the result with HDBSCAN, and attributes features to the latent dimensions with SHAP values from XGBoost.

Reference implementation: the gaudi R package (GPL-3), Castellano-Escuder et al., Nat Commun 16:5771 (2025), doi:10.1038/s41467-025-60822-1.

Use it when you have two or more omics layers measured on the same samples and want sample clusters plus an interpretable account of which features drive them. It is a clustering and stratification method, not a supervised predictor and not an imputation method.

Everything documented below was verified against gaudi 0.1.18 — the Python route by execution, the R route by executing the install, load and availability blocks and by reading the 0.1.18 sources for the parameter defaults. The failure modes in Sharp edges are real and reproducible — read that section before running anything on real data. Several of them fail silently, and one of them silently deletes samples.

No R available, or Python-only? GAUDI is a composition of UMAP, HDBSCAN and SHAP — the R package is one implementation of the method, not the method itself. references/python-implementation.md states the algorithm language-agnostically and gives a tested Python version that reproduces the R package's clusters exactly on separable data. Do not install R just to run this. Everything in this document other than the R call syntax — the input contract, the parameter choices, the sharp edges, the reporting standards — applies to either route.

Why two stages of UMAP

The first UMAP normalises each layer's contribution: a 25,000-feature transcriptome and a 700-feature metabolome both become an n×4 embedding, so the integration is not dominated by whichever layer happens to be highest-dimensional. The second UMAP integrates the standardised representations. The paper validated this against nine alternative integration strategies (intersection, union, subtraction, joint matrix factorisation, UMAP on raw concatenated data) and selected concatenate-then-UMAP on cluster purity and silhouette.

HDBSCAN is used rather than k-means because the second UMAP's latent space is non-linear: HDBSCAN handles clusters of varying density and irregular shape, and does not require the cluster count in advance. The cost is a dedicated noise label (cluster 0) that you must handle deliberately.

Install and load

gaudi is not on CRAN or Bioconductor. If R is unavailable, skip this section entirely and use references/python-implementation.md.

# install.packages("devtools")
devtools::install_github("hirscheylab/gaudi")

Some functions need packages attached, not merely installed (see Sharp edges). The safe preamble:

library(gaudi)
library(dplyr)       # required by the survival helpers
library(patchwork)   # required by plot_gaudi_grid()

Check what is present before writing a script:

Rscript -e 'for (p in c("gaudi","uwot","dbscan","xgboost","SHAPforxgboost","dplyr","patchwork","sva","vegan","survminer","clusterProfiler","org.Hs.eg.db")) cat(sprintf("%-16s %s\n", p, ifelse(requireNamespace(p, quietly=TRUE), as.character(packageVersion(p)), "MISSING")))'

Input contract

A named list of matrices or data frames, samples in rows, features in columns, with row names set to sample IDs:

omics_list <- list(
  expression  = expr_mat,    # n samples x p1 genes
  methylation = meth_mat,    # n samples x p2 probes
  protein     = prot_mat     # n samples x p3 proteins
)

Name the list. Element names carry through to individual_factors and make on_omics indices legible — the metagenes slot is unnamed, so the list order is your only record of which layer is which.

Layers need not share features or dimensionality, but they must overlap in samples. GAUDI silently intersects on sample names and returns only the intersection, so a naming mismatch shows up as a quietly smaller result rather than an error. Check it yourself:

Reduce(intersect, lapply(omics_list, rownames))   # confirm before running

Pass samples_in_rows = FALSE if your matrices are features×samples. Row names are mandatory in either orientation — unnamed input fails with a confusing n_neighbors must be smaller than the dataset size.

Preprocess before GAUDI

GAUDI does no normalisation. Do it per layer, upstream:

  • Impute missing values. GAUDI cannot tolerate NA — see Sharp edges.
  • Scale/standardise per layer (z-score is the usual choice).
  • Apply the layer-appropriate transform first: log-CPM or VST for counts, M-values for methylation betas, log for intensity-based metabolomics.
  • Drop features that are constant or near-constant.

Running it

res <- gaudi(
  omics_list,
  umap_params      = list(n_neighbors = 15, n_components = 4),  # per-layer
  umap_params_conc = list(n_neighbors = 15, n_components = 2),  # integrated
  min_pts          = NULL,       # NULL -> max(floor(0.03 * n), 2)  <- SET THIS
  xgboost_params   = list(lambda = 0, eta = 0.5, gamma = 50,
                          max_depth = 10, subsample = 0.95),
  compute_features = TRUE,
  combine_omics    = FALSE,
  samples_in_rows  = TRUE,
  reassign_cluster_zero = FALSE,
  method           = "xgboost"   # or "rf"
)

metric = "euclidean" and min_dist = 0.01 are not arguments — they are uwot::umap defaults that GAUDI inherits. To change them, add them to umap_params.

Published defaults and their stated rationale:

parameter default why
n_neighbors 15 uwot default; balances local/global structure
n_components (stage 1) 4 enough for primary variance components per layer
n_components (stage 2) 2 visualisation and interpretation
min_dist 0.01 tight local packing, favours compact clusters
lambda 0 no regularisation — sparse biological signal
eta 0.5 moderate learning rate
gamma 50 sparse trees; few features retained
max_depth 10 admits feature interactions
subsample 0.95 stable feature selection across runs

Sizing guidance from the package: reduce n_neighbors below 100 samples, consider raising it above 1000. Scaling is roughly linear; 1000 samples × 25,000 features × 5 layers runs in under 10 minutes on standard hardware.

Result object

GAUDIObject, an S4 object with six slots:

slot contents
@factors data.frame, n × 3: UMAP1, UMAP2, clust; row names = sample IDs
@clusters numeric vector of cluster labels; 0 = HDBSCAN noise
@silhouette_score mean silhouette over the 2-D embedding; 0 if one cluster
@individual_factors named list of per-layer n×4 embeddings
@metagenes unnamed list, one data.frame per layer: contrib1, contrib2

metagenes[[i]] holds mean absolute SHAP contributions, so values are non-negative — magnitude only, no direction. Rows are sorted by |contrib1| descending. Because gamma = 50 forces sparse trees, most features are exactly 0; that is the intended behaviour, not a failure. To recover direction, compare the feature across clusters directly (e.g. limma).

res@silhouette_score
table(res@clusters)                       # inspect the size of cluster 0
head(res@metagenes[[1]], 20)              # top features for layer 1
nonzero <- sum(res@metagenes[[1]]$contrib1 > 0)

Consensus variant

c_gaudi() runs gaudi() n_max times, picks the highest-silhouette run as a reference, Procrustes-aligns the others onto it (vegan::procrustes), takes the per-sample median embedding, and clusters that. Use it when UMAP run-to-run variability matters — reporting a single stochastic embedding as the result is not defensible.

res_c <- c_gaudi(omics_list, n_max = 10, min_pts = 12)

Two differences to know: c_gaudi defaults to different XGBoost parameters (eta = 1, gamma = 100 — sparser than gaudi()), and it returns individual_factors empty. Pass xgboost_params explicitly if you want metagenes comparable to gaudi().

Sharp edges

These are verified defects and traps in gaudi 0.1.18, not hypotheticals.

compute_features = TRUE — the default — fails on xgboost R 3.x

The default path dies before it returns, on any current R toolchain:

Error in `colnames<-`(`*tmp*`, value = c(colnames(X_train), "(Intercept)")) :
  attempt to set 'colnames' on an object with less than two dimensions
Calls: gaudi -> xgboost_model -> <Anonymous> -> colnames<-

Nothing in gaudi changed. gaudi::xgboost_model calls the pre-3.0 R interface (xgboost(data =, label =, params =, verbose =)); xgboost 3.x renamed data to x and label to y and dropped params, and its deprecation shims let training complete with the wrong arguments rather than stopping. predict(mod, x, predcontrib = TRUE) then hands back a plain numeric vector instead of an n x (p+1) matrix, and SHAPforxgboost::shap.values fails setting column names on a one-dimensional object.

Everything downstream of the metagenes is unreachable while this holds: the metagenes slot, plot_metagenes(), plot_gaudi_grid() and gaudi_enrichment(). compute_features = FALSE runs clean, and gives you clusters and factors without the SHAP attribution.

This is the failure a source read cannot find. The gaudi sources are unchanged and self-consistent; what moved was a dependency underneath them. Check it against your own installed xgboost rather than trusting either the package or this page.

NA silently deletes samples, then crashes

align_omics() computes apply(x, 2, sd) > 0 on the features×samples matrix and keeps columns passing the test. A single NA anywhere in a sample's column makes its sd return NA, the test yields NA, and that sample is dropped. Because layers are filtered independently, they end up with different sample counts and the run dies in dplyr::bind_cols with:

Can't recycle `..1` (size 59) to match `..2` (size 60).

That message names neither NA nor the sample. Impute before calling GAUDI and assert completeness:

stopifnot(!any(sapply(omics_list, function(x) anyNA(x))))

For >20% missingness, drop the affected samples or features rather than imputing.

The zero-SD filter removes constant samples, not constant features

The same line intends to drop zero-variance features but, given the matrix orientation at that point, tests samples. Consequences: a constant/flatlined sample is dropped without warning, and genuinely constant features survive into the model. So filter zero-variance features yourself upstream, and verify the returned sample count:

stopifnot(nrow(res@factors) == length(Reduce(intersect, lapply(omics_list, rownames))))

The default min_pts fragments small datasets

min_pts defaults to max(floor(0.03 * n), 2). Below ~100 samples this collapses to 2 or 3, and HDBSCAN returns a shower of doubleton clusters. On a 60-sample, 3-group synthetic set with clean separation:

min_pts clusters found noise silhouette
2 (the default at n=60) 19 4 0.33
5 3 0 0.96
10 3 0 0.96
18 3 0 0.94

The truth was 3. Always set min_pts explicitly below ~150 samples, and sweep it rather than accepting the first answer:

for (mp in c(5, 8, 10, 15, 20)) {
  r <- gaudi(omics_list, min_pts = mp, compute_features = FALSE)
  cat(sprintf("min_pts=%2d  k=%2d  noise=%3d  sil=%.3f\n",
      mp, length(setdiff(unique(r@clusters), 0)), sum(r@clusters == 0),
      r@silhouette_score))
}

Report the sweep, not just the chosen value. A cluster count that swings wildly across min_pts is a finding about the data, and silhouette alone will not tell you the count is wrong — it is computed on the same embedding that produced the clusters, so it rewards the geometry UMAP already imposed. Anchor cluster counts to something external (survival, known subtypes, held-out annotation) before claiming they are biological.

gaudi_enrichment() only works on factor 1

Metagene rows are sorted by |contrib1|, and clusterProfiler::gseGO requires a decreasing-sorted input. So on_factor = 2 fails with geneList should be a decreasing sorted vector. Sort it yourself for dimension 2:

gl <- res@metagenes[[1]][["contrib2"]]
names(gl) <- rownames(res@metagenes[[1]])
gl <- sort(gl[gl != 0], decreasing = TRUE)
gse <- clusterProfiler::gseGO(gl, ont = "ALL", OrgDb = org.Hs.eg.db,
                              keyType = "SYMBOL", minGSSize = 2, maxGSSize = 200)

Also: gaudi_enrichment() hardcodes pvalueCutoff = 1 and pAdjustMethod = "none", so it returns everything unfiltered and uncorrected. You must correct and filter yourself:

e <- gaudi_enrichment(res, on_omics = 1, organism = org.Hs.eg.db)
e$padj <- p.adjust(e$pvalue, method = "BH")
e <- e[e$padj < 0.05, ]

Never report rows straight out of that function as significant. Note also that SHAP scores are all non-negative, so GSEA is running on a one-sided statistic — fgsea will warn about scoreType. Treat the result as ranked-magnitude enrichment, not directional enrichment.

Row names of the chosen layer must be gene symbols for this to work.

Survival helpers need data frames and dplyr attached

get_pairwise_survival_data() and plot_survival() are documented as accepting named vectors. They do not — the vector path hits an unquoted := and dies with could not find function ":=". Pass data frames with a time / censor column and sample row names. Both functions also call bare mutate() / across(), so library(dplyr) must be attached.

library(dplyr)
time_df   <- data.frame(time   = clin$os_days,   row.names = clin$sample_id)
censor_df <- data.frame(censor = clin$os_event,  row.names = clin$sample_id)

pairs <- get_pairwise_survival_data(res, time_df, censor_df)
head(pairs)   # clusters, pval, distance, quality

plot_survival(res, time_df, censor_df, pairs[["clusters"]][[1]])

quality is 0.8 * normalised(-log p) + 0.2 * scaled centroid distance — an uninterpretable triage heuristic for ranking which cluster pair to examine, not a statistic. The pval column is a raw pairwise log-rank p-value across all choose(k, 2) comparisons; correct it before reporting:

pairs$padj <- p.adjust(pairs$pval, method = "BH")

plot_survival() drops cluster 0 and any NA rows automatically.

plot_gaudi_grid() needs patchwork attached

It composes panels with | and /. Without library(patchwork) it fails with non-numeric argument to binary operator. plot_factors() and plot_metagenes() are unaffected.

combine_omics = TRUE shifts every metagene index

With combine_omics = TRUE, GAUDI ComBat-corrects the concatenated layers (treating layer identity as batch) and prepends that combined result to metagenes. For a two-layer input, length(metagenes) becomes 3: [[1]] = combined, [[2]] = layer 1, [[3]] = layer 2. Every on_omics index shifts by one. Verify by row count:

sapply(res@metagenes, nrow)

Cluster 0 is noise, and dropping it is a decision

drop_clusters(res, clusters = 0) removes noise samples; reassign_cluster_zero = TRUE instead assigns each to its nearest neighbour's cluster. Both change downstream n. State which you used and why — silently absorbing noise into clusters inflates apparent separation. Report the noise fraction either way.

Note that reassign_cluster_zero = TRUE reorders @factors rows (noise samples move to the end), so re-match by row name rather than position when joining metadata.

A defensible run

library(gaudi); library(dplyr); library(patchwork)
set.seed(42)   # UMAP and XGBoost are both stochastic

stopifnot(!any(sapply(omics_list, anyNA)))
common <- Reduce(intersect, lapply(omics_list, rownames))
message(length(common), " samples shared across ", length(omics_list), " layers")

res <- gaudi(omics_list, min_pts = 12, samples_in_rows = TRUE)

cat("k =", length(setdiff(unique(res@clusters), 0)),
    "| noise =", sum(res@clusters == 0),
    "| silhouette =", round(res@silhouette_score, 3), "\n")
stopifnot(nrow(res@factors) == length(common))

ggplot2::ggsave("results/gaudi_embedding.png", plot_factors(res, label_size = 4),
                width = 7, height = 6, dpi = 300)
ggplot2::ggsave("results/gaudi_grid.png", plot_gaudi_grid(res, top_features = 15),
                width = 14, height = 7, dpi = 300)

clusters <- tibble::rownames_to_column(res@factors, "sample_id")
readr::write_csv(clusters, "results/gaudi_clusters.csv")

Record for provenance — min_pts (and the sweep), the seed, the per-layer preprocessing, gaudi version (res@gaudiVersion), noise fraction, silhouette, and how cluster 0 was handled.

For a guarded end-to-end version of this — the same run with every failure mode in Sharp edges turned into a loud error, plus a min_pts sweep and a provenance file — see references/analysis-template.md. Copy the block into gaudi_analysis.R, fill in the CONFIG section, and run it with Rscript.

Reporting

State plainly:

  • number of layers, features per layer, and the final sample n after intersection
  • preprocessing and imputation per layer
  • min_pts and how it was chosen; the cluster-count sweep
  • clusters found, noise fraction, treatment of cluster 0
  • silhouette — and that it is computed on the same embedding it evaluates
  • random seed; whether gaudi() or c_gaudi() was used
  • for metagenes: that SHAP values are unsigned importance
  • multiple-testing correction for pairwise survival and enrichment

Do not claim clusters are biological subtypes on silhouette alone. Do not read UMAP inter-cluster distances as quantitative — UMAP preserves local neighbourhoods, not global geometry, so the layout is topology, not a metric space. Cluster assignments are the output; the coordinates are a picture of them.

Limits

  • Unsupervised clustering, not classification or prediction. No predict(); new samples require re-running the whole integration.
  • Stochastic. Set a seed; prefer c_gaudi() when stability matters.
  • Complete data required. No NA tolerance — impute upstream.
  • Requires matched samples across layers. For unmatched single-cell layers, the paper matches cells first (PCA + greedy matching, validated with Procrustes) before GAUDI; that matching is a separate analysis with its own assumptions.
  • Unsigned feature attribution. Magnitude only.
  • Silhouette is internal. It cannot validate that the cluster count is biologically real.

Try it

A self-contained check that this skill still works. No account, no key, no network beyond installing the packages.

Data — generated inline, which is why datasets: is empty. The method needs a known partition to be checked against, and no public multi-omics dataset supplies one: three groups of twenty samples across two layers, with the signal in the first 30 features of each, is a case where the right answer is three clusters and anything else is a failure. The generator is pinned here so the figures below can be reproduced exactly rather than approximately.

What this does and does not exercise. It runs the Python route from references/python-implementation.md, inlined so the block runs cold in an empty directory. The R package is not exercised — installing it needs an R toolchain — so the Sharp edges above remain source-confirmed rather than executed.

Run — pip install numpy pandas scikit-learn umap-learn hdbscan xgboost shap, then:

import inspect
import numpy as np, pandas as pd, umap, hdbscan
from sklearn.metrics import silhouette_score, adjusted_rand_score

def gaudi(omics, min_pts=None, n_components_layer=4, n_components_conc=2,
          n_neighbors=15, min_dist=0.01, compute_features=False, seed=42):
    """The implementation from references/python-implementation.md, metagenes off."""
    common = sorted(set.intersection(*[set(df.index) for df in omics.values()]))
    if len(common) < 10:
        raise ValueError(f"only {len(common)} samples shared across layers")
    omics = {k: v.loc[common] for k, v in omics.items()}
    for k, v in omics.items():
        if v.isna().any().any():
            raise ValueError(f"NA in layer '{k}' — impute upstream")
        omics[k] = v.loc[:, v.std(axis=0, ddof=1) > 0]
    emb = [umap.UMAP(n_neighbors=n_neighbors, n_components=n_components_layer,
                     min_dist=min_dist, metric="euclidean",
                     random_state=seed).fit_transform(v.values) for v in omics.values()]
    integ = umap.UMAP(n_neighbors=n_neighbors, n_components=n_components_conc,
                      min_dist=min_dist, metric="euclidean",
                      random_state=seed).fit_transform(np.hstack(emb))
    n = len(common)
    if min_pts is None:
        min_pts = int(np.floor(0.03 * n))
    min_pts = max(min_pts, 2)
    lab = hdbscan.HDBSCAN(min_cluster_size=min_pts,
                          min_samples=min_pts).fit_predict(integ) + 1
    sil = float(silhouette_score(integ, lab)) if len(set(lab)) > 1 else 0.0
    return dict(clusters=lab, silhouette_score=sil, n=n)

def layers(shift, n_per=20, seed=0):
    """Three groups; signal in the first 30 features of each layer."""
    rng = np.random.default_rng(seed)
    truth = np.repeat([0, 1, 2], n_per)
    ids = [f"S{i:02d}" for i in range(n_per * 3)]
    out = {}
    for name, p in (("expression", 300), ("methylation", 150)):
        X = rng.normal(size=(n_per * 3, p))
        for g in range(3):
            X[truth == g, :30] += shift * (g - 1)
        out[name] = pd.DataFrame(X, index=ids, columns=[f"{name[:4]}{j}" for j in range(p)])
    return out, truth

k = lambda lab: len(set(lab) - {0})

# --- INVARIANTS: a failure here means the skill is wrong ----------------------
clean, truth = layers(shift=3.0)
r = gaudi(clean, min_pts=5)
assert k(r["clusters"]) == 3, k(r["clusters"])          # separable groups are recovered
assert (r["clusters"] == 0).sum() == 0                  # nothing left as noise
assert adjusted_rand_score(truth, r["clusters"]) == 1.0 # exactly the true partition
assert r["clusters"].min() >= 0                         # R's convention: 0 = noise, never -1
print(f"clean (shift 3.0), min_pts=5 : k={k(r['clusters'])} noise=0 ARI=1.000 "
      f"sil={r['silhouette_score']:.3f}")

# umap-learn's min_dist default is 10x uwot's, which GAUDI inherits. The
# implementation passes 0.01 explicitly; if this ever fails, the porting note changed.
assert inspect.signature(umap.UMAP).parameters["min_dist"].default == 0.1
print("umap-learn min_dist default  : 0.1 (uwot's is 0.01 — pass it explicitly)")

# Error paths: the Python route rejects what the R package silently drops.
na = {n: d.copy() for n, d in clean.items()}
na["expression"].iloc[3, 7] = np.nan
try:
    gaudi(na, min_pts=5); raise AssertionError("NA should have raised")
except ValueError as e:
    print("NA in a layer                :", e)
try:
    gaudi({n: d.iloc[:5] for n, d in clean.items()}, min_pts=2)
    raise AssertionError("5 samples should have raised")
except ValueError as e:
    print("too few shared samples       :", e)

# --- OBSERVED, 2026-08-28: drift to investigate, not a bug -------------------
print("\nmin_pts x separation (truth k=3, n=60, default min_pts = 2):")
print(f"{'shift':>6} {'default':>18} {'min_pts=5':>18} {'min_pts=10':>18}")
for shift in (0.8, 3.0):
    om, tr = layers(shift=shift)
    cells = []
    for mp in (None, 5, 10):
        rr = gaudi(om, min_pts=mp)
        cells.append(f"k={k(rr['clusters']):2d} ARI={adjusted_rand_score(tr, rr['clusters']):.2f}")
    print(f"{shift:>6} {cells[0]:>18} {cells[1]:>18} {cells[2]:>18}")

# min_pts at or above the group size returns NOTHING, silently: no exception,
# every sample labelled noise, and silhouette reported as 0.0.
big = gaudi(clean, min_pts=20)
print(f"\nmin_pts=20 (group size 20)   : k={k(big['clusters'])} "
      f"noise={(big['clusters'] == 0).sum()} sil={big['silhouette_score']:.3f}")
if k(big["clusters"]) == 0:
    assert big["silhouette_score"] == 0.0   # an empty result, reported as a number

Expect — the assertions above are invariants: on separable data the method returns the true partition, noise is labelled 0 and never -1, and both error paths raise rather than returning a quietly smaller result. The numbers below are observed on the versions in verified: and dated — a mismatch is drift to investigate, not a bug. UMAP prints UserWarning: n_jobs value 1 overridden to 1 by setting random_state once per embedding; that is expected, and it is the seed doing its job.

clean (shift 3.0), min_pts=5 : k=3 noise=0 ARI=1.000 sil=0.941
umap-learn min_dist default  : 0.1 (uwot's is 0.01 — pass it explicitly)
NA in a layer                : NA in layer 'expression' — impute upstream
too few shared samples       : only 5 samples shared across layers

min_pts x separation (truth k=3, n=60, default min_pts = 2):
 shift            default          min_pts=5         min_pts=10
   0.8      k= 5 ARI=0.70      k= 3 ARI=0.75      k= 3 ARI=0.78
   3.0      k= 3 ARI=1.00      k= 3 ARI=1.00      k= 3 ARI=1.00

min_pts=20 (group size 20)   : k=0 noise=60 sil=0.000

Two things in that table are the point of running it. The default min_pts is not safe merely because a clean dataset survives it — at shift 3.0 it returns the truth, and at shift 0.8 the same default splits three groups into five. Separation decides, and you do not know your separation in advance, which is why the sweep is not optional. And min_pts=20, at the group size, returns every sample as noise with a silhouette of 0.0 and no error at all — an empty result that reads like a computed one. This collapse is specific to the Python route. The same sweep run in R against gaudi 0.1.18 returns k=3, noise=0 at every value from 5 to 20 on the same shape of data, so the exposure is shared but the size of the effect is not. Bound the sweep from above as well as below on the Python route, and read the noise count before the silhouette on either.

Reference files

  • references/api-reference.md — every exported signature, the full default set for gaudi() and c_gaudi(), the GAUDIObject slot table, and a behaviour matrix of what was observed empirically. Read it when you need an argument you cannot find above, or to confirm a slot's shape before indexing.
  • references/analysis-template.md — a complete guarded analysis script. Copy the fenced block to gaudi_analysis.R and edit the CONFIG section.
  • references/python-implementation.md — the method stated language-agnostically, plus a tested Python implementation for when R is unavailable. Includes the porting traps that make a naive translation silently wrong (min_dist defaults differ 10×; the noise label differs) and a verified account of what reproduces across implementations and what does not.

Sources