Skip to content

Make concurrent one_off_* and ordination computes safe in one process - #88

Merged
sfiligoi merged 27 commits into
mainfrom
concurrency-safety
Aug 7, 2026
Merged

Make concurrent one_off_* and ordination computes safe in one process#88
sfiligoi merged 27 commits into
mainfrom
concurrency-safety

Conversation

@wasade

@wasade wasade commented Jul 31, 2026

Copy link
Copy Markdown
Member

libssu can be loaded into a host process — a server, a notebook, an embedding
extension — that wants several independent computes running at once. Today that
does not work, for two separate reasons.

It crashes. Two threads in one_off_matrix_inmem_fp32_v3, under ASan:
SEGV on 0x0 READ at unifrac_internal.cpp:74. register_report_status()
callocs the progress-flag array per compute and remove_report_status() frees
and NULLs it at the end, so whichever compute finishes first pulls the array out
from under the others. Serially the same binary is clean.

And where it doesn't crash, it returns wrong answers. Subsampling, PCoA and
PERMANOVA all draw from process-global std::mt19937 state — ours in
skbio_alt.cpp, and scikit-bio-binaries' own. Concurrent callers interleave
their draws. Measured on the 6-sample fixture: 4 threads each seeding 42, and
52 of 161 results disagreed with the serial seeded answer. This is why the
fix adds API surface rather than just documenting "don't do that".

What changed

  • report_status is a static std::atomic<bool>[]. No per-compute allocation,
    and remove_report_status() is deleted rather than refcounted — the state is
    process-wide and genuinely has no owner.
  • Per-call seeds: one_off_matrix_inmem{,_fp32}_v4, pcoa{,_fp32,_mixed}_seeded,
    compute_permanova_inmem_fp{64,32}_seeded. At seed >= 0 these build a
    generator local to the call and touch no shared state. Needs no
    scikit-bio-binaries change — those entry points already take a seed, we were
    hardcoding -1.
  • Every existing entry point keeps its behaviour, defined as its seeded form at
    seed = -1.
  • Two crashes found in the same code: n_substeps = 0 was a SIGFPE and
    n_substeps above the stripe count was a heap-buffer-overflow, both reachable
    from the in-memory entry points.
  • api.hpp gets an include guard.
  • New CI job building libssu_inmem.a and running a concurrency suite against
    it. Nothing in CI built that archive, and it is a distinct configuration —
    UNIFRAC_WASM semantics on a multithreaded target.

Judgement calls

  • No n_threads argument. nthreads-var is a per-task ICV, so
    omp_set_num_threads() on a calling thread already gives that thread its own
    team width without affecting others (verified against libgomp: two threads
    setting 2 and 8 each got exactly that). Adding a parameter would mean
    threading it through ~40 #pragma omp sites in unifrac_task_impl.hpp.
  • Ordination reproduces to a tolerance, not bit-exactly. A concurrent
    unifrac matrix is bit-identical to the serial one; a concurrent
    pcoa*_seeded is not, because the randomized SVD's parallel reductions are
    not order-pinned — 2.8e-16 measured, against ~0.5 for a genuinely different
    seed. The tests assert that bound and the docs say so rather than implying a
    bitwise-stable ordination.
  • pcoa*_seeded are EXTERN, the older pcoa* are not. The old names have
    C++ linkage and no combined/libssu.c wrapper, so they are missing from the
    dispatcher that gets installed. New names have no ABI history, so they are
    declared EXTERN and wrapped; changing the old ones would break anyone
    linking the mangled symbols.

Testing

test_su 2397, test_su_api 512, test_ska 435, test_api 403 assertions,
zero failures. capi_test, capi_inmem_test, ci/crawford_test.sh, and the
in-memory suite plain and under ASan. Each fix has a test that failed before it;
the n_substeps cases were rebuilt against the pre-fix commit to confirm they
still fault there.

Not covered

GPU/ACC builds — concurrent computes there share one device and one queue.
find_eigens_fast still passes seed = -1 internally; nothing calls it, and
its declared fp32 name (find_eigens_fast_p32) has never matched the exported
symbol (find_eigens_fast_fp32), which is a pre-existing bug worth its own fix.
Accelerator-detection caches race benignly on first call and a sanitizer will
flag them.

wasade and others added 12 commits July 30, 2026 15:40
Two or more one_off_matrix_inmem_fp32_v3 computes in flight in one
process crash today: su::process_stripes brackets every compute with
register_report_status()/remove_report_status(), and the latter frees
and NULLs the process-global report_status array while another compute
is still dereferencing report_status[tid] in su::try_report.

Confirmed under ASan and in test_su itself:

  Thread 7 received signal SIGSEGV
  #0 su::try_report(...)            src/unifrac_internal.cpp:74
  #1 unifracTT<...>                 src/unifrac_cmp.cpp:143
  #2 su::process_stripes(...)       src/unifrac.cpp:544

Reported by duckdb-miint, which today wraps every libssu call in a
process-wide mutex to avoid this, serializing a stage of their pipeline
that is otherwise embarrassingly parallel.

test_concurrent_matrix_inmem runs 4 threads x 25 computes against a
serial reference and requires bit-exact agreement (the unweighted
compute is deterministic). test_concurrent_faith_pd_inmem pins the same
contract for faith_pd_inmem, which does not register progress state and
is expected to already be safe.

Worker threads never call ASSERT -- the harness counters in
test_helper.hpp are not thread-safe -- so each records into its own
slot and the main thread checks every slot after joining.

Also add -lpthread to the test_su_api link line; test_su already had it.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
report_status becomes a statically allocated std::atomic<bool> array
instead of a per-compute calloc/free, so a compute can no longer free or
NULL state another concurrent compute is reading. This is the fix for the
SIGSEGV added as a failing test in the previous commit.

  static std::atomic<bool> report_status[CPU_SETSIZE];

Static storage is zero-initialized before any dynamic initialization and
that alone is enough for any thread to use the atomics, so nothing has to
be set up at runtime. Access is memory_order_relaxed: each flag is a
standalone notification that orders nothing else, so this compiles to the
same plain load/store as the bool array it replaces.

register_report_status() now only installs the SIGUSR1 handler, exactly
once however many computes come and go, and remove_report_status() is
deleted outright along with its only call site in process_stripes. It was
never part of the installed public surface (install_lib/install_inmem
export api.hpp, task_parameters.hpp and status_enum.hpp only), so there
is no API or ABI break.

Two incidental fixes fall out:

  - sig_handler no longer calls fprintf. That was reached whenever a
    SIGUSR1 arrived between computes, and fprintf is not
    async-signal-safe, so the old handler was itself undefined behavior.
  - try_report bounds task_p->tid against CPU_SETSIZE. tid is the task
    index, bounded by the stripe count rather than by any CPU count, so a
    caller asking for more substeps than CPU_SETSIZE used to index past
    the allocation. Those tasks now simply do not report progress.

pthread_mutex_init/destroy on printf_mutex are gone; it keeps its static
PTHREAD_MUTEX_INITIALIZER. Destroying a mutex another compute might be
locking was undefined behavior, and the mutex now has a real job, since
concurrent computes can reach sync_printf at the same time.

Behavior change worth noting for the ssu/faithpd tools: a SIGUSR1 sent
while no compute is running used to print "Cannot report status." Flags
are no longer torn down between computes, so it now sets them silently
and the next compute reports at its first checkpoint.

test_concurrent_matrix_inmem_reporting covers the signal path itself
under concurrency, raising SIGUSR1 before the workers start so the
coverage is deterministic rather than racing delivery against a compute
that may already have finished.

Verified: test_su, test_su_api, test_ska and test_api all pass (3542
assertions total); the previously segfaulting concurrent matrix compute
is ASan-clean over 801 computes on 4 threads, and TSan reports no races
for concurrent non-subsampled computes with bit-identical results. The
remaining TSan findings are the global-RNG races in skbio_alt.cpp, which
a later commit addresses.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
n_substeps says how many tasks to split a stripe range across and comes
straight from the caller. Two values were unsafe:

  0            set_tasks computes ceil(n_range / n_tasks), so this is a
               division by zero -- SIGFPE, reproduced as exit 136 in
               test_su.
  > n_stripes  set_tasks hands trailing tasks start == stop == n_stripes,
               and UnifracTaskVector's constructor reads
               dm_stripes[task_p->start] unconditionally
               (unifrac_task.hpp:60), one past the end. ASan:
               heap-buffer-overflow, reachable with a 6-sample table and
               n_substeps=4.

The overflow reads a pointer that an empty stripe range never
dereferences, which is why it needs a sanitizer to see and why it went
unnoticed.

Two of the three set_tasks call sites already clamped the second case
with duplicated code; one_off_matrix_T had no guard at all, and that is
the path behind every one_off_matrix* entry point. clamp_substeps now
holds the rule for all three, so the two that were already clamping also
pick up the zero guard.

The bound is the number of stripes the tasks will actually divide, which
is not the same as dm_stripes.size(): partial_v3 deliberately allocates
the full stripe vector while computing only the caller's sub-range, so
clamping it against the total left the same overflow reachable for any
sub-range that does not start at stripe 0 --

  ssu --mode partial --start 2 --stop 10 --n-substeps 9

on a 20-sample table, or the equivalent through the public partial() and
partial_v3() entry points. effective_stripe_stop() now owns the
"stripe_stop <= stripe_start means through the last stripe" rule that
set_tasks used to apply privately, so a caller can compute the same
range set_tasks will use rather than guessing.

Note for anyone reading the new tests: results are asserted bit-identical
across n_substeps values. That holds because the tree-traversal chunk
size is a compile-time constant and every accumulator is written by
exactly one task, so splitting the same stripes differently cannot change
the summation order for any output cell. It would stop holding if the
kernel's chunk size ever became a function of task size.

Still unvalidated and out of scope here: stripe_start is never checked
against the stripe count, so a caller passing one past the end underflows
the range arithmetic in set_tasks. That is orthogonal to n_substeps and
reachable even at n_substeps=1.

Verified: the whole test_su suite is ASan-clean (2360 assertions, exit 0),
including the concurrency suites; test_su, test_ska, test_api and
test_su_api all pass; ci/crawford_test.sh real-data regression passes.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
A subsampled compute draws from the process-global std::mt19937 in
skbio_alt.cpp. Concurrent callers therefore both race on it -- TSan
reports data races on mt19937::operator() (skbio_alt.cpp:240) and
mt19937::seed() (:35) -- and cannot get a reproducible result without
holding a process-wide lock across ssu_set_random_seed()-then-compute.
Measured: 4 threads x 40 concurrent seeded computes gave 15 races and 4
results that disagreed with the serial answer.

one_off_matrix_inmem_v4 and one_off_matrix_inmem_fp32_v4 add an int seed
next to the subsample arguments it governs:

  seed >= 0   deterministic draw from the seed alone, touching no shared
              state, so concurrent callers need no lock
  seed <  0   draw from the global RNG, which is what v3 always did

The rule and its wording are lifted from subsample_table_inmem_seeded,
which established this convention for the subsample entry point; v4 is
the same idea applied where the compute itself does the subsampling.
Same harness on the v4 path: 0 races, 0 mismatches.

v3 keeps its exact signature and behaviour and moves to api_compat.hpp as
a wrapper passing seed = -1, which is how v2 is already handled there.
That also lets combined/libssu.c dispatch only v4 and pick up v3 from
api_compat.hpp for free -- the same layering it already had for v2 over
v3. Without that wrapper the new entry points would have been unresolved
symbols for everyone consuming the dlopen-dispatched libssu.so rather
than libssu_inmem.a.

The subsampled table is now heap-allocated behind a unique_ptr because
the two branches differ in type; biom_inmem's destructor is virtual, so
delete through the base frees what the old stack local did.

Deliberately not covered: pcoa* and compute_permanova_inmem_* still pass
-1 to the dependency, so they use its global generator and remain unsafe
to call concurrently. Threading a seed to them is a separate change.

Note the bit-exactness the new tests assert relies on every call seeing
the same OpenMP width, since the subsample draw is distributed across the
team (biom_subsampled.cpp:249). That holds within a process because
nthreads-var is a per-thread ICV that nothing here changes, and OMP_DYNAMIC
defaults off. It is not a claim of reproducibility across thread counts.

Verified: test_su, test_ska, test_api, test_su_api (2388/426/403/503) and
test/capi_test, test/capi_inmem_test all pass; test_su is ASan-clean;
make inmem_static exports v3 and v4; combined/libssu.so exports v2, v3
and v4; ci/crawford_test.sh passes.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The honest answer to "which entry points may I call concurrently, and
what state is global?" could until now only be derived by reading
unifrac_internal.cpp and skbio_alt.cpp. README.md gains a "Calling the
library concurrently" section stating it, and api.hpp a short version
pointing at it.

Every claim was traced to the implementation. The ones worth surfacing
here because they are easy to get wrong:

  - There are two chained process-global generators, not one:
    ssu_set_random_seed reseeds this library's mt19937 and then consumes
    a draw from it to derive skbb's seed. Any other RNG-consuming call
    landing in between desynchronizes the sequence, which breaks
    reproducibility even single-threaded.
  - skbb_pcoa_fsvd_* with a non-negative seed touches no shared state, so
    it is the direct route for callers who need concurrent ordination
    while libssu's own pcoa* still hard-codes -1. skbb_permanova_* avoids
    the generator the same way but still consults a cached CPU-dispatch
    flag on x86_64, so it gets a weaker guarantee.
  - The accelerator detection caches (proc_use_acc, skbio_use_acc, and
    the dependency's equivalent) race on first call from every entry
    point documented as safe. Benign -- detection is a pure function of
    the environment, so racing writers store the same value -- but a
    sanitizer will flag it and concurrent info logging can interleave.
  - GPU/ACC builds are excluded, and not for lack of trying: the task
    kernels use bare `#pragma acc ... async` and `#pragma acc wait`,
    which share one default queue on one device, so a wait in one thread
    would block on another thread's kernels.

api.hpp also gets the include guard it never had. Without it a
translation unit that reaches the header twice redefines every typedef;
verified that a doubly-including C and C++ TU fails to compile before
this change and succeeds after. Downstream carries a wrapper header
solely to work around this, which can now go away.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
subsample_depth > 0 distributes the draw across the OpenMP team: one
generator per thread seeded in turn from the caller's seed, with
observations assigned to threads by the schedule. Both the number of
generators and the observation-to-generator mapping depend on team size,
so the same seed can produce a different subsampled matrix at a different
width.

Measured on the src/test.biom fixture at seed 42: widths 1, 2 and 4 agree,
width 8 differs, and each width is stable across repeats.

This does not affect concurrent callers in one process, which each get
their own team of the same size -- that is what the new tests rely on --
but it does affect any caller whose width varies per call, for instance
one taking it from a host thread-pool setting. Such a caller gets
reproducibility only for a fixed (seed, width) pair, which is worth
knowing before it turns into a confusing bug report.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
libssu_inmem.a is what downstream projects embed, and nothing in CI built
it. It is also its own configuration: UNIFRAC_WASM is defined, so no SIGUSR1
handler is installed and CPU_SETSIZE takes its fallback, yet unlike the WASM
build it is multi-threaded. The concurrency suite in test_su.cpp covers the
same entry points, but only as compiled for libssu.so.

Adds a native concurrency test that links the archive, plus inmem_test and
inmem_test_asan targets and a CI job that runs both. The job also checks the
archive exports the v3 and v4 in-memory entry points, since a dropped
translation unit would otherwise only surface downstream.

Confirmed the gate has teeth in this configuration: built at the commit
before the n_substeps clamp, the test's substep cases fault under ASan --
SIGFPE at api.cpp:477 for 0, heap-buffer-overflow at unifrac_task.hpp:60
for 4 -- both through one_off_matrix_inmem_fp32_v3.

skbb comes from conda-forge, headers and library from the same package, so
the job needs no second checkout; INMEM_SKBB_EXTERN became overridable to
allow that.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
su::pcoa, su::pcoa_inplace and su::permanova passed seed = -1 into
scikit-bio-binaries unconditionally, and -1 means "draw from skbb's
process-global generator". So pcoa, pcoa_fp32, pcoa_mixed and
compute_permanova_inmem_fp64/fp32 could not be called concurrently, and could
not be made reproducible by any caller: a concurrent PCoA would consume draws
another one was relying on.

This needs no change in the dependency. skbb_pcoa_fsvd_* and skbb_permanova_*
already take an int seed, and at seed >= 0 both build a generator local to the
call. The seed is now threaded through the su:: layer as a trailing parameter
defaulting to -1, so every existing call site keeps its behavior, and reaches
the new pcoa_seeded / pcoa_fp32_seeded / pcoa_mixed_seeded and
compute_permanova_inmem_fp64_seeded / _fp32_seeded entry points. The non-seeded
forms are now defined as their seeded form at -1, the same shape
subsample_table_inmem already had.

Concurrent ordination is reproducible but not bit-exact, unlike a concurrent
unifrac matrix: the randomized SVD accumulates through parallel reductions
whose order is not pinned, so a compute competing with others drifts by ULPs --
2.8e-16 measured across 4 threads, against ~0.5 for a changed seed. The tests
assert that tolerance rather than equality, and README.md and api.hpp say so
rather than promising bit-exactness.

find_eigens_fast keeps its hardcoded -1; README records that as deliberate.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Cleanup pass over the branch. Net -157 lines, no coverage change.

The one functional fix: pcoa_seeded / _fp32_seeded / _mixed_seeded are now
EXTERN with wrappers in combined/libssu.c. They had inherited the C++ linkage
of the older pcoa*, which has no wrapper -- so they were absent from the
dispatcher shim that `make all` installs, and README was listing an entry point
as concurrency-safe that most consumers could not call at all. The older three
keep their linkage; changing it would break anyone linking the mangled names.

Test scaffolding: test_su.cpp rebuilt the fixture structs eight times, carried
twin 35-line matrix workers differing only in the entry point, and repeated the
spawn/join block six times. The inmem suite already had the factored shape, so
that shape moved back: inmem_fixture::{make_table,make_tree,run_matrix} and
concurrency_fixture::run_workers. Two hand-rolled vector comparisons became
operator==, and the inmem suite now uses almost_equal from the header it
already included.

Comments: the seed-semantics paragraph appeared eleven times and the
ordination-tolerance explanation nine, with the same measurement quoted as
2.8e-16, ~3e-16 and ~1e-15 in three files. README stays normative and the
copies point at it. Every measured number, trap warning and non-obvious
constraint is kept.

Also: the ASan objects get their own suffix and archive, so the plain and
instrumented builds no longer invalidate each other and the clean-on-both-sides
dance is gone (17s -> 7s locally, and CI passes -j); the three su:: call sites
that take the global-RNG path now say /*seed*/ -1 like the C layer does; and
api_compat.hpp states the criterion for what belongs there, which has a real
consequence -- anything defined there needs no dlsym stub.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Every one of these was untracked and unignored, which is how a stray `git add
-A` could sweep them into a commit: the built binaries (none have a suffix, so
they need naming individually), the generated accelerator translation units
(the cpu ones were already listed, the acc ones were not), emscripten's default
a.out.js, and the outputs ci/crawford_test.sh writes into its own directory.
Same contract, 87 lines down to 59. Cuts the derivations and keeps the
conclusions: every measured number, every trap, and the per-symbol lists stay,
while the mechanism behind each one goes to the code comment that owns it.
test_concurrent_compute_permanova_inmem_fp64_seeded was the only red leg on
PR #88: (linux-gpu-cuda, ompgpu), in test_su_api, 1 of 512. Not a tolerance
problem and not ours -- it is a real defect in the dependency that our test is
the first thing to reach.

scikit-bio-binaries' permanova_perm_fp_sW_T sizes the device buffer for
permutted_sWs at n_perm and copies back n_perm elements, but launches the
kernel over n_perm+1 groupings against a host array of n_perm+1
(src/distance/permanova.cpp). So the last permutation's pseudo-F is never
written on the host, the device write runs one element past its allocation, and
permanova_T reads permutted_fstats[n_perm] while counting. Every GPU p-value
carries one count of heap garbage. Serially that garbage is stable and the
answer looks reproducible; with four callers churning the heap it is not, and
one count is 1/(n_perm+1) -- exactly the test's PVALUE_TOL.

Only test_su_api reaches it. test_su and test_ska link the CPU-only skbio_alt,
which forces skbb to the CPU, and the pre-existing test_permanova_inmem uses
999 permutations, so one count is 1/1000 and stays inside the same tolerance.

fstat is index 0 and is copied back correctly, so it stays asserted everywhere;
the p-value check is now gated on skbb actually running on the CPU, queried via
skbb_get_acc_mode(). test_su_api gains -lskbb for that query alone -- the code
under test is still reached only through libssu.

Also fixes an adjacent doc bug found while writing this up: README documented
UNIFRAC_SKBIO_USE_GPU for disabling GPU offload in the ordination path, and
nothing reads that variable. The one that works is skbb's own SKBB_USE_GPU.

Suites unchanged: test_su 2397, test_su_api 512, test_ska 435, test_api 403,
zero failures, plus inmem_test. The skip branch was exercised by preloading a
stub skbb_get_acc_mode() that reports an NVIDIA GPU.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
wasade and others added 6 commits July 31, 2026 12:51
scikit-bio/scikit-bio-binaries#15, filed with the full diagnosis. Comment and
README text only; no behaviour change. Suites unchanged at 2397 and 512.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Run 30656219910 was green on (linux-gpu-cuda, ompgpu), but not because the new
gate fired -- test_su_api reported "NVIDIA GPU not detected" that time and took
the CPU path, so the p-value was still asserted. The same binary on the same
runner had reported "NVIDIA GPU detected" on run 30652784446. Detection is
intermittent between processes within a single job, so a green run says nothing
by itself about which path the suite took.

Print the mode unconditionally, so the log records it either way.

Note this also corrects the reasoning in b6e00b3's message: test_su and test_ska
do have the ACC code compiled in -- their "NVIDIA GPU not detected" lines are
inside #if defined(UNIFRAC_ENABLE_ACC_NV) -- so what keeps them off the device
is failed detection in that process, not a CPU-only build.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Replaces the approach in b6e00b3..d2a2897, which was wrong twice over.

First, the link change had a side effect I did not anticipate: adding $(SKBBLIB)
to test_su_api stopped unifrac's NVIDIA detection from finding the device, so
the whole 512-assertion suite quietly stopped running on the GPU. That, not the
new gate, is why run 30656219910 was green -- the gate never fired, and the
p-value was still being asserted on the CPU path. Across three runs the
correlation is exact: test_api (no -lskbb) detects the GPU every time, test_su
(always -lskbb) never does, and test_su_api flipped from detecting to not
detecting exactly when the flag was added. Likely a link-order effect, since
skbb's loader uses dlmopen and two CUDA runtimes in separate namespaces is the
obvious suspect, but that part is unproven -- there is no GPU on the box where
this was diagnosed. Either way the flag is reverted, with a comment so nobody
adds it back.

Second, gating the assertion was the weaker fix. A p-value is a rank over
n_perm+1 values, and README already promises only that it holds or steps by
1/n_perm, so asserting equality asserted more than the library guarantees --
on CPU too, where ULP drift in the s_T reduction can flip a count near a tail
boundary. Allowing exactly one step keeps a real p-value check on the GPU,
where the gate would have checked nothing at all, and needs no new public API.

scikit-bio/scikit-bio-binaries#15 is filed with the full diagnosis and is
documented in README; unifrac's suite is not the right place to gate it.

Measured rather than assumed, on this fixture at ORD_SEED = 7: of 59 other
seeds, 49 move the p-value past the new bound and 10 do not, so one comparison
catches a leaked seed about five times in six and a leak has to survive a
hundred of them. Widening from 1e-2 costs 2 of those 10. fstat turns out not to
depend on the seed at all -- identical for all 59 -- so it is the corruption
check, not the seed check; the comment now says so.

Suites unchanged: test_su 2397, test_su_api 512, test_ska 435, test_api 403,
zero failures, plus inmem_test. Net diff against 124d7ab is now one tolerance
constant and two comments.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Review of 124d7ab..c249bce found five things worth fixing and two of my
statements to be false. Verified each against the source rather than taking
them on trust; all five held.

Blocking, now fixed:

- tests/inmem/test_concurrency_inmem.cpp was left on 1e-2 with the same
  fixture, seed and permutation count. Both bounds are now named constants
  derived the same way, so they cannot drift apart.
- The PVALUE_TOL comment cited "the permanova case in tests/inmem" for
  seed consumption. There is no such case -- the seed check there is pcoa only.
  Rather than just delete the reference, added the missing permanova one; the
  gap was real and the false citation is what surfaced it.
- It also cited test_pcoa_seeded without saying it lives in test_ska.cpp.
- README and test_ska.cpp both said a p-value steps by 1/n_perm. It is
  1/(n_perm+1) -- the unpermuted statistic is in the distribution. Anyone
  deriving a tolerance from the documented contract got 1.5/99.
- The comment justified the bound with two mechanisms and budgeted one count,
  which does not add up. Settled it by measuring instead of arguing: at a fixed
  OpenMP width, over 500 repeats, fstat takes two values ~9e-16 apart and the
  p-value takes exactly one. The CPU clause was speculation and is dead on this
  fixture, leaving skbb#15 -- provably at most one count -- as the only live
  mechanism. One step is the right budget; the comment now says why.

That measurement turned up something the review did not ask about, and it is
the more useful finding: a seeded PERMANOVA p-value depends on the caller's
OpenMP width. Not ULP drift -- skbb sizes its permutation chunk as
2*omp_get_max_threads()*16, so the width selects the chunking and therefore the
whole permutation set. Width 1 gives 0.55 on this fixture, every width >= 2
gives 0.49. Documented in README next to the identical caveat for a seeded
subsample. It does not affect the test, whose threads all share one width.

Also from the review:

- The 59/49/10 figures were prose that would rot silently the moment anyone
  touched the fixture or the seed. Replaced with an assertion over a spread of
  seeds, the pattern test_ska.cpp already uses, so the tolerance's
  discriminating power is checked on every run. Hedged the overclaim that a
  leaked seed cannot survive N_THREADS*N_ITERS trials -- true for a fresh
  stream per call, not for a shared generator returning one wrong answer.
- The Makefile comment stated an unproven dlmopen theory as fact. It now
  separates what was observed from what is suspected, records that test_su has
  the same problem and cannot drop the flag, and says why nobody should need
  skbb symbols in an API_ONLY binary.
- README's skbb#15 paragraph moved to a Known issues section under GPU
  support. It corrupts a single serial GPU PERMANOVA just as well, so burying
  it under "Calling the library concurrently" hid it from the people it
  affects. Fixed two overstatements there too: the last permutation is
  computed, just written past the end of an undersized buffer and never copied
  back, and one garbage count only moves the p-value when it compares
  differently. Dropped the changelog parenthetical; noted both GPU env vars
  latch on first compute.

Acted on one suggestion beyond the diff, because it is what made this hard to
diagnose: ssu_load_check() in ssu_ld.c discarded dlerror(), so "no such
library" and "library present but would not load" both surfaced as "NVIDIA GPU
not detected". It now prints the reason under UNIFRAC_CPU_INFO or
UNIFRAC_GPU_INFO. Three CI runs of correlation would have been one.

Declined: the review's scoped skbb_set_acc_mode() save/force/restore as an
alternative to widening. It needs the skbb symbol in test_su_api, which is the
link change this PR just reverted for taking the suite off the GPU.

test_su 2403, test_su_api 518 (+6 each, the new seed assertions), test_ska 435,
test_api 403, zero failures; inmem_test and inmem_test_asan both clean.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The dlerror() reporting added in 28455dc fires from a detection path that
re-probes on every call, so run 30717704778's GPU leg carried a few hundred
copies of "Could not load shared library libssu_acc_amd.so: libelf.so.1" --
in a log whose whole purpose is diagnosing this class of problem. Latch it per
variant per process.

The message itself is doing its job and stays: it is how we now know that
runner's AMD variant is present but unloadable, which previously showed up only
as "AMD GPU not detected".

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Unrelated to the rest of this branch; splittable if preferred. It is here
because it is what is currently keeping the PR red.

rapi_test began failing with "there is no package called 'Rcpp'" on run
30717704778 (linux-gpu-cuda) and again on 30718905872 (ubuntu-24.04-arm). Not a
runner or an R-version problem: in the first of those runs, ubuntu-latest/cpu
fetched Rcpp_1.1.2.tar.gz and built it against the same R 4.6.1, minutes apart
from the leg that could not.

The pinned mirror is the cause. Fetching http://lib.stat.cmu.edu/R/CRAN's
PACKAGES.gz returns nothing readable, and R renders an unreadable index as
"package is not available for this version of R", which sends you looking at
the R version instead of the transport. cloud.r-project.org serves Rcpp 1.1.2 --
the same version the passing leg built -- over https and is a CDN rather than a
single host, so it does not have this failure mode.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
@sfiligoi

sfiligoi commented Aug 3, 2026

Copy link
Copy Markdown
Collaborator

@wasade Since you are adding a new set of _v4 functions, do you mind adding another parameter to them:

  • device_id

with <0 meaning host, and >=0 meaning that input (tree and samples) and output data (DM) is in GPU memory.

It is fine if it >=0 remains unsupported in your PR... Just return an error code.
I can add full support later.

@wasade

wasade commented Aug 3, 2026

Copy link
Copy Markdown
Member Author

Certainly, will add soon

Requested by @sfiligoi on #88: since v4 is new API anyway, take a device_id now
so the signature does not have to change again when device-resident compute
lands. < 0 is host memory; >= 0 names an accelerator and returns an error for
the moment, to be implemented later.

Placed after seed and before mmap_dir, grouping it with the other data-placement
argument and keeping result last. Easy to move if you would rather it sat
elsewhere.

>= 0 returns a new unsupported_device rather than an existing code:
invalid_method is about the method string, and whoever implements this wants a
code they can distinguish and then stop returning. Appended to the enum, since
the numeric values are ABI; there are no exhaustive switches over ComputeStatus,
so nothing starts warning.

The check is the first thing both functions do, before any allocation, so a
caller who asks for a device cannot be handed a host-computed answer and
*result is left untouched. v3 delegates at device_id = -1, alongside the seed
= -1 it already passed.

Verified by disabling the guard and confirming the new assertions fail (8/526),
then restoring it. test_su 2411, test_su_api 526, test_ska 435, test_api 403,
zero failures; combined/ dispatcher compiles, inmem_test, capi_test,
capi_inmem_test and crawford_test.sh all pass.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
@wasade

wasade commented Aug 3, 2026

Copy link
Copy Markdown
Member Author

…added

Comment thread src/Makefile Outdated
Comment thread src/inmem_build.mk
# under concurrency, not allocation hygiene in whatever skbb is installed.
# The archive keeps its shipped -O3; the n_substeps overflow was confirmed to
# report at that level.
INMEM_ASAN_FLAGS ?= -fsanitize=address -g

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Does this work with all compilers?
At least the NVIDIA HPC SDK and AMD GPU-aware compilers? (in addition to vanilla gcc and clang)

Comment thread src/api.hpp Outdated
Comment thread combined/libssu.c
Comment thread src/api.cpp
}
#endif // UNIFRAC_WASM (file-based one_off_matrix wrappers)

compute_status one_off_matrix_inmem_v3(const support_biom_t *table_data, const support_bptree_t *tree_data,

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We need to have the implementation of both v3 and v4, for backwards compatibility.

Note: Just renaming the templated versions is OK, as those are internal. But we must preserve the external API interface.

wasade and others added 2 commits August 3, 2026 17:55
Four review comments from @sfiligoi, all of which held up when checked
against the code.

Keep dispatching v3 (libssu.c). This was a real break, not a style point.
combined/libssu.c dlopens a variant library versioned independently of itself,
so an entry point it resolves by dlsym is answered by whatever variant is
installed. Forwarding v3 to v4 inside the dispatcher instead made every v3 call
require a v4 symbol in the variant, and ssu_load() exits the process when a
symbol is missing -- so a new dispatcher over a pre-v4 variant would have killed
the caller. v3 and v4 are both dispatched again; api_compat.hpp's v3 forwarders
are guarded by UNIFRAC_COMPAT_SKIP_INMEM_V3 so src/libssu.so still gets them and
the two definitions do not collide.

The rule at the top of api_compat.hpp is what produced this, so it is rewritten:
it said only that anything defined there needs no dlsym stub, and omitted that
anything defined there is answered inside the dispatcher and therefore pins
every call to the newest variant ABI.

Drop EXTERN from pcoa*_seeded (api.hpp, libssu.c). Matches the non-seeded pcoa*
they extend, per the convention that non-EXTERN names are internal. Checked the
consumer first, since reachability is why they were EXTERN: duckdb-miint links
libssu_inmem.a / libssu_wasm.a and calls skbb's pcoa directly, so nothing needs
them through the dispatcher. compute_permanova_inmem_*_seeded stay EXTERN --
their non-seeded peers are. Consequence: test_concurrent_pcoa can no longer link
under API_ONLY, so it is guarded like the other internal-only cases in that file
and still runs in test_su and in tests/inmem.

Drop -lpthread from test_su_api (Makefile). The point of that binary is that a
consumer needs nothing but -l$(SSU); the flag was added on the assumption that
macOS/arm required it, which looks never to have been tested. It links clean
without it. Folded the recipe's two stacked comments into one.

Gate the ASan targets on gcc/clang (inmem_build.mk). Both -fsanitize=address and
the -print-file-name=libasan.so lookup are gcc/clang spellings, and NVHPC/AOMP
take neither. Nothing builds this target with those today -- the archive is
CPU-only and INMEM_CXX defaults to $(CXX) -- but it now fails with a clear
message rather than producing an uninstrumented binary. Verified with a stub
reporting itself as nvc++.

Not changed: the fifth comment, on api.cpp:797. The external v3 entry points are
unchanged and still exported -- diffing the EXTERN surface against main shows
additions only -- and the rename there was of the internal template, which the
comment explicitly allows.

test_su 2411, test_su_api 511 (pcoa suite moved out of the API_ONLY build),
test_ska 435, test_api 403, zero failures; dispatcher builds and dlsyms v3 and
v4, inmem_test, inmem_test_asan, capi_test, capi_inmem_test and crawford_test.sh
all pass.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Fallout from 40827bf, caught by build-and-test-inmem: dropping EXTERN from
pcoa_seeded / _fp32_seeded / _mixed_seeded gave them C++ linkage, so the archive
holds mangled names and the exact-match check reported MISSING SYMBOL for all
three. The symbols are present -- nm -C shows pcoa_seeded(double const*, ...).

Split the loop: the EXTERN "C" entries keep the exact-name match, the three
pcoa* are matched on the demangled signature. Still a real check that the TU
made it into the archive, just spelled for the linkage each one has.

Verified against the local archive, including that a name which is not there
does not match.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Comment thread src/api_compat.hpp Outdated
Comment thread combined/libssu.c Outdated
Comment thread combined/libssu.c
Comment thread src/api.hpp Outdated
Three more comments from @sfiligoi. The first two turn out to be one point, and
answering it properly shows my earlier triage of the api.cpp:797 comment was
wrong.

"Why do we need this??? We have the older (e.g. _v2) in this file without any
guardrails!" (api_compat.hpp) and "I don't think we want any conditional code
guardrails. All the _v3 functions should always compile." (libssu.c) are both
answered by moving one_off_matrix_inmem_v3 / _fp32_v3 out of api_compat.hpp and
back into api.cpp, which is where main had them.

The convention, stated plainly: api.cpp holds a real implementation of every
version the dispatcher resolves by dlsym; combined/libssu.c has an entry for
each; api_compat.hpp holds only the versions that are never dispatched, v2 and
older. Moving v3 into api_compat.hpp is what produced both earlier problems --
the dispatcher regression fixed in 40827bf, and then the guard I invented to
stop the duplicate definition. So the guard was a symptom of the wrong file, not
something that needed a mechanism. It is gone; the only #ifndef left in libssu.c
is the pre-existing BASIC_ONLY. api_compat.hpp is now byte-identical to main.

That also means the api.cpp:797 comment was not a misreading, as I had it in
triage: "we need to have the implementation of both v3 and v4" was saying where
v3 belongs, not asking whether the symbol still existed. v3 is now implemented
next to v4 in api.cpp, as a thin forwarder at seed = -1, device_id = -1 --
mirroring the file-based one_off_matrix_v3, which already forwards to the v4
template rather than duplicating its validation.

"Could you avoid making changes that are just eye candy": fair, and worse than
the diff showed. I had reflowed both v3 wrapper signatures onto two lines and
stripped a trailing space, putting formatting noise directly beside a real
semantic change. The v3 block in libssu.c is restored verbatim from main and the
v4 block is purely additive.

"Do we need this comment here? If you think it is important, it should be present
above, as it covers (the existing) pcoa, too": moved. The linkage note now sits
above the whole pcoa group, where it applies to all six names, and the seeded
block keeps only its one-line pointer to the seed semantics. Keeping it because
it is the thing that led to EXTERN being added in the first place.

test_su 2411, test_su_api 511, test_ska 435, test_api 403, zero failures. Both
libraries build with no duplicate-symbol error, which is itself the check that
the guard is unnecessary; the variant exports v3 and v4 and the dispatcher
resolves all four. inmem_test, inmem_test_asan, capi_test, capi_inmem_test and
crawford_test.sh pass.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Comment thread src/unifrac.cpp
Comment thread src/unifrac_internal.hpp
@sfiligoi

sfiligoi commented Aug 4, 2026

Copy link
Copy Markdown
Collaborator

@wasade
At this point, SIGUSR1 code handling looks like a major hack.
It does not belong to a multi-threaded library.
Do we still need it?
Should we just get rid of it, instead?
Or move it to the CLI ssu code?

Comment thread README.md Outdated
concurrently and they stop agreeing.

This is
[scikit-bio-binaries#15](https://github.com/scikit-bio/scikit-bio-binaries/issues/15).

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think this is the right place to apply a patch.
Let's just fix it as a separate PR.

wasade and others added 4 commits August 4, 2026 10:18
@sfiligoi: "I don't think this is the right place to apply a patch. Let's just
fix it as a separate PR."

The "Known issues" section documented scikit-bio-binaries#15 -- a device buffer
sized n_perm while the kernel runs over n_perm + 1 -- and recommended
SKBB_USE_GPU=N as a workaround. That belongs upstream, not in this README, and
not in this PR.

The concurrency contract keeps a one-clause factual note, because a caller
comparing p-values across concurrent computes needs to know one of the
1/(n_perm + 1) steps is already spent on a GPU. It links the upstream issue and
recommends nothing. src/test_su.cpp and the inmem suite already cite the issue
directly, so nothing dangles.

No code change: PVALUE_TOL is unchanged and still allows exactly the one rank
step README has always promised.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
@sfiligoi: "At this point, SIGUSR1 code handling looks like a major hack. It
does not belong to a multi-threaded library. Do we still need it? Should we just
get rid of it, instead? Or move it to the CLI ssu code?"

Removed rather than moved, per @wasade.

Installing a process-wide signal handler was a side effect of computing a
distance matrix: su::process_stripes called register_report_status(), so any
embedder -- duckdb-miint, for one -- had its SIGUSR1 disposition replaced
without asking. Nothing in the library restored it.

What went, in full:

  unifrac_internal.{hpp,cpp}  report_status[CPU_SETSIZE] and its static_assert,
                              sig_handler, register_report_status, try_report,
                              sync_printf and its pthread mutex, and the
                              CPU_SETSIZE fallback that existed only to size
                              that array
  unifrac.cpp                 the register_report_status() call
  unifrac_cmp.cpp             both su::try_report() call sites
  su.cpp                      ssu_sig_handler, the signal() call in main(), and
                              the usage() text documenting the feature
  test_su.cpp                 test_concurrent_matrix_inmem_reporting
  README.md                   the SIGUSR1 caveat in the concurrency contract

The CLI's own handler went too. It only ever printed "Status cannot be
reported." -- a placeholder for before the library took over -- so keeping it
would have advertised a feature that reports nothing, and usage() promised real
"tid:... k:... total:..." output. SIGUSR1 now takes its default disposition.
That is a user-visible removal: `kill -USR1 <ssu pid>` terminates ssu instead of
printing a status line.

Two notes on what this incidentally fixes and does not touch:

- On main, remove_report_status() restored SIG_DFL rather than the CLI's
  handler, so after the first compute a SIGUSR1 killed ssu. That is gone with
  the feature.
- src/unifrac.cpp's #include <signal.h> was already unused on main and is
  dropped as the last signal-related artifact in the compute path. Its
  neighbouring dead <stdarg.h>/<pthread.h> are left alone -- unrelated.

su::progressbar (unifrac.cpp) is untouched: also dead code, zero callers on main
and here, but a \r progress bar rather than anything signal-driven.

Verified: test_su 2397, test_su_api 497, test_ska 435, test_api 403, zero
failures; the counts drop by exactly the 14 assertions and one suite removed,
and no stray "tid:" lines remain. inmem_test and inmem_test_asan OK, dispatcher
builds, capi_test/capi_inmem_test and crawford_test.sh pass. Zero new warnings.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Adding unsupported_device to ComputeStatus widened the enum to 10 values but
left su.cpp's compute_status_messages at 9 entries, indices 0 through 8. Any
"Compute failed in ...: %s" report for that status would have indexed one past
the end of the array. Not reachable from ssu today -- only
one_off_matrix_inmem{,_fp32}_v4 return it, and the CLI has no device_id -- but
the table is indexed by whatever status the library hands back, so this was a
latent out-of-bounds read on an error path, which is the worst place for one.

Adds the missing message, drops the hard-coded bound, and asserts the count
against the enum so the next appended status cannot reintroduce the gap.
Verified the assert fires: removing the new entry gives "static assertion
failed: compute_status_messages needs one entry per ComputeStatus".

src/faithpd.cpp has its own 7-entry copy and is deliberately not touched here.
It has been short since grouping_missing was added, well before this branch, and
faith_pd_one_off returns nothing above output_error (index 6), so it is not
currently reachable either. Worth its own fix.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Zero callers. Verified against the whole working tree, not just tracked files --
the build generates unifrac_task_noclass_*.cpp / unifrac_accapi_*.cpp, which
git ls-files does not see -- and the only occurrence anywhere was the definition
itself. No header declares it, tracked or generated, so nothing could call it
without hand-declaring the symbol.

Not part of any API: global scope with C++ linkage (_Z11progressbarf), not
EXTERN, so it was never wrapped in ../combined/libssu.c and never appeared in
the dispatcher's exports -- only in src/libssu.so, and now in neither. Unrelated
to SIGUSR1; it wrote a \r bar to stdout and predates this branch.

Verified: test_su 2397, test_su_api 497, test_ska 435, test_api 403, zero
failures; inmem_test and inmem_test_asan OK, dispatcher builds,
capi_test/capi_inmem_test and crawford_test.sh pass. No new warnings.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Comment thread src/Makefile Outdated
@sfiligoi

sfiligoi commented Aug 7, 2026

Copy link
Copy Markdown
Collaborator

The GPU runners are indeed not detecting the GPU.
Need to diagnose if it is a CI-setup problem or a real problem.

Never mind.
Only test_su and test_ska are reporting NVIDIA GPU not detected and apparently this was always the case.

Need to investigate and document if this is indeed expected.
But unrelated to this PR.

@sfiligoi on src/Makefile:214: "I don't think this belongs in here. This is the
Makefile, not a build log."

Correct, and the same thing had crept into several other files, so this is a
sweep rather than a one-line fix. The rule was: keep what tells you why the code
is the way it is, drop what records how we found out.

Removed:

  src/Makefile        the whole GPU-detection incident writeup -- dated run
                      counts, which binary flipped when a flag was briefly
                      added, an unresolved dlmopen theory. The link invariant
                      and what an undefined reference means are kept, in three
                      lines instead of seventeen.
  main.yml            the Rcpp mirror outage narrative: date, the failing
                      index, and two CI run numbers. Kept why a CDN and not one
                      university mirror.
  ssu_ld.c            "takes several CI runs of correlation" and "added a few
                      hundred lines to a CI log" -- the second also carried a
                      figure corrected later. Kept both rationales.
  inmem_build.mk      "the n_substeps cases in particular corrupted the heap
                      silently before they were clamped" and "was confirmed to
                      report at that level".
  test_su.cpp         "over 500 repeats per width", "49 of 59 alternative seeds
                      clear the bound", "bit-identical for all 59", and the
                      framing of 1e-2 as the bound this branch replaced. The
                      conclusions those measurements support are kept; a
                      procedure nobody will rerun is not.
  inmem test          "nothing proved they link and run in it until this ran
                      here", and how the missing permanova seed check came to
                      light. Both describe writing the test, not the test.

Also made a few phrasings timeless rather than historical ("predates the
per-call seed" -> "has no per-call seed"; "what every caller did before the
parameter existed"), trimmed a speculative note about an entry point nobody has
asked for, and fixed test_ska.cpp saying a p-value is a rank out of n_perm where
the same file says n_perm+1 twenty lines earlier.

No code change. test_su 2397, test_su_api 497, test_ska 435, test_api 403, zero
failures; inmem_test and inmem_test_asan OK; workflow YAML still parses to the
same three jobs.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

@sfiligoi sfiligoi left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good now.

@sfiligoi
sfiligoi merged commit a9cea63 into main Aug 7, 2026
13 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants