[med-svn] [Git][med-team/scikit-bio-binaries][upstream] New upstream version 1.1.0
Andreas Tille (@tille)
gitlab at salsa.debian.org
Mon Sep 7 03:39:28 BST 2026
Andreas Tille pushed to branch upstream at Debian Med / scikit-bio-binaries
Commits:
ccba9dc0 by Andreas Tille at 2026-09-07T04:33:41+02:00
New upstream version 1.1.0
- - - - -
28 changed files:
- .github/workflows/main.yml
- + .gitignore
- Makefile
- README.rst
- + api_tests/wasm/Makefile
- + scripts/fetch_eigen.sh
- + scripts/test_eigen_wasm.sh
- src/Makefile
- src/distance/permanova.cpp
- src/distance/permanova_dyn_impl.hpp
- + src/inmem_build.mk
- + src/ordination/linalg_backend.hpp
- + src/ordination/linalg_backend_eigen.cpp
- + src/ordination/linalg_backend_lapacke.cpp
- src/ordination/principal_coordinate_analysis.cpp
- + src/tests/wasm/generate_pcoa_expected.cpp
- + src/tests/wasm/generate_permanova_expected.cpp
- + src/tests/wasm/pcoa_inputs.hpp
- + src/tests/wasm/permanova_inputs.hpp
- + src/tests/wasm/test_center_wasm.cpp
- + src/tests/wasm/test_pcoa_wasm.cpp
- + src/tests/wasm/test_permanova_wasm.cpp
- + src/tests/wasm/test_smoke.cpp
- src/tools/skbb_generate_helper.py
- + src/util/portable_shuffle.hpp
- src/util/skbb_detect_acc.cpp
- src/util/skbb_dgb_info.hpp
- + src/wasm/emscripten_build.mk
Changes:
=====================================
.github/workflows/main.yml
=====================================
@@ -13,7 +13,10 @@ jobs:
build-and-test:
strategy:
matrix:
- os: [ubuntu-latest, macos-13, macos-latest, ubuntu-24.04-arm]
+ # macos-13 dropped: GitHub-hosted macos-13 runners were retired
+ # in early 2026; jobs queue for 24h and then time out. macos-latest
+ # is the current supported macOS runner and stays in the matrix.
+ os: [ubuntu-latest, macos-latest, ubuntu-24.04-arm]
runs-on: ${{ matrix.os }}
steps:
- uses: actions/checkout at v3
@@ -68,3 +71,69 @@ jobs:
# and a weird number to potentially catch potential bugs
export OMP_NUM_THREADS=3
make test
+
+ build-and-test-wasm:
+ # Emscripten WASM build: produces libskbb_wasm.a (single-threaded,
+ # CPU-only, Eigen-backed) and runs all WASM tests under node. Runs
+ # in parallel with the native matrix above.
+ #
+ # The WASM tests compare against expected values produced by the
+ # NATIVE build (see src/tests/wasm/generate_*_expected.cpp). The
+ # generated headers are not committed, so this job needs both the
+ # emsdk toolchain and a native gcc + cblas + lapacke to build the
+ # generators.
+ runs-on: ubuntu-latest
+ steps:
+ - uses: actions/checkout at v3
+ # Canonical setup-emsdk action (moved from mymindstorm/ to
+ # emscripten-core/ in v16). Pin major for reproducibility.
+ - uses: emscripten-core/setup-emsdk at v16
+ with:
+ version: 5.0.3
+ actions-cache-folder: emsdk-cache
+ - uses: actions/setup-node at v4
+ with:
+ node-version: '20'
+ - name: Install native build deps for expected-value generators
+ run: |
+ # Headers vs link symbols, on Ubuntu, are split across packages:
+ # - libopenblas-dev provides libopenblas.so (cblas + LAPACK
+ # Fortran symbols; USE_LAPACK=1 in the
+ # Ubuntu build). It does NOT provide
+ # cblas.h at the default include path.
+ # - libblas-dev provides cblas.h at
+ # /usr/include/x86_64-linux-gnu/cblas.h
+ # via Debian alternatives.
+ # - liblapacke-dev provides lapacke.h plus liblapacke.so
+ # (the LAPACKE C wrappers). The wrappers
+ # forward to whichever BLAS implementation
+ # `-lopenblas` or alternates is linked.
+ # We override BLASLIB to "-llapacke -lopenblas" on the test
+ # steps below — the default "-llapacke -lcblas" would need
+ # libatlas-base-dev for libcblas.so, which we don't install.
+ sudo apt-get update -qq
+ sudo apt-get install -y --no-install-recommends \
+ g++ make libopenblas-dev libblas-dev liblapacke-dev
+ - name: Cache Eigen headers
+ uses: actions/cache at v4
+ with:
+ path: .wasm-cache/eigen
+ key: eigen-3.4.0
+ - name: Fetch Eigen (cached)
+ run: scripts/fetch_eigen.sh
+ - name: Verify Eigen cache
+ run: scripts/test_eigen_wasm.sh
+ - name: Build libskbb_wasm.a
+ run: make wasm
+ - name: Run WASM unit tests
+ # libskbb_cpu.a and the expected-value generators are built as
+ # transitive prereqs of tests/wasm/expected/*.h. NOGPU=1 keeps the
+ # native portion CPU-only. BLASLIB substitutes the conda-flavoured
+ # default `-llapacke -lcblas` with `-llapacke -lopenblas`:
+ # liblapacke-dev provides the LAPACKE C wrappers, and
+ # libopenblas-dev provides the cblas + LAPACK Fortran symbols
+ # those wrappers call. Ubuntu does not ship a libcblas.so under
+ # that exact name (ATLAS does, but we don't install it).
+ run: env NOGPU=1 BLASLIB="-llapacke -lopenblas" make wasm_test
+ - name: Run WASM public-API tests
+ run: env NOGPU=1 BLASLIB="-llapacke -lopenblas" make wasm_api_test
=====================================
.gitignore
=====================================
@@ -0,0 +1,7 @@
+.wasm-cache/
+src/tools/__pycache__/
+
+# Generated by `make -C src wasm_test` from the native build's output;
+# never committed (see https://github.com/scikit-bio/scikit-bio-binaries/pull/12
+# discussion).
+src/tests/wasm/expected/
=====================================
Makefile
=====================================
@@ -1,10 +1,33 @@
-.PHONY: all api install test_bins test clean clean_install
+.PHONY: all api install test_bins test clean clean_install wasm wasm_clean
+
+# The native api/install paths and the wasm path both depend on generated
+# source files under src/ (permanova_cpu.cpp, skbb_accapi_cpu.cpp, produced
+# by python scripts). Each sub-make has its own DAG, so two concurrent
+# top-level targets (e.g. `make -j all wasm`) could fire the generators
+# in parallel and truncate the same output file. Force top-level targets
+# to serialize; sub-makes retain their internal parallelism.
+.NOTPARALLEL:
all:
$(MAKE) api
$(MAKE) install
$(MAKE) test_bins
+wasm:
+ scripts/fetch_eigen.sh
+ cd src && $(MAKE) wasm
+
+wasm_test:
+ scripts/fetch_eigen.sh
+ cd src && $(MAKE) wasm_test
+
+wasm_api_test: wasm
+ cd api_tests/wasm && $(MAKE) wasm_test
+
+wasm_clean:
+ cd src && $(MAKE) wasm_clean
+ cd api_tests/wasm && $(MAKE) wasm_clean
+
api:
cd src && $(MAKE) api
=====================================
README.rst
=====================================
@@ -11,7 +11,7 @@
Installation
------------
-You can install the latest release of scikit-bio-binaries using ``conda`` (`<https://www.anaconda.com/docs/getting-started/miniconda/main>`_)::
+You can install the latest release of scikit-bio-binaries using ``conda`` (`<https://github.com/conda-forge/miniforge>`_)::
conda install -c conda-forge scikit-bio-binaries
@@ -25,7 +25,7 @@ Local build instructions
The package can be build from source on your local machine.
We recommend using the conda-provided compilers and libraries, but system-installed ones should work as well.
-If you decide to create a dedicated build environment in ``conda`` (`<https://www.anaconda.com/docs/getting-started/miniconda/main>`_)::
+If you decide to create a dedicated build environment in ``conda`` (`<https://github.com/conda-forge/miniforge>`_)::
cplatform=`conda info |awk '/platform/{print $3}'`
if [[ "$(uname -s)" == "Linux" ]];
@@ -42,6 +42,69 @@ To test that the build succeeded, run::
make test
+WebAssembly (WASM) build
+~~~~~~~~~~~~~~~~~~~~~~~~
+
+scikit-bio-binaries can also be compiled to WebAssembly for use in
+browser-targeted or ``duckdb-wasm``-based projects. This build produces
+a static archive, ``libskbb_wasm.a``, that downstream emscripten
+projects can link into a final ``.wasm`` module.
+
+The WASM variant is single-threaded, CPU-only (no OpenMP, no GPU, no
+pthread) and uses `Eigen 3.4.0 <https://eigen.tuxfamily.org/>`_ as the
+linear-algebra backend instead of ``cblas``/``LAPACKE``. Eigen is
+header-only, so no external numerical library needs to be built under
+emscripten. The public ``skbb_*`` C symbol set is identical to the
+native build.
+
+Prerequisites:
+
+- An activated `emsdk <https://emscripten.org/docs/getting_started/downloads.html>`_
+ (emsdk 5.0.3 or newer recommended).
+- ``node`` 18+ available on ``PATH`` (only needed to run the WASM tests).
+
+Build::
+
+ scripts/fetch_eigen.sh # downloads Eigen 3.4.0 into .wasm-cache/eigen
+ make wasm # produces src/libskbb_wasm.a
+
+Run the WASM test suites::
+
+ make wasm_test # smoke, PERMANOVA, centering, PCoA
+ make wasm_api_test # public C API parity (api_tests/wasm)
+
+Expected tolerances (native LAPACK vs WASM Eigen, at the test matrices
+currently shipped):
+
+- ``mat_to_centered``: 1e-6 absolute per element.
+- PCoA / FSVD: 1e-6 absolute for eigenvalues and proportion explained;
+ 1e-3 absolute for sample coordinates (sign-adjusted per axis, since
+ eigenvectors are unique only up to sign). In practice the observed
+ drift is at machine epsilon (~1e-15), but the headroom is kept for
+ larger or more ill-conditioned inputs.
+- PERMANOVA: ``fstat`` is bit-identical native vs WASM at a fixed seed
+ (same arithmetic, same ``std::mt19937``, and the shuffle itself is a
+ portable Fisher-Yates in ``src/util/portable_shuffle.hpp`` rather
+ than ``std::shuffle``, so the permutation sequence is also
+ reproducible across toolchains). ``pvalue`` is therefore also
+ expected to be bit-identical at a fixed seed under single-thread
+ native and WASM. The WASM build is internally deterministic: the
+ same inputs and seed always produce the same outputs.
+
+Downstream linking example::
+
+ emcc my_code.c src/libskbb_wasm.a -sEXIT_RUNTIME=1 \\
+ -I src/extern -o my_code.js
+
+Limitations / out of scope:
+
+- No GPU offload (NVIDIA/AMD).
+- No multithreading. A future pthread-enabled variant would require
+ rebuilding with ``-pthread`` and serving the hosting page with the
+ appropriate COOP/COEP headers.
+- No ``SKBB_ENABLE_CPU_X86_LEVELS`` dispatch (wasm32 has its own
+ SIMD story; not wired up yet).
+
GPU support
~~~~~~~~~~~
=====================================
api_tests/wasm/Makefile
=====================================
@@ -0,0 +1,53 @@
+# WASM variant of the public-API tests.
+#
+# Compiles the existing api_tests/test_*.c files under emcc and runs them
+# under node, linking against the static libskbb_wasm.a produced by
+# `make -C src wasm`. This proves the skbb_* C surface behaves identically
+# to the native build, within documented tolerances.
+
+.PHONY: all wasm_test wasm_clean
+
+REPO_ROOT := $(abspath $(dir $(lastword $(MAKEFILE_LIST)))/../..)
+SRC_DIR := $(REPO_ROOT)/src
+WASM_LIB := $(SRC_DIR)/libskbb_wasm.a
+EIGEN_INC := $(REPO_ROOT)/.wasm-cache/eigen/include
+
+WASM_CC := emcc
+WASM_CFLAGS := -std=c99 -O2 -Wall \
+ -I$(REPO_ROOT)/api_tests/wasm/include
+WASM_LDFLAGS := -sEXIT_RUNTIME=1 -sALLOW_MEMORY_GROWTH=1 \
+ -sENVIRONMENT=node -sNODERAWFS=0
+
+# Stage the skbb public headers under scikit-bio-binaries/ so the .c
+# files' `#include "scikit-bio-binaries/..."` lines resolve under emcc.
+include/scikit-bio-binaries/util.h: $(SRC_DIR)/extern/util.h
+ @mkdir -p $(@D)
+ cp $< $@
+include/scikit-bio-binaries/distance.h: $(SRC_DIR)/extern/distance.h
+ @mkdir -p $(@D)
+ cp $< $@
+include/scikit-bio-binaries/ordination.h: $(SRC_DIR)/extern/ordination.h
+ @mkdir -p $(@D)
+ cp $< $@
+
+STAGED_HEADERS := include/scikit-bio-binaries/util.h \
+ include/scikit-bio-binaries/distance.h \
+ include/scikit-bio-binaries/ordination.h
+
+test_distance_wasm.js: ../test_distance.c $(STAGED_HEADERS) $(WASM_LIB)
+ $(WASM_CC) $(WASM_CFLAGS) $< $(WASM_LIB) $(WASM_LDFLAGS) -o $@
+
+test_ordination_wasm.js: ../test_ordination.c $(STAGED_HEADERS) $(WASM_LIB)
+ $(WASM_CC) $(WASM_CFLAGS) $< $(WASM_LIB) $(WASM_LDFLAGS) -o $@
+
+wasm_test: test_distance_wasm.js test_ordination_wasm.js
+ @echo "--- api/distance ---"
+ node test_distance_wasm.js
+ @echo "--- api/ordination ---"
+ node test_ordination_wasm.js
+
+all: wasm_test
+
+wasm_clean:
+ rm -f *.js *.wasm
+ rm -rf include
=====================================
scripts/fetch_eigen.sh
=====================================
@@ -0,0 +1,84 @@
+#!/usr/bin/env bash
+# Fetch Eigen headers for the WASM build of scikit-bio-binaries.
+#
+# Eigen is header-only: there is no compile step. We simply download a
+# pinned release tarball and extract the Eigen/ subtree into
+# .wasm-cache/eigen/include/. Idempotent: cache hit returns immediately.
+#
+# Pinned to Eigen 3.4.0 (stable release, ABI-compatible with skbb usage).
+#
+# Inputs (env, all optional):
+# EIGEN_VERSION release tag (default: 3.4.0)
+# FORCE_REFRESH=1 wipe the cache before fetching
+
+set -euo pipefail
+
+EIGEN_VERSION="${EIGEN_VERSION:-3.4.0}"
+EIGEN_SHA256="${EIGEN_SHA256:-8586084f71f9bde545ee7fa6d00288b264a2b7ac3607b974e54d13e7162c1c72}"
+EIGEN_URL="${EIGEN_URL:-https://gitlab.com/libeigen/eigen/-/archive/${EIGEN_VERSION}/eigen-${EIGEN_VERSION}.tar.gz}"
+
+REPO_ROOT="$(cd -- "$(dirname -- "${BASH_SOURCE[0]}")/.." && pwd)"
+CACHE="${REPO_ROOT}/.wasm-cache/eigen"
+INCLUDE_DIR="${CACHE}/include"
+
+log() { echo "[fetch_eigen] $*"; }
+
+if [[ "${FORCE_REFRESH:-0}" == "1" ]]; then
+ log "FORCE_REFRESH=1 — clearing cache"
+ rm -rf "${CACHE}"
+fi
+
+if [[ -f "${INCLUDE_DIR}/Eigen/Dense" ]]; then
+ log "cache hit: ${INCLUDE_DIR}/Eigen/Dense already present — skipping"
+ exit 0
+fi
+
+mkdir -p "${CACHE}"
+TARBALL="${CACHE}/eigen-${EIGEN_VERSION}.tar.gz"
+
+if [[ ! -f "${TARBALL}" ]]; then
+ log "downloading ${EIGEN_URL}"
+ if command -v curl >/dev/null 2>&1; then
+ curl -fL --retry 3 --retry-delay 2 -o "${TARBALL}" "${EIGEN_URL}"
+ elif command -v wget >/dev/null 2>&1; then
+ wget -q -O "${TARBALL}" "${EIGEN_URL}"
+ else
+ echo "ERROR: neither curl nor wget available" >&2
+ exit 1
+ fi
+fi
+
+# Verify checksum. Skip verification if EIGEN_SHA256 is explicitly empty.
+if [[ -n "${EIGEN_SHA256}" ]]; then
+ log "verifying checksum"
+ ACTUAL="$(sha256sum "${TARBALL}" | awk '{print $1}')"
+ if [[ "${ACTUAL}" != "${EIGEN_SHA256}" ]]; then
+ echo "ERROR: sha256 mismatch on ${TARBALL}" >&2
+ echo " expected: ${EIGEN_SHA256}" >&2
+ echo " actual: ${ACTUAL}" >&2
+ rm -f "${TARBALL}"
+ exit 1
+ fi
+fi
+
+log "extracting into ${INCLUDE_DIR}"
+mkdir -p "${INCLUDE_DIR}"
+TMPDIR_EXTRACT="$(mktemp -d)"
+trap 'rm -rf "${TMPDIR_EXTRACT}"' EXIT
+tar -xzf "${TARBALL}" -C "${TMPDIR_EXTRACT}"
+
+# Tarball layout: eigen-<version>/Eigen/..., eigen-<version>/unsupported/...
+EXTRACTED="${TMPDIR_EXTRACT}/eigen-${EIGEN_VERSION}"
+if [[ ! -d "${EXTRACTED}/Eigen" ]]; then
+ echo "ERROR: unexpected tarball layout under ${EXTRACTED}" >&2
+ ls -la "${TMPDIR_EXTRACT}" >&2
+ exit 1
+fi
+
+rm -rf "${INCLUDE_DIR}/Eigen" "${INCLUDE_DIR}/unsupported"
+cp -r "${EXTRACTED}/Eigen" "${INCLUDE_DIR}/Eigen"
+if [[ -d "${EXTRACTED}/unsupported" ]]; then
+ cp -r "${EXTRACTED}/unsupported" "${INCLUDE_DIR}/unsupported"
+fi
+
+log "done: ${INCLUDE_DIR}/Eigen"
=====================================
scripts/test_eigen_wasm.sh
=====================================
@@ -0,0 +1,17 @@
+#!/usr/bin/env bash
+# Smoke test for the Eigen fetch. Header-only, so all we check is that
+# fetch_eigen.sh populated the expected include tree.
+
+set -euo pipefail
+
+REPO_ROOT="$(cd -- "$(dirname -- "${BASH_SOURCE[0]}")/.." && pwd)"
+CACHE="${REPO_ROOT}/.wasm-cache/eigen"
+
+fail() { echo "FAIL: $1" >&2; exit 1; }
+
+[[ -f "${CACHE}/include/Eigen/Dense" ]] || fail "missing ${CACHE}/include/Eigen/Dense"
+[[ -f "${CACHE}/include/Eigen/QR" ]] || fail "missing ${CACHE}/include/Eigen/QR"
+[[ -f "${CACHE}/include/Eigen/SVD" ]] || fail "missing ${CACHE}/include/Eigen/SVD"
+[[ -f "${CACHE}/include/Eigen/Core" ]] || fail "missing ${CACHE}/include/Eigen/Core"
+
+echo "OK: Eigen cache looks valid at ${CACHE}"
=====================================
src/Makefile
=====================================
@@ -1,4 +1,4 @@
-.PHONY: all api install test_bins test clean clean_install
+.PHONY: all api install test_bins test clean clean_install wasm wasm_test wasm_clean
all:
$(MAKE) api
@@ -27,7 +27,11 @@ ifndef NOGPU
LDFLAGS += -ldl
endif
-BLASLIB=-llapacke -lcblas
+# Conditional assignment so callers can override via env or command line
+# (e.g. BLASLIB=-lopenblas on Ubuntu where libopenblas-dev does not ship
+# the libcblas.so / liblapacke.so individual link names). Without `?=`,
+# this assignment would shadow the env variable inside the Makefile.
+BLASLIB ?= -llapacke -lcblas
ifeq ($(PLATFORM),Darwin)
SO_LDDFLAGS = -L$(PREFIX)/lib -fopenmp -dynamiclib -install_name @rpath/libskbb.so
@@ -96,24 +100,59 @@ ifneq ($(SKBB_ENABLE_ACC_AMD),)
endif
endif
+# WebAssembly build (emsdk required; see scripts/fetch_eigen.sh).
+# Provides WASM_CXX / WASM_CXXFLAGS used by the skbb_cpu_tu macro below,
+# plus the libskbb_wasm.a archive rule and WASM-only objects/tests.
+include wasm/emscripten_build.mk
+
+# Native in-memory static-archive build. Eigen backend, OpenMP enabled,
+# no LAPACKE/cblas/GPU. Provides INMEM_CXX / INMEM_CXXFLAGS used by the
+# skbb_cpu_tu macro below, plus the libskbb_inmem.a archive rule.
+include inmem_build.mk
+
+wasm: libskbb_wasm.a
clean:
# all source files are in subdirectories
rm -f *.cpp *.hpp *.h *.cu
rm -f *.o *.so *.a *.exe
+##
+# Shared per-translation-unit recipe.
+#
+# A handful of CPU-only sources are compiled into BOTH the native libskbb
+# (via $(CXX) $(CXXFLAGS)) and the WASM libskbb_wasm (via $(WASM_CXX)
+# $(WASM_CXXFLAGS)). The two builds previously had parallel per-file rules
+# in src/Makefile and src/wasm/emscripten_build.mk; this macro emits both
+# from a single point of definition so adding a new TU only needs one line.
+#
+# Native-only rules (x86 dispatch, GPU dispatch) and WASM-only rules
+# (Eigen linalg backend) remain explicit further down — the macro is for
+# the symmetric case only.
+#
+# Args:
+# $(1) — object stem (e.g. `util_rand` => `util_rand.o`, `util_rand.wasm.o`)
+# $(2) — source path
+# $(3) — header dependencies (space separated)
+# $(4) — extra compile flags applied to BOTH builds (e.g. -DSKBB_ACC_NM=skbb_cpu)
+
+define skbb_cpu_tu
+$(1).o: $(2) $(3)
+ $$(CXX) $$(CXXFLAGS) $(4) -c $$< -o $$@
+$(1).wasm.o: $(2) $(3)
+ $$(WASM_CXX) $$(WASM_CXXFLAGS) $(4) -c $$< -o $$@
+$(1).inmem.o: $(2) $(3)
+ $$(INMEM_CXX) $$(INMEM_CXXFLAGS) $(4) -c $$< -o $$@
+endef
+
##
# Utility/helper modules
SKBB_OBJS := util_rand.o
-
-util_rand.o: util/rand.cpp util/rand.hpp
- $(CXX) $(CXXFLAGS) -c $< -o $@
+$(eval $(call skbb_cpu_tu,util_rand,util/rand.cpp,util/rand.hpp,))
SKBB_OBJS += skbb_detect_acc.o
-
-skbb_detect_acc.o: util/skbb_detect_acc.cpp util/skbb_detect_acc.hpp util/skbb_accapi.hpp
- $(CXX) $(CXXFLAGS) -c $< -o $@
+$(eval $(call skbb_cpu_tu,skbb_detect_acc,util/skbb_detect_acc.cpp,util/skbb_detect_acc.hpp util/skbb_accapi.hpp,))
##
# skbb_accapi dynamic code wrappers generation
@@ -128,8 +167,7 @@ SKBB_OBJS += skbb_accapi_cpu.o
skbb_accapi_cpu.cpp: util/skbb_accapi.hpp
./tools/generate_skbb_accapi.py cpu direct > $@
-skbb_accapi_cpu.o: skbb_accapi_cpu.cpp util/skbb_accapi.hpp util/skbb_accapi_impl.hpp
- $(CXX) $(CXXFLAGS) -DSKBB_ACC_NM=skbb_cpu -c $< -o $@
+$(eval $(call skbb_cpu_tu,skbb_accapi_cpu,skbb_accapi_cpu.cpp,util/skbb_accapi.hpp util/skbb_accapi_impl.hpp,-DSKBB_ACC_NM=skbb_cpu))
#
# True accelerated variants will use a separate compiler
@@ -176,10 +214,7 @@ endif
# Distance methods
SKBB_OBJS += dist_permanova.o
-
-dist_permanova.o: distance/permanova.cpp distance/permanova.hpp distance/permanova_dyn.hpp util/skbb_accapi.hpp util/skbb_detect_acc.hpp util/rand.hpp util/skbb_dgb_info.hpp
- $(CXX) $(CXXFLAGS) -c $< -o $@
-
+$(eval $(call skbb_cpu_tu,dist_permanova,distance/permanova.cpp,distance/permanova.hpp distance/permanova_dyn.hpp util/skbb_accapi.hpp util/skbb_detect_acc.hpp util/rand.hpp util/portable_shuffle.hpp util/skbb_dgb_info.hpp,))
#
# permanova dynamic code wrappers generation
@@ -191,8 +226,7 @@ SKBB_OBJS += permanova_cpu.o
permanova_cpu.cpp: distance/permanova_dyn.hpp
./tools/generate_permanova_dyn.py cpu direct > $@
-permanova_cpu.o: permanova_cpu.cpp distance/permanova_dyn.hpp distance/permanova_dyn_impl.hpp
- $(CXX) $(CXXFLAGS) -DSKBB_ACC_NM=skbb_cpu -c $< -o $@
+$(eval $(call skbb_cpu_tu,permanova_cpu,permanova_cpu.cpp,distance/permanova_dyn.hpp distance/permanova_dyn_impl.hpp,-DSKBB_ACC_NM=skbb_cpu))
ifeq ($(SKBB_ENABLE_CPU_X86V3),1)
SKBB_OBJS += permanova_cpu_x86_v3.o
@@ -289,8 +323,17 @@ endif
# Ordination methods
SKBB_OBJS += ord_pcoa.o
+$(eval $(call skbb_cpu_tu,ord_pcoa,ordination/principal_coordinate_analysis.cpp,ordination/principal_coordinate_analysis.hpp ordination/linalg_backend.hpp util/skbb_dgb_info.hpp,))
+
+# Linalg backend split: native uses cblas+LAPACKE, WASM uses Eigen.
+# These are different source files, so the symmetric skbb_cpu_tu macro
+# doesn't apply — keep two explicit rules. The WASM rule lives in
+# wasm/emscripten_build.mk to keep the Eigen include path local to that
+# file.
+
+SKBB_OBJS += ord_linalg_backend.o
-ord_pcoa.o: ordination/principal_coordinate_analysis.cpp ordination/principal_coordinate_analysis.hpp util/skbb_dgb_info.hpp
+ord_linalg_backend.o: ordination/linalg_backend_lapacke.cpp ordination/linalg_backend.hpp
$(CXX) $(CXXFLAGS) -c $< -o $@
##
@@ -298,14 +341,9 @@ ord_pcoa.o: ordination/principal_coordinate_analysis.cpp ordination/principal_co
SHBB_EXTERN_OBJS := skbb_extern_distance.o skbb_extern_ordination.o skbb_extern_util.o
SHBB_EXTERN_HS := distance.h ordination.h util.h
-skbb_extern_util.o: extern/skbb_util.cpp extern/util.h util/rand.hpp util/skbb_detect_acc.hpp
- $(CXX) $(CXXFLAGS) -c $< -o $@
-
-skbb_extern_distance.o: extern/skbb_distance.cpp extern/distance.h distance/permanova.hpp
- $(CXX) $(CXXFLAGS) -c $< -o $@
-
-skbb_extern_ordination.o: extern/skbb_ordination.cpp extern/ordination.h ordination/principal_coordinate_analysis.hpp
- $(CXX) $(CXXFLAGS) -c $< -o $@
+$(eval $(call skbb_cpu_tu,skbb_extern_util,extern/skbb_util.cpp,extern/util.h util/rand.hpp util/skbb_detect_acc.hpp,))
+$(eval $(call skbb_cpu_tu,skbb_extern_distance,extern/skbb_distance.cpp,extern/distance.h distance/permanova.hpp,))
+$(eval $(call skbb_cpu_tu,skbb_extern_ordination,extern/skbb_ordination.cpp,extern/ordination.h ordination/principal_coordinate_analysis.hpp,))
##
=====================================
src/distance/permanova.cpp
=====================================
@@ -48,8 +48,9 @@
#include "util/rand.hpp"
+#include "util/portable_shuffle.hpp"
-#include <stdlib.h>
+#include <stdlib.h>
#include <algorithm>
@@ -158,7 +159,12 @@ static inline void permanova_perm_fp_sW_T(const uint32_t n_dims,
if (p!=0) { // do not permute the first one
const uint32_t grouping_el = p-tp;
uint32_t *my_grouping = permutted_groupings + uint64_t(grouping_el)*uint64_t(n_dims);
- std::shuffle(my_grouping, my_grouping+n_dims, randomGenerators.get_random_generator(grouping_el));
+ // Portable Fisher-Yates: std::shuffle delegates to
+ // std::uniform_int_distribution, whose output is implementation
+ // defined across libstdc++/libc++ and breaks cross-toolchain
+ // reproducibility. See util/portable_shuffle.hpp for rationale.
+ skbb::portable_shuffle(my_grouping, n_dims,
+ randomGenerators.get_random_generator(grouping_el));
}
}
// now call the actual permanova
=====================================
src/distance/permanova_dyn_impl.hpp
=====================================
@@ -40,7 +40,9 @@
#elif !(defined(_OPENACC) || defined(OMPGPU))
+#if defined(_OPENMP)
#include <omp.h>
+#endif
#define SKBB_CPU Y
@@ -51,7 +53,22 @@ static inline int pmn_get_max_parallelism_T() {
// No good reason to do more than max threads
// (but use 2x to reduce thread spawning overhead)
// but we do use 16x blocking, so account for that, too
+#if defined(_OPENMP)
return 2*omp_get_max_threads()*16;
+#elif defined(SKBB_WASM)
+ // WASM build: single-threaded, no OpenMP. Use 2*1*16 = 32 so the
+ // RNG-chunking scheme matches what a native OMP_NUM_THREADS=1 run
+ // would produce, keeping PERMANOVA results reproducible across the
+ // two toolchains at a fixed seed.
+ return 32;
+#else
+ // Reject other combinations at compile time rather than silently
+ // choosing an arbitrary chunk size. A non-OpenMP non-WASM native CPU
+ // build was never supported (pre-WASM code required <omp.h>), and
+ // changing its PERMANOVA chunk size would silently change seeded
+ // pvalues. If you hit this, set -fopenmp or define SKBB_WASM.
+#error "pmn_get_max_parallelism_T: no chunk size for non-OpenMP non-WASM CPU builds"
+#endif
#elif defined(SKBB_CUDA)
int deviceID;
=====================================
src/inmem_build.mk
=====================================
@@ -0,0 +1,79 @@
+# Native in-memory static-archive build for scikit-bio-binaries.
+#
+# Activated by the top-level `inmem_static` target. Produces
+# libskbb_inmem.a — the same in-memory subset shipped as libskbb_wasm.a
+# (Eigen backend, no LAPACKE/cblas, no GPU, no CPU-arch dispatch),
+# but compiled with the host toolchain so OpenMP is enabled. Intended
+# for downstream native projects that want to embed the in-memory
+# scikit-bio-binaries surface without dragging in BLAS/LAPACK.
+#
+# Per-translation-unit rules for sources shared with the native and
+# WASM builds are defined in src/Makefile via the `skbb_cpu_tu` canned
+# recipe (one definition emits .o, .wasm.o, AND .inmem.o rules). This
+# file owns:
+# - the inmem compiler vars (used by skbb_cpu_tu)
+# - the inmem-only TU (Eigen linalg backend, native CXX)
+# - the libskbb_inmem.a archive rule
+# - install / clean targets
+#
+# Invocation: this Makefile fragment is included from src/Makefile and
+# must be run with src/ as the working directory.
+
+# Reuse the Eigen drop fetched by scripts/fetch_eigen.sh for the WASM
+# build — the headers don't care which compiler is consuming them.
+INMEM_REPO_ROOT := $(abspath $(dir $(lastword $(MAKEFILE_LIST)))/../)
+INMEM_EIGEN_INC := $(INMEM_REPO_ROOT)/.wasm-cache/eigen/include
+
+INMEM_CXX ?= $(CXX)
+INMEM_AR ?= ar
+INMEM_MPFLAG ?= -fopenmp
+INMEM_CXXFLAGS := -std=c++17 -O3 -Wall -fPIC -I. \
+ -I$(INMEM_EIGEN_INC) \
+ -DSKBB_BLAS_BACKEND_EIGEN=1 \
+ -DEIGEN_DONT_PARALLELIZE \
+ -DNOGPU=1 \
+ $(INMEM_MPFLAG) \
+ -Wno-unknown-pragmas
+
+# Object list for the inmem archive. Stems must match the .inmem.o
+# targets emitted by skbb_cpu_tu in src/Makefile, plus the Eigen-backed
+# linalg implementation defined below. Mirrors WASM_OBJS — same TU set.
+INMEM_OBJS := \
+ util_rand.inmem.o \
+ skbb_detect_acc.inmem.o \
+ skbb_accapi_cpu.inmem.o \
+ dist_permanova.inmem.o \
+ permanova_cpu.inmem.o \
+ ord_pcoa.inmem.o \
+ ord_linalg_backend_eigen.inmem.o \
+ skbb_extern_util.inmem.o \
+ skbb_extern_distance.inmem.o \
+ skbb_extern_ordination.inmem.o
+
+# Inmem-only TU: the Eigen-backed linalg backend has no symmetric-macro
+# equivalent because the native build uses linalg_backend_lapacke.cpp
+# under the same `ord_linalg_backend.o` stem. Same source as the WASM
+# rule, just compiled with $(INMEM_CXX) instead of em++.
+ord_linalg_backend_eigen.inmem.o: ordination/linalg_backend_eigen.cpp ordination/linalg_backend.hpp
+ $(INMEM_CXX) $(INMEM_CXXFLAGS) -c $< -o $@
+
+libskbb_inmem.a: $(INMEM_OBJS)
+ rm -f $@
+ $(INMEM_AR) rcs $@ $(INMEM_OBJS)
+
+inmem_static: libskbb_inmem.a
+
+# Install (archive + public headers under a stable prefix layout).
+# Guards against the empty-PREFIX footgun: the top-level Makefile falls
+# back to CONDA_PREFIX, but if both are unset `mkdir -p /lib` would
+# silently target the root filesystem.
+install_inmem: libskbb_inmem.a
+ @test -n "$(PREFIX)" || { echo "ERROR: PREFIX is unset (and CONDA_PREFIX is unset). Pass PREFIX=/path or activate a conda env."; exit 1; }
+ mkdir -p "${PREFIX}/lib" "${PREFIX}/include/scikit-bio-binaries"
+ rm -f "${PREFIX}/lib/libskbb_inmem.a"; cp libskbb_inmem.a "${PREFIX}/lib/"
+ for f in $(SHBB_EXTERN_HS); do rm -f "${PREFIX}/include/scikit-bio-binaries/$${f}"; cp "extern/$${f}" "${PREFIX}/include/scikit-bio-binaries/"; done
+
+inmem_clean:
+ rm -f libskbb_inmem.a *.inmem.o
+
+.PHONY: inmem_static install_inmem inmem_clean
=====================================
src/ordination/linalg_backend.hpp
=====================================
@@ -0,0 +1,77 @@
+/*
+ * BSD 3-Clause License
+ *
+ * Copyright (c) 2025--, scikit-bio development team.
+ * All rights reserved.
+ *
+ * See LICENSE file for more details
+ */
+
+/*
+ * Narrow linear-algebra backend used by PCoA / FSVD.
+ *
+ * The native build links a BLAS/LAPACKE implementation (OpenBLAS / Netlib).
+ * The WebAssembly build ships without one and uses Eigen instead.
+ *
+ * Backend selection (compile-time):
+ * SKBB_BLAS_BACKEND_EIGEN=1 -> Eigen header-only
+ * otherwise -> cblas + LAPACKE
+ *
+ * The only entry points are the three primitives our PCoA/FSVD pipeline
+ * actually needs. We do not attempt a general BLAS facade.
+ */
+
+#ifndef SKBB_LINALG_BACKEND_HPP
+#define SKBB_LINALG_BACKEND_HPP
+
+#include <stdint.h>
+
+namespace skbb {
+namespace linalg {
+
+/*
+ * Column-major general matrix multiply: C = A * B.
+ * A is (m x k), B is (k x n), C is (m x n).
+ * Alpha and beta are fixed at 1.0 and 0.0 respectively (PCoA's only use).
+ * All strides equal the leading dimension (no sub-views).
+ */
+void gemm_nn(uint32_t m, uint32_t n, uint32_t k,
+ const double *A, const double *B, double *C);
+void gemm_nn(uint32_t m, uint32_t n, uint32_t k,
+ const float *A, const float *B, float *C);
+
+/*
+ * In-place QR: on input, H is (rows x cols); on output, H holds Q
+ * (rows x qcols) where qcols = min(rows, cols). Returns 0 on success.
+ *
+ * This replaces the {LAPACKE_dgeqrf + LAPACKE_dorgqr} pair used in
+ * principal_coordinate_analysis.cpp to build Q from H in-place.
+ */
+int qr_inplace(uint32_t rows, uint32_t cols, double *H, uint32_t &qcols);
+int qr_inplace(uint32_t rows, uint32_t cols, float *H, uint32_t &qcols);
+
+/*
+ * SVD "N, O" variant (matches LAPACKE_dgesvd flags 'N','O'):
+ * - jobu = 'N' : do not return U
+ * - jobvt = 'O' : overwrite T with V^T
+ *
+ * Layout contract:
+ * T is a column-major (rows x cols) buffer. On output, V^T occupies
+ * the first k = min(rows,cols) rows of T, i.e. for all (r,c) with
+ * 0 <= r < k and 0 <= c < cols, T[r + c*rows] = V^T[r,c]. The
+ * trailing rows in [k, rows) are left undefined by LAPACK when
+ * rows > cols (jobvt='O' only writes the k V^T rows). Consumers MUST
+ * use stride `rows` (the full buffer leading dimension); reading k
+ * packed rows as a (k x cols) contiguous submatrix is invalid and
+ * will observe garbage in the gaps.
+ *
+ * S (length >= min(rows,cols)) receives the singular values in
+ * descending order. Returns 0 on success.
+ */
+int svd_no(uint32_t rows, uint32_t cols, double *T, double *S);
+int svd_no(uint32_t rows, uint32_t cols, float *T, float *S);
+
+} // namespace linalg
+} // namespace skbb
+
+#endif
=====================================
src/ordination/linalg_backend_eigen.cpp
=====================================
@@ -0,0 +1,122 @@
+/*
+ * BSD 3-Clause License
+ *
+ * Copyright (c) 2025--, scikit-bio development team.
+ * All rights reserved.
+ *
+ * See LICENSE file for more details
+ */
+
+/*
+ * Header-only Eigen backend used by the WebAssembly build of
+ * scikit-bio-binaries, where a system BLAS/LAPACKE is not available.
+ *
+ * Each routine implements exactly the LAPACK-shaped contract described in
+ * linalg_backend.hpp: same dimensions, same in/out buffers, column-major
+ * layout. Eigen does the heavy lifting; we just own the storage shuffle.
+ */
+
+#include "linalg_backend.hpp"
+
+#include <algorithm>
+
+// Eigen is sensitive to assert triggers in optimized builds; keep them on
+// in debug builds only. The WASM build compiles with -O3 -fno-exceptions.
+#ifdef NDEBUG
+# define EIGEN_NO_DEBUG
+#endif
+// We never allocate temporaries inside real-time paths, but our PCoA
+// routine is called once per distance matrix; default policy is fine.
+
+#include <Eigen/Core>
+#include <Eigen/QR>
+#include <Eigen/SVD>
+
+namespace skbb {
+namespace linalg {
+
+template <class T>
+using ColMat = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>;
+template <class T>
+using CMapConst= Eigen::Map<const ColMat<T>>;
+template <class T>
+using CMap = Eigen::Map<ColMat<T>>;
+
+template <class T>
+static inline void gemm_nn_T(uint32_t m, uint32_t n, uint32_t k,
+ const T *A, const T *B, T *C) {
+ CMapConst<T> a(A, m, k);
+ CMapConst<T> b(B, k, n);
+ CMap<T> c(C, m, n);
+ c.noalias() = a * b;
+}
+
+void gemm_nn(uint32_t m, uint32_t n, uint32_t k,
+ const double *A, const double *B, double *C) {
+ gemm_nn_T<double>(m, n, k, A, B, C);
+}
+void gemm_nn(uint32_t m, uint32_t n, uint32_t k,
+ const float *A, const float *B, float *C) {
+ gemm_nn_T<float>(m, n, k, A, B, C);
+}
+
+template <class T>
+static inline int qr_inplace_T(uint32_t rows, uint32_t cols, T *H, uint32_t &qcols) {
+ qcols = std::min(rows, cols);
+
+ // HouseholderQR takes its input by value; we copy H into an owning
+ // matrix, factor, then extract the thin Q (rows x qcols) by applying
+ // householderQ() to an identity of the desired width.
+ ColMat<T> A = CMapConst<T>(H, rows, cols);
+ Eigen::HouseholderQR<ColMat<T>> qr(A);
+
+ ColMat<T> Q = ColMat<T>::Identity(rows, qcols);
+ Q = qr.householderQ() * Q;
+
+ CMap<T>(H, rows, qcols) = Q;
+ return 0;
+}
+
+int qr_inplace(uint32_t rows, uint32_t cols, double *H, uint32_t &qcols) {
+ return qr_inplace_T<double>(rows, cols, H, qcols);
+}
+int qr_inplace(uint32_t rows, uint32_t cols, float *H, uint32_t &qcols) {
+ return qr_inplace_T<float>(rows, cols, H, qcols);
+}
+
+template <class T>
+static inline int svd_no_T(uint32_t rows, uint32_t cols, T *T_mat, T *S) {
+ // Match LAPACKE_gesvd(jobu='N', jobvt='O'): compute V^T, no U.
+ // See linalg_backend.hpp for the output layout contract — we write
+ // V^T into the first k rows of T_mat with stride `rows`; consumers
+ // must read with the same stride.
+ ColMat<T> A = CMapConst<T>(T_mat, rows, cols);
+ Eigen::BDCSVD<ColMat<T>> svd(A, Eigen::ComputeThinV);
+
+ const uint32_t k = std::min(rows, cols);
+ for (uint32_t i = 0; i < k; ++i) {
+ S[i] = svd.singularValues()(i);
+ }
+
+ // Eigen's V is (cols x k); its transpose is (k x cols). Write
+ // Vt into rows [0, k) of the (rows x cols) T_mat buffer using
+ // the same stride `rows` the caller expects from LAPACKE.
+ ColMat<T> Vt = svd.matrixV().transpose();
+ CMap<T> dst(T_mat, rows, cols);
+ for (uint32_t c = 0; c < cols; ++c) {
+ for (uint32_t r = 0; r < k; ++r) {
+ dst(r, c) = Vt(r, c);
+ }
+ }
+ return 0;
+}
+
+int svd_no(uint32_t rows, uint32_t cols, double *T, double *S) {
+ return svd_no_T<double>(rows, cols, T, S);
+}
+int svd_no(uint32_t rows, uint32_t cols, float *T, float *S) {
+ return svd_no_T<float>(rows, cols, T, S);
+}
+
+} // namespace linalg
+} // namespace skbb
=====================================
src/ordination/linalg_backend_lapacke.cpp
=====================================
@@ -0,0 +1,81 @@
+/*
+ * BSD 3-Clause License
+ *
+ * Copyright (c) 2025--, scikit-bio development team.
+ * All rights reserved.
+ *
+ * See LICENSE file for more details
+ */
+
+/*
+ * Native linear-algebra backend: forwards to cblas + LAPACKE. Exact
+ * behaviour of the pre-refactor principal_coordinate_analysis.cpp.
+ */
+
+#include "linalg_backend.hpp"
+
+#include <algorithm>
+#include <cstdlib>
+#include <cblas.h>
+#include <lapacke.h>
+
+namespace skbb {
+namespace linalg {
+
+void gemm_nn(uint32_t m, uint32_t n, uint32_t k,
+ const double *A, const double *B, double *C) {
+ cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
+ m, n, k, 1.0, A, m, B, k, 0.0, C, m);
+}
+
+void gemm_nn(uint32_t m, uint32_t n, uint32_t k,
+ const float *A, const float *B, float *C) {
+ cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
+ m, n, k, 1.0f, A, m, B, k, 0.0f, C, m);
+}
+
+int qr_inplace(uint32_t rows, uint32_t cols, double *H, uint32_t &qcols) {
+ qcols = std::min(rows, cols);
+ double *tau = new double[qcols];
+ int rc = LAPACKE_dgeqrf(LAPACK_COL_MAJOR, rows, cols, H, rows, tau);
+ if (rc == 0) {
+ rc = LAPACKE_dorgqr(LAPACK_COL_MAJOR, rows, qcols, qcols, H, rows, tau);
+ }
+ delete[] tau;
+ return rc;
+}
+
+int qr_inplace(uint32_t rows, uint32_t cols, float *H, uint32_t &qcols) {
+ qcols = std::min(rows, cols);
+ float *tau = new float[qcols];
+ int rc = LAPACKE_sgeqrf(LAPACK_COL_MAJOR, rows, cols, H, rows, tau);
+ if (rc == 0) {
+ rc = LAPACKE_sorgqr(LAPACK_COL_MAJOR, rows, qcols, qcols, H, rows, tau);
+ }
+ delete[] tau;
+ return rc;
+}
+
+int svd_no(uint32_t rows, uint32_t cols, double *T, double *S) {
+ // LAPACKE contract: superb length >= min(rows,cols) - 1. Use
+ // min(rows,cols) so a 1x1 matrix (min==1) still has a non-zero
+ // allocation; LAPACK won't read past min-1 anyway.
+ const uint32_t superb_len = std::min(rows, cols);
+ double *superb = (double *) malloc(sizeof(double) * superb_len);
+ int rc = LAPACKE_dgesvd(LAPACK_COL_MAJOR, 'N', 'O', rows, cols,
+ T, rows, S, NULL, rows, NULL, cols, superb);
+ free(superb);
+ return rc;
+}
+
+int svd_no(uint32_t rows, uint32_t cols, float *T, float *S) {
+ const uint32_t superb_len = std::min(rows, cols);
+ float *superb = (float *) malloc(sizeof(float) * superb_len);
+ int rc = LAPACKE_sgesvd(LAPACK_COL_MAJOR, 'N', 'O', rows, cols,
+ T, rows, S, NULL, rows, NULL, cols, superb);
+ free(superb);
+ return rc;
+}
+
+} // namespace linalg
+} // namespace skbb
=====================================
src/ordination/principal_coordinate_analysis.cpp
=====================================
@@ -15,14 +15,14 @@
#include "principal_coordinate_analysis.hpp"
#include "util/skbb_dgb_info.hpp"
#include "util/rand.hpp"
+#include "linalg_backend.hpp"
-#include <stdlib.h>
+#include <stdlib.h>
+#ifdef _OPENMP
#include <omp.h>
+#endif
#include <algorithm>
-#include <cblas.h>
-#include <lapacke.h>
-
//
// ======================= PCoA ========================
//
@@ -195,18 +195,8 @@ void skbb::mat_to_centered(const uint32_t n_dims, const double mat[], float cen
// mat must be cols x rows
// other must be cols x rows (ColOrder... rows elements together)
template<class TReal>
-static inline void mat_dot_T(const TReal *mat, const TReal *other, const uint32_t rows, const uint32_t cols, TReal *out);
-
-template<>
-inline void mat_dot_T<double>(const double *mat, const double *other, const uint32_t rows, const uint32_t cols, double *out)
-{
- cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans, rows , cols, rows, 1.0, mat, rows, other, rows, 0.0, out, rows);
-}
-
-template<>
-inline void mat_dot_T<float>(const float *mat, const float *other, const uint32_t rows, const uint32_t cols, float *out)
-{
- cblas_sgemm(CblasColMajor,CblasNoTrans,CblasNoTrans, rows , cols, rows, 1.0, mat, rows, other, rows, 0.0, out, rows);
+static inline void mat_dot_T(const TReal *mat, const TReal *other, const uint32_t rows, const uint32_t cols, TReal *out) {
+ skbb::linalg::gemm_nn(rows, cols, rows, mat, other, out);
}
// Expects FORTRAN-style ColOrder
@@ -249,38 +239,12 @@ static inline void centered_randomize_T(const TReal centered[], const uint32_t n
free(tmp);
}
-// templated LAPACKE wrapper
-
// Compute QR
// H is in,overwritten by Q on out
// H is (r x c), Q is (r x qc), with rc<=c
template<class TReal>
-static inline int qr_i_T(const uint32_t rows, const uint32_t cols, TReal *H, uint32_t &qcols);
-
-template<>
-inline int qr_i_T<double>(const uint32_t rows, const uint32_t cols, double *H, uint32_t &qcols) {
- qcols= std::min(rows,cols);
- double *tau= new double[qcols];
- int rc = LAPACKE_dgeqrf(LAPACK_COL_MAJOR, rows, cols, H, rows, tau);
- if (rc==0) {
- qcols= std::min(rows,cols);
- rc = LAPACKE_dorgqr(LAPACK_COL_MAJOR, rows, qcols, qcols, H, rows, tau);
- }
- delete[] tau;
- return rc;
-}
-
-template<>
-inline int qr_i_T<float>(const uint32_t rows, const uint32_t cols, float *H, uint32_t &qcols) {
- qcols= std::min(rows,cols);
- float *tau= new float[qcols];
- int rc = LAPACKE_sgeqrf(LAPACK_COL_MAJOR, rows, cols, H, rows, tau);
- if (rc==0) {
- qcols= std::min(rows,cols);
- rc = LAPACKE_sorgqr(LAPACK_COL_MAJOR, rows, qcols, qcols, H, rows, tau);
- }
- delete[] tau;
- return rc;
+static inline int qr_i_T(const uint32_t rows, const uint32_t cols, TReal *H, uint32_t &qcols) {
+ return skbb::linalg::qr_inplace(rows, cols, H, qcols);
}
namespace skbb {
@@ -326,47 +290,31 @@ class QR {
template<>
inline void skbb::QR<double>::qdot_r_sq(const double *mat, double *res) {
- cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans, rows , cols, rows, 1.0, mat, rows, Q, rows, 0.0, res, rows);
+ skbb::linalg::gemm_nn(rows, cols, rows, mat, Q, res);
}
template<>
inline void skbb::QR<float>::qdot_r_sq(const float *mat, float *res) {
- cblas_sgemm(CblasColMajor,CblasNoTrans,CblasNoTrans, rows , cols, rows, 1.0, mat, rows, Q, rows, 0.0, res, rows);
+ skbb::linalg::gemm_nn(rows, cols, rows, mat, Q, res);
}
template<>
inline void skbb::QR<double>::qdot_l_sq(const double *mat, double *res) {
- cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans, rows , cols, cols, 1.0, Q, rows, mat, cols, 0.0, res, rows);
+ skbb::linalg::gemm_nn(rows, cols, cols, Q, mat, res);
}
template<>
void skbb::QR<float>::qdot_l_sq(const float *mat, float *res) {
- cblas_sgemm(CblasColMajor,CblasNoTrans,CblasNoTrans, rows , cols, cols, 1.0, Q, rows, mat, cols, 0.0, res, rows);
+ skbb::linalg::gemm_nn(rows, cols, cols, Q, mat, res);
}
// compute svd, and return S and V
// T = input
// S output
-// T is Vt on output
+// T is Vt on output (first min(rows,cols) rows hold V^T)
template<class TReal>
-inline int svd_it_T(const uint32_t rows, const uint32_t cols, TReal *T, TReal *S);
-
-template<>
-inline int svd_it_T<double>(const uint32_t rows, const uint32_t cols, double *T, double *S) {
- double *superb = (double *) malloc(sizeof(double)*rows);
- int res =LAPACKE_dgesvd(LAPACK_COL_MAJOR, 'N', 'O', rows, cols, T, rows, S, NULL, rows, NULL, cols, superb);
- free(superb);
-
- return res;
-}
-
-template<>
-inline int svd_it_T<float>(const uint32_t rows, const uint32_t cols, float *T, float *S) {
- float *superb = (float *) malloc(sizeof(float)*rows);
- int res =LAPACKE_sgesvd(LAPACK_COL_MAJOR, 'N', 'O', rows, cols, T, rows, S, NULL, rows, NULL, cols, superb);
- free(superb);
-
- return res;
+inline int svd_it_T(const uint32_t rows, const uint32_t cols, TReal *T, TReal *S) {
+ return skbb::linalg::svd_no(rows, cols, T, S);
}
// square matrix transpose, with org not alingned
=====================================
src/tests/wasm/generate_pcoa_expected.cpp
=====================================
@@ -0,0 +1,117 @@
+/*
+ * BSD 3-Clause License
+ *
+ * Copyright (c) 2025--, scikit-bio development team.
+ * All rights reserved.
+ *
+ * See LICENSE file for more details
+ */
+
+/*
+ * Native binary that runs each PCoACase through skbb::pcoa_fsvd and
+ * writes the observed (eigenvalues, samples, proportion_explained) as a
+ * C header consumed by test_pcoa_wasm.cpp.
+ *
+ * Unlike PERMANOVA, PCoA results are allowed to differ between
+ * native-LAPACK and WASM-Eigen paths (different QR/SVD algorithms can
+ * produce different rotations / sign conventions). The test uses
+ * absolute-value comparison per axis and a generous tolerance; the
+ * expected values captured here are the LAPACK reference.
+ *
+ * MUST be invoked with OMP_NUM_THREADS=1. PCoA doesn't use the RNG-
+ * chunking heuristic directly, but skbb::set_random_seed(s) is the only
+ * input to Halko's Gaussian matrix: if something changes upstream that
+ * makes this sensitive to thread count we want a deterministic baseline.
+ */
+
+#include "ordination/principal_coordinate_analysis.hpp"
+#include "util/rand.hpp"
+
+#include "tests/wasm/pcoa_inputs.hpp"
+
+#include <cstdio>
+#include <cstdlib>
+#include <string>
+
+static void emit_double_array(const char *name, const double *buf, unsigned int n) {
+ std::printf("static const double %s[] = {\n", name);
+ for (unsigned int i = 0; i < n; ++i) {
+ std::printf(" %a%s\n", buf[i], (i + 1 == n) ? "" : ",");
+ }
+ std::printf("};\n\n");
+}
+
+int main() {
+ const char *nt = std::getenv("OMP_NUM_THREADS");
+ if (nt == nullptr || std::string(nt) != "1") {
+ std::fprintf(stderr,
+ "ERROR: generate_pcoa_expected must run with "
+ "OMP_NUM_THREADS=1 (got '%s').\n",
+ nt ? nt : "(unset)");
+ return 2;
+ }
+
+ using namespace skbb_wasm_test;
+
+ std::printf("/* GENERATED by generate_pcoa_expected.cpp — do not edit. */\n");
+ std::printf("#ifndef SKBB_WASM_PCOA_EXPECTED_H\n");
+ std::printf("#define SKBB_WASM_PCOA_EXPECTED_H\n\n");
+ std::printf("namespace skbb_wasm_test {\n\n");
+
+ for (unsigned int i = 0; i < kPCoACaseCount; ++i) {
+ const PCoACase &c = kPCoACases[i];
+ const unsigned int n = c.n;
+ const unsigned int k = c.n_eighs;
+
+ double *eigenvalues = new double[k];
+ double *samples = new double[uint64_t(n) * k];
+ double *prop = new double[k];
+
+ // Seed skbb's global RNG as well, for parity with the WASM path
+ // even though pcoa_fsvd itself accepts a seed parameter.
+ skbb::set_random_seed(static_cast<unsigned int>(c.seed));
+ skbb::pcoa_fsvd(n, c.mat, k, c.seed,
+ eigenvalues, samples, prop);
+
+ char base[128];
+ std::snprintf(base, sizeof(base), "pcoa_%s", c.name);
+
+ {
+ char buf[160];
+ std::snprintf(buf, sizeof(buf), "%s_eigenvalues", base);
+ emit_double_array(buf, eigenvalues, k);
+ }
+ {
+ char buf[160];
+ std::snprintf(buf, sizeof(buf), "%s_samples", base);
+ emit_double_array(buf, samples, uint64_t(n) * k);
+ }
+ {
+ char buf[160];
+ std::snprintf(buf, sizeof(buf), "%s_proportion_explained", base);
+ emit_double_array(buf, prop, k);
+ }
+
+ delete[] eigenvalues;
+ delete[] samples;
+ delete[] prop;
+ }
+
+ // Emit a dispatch table so the test can look up expected arrays by index.
+ std::printf("struct PCoAExpected {\n"
+ " const double *eigenvalues;\n"
+ " const double *samples;\n"
+ " const double *proportion_explained;\n"
+ "};\n\n");
+ std::printf("static const PCoAExpected kPCoAExpected[] = {\n");
+ for (unsigned int i = 0; i < kPCoACaseCount; ++i) {
+ const PCoACase &c = kPCoACases[i];
+ std::printf(" { pcoa_%s_eigenvalues, pcoa_%s_samples, pcoa_%s_proportion_explained },\n",
+ c.name, c.name, c.name);
+ }
+ std::printf("};\n\n");
+
+ std::printf("} // namespace skbb_wasm_test\n\n");
+ std::printf("#endif\n");
+ return 0;
+}
=====================================
src/tests/wasm/generate_permanova_expected.cpp
=====================================
@@ -0,0 +1,70 @@
+/*
+ * BSD 3-Clause License
+ *
+ * Copyright (c) 2025--, scikit-bio development team.
+ * All rights reserved.
+ *
+ * See LICENSE file for more details
+ */
+
+/*
+ * Native binary that runs each PermanovaCase through skbb::permanova and
+ * emits the observed (fstat, pvalue) values as a C header consumed by
+ * test_permanova_wasm.cpp.
+ *
+ * The header is checked into the repo so the WASM test is hermetic:
+ * it doesn't need skbb's native build at test time, only the precomputed
+ * expected values.
+ *
+ * MUST be invoked with OMP_NUM_THREADS=1. The WASM build is
+ * single-threaded, and pmn_get_max_parallelism() returns 32 there
+ * (matching 2*1*16 on native). Any other thread count changes the RNG
+ * chunking and thus the permutation sequence, breaking bit-equality.
+ * The Makefile rule enforces this.
+ *
+ * Regenerate with: make -C src wasm_regen_permanova_expected
+ */
+
+#include "distance/permanova.hpp"
+#include "tests/wasm/permanova_inputs.hpp"
+
+#include <cstdio>
+#include <cstdlib>
+#include <string>
+
+int main() {
+ // Belt-and-suspenders: abort if the operator forgot OMP_NUM_THREADS=1.
+ const char *nt = std::getenv("OMP_NUM_THREADS");
+ if (nt == nullptr || std::string(nt) != "1") {
+ std::fprintf(stderr,
+ "ERROR: generate_permanova_expected must run with "
+ "OMP_NUM_THREADS=1 (got '%s').\n",
+ nt ? nt : "(unset)");
+ return 2;
+ }
+
+ using namespace skbb_wasm_test;
+
+ std::printf("/* GENERATED by generate_permanova_expected.cpp — do not edit. */\n");
+ std::printf("#ifndef SKBB_WASM_PERMANOVA_EXPECTED_H\n");
+ std::printf("#define SKBB_WASM_PERMANOVA_EXPECTED_H\n\n");
+ std::printf("namespace skbb_wasm_test {\n\n");
+ std::printf("struct PermanovaExpected { double fstat; double pvalue; };\n\n");
+ std::printf("static const PermanovaExpected kPermanovaExpected[] = {\n");
+
+ for (unsigned int i = 0; i < kPermanovaCaseCount; ++i) {
+ const PermanovaCase &c = kPermanovaCases[i];
+ double fstat = 0.0, pvalue = 0.0;
+ skbb::permanova(c.n, c.mat, c.grouping, c.n_perm, c.seed, fstat, pvalue);
+ // %a emits the exact hex representation of the double, so no
+ // precision is lost in the round-trip. The decimal comment is for
+ // human readability when the header shows up in diffs.
+ std::printf(" { %a, %a }, /* %s: fstat=%.17g pvalue=%.17g */\n",
+ fstat, pvalue, c.name, fstat, pvalue);
+ }
+
+ std::printf("};\n\n");
+ std::printf("} // namespace skbb_wasm_test\n\n");
+ std::printf("#endif\n");
+ return 0;
+}
=====================================
src/tests/wasm/pcoa_inputs.hpp
=====================================
@@ -0,0 +1,59 @@
+/*
+ * BSD 3-Clause License
+ *
+ * Copyright (c) 2025--, scikit-bio development team.
+ * All rights reserved.
+ *
+ * See LICENSE file for more details
+ */
+
+/*
+ * Fixed PCoA inputs for native-vs-WASM correctness comparison.
+ *
+ * For each case we pin a deterministic seed so FSVD's Halko-randomized
+ * initial matrix is reproducible. The distance matrices themselves are
+ * checked in so the generator and the WASM test read from the same data.
+ *
+ * The synthetic matrices are built to have a well-separated eigenvalue
+ * spectrum (via a dominant rank-k contribution), so per-axis sample
+ * coordinates can be compared up to sign with a tight tolerance instead
+ * of collapsing into a subspace-only check.
+ */
+
+#ifndef SKBB_WASM_PCOA_INPUTS_HPP
+#define SKBB_WASM_PCOA_INPUTS_HPP
+
+#include <stdint.h>
+
+namespace skbb_wasm_test {
+
+struct PCoACase {
+ const char *name;
+ unsigned int n; // matrix dimension (n_samples)
+ unsigned int n_eighs; // number of eigenvalues requested
+ int seed;
+ const double *mat; // row-major, n*n doubles
+};
+
+// ---------------------------------------------------------------------
+// Case P1: uuf_6x6 — same 6x6 matrix as the native test_center_mat test.
+// ---------------------------------------------------------------------
+static const double pcoa_uuf_6x6[] = {
+ 0.0000000000, 0.2000000000, 0.5714285714, 0.6000000000, 0.5000000000, 0.2000000000,
+ 0.2000000000, 0.0000000000, 0.4285714286, 0.6666666667, 0.6000000000, 0.3333333333,
+ 0.5714285714, 0.4285714286, 0.0000000000, 0.7142857143, 0.8571428571, 0.4285714286,
+ 0.6000000000, 0.6666666667, 0.7142857143, 0.0000000000, 0.3333333333, 0.4000000000,
+ 0.5000000000, 0.6000000000, 0.8571428571, 0.3333333333, 0.0000000000, 0.6000000000,
+ 0.2000000000, 0.3333333333, 0.4285714286, 0.4000000000, 0.6000000000, 0.0000000000,
+};
+
+static const PCoACase kPCoACases[] = {
+ { "uuf_6x6_k3", 6, 3, 12345, pcoa_uuf_6x6 },
+ { "uuf_6x6_k2", 6, 2, 12345, pcoa_uuf_6x6 },
+};
+static const unsigned int kPCoACaseCount =
+ sizeof(kPCoACases) / sizeof(kPCoACases[0]);
+
+} // namespace skbb_wasm_test
+
+#endif
=====================================
src/tests/wasm/permanova_inputs.hpp
=====================================
@@ -0,0 +1,80 @@
+/*
+ * BSD 3-Clause License
+ *
+ * Copyright (c) 2025--, scikit-bio development team.
+ * All rights reserved.
+ *
+ * See LICENSE file for more details
+ */
+
+/*
+ * Shared fixed inputs for the WASM PERMANOVA correctness test.
+ * Consumed by both generate_permanova_expected (native) and
+ * test_permanova_wasm (WASM), so the two programs cannot drift.
+ *
+ * Every case pins a seed so std::mt19937 + std::shuffle produce
+ * identical permutation sequences across native and WASM, yielding
+ * bit-identical fstat and pvalue.
+ */
+
+#ifndef SKBB_WASM_PERMANOVA_INPUTS_HPP
+#define SKBB_WASM_PERMANOVA_INPUTS_HPP
+
+#include <stdint.h>
+
+namespace skbb_wasm_test {
+
+struct PermanovaCase {
+ const char *name;
+ unsigned int n;
+ const double *mat;
+ const unsigned int *grouping;
+ unsigned int n_perm;
+ int seed;
+};
+
+// Case 1: 4x4, two obvious groups; matches the smoke test.
+static const double case_smoke_mat[] = {
+ 0.0, 0.2, 0.4, 0.1,
+ 0.2, 0.0, 0.3, 0.5,
+ 0.4, 0.3, 0.0, 0.2,
+ 0.1, 0.5, 0.2, 0.0,
+};
+static const unsigned int case_smoke_grp[] = {0, 0, 1, 1};
+
+// Case 2: 6x6 unweighted-unifrac-shaped matrix from the existing native
+// test suite (src/tests/test_pcoa.cpp:test_center_mat); three groups.
+static const double case_uuf_mat[] = {
+ 0.0000000000, 0.2000000000, 0.5714285714, 0.6000000000, 0.5000000000, 0.2000000000,
+ 0.2000000000, 0.0000000000, 0.4285714286, 0.6666666667, 0.6000000000, 0.3333333333,
+ 0.5714285714, 0.4285714286, 0.0000000000, 0.7142857143, 0.8571428571, 0.4285714286,
+ 0.6000000000, 0.6666666667, 0.7142857143, 0.0000000000, 0.3333333333, 0.4000000000,
+ 0.5000000000, 0.6000000000, 0.8571428571, 0.3333333333, 0.0000000000, 0.6000000000,
+ 0.2000000000, 0.3333333333, 0.4285714286, 0.4000000000, 0.6000000000, 0.0000000000,
+};
+static const unsigned int case_uuf_grp[] = {0, 0, 1, 1, 2, 2};
+
+// Case 3: 8x8 equal-group no-signal (ties), from tests/test_permanova.cpp style.
+static const double case_ties_mat[] = {
+ 0.0, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
+ 0.5, 0.0, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
+ 0.5, 0.5, 0.0, 0.5, 0.5, 0.5, 0.5, 0.5,
+ 0.5, 0.5, 0.5, 0.0, 0.5, 0.5, 0.5, 0.5,
+ 0.5, 0.5, 0.5, 0.5, 0.0, 0.5, 0.5, 0.5,
+ 0.5, 0.5, 0.5, 0.5, 0.5, 0.0, 0.5, 0.5,
+ 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.0, 0.5,
+ 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.0,
+};
+static const unsigned int case_ties_grp[] = {0, 0, 0, 0, 1, 1, 1, 1};
+
+static const PermanovaCase kPermanovaCases[] = {
+ { "smoke_4x4", 4, case_smoke_mat, case_smoke_grp, 99, 42 },
+ { "uuf_6x6", 6, case_uuf_mat, case_uuf_grp, 199, 1 },
+ { "ties_8x8", 8, case_ties_mat, case_ties_grp, 499, 7 },
+};
+static const unsigned int kPermanovaCaseCount =
+ sizeof(kPermanovaCases) / sizeof(kPermanovaCases[0]);
+
+} // namespace skbb_wasm_test
+
+#endif
=====================================
src/tests/wasm/test_center_wasm.cpp
=====================================
@@ -0,0 +1,105 @@
+/*
+ * BSD 3-Clause License
+ *
+ * Copyright (c) 2025--, scikit-bio development team.
+ * All rights reserved.
+ *
+ * See LICENSE file for more details
+ */
+
+/*
+ * Stage 5: mat_to_centered correctness under WASM.
+ *
+ * Uses the same 6x6 unweighted-unifrac distance matrix and the same
+ * hard-coded expected centered matrix as the existing native test
+ * (src/tests/test_pcoa.cpp:test_center_mat). Tolerance 1e-6 matches
+ * the native test's vec_almost_equal.
+ *
+ * Centering is deterministic, platform-independent arithmetic (no RNG,
+ * no BLAS — it's a hand-rolled O(n^2) double-loop in
+ * principal_coordinate_analysis.cpp). This test should always pass
+ * bit-exactly modulo the accumulated 1e-6 rounding in row_sum.
+ */
+
+#include "scikit-bio-binaries/ordination.h"
+
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+
+static int failures = 0;
+#define CHECK(cond) do { \
+ if (!(cond)) { \
+ std::fprintf(stderr, "FAIL [%s:%d] %s\n", __FILE__, __LINE__, #cond); \
+ ++failures; \
+ } \
+} while (0)
+
+int main() {
+ // Input: unweighted-UniFrac distance matrix from test.biom, per
+ // src/tests/test_pcoa.cpp:16-26.
+ double matrix[] = {
+ 0.0000000000, 0.2000000000, 0.5714285714, 0.6000000000, 0.5000000000, 0.2000000000,
+ 0.2000000000, 0.0000000000, 0.4285714286, 0.6666666667, 0.6000000000, 0.3333333333,
+ 0.5714285714, 0.4285714286, 0.0000000000, 0.7142857143, 0.8571428571, 0.4285714286,
+ 0.6000000000, 0.6666666667, 0.7142857143, 0.0000000000, 0.3333333333, 0.4000000000,
+ 0.5000000000, 0.6000000000, 0.8571428571, 0.3333333333, 0.0000000000, 0.6000000000,
+ 0.2000000000, 0.3333333333, 0.4285714286, 0.4000000000, 0.6000000000, 0.0000000000,
+ };
+ const unsigned int n = 6;
+
+ // Expected centered matrix from the same native test source.
+ const double exp[] = {
+ 0.05343726, 0.04366213, -0.0329743, -0.07912698, -0.00495654, 0.01995843,
+ 0.04366213, 0.073887, 0.04867914, -0.11112434, -0.04973167, -0.00537226,
+ -0.0329743, 0.04867914, 0.20714475, -0.07737528, -0.17044974, 0.02497543,
+ -0.07912698, -0.11112434, -0.07737528, 0.14830877, 0.11192366, 0.00739418,
+ -0.00495654, -0.04973167, -0.17044974, 0.11192366, 0.18664966, -0.07343537,
+ 0.01995843, -0.00537226, 0.02497543, 0.00739418, -0.07343537, 0.02647959,
+ };
+
+ const double tol = 1e-6;
+
+ // fp64 (out-of-place)
+ {
+ double centered[36];
+ skbb_center_distance_matrix_fp64(n, matrix, centered);
+ for (int i = 0; i < 36; ++i) {
+ CHECK(std::fabs(centered[i] - exp[i]) < tol);
+ }
+ }
+
+ // fp64 -> fp32 mixed
+ {
+ float centered_fp32[36];
+ skbb_center_distance_matrix_fp64_to_fp32(n, matrix, centered_fp32);
+ for (int i = 0; i < 36; ++i) {
+ CHECK(std::fabs(static_cast<double>(centered_fp32[i]) - exp[i]) < tol);
+ }
+ }
+
+ // fp32 out-of-place
+ {
+ float matrix_fp32[36];
+ for (int i = 0; i < 36; ++i) matrix_fp32[i] = static_cast<float>(matrix[i]);
+ float centered_fp32[36];
+ skbb_center_distance_matrix_fp32(n, matrix_fp32, centered_fp32);
+ for (int i = 0; i < 36; ++i) {
+ CHECK(std::fabs(static_cast<double>(centered_fp32[i]) - exp[i]) < tol);
+ }
+ }
+
+ // fp64 in-place
+ {
+ double buf[36];
+ for (int i = 0; i < 36; ++i) buf[i] = matrix[i];
+ skbb_center_distance_matrix_fp64(n, buf, buf);
+ for (int i = 0; i < 36; ++i) {
+ CHECK(std::fabs(buf[i] - exp[i]) < tol);
+ }
+ }
+
+ std::printf("%s: centering checks, %d failures\n",
+ failures == 0 ? "PASS" : "FAIL", failures);
+ return failures == 0 ? 0 : 1;
+}
=====================================
src/tests/wasm/test_pcoa_wasm.cpp
=====================================
@@ -0,0 +1,161 @@
+/*
+ * BSD 3-Clause License
+ *
+ * Copyright (c) 2025--, scikit-bio development team.
+ * All rights reserved.
+ *
+ * See LICENSE file for more details
+ */
+
+/*
+ * Stage 6: PCoA (FSVD) correctness under WASM.
+ *
+ * Replays each PCoACase through the WASM build of skbb_pcoa_fsvd_fp64 and
+ * compares against the native (LAPACK) expected values.
+ *
+ * Tolerance policy:
+ * - eigenvalues: 1e-6 absolute. These are a deterministic function of
+ * the centered distance matrix's dominant eigenspace and should
+ * agree closely across QR/SVD implementations, even though Eigen's
+ * BDCSVD and LAPACK's dgesvd use different bidiagonalization paths.
+ * - proportion_explained: same tolerance (it's eigenvalues / trace).
+ * - samples (coordinates): sign-adjusted per axis, 1e-3 absolute.
+ * Eigenvectors are unique only up to sign (and, in degenerate
+ * eigenspaces, up to rotation). We compare axes independently and
+ * pick the sign that minimizes per-case L2 distance.
+ *
+ * If a future test case has near-degenerate eigenvalues the sign-only
+ * fixup is insufficient and the test would need a subspace-distance
+ * helper. The current test cases (6x6 unweighted-UniFrac) have a
+ * well-separated spectrum and don't need that.
+ */
+
+#include "scikit-bio-binaries/ordination.h"
+#include "scikit-bio-binaries/util.h"
+
+#include "tests/wasm/pcoa_inputs.hpp"
+#include "tests/wasm/expected/pcoa_expected.h"
+
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+
+static int failures = 0;
+
+static bool almost_equal_vec(const char *label, const double *got,
+ const double *want, unsigned int len,
+ double tol) {
+ double max_abs_err = 0.0;
+ for (unsigned int i = 0; i < len; ++i) {
+ double e = std::fabs(got[i] - want[i]);
+ if (e > max_abs_err) max_abs_err = e;
+ }
+ bool ok = max_abs_err <= tol;
+ std::printf(" %s %s max|delta|=%.3e tol=%.3e len=%u\n",
+ label, ok ? "OK" : "FAIL", max_abs_err, tol, len);
+ return ok;
+}
+
+// Compare two same-shape coordinate matrices (rows x cols, row-major) where
+// each column is an eigenvector axis up to sign. For every axis, try both
+// sign assignments and accept if either fits within tol.
+static bool almost_equal_samples_up_to_sign(const double *got,
+ const double *want,
+ unsigned int rows,
+ unsigned int cols,
+ double tol) {
+ double worst = 0.0;
+ unsigned int worst_col = 0;
+ int worst_sign = 0;
+ for (unsigned int col = 0; col < cols; ++col) {
+ double err_pos = 0.0, err_neg = 0.0;
+ for (unsigned int row = 0; row < rows; ++row) {
+ double g = got[row * cols + col];
+ double w = want[row * cols + col];
+ double ep = std::fabs(g - w);
+ double en = std::fabs(g + w);
+ if (ep > err_pos) err_pos = ep;
+ if (en > err_neg) err_neg = en;
+ }
+ double ax_err = std::fmin(err_pos, err_neg);
+ int ax_sign = (err_pos <= err_neg) ? +1 : -1;
+ if (ax_err > worst) {
+ worst = ax_err;
+ worst_col = col;
+ worst_sign = ax_sign;
+ }
+ }
+ bool ok = worst <= tol;
+ std::printf(" samples %s worst-axis=%u sign=%+d max|delta|=%.3e tol=%.3e\n",
+ ok ? "OK" : "FAIL", worst_col, worst_sign, worst, tol);
+ return ok;
+}
+
+int main() {
+ using namespace skbb_wasm_test;
+
+ static_assert(sizeof(kPCoAExpected) / sizeof(kPCoAExpected[0])
+ == kPCoACaseCount,
+ "expected-table length must match input-table length");
+
+ const double eigen_tol = 1e-6;
+ const double prop_tol = 1e-6;
+ const double sample_tol = 1e-3;
+
+ for (unsigned int i = 0; i < kPCoACaseCount; ++i) {
+ const PCoACase &c = kPCoACases[i];
+ const PCoAExpected &e = kPCoAExpected[i];
+
+ const unsigned int n = c.n;
+ const unsigned int k = c.n_eighs;
+
+ double *eigenvalues = new double[k];
+ double *samples = new double[static_cast<uint64_t>(n) * k];
+ double *prop = new double[k];
+
+ // Seed skbb's global RNG for Halko's Gaussian matrix, matching
+ // the native generator.
+ skbb_set_random_seed(static_cast<unsigned int>(c.seed));
+ skbb_pcoa_fsvd_fp64(n, c.mat, k, c.seed,
+ eigenvalues, samples, prop);
+
+ std::printf(" %s\n", c.name);
+ bool ok = true;
+ ok &= almost_equal_vec("eigenvalues", eigenvalues, e.eigenvalues, k, eigen_tol);
+ ok &= almost_equal_vec("prop_expl ", prop, e.proportion_explained, k, prop_tol);
+ ok &= almost_equal_samples_up_to_sign(samples, e.samples, n, k, sample_tol);
+
+ // Determinism: same seed, same inputs -> identical output on this build.
+ double *eigenvalues2 = new double[k];
+ double *samples2 = new double[static_cast<uint64_t>(n) * k];
+ double *prop2 = new double[k];
+ skbb_set_random_seed(static_cast<unsigned int>(c.seed));
+ skbb_pcoa_fsvd_fp64(n, c.mat, k, c.seed,
+ eigenvalues2, samples2, prop2);
+
+ bool det_eigen = true, det_samp = true, det_prop = true;
+ for (unsigned int j = 0; j < k; ++j) {
+ if (eigenvalues[j] != eigenvalues2[j]) det_eigen = false;
+ if (prop[j] != prop2[j]) det_prop = false;
+ }
+ for (uint64_t j = 0; j < static_cast<uint64_t>(n) * k; ++j) {
+ if (samples[j] != samples2[j]) det_samp = false;
+ }
+ std::printf(" determinism %s (eigen=%s samples=%s prop=%s)\n",
+ (det_eigen && det_samp && det_prop) ? "OK" : "FAIL",
+ det_eigen ? "OK" : "FAIL",
+ det_samp ? "OK" : "FAIL",
+ det_prop ? "OK" : "FAIL");
+ if (!(det_eigen && det_samp && det_prop)) ok = false;
+
+ if (!ok) ++failures;
+
+ delete[] eigenvalues; delete[] samples; delete[] prop;
+ delete[] eigenvalues2; delete[] samples2; delete[] prop2;
+ }
+
+ std::printf("%s: %u cases, %d failed\n",
+ failures == 0 ? "PASS" : "FAIL",
+ kPCoACaseCount, failures);
+ return failures == 0 ? 0 : 1;
+}
=====================================
src/tests/wasm/test_permanova_wasm.cpp
=====================================
@@ -0,0 +1,86 @@
+/*
+ * BSD 3-Clause License
+ *
+ * Copyright (c) 2025--, scikit-bio development team.
+ * All rights reserved.
+ *
+ * See LICENSE file for more details
+ */
+
+/*
+ * Stage 3: PERMANOVA correctness under WASM.
+ *
+ * Replays the same fixed inputs (permanova_inputs.hpp) through the WASM
+ * build of skbb_permanova_fp64 and compares against the native expected
+ * values captured in expected/permanova_expected.h. Expected values are
+ * generated by a single-threaded native build — see
+ * generate_permanova_expected.cpp.
+ *
+ * Tolerance policy (both bit-identical):
+ * - fstat: the unpermuted pseudo-F is a deterministic function of
+ * the distance matrix and grouping with no RNG dependence.
+ * - pvalue: skbb uses a portable Fisher-Yates (see
+ * util/portable_shuffle.hpp) built on the raw std::mt19937
+ * 32-bit output. mt19937 is fully specified by the C++
+ * standard, so the permutation sequence is reproducible
+ * across libstdc++ (native gcc) and libc++ (emcc).
+ *
+ * Determinism check: the WASM permanova is run twice with identical
+ * inputs and the results must be bit-identical — a sanity check against
+ * RNG state leakage between calls.
+ *
+ * Exit 0 on success, nonzero on any mismatch.
+ */
+
+#include "scikit-bio-binaries/distance.h"
+
+#include "tests/wasm/permanova_inputs.hpp"
+#include "tests/wasm/expected/permanova_expected.h"
+
+#include <cstdio>
+#include <cstdlib>
+
+int main() {
+ using namespace skbb_wasm_test;
+
+ static_assert(sizeof(kPermanovaExpected) / sizeof(kPermanovaExpected[0])
+ == kPermanovaCaseCount,
+ "expected-table length must match input-table length");
+
+ int failures = 0;
+
+ for (unsigned int i = 0; i < kPermanovaCaseCount; ++i) {
+ const PermanovaCase &c = kPermanovaCases[i];
+ const PermanovaExpected &e = kPermanovaExpected[i];
+
+ double fstat = 0.0, pvalue = 0.0;
+ skbb_permanova_fp64(c.n, c.mat, c.grouping,
+ c.n_perm, c.seed, &fstat, &pvalue);
+
+ // Determinism: second run with identical inputs, same seed.
+ double fstat2 = 0.0, pvalue2 = 0.0;
+ skbb_permanova_fp64(c.n, c.mat, c.grouping,
+ c.n_perm, c.seed, &fstat2, &pvalue2);
+
+ const bool fstat_ok = (fstat == e.fstat);
+ const bool pvalue_ok = (pvalue == e.pvalue);
+ const bool det_ok = (fstat == fstat2) && (pvalue == pvalue2);
+
+ std::printf(" %-12s fstat %s (got %.17g, want %.17g) "
+ "pvalue %s (got %.17g, want %.17g) "
+ "det %s\n",
+ c.name,
+ fstat_ok ? "OK" : "FAIL", fstat, e.fstat,
+ pvalue_ok ? "OK" : "FAIL", pvalue, e.pvalue,
+ det_ok ? "OK" : "FAIL");
+
+ if (!fstat_ok) ++failures;
+ if (!pvalue_ok) ++failures;
+ if (!det_ok) ++failures;
+ }
+
+ std::printf("%s: %u cases, %d failures\n",
+ failures == 0 ? "PASS" : "FAIL",
+ kPermanovaCaseCount, failures);
+ return failures == 0 ? 0 : 1;
+}
=====================================
src/tests/wasm/test_smoke.cpp
=====================================
@@ -0,0 +1,60 @@
+/*
+ * BSD 3-Clause License
+ *
+ * Copyright (c) 2025--, scikit-bio development team.
+ * All rights reserved.
+ *
+ * See LICENSE file for more details
+ */
+
+/*
+ * Stage 2 smoke test: confirms libskbb_wasm.a links and runs under node.
+ * Exercises only the non-BLAS surface (get_api_version + permanova_fp64).
+ * Full numerical correctness is validated in test_permanova_wasm.cpp
+ * (Stage 3) and test_pcoa_wasm.cpp (Stage 6).
+ *
+ * Exit 0 on success, nonzero on any assertion.
+ */
+
+#include "scikit-bio-binaries/util.h"
+#include "scikit-bio-binaries/distance.h"
+
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+
+static int failures = 0;
+#define CHECK(cond) do { \
+ if (!(cond)) { \
+ std::fprintf(stderr, "FAIL [%s:%d] %s\n", __FILE__, __LINE__, #cond); \
+ ++failures; \
+ } \
+} while (0)
+
+int main() {
+ CHECK(skbb_get_api_version() >= 1);
+
+ // Minimal PERMANOVA smoke: 4 samples, 2 groups, no separation.
+ // Expected: f-stat finite and non-negative; pvalue in [0, 1].
+ const unsigned int n = 4;
+ const double mat[] = {
+ 0.0, 0.2, 0.4, 0.1,
+ 0.2, 0.0, 0.3, 0.5,
+ 0.4, 0.3, 0.0, 0.2,
+ 0.1, 0.5, 0.2, 0.0,
+ };
+ const unsigned int grouping[] = {0, 0, 1, 1};
+
+ double fstat = -1.0, pvalue = -1.0;
+ skbb_permanova_fp64(n, mat, grouping, /*n_perm=*/99, /*seed=*/42,
+ &fstat, &pvalue);
+
+ CHECK(std::isfinite(fstat));
+ CHECK(fstat >= 0.0);
+ CHECK(pvalue >= 0.0);
+ CHECK(pvalue <= 1.0);
+
+ std::printf("smoke: api_version=%u fstat=%.6f pvalue=%.6f failures=%d\n",
+ skbb_get_api_version(), fstat, pvalue, failures);
+ return failures == 0 ? 0 : 1;
+}
=====================================
src/tools/skbb_generate_helper.py
=====================================
@@ -234,7 +234,7 @@ def print_body(method,lines,nmspace):
line = lines[i]
i+=1
- for ft in ftypes:
+ for ft in sorted(ftypes):
print_func_args(method,ftype,nmspace,fname,ft,fargs)
print('');
=====================================
src/util/portable_shuffle.hpp
=====================================
@@ -0,0 +1,70 @@
+/*
+ * BSD 3-Clause License
+ *
+ * Copyright (c) 2025--, scikit-bio development team.
+ * All rights reserved.
+ *
+ * See LICENSE file for more details
+ */
+
+/*
+ * Portable Fisher-Yates shuffle built directly on top of std::mt19937.
+ *
+ * We cannot use std::shuffle: it delegates to std::uniform_int_distribution,
+ * whose exact mapping from RNG output to bounded integer is implementation
+ * defined (libstdc++, libc++, and MSVC each make a different choice). That
+ * means seeded results from std::shuffle are not reproducible across
+ * toolchains — which breaks bit-identical PERMANOVA p-values between the
+ * native gcc/libstdc++ build and the emscripten/libc++ WASM build.
+ *
+ * std::mt19937 itself is fully specified by the C++ standard (period,
+ * state, output sequence), so building the bounded reduction ourselves
+ * removes the only source of cross-toolchain drift in PERMANOVA.
+ *
+ * The reduction here is the textbook rejection method:
+ * threshold = (2^32 - bound) % bound
+ * loop: r = rng(); if (r < threshold) repeat; else return r % bound
+ * which is unbiased and deterministic across platforms.
+ */
+
+#ifndef SKBB_PORTABLE_SHUFFLE_HPP
+#define SKBB_PORTABLE_SHUFFLE_HPP
+
+#include <stdint.h>
+
+namespace skbb {
+
+// Uniform integer in [0, bound) using rejection on the raw 32-bit output of
+// the supplied RNG. Requires bound >= 1.
+template <class RNG>
+static inline uint32_t portable_uniform_u32(uint32_t bound, RNG &rng) {
+ // mt19937 returns values in [0, 2^32). `(uint32_t)(-bound)` is
+ // 2^32 - bound in 32-bit wrap-around arithmetic, and modding that by
+ // bound gives the smallest rejection threshold that yields a uniform
+ // distribution in [0, bound).
+ const uint32_t threshold =
+ static_cast<uint32_t>(0u - bound) % bound;
+ uint32_t r;
+ do {
+ r = static_cast<uint32_t>(rng());
+ } while (r < threshold);
+ return r % bound;
+}
+
+// In-place Fisher-Yates shuffle of `first[0..n)` driven by `rng`. Matches
+// the element order std::shuffle produces on libstdc++ only by coincidence
+// — use this explicitly when you want cross-toolchain reproducibility.
+template <class T, class RNG>
+static inline void portable_shuffle(T *first, uint32_t n, RNG &rng) {
+ for (uint32_t i = n; i > 1; --i) {
+ const uint32_t j = portable_uniform_u32(i, rng);
+ // swap first[i-1], first[j]
+ T tmp = first[i - 1];
+ first[i - 1] = first[j];
+ first[j] = tmp;
+ }
+}
+
+} // namespace skbb
+
+#endif
=====================================
src/util/skbb_detect_acc.cpp
=====================================
@@ -27,6 +27,7 @@
#endif
#include <stdlib.h>
+#include <cstdlib>
#include <string>
#include <string.h>
=====================================
src/util/skbb_dgb_info.hpp
=====================================
@@ -11,6 +11,8 @@
#define SKBB_DBG_INFO_HPP
#include <chrono>
+#include <cstdlib>
+#include <string>
// To be used once per function
#define SETUP_TDBG(method) const char *tdbg_method=method; \
=====================================
src/wasm/emscripten_build.mk
=====================================
@@ -0,0 +1,149 @@
+# WebAssembly build rules for scikit-bio-binaries.
+#
+# Activated by the top-level `wasm` target. Produces libskbb_wasm.a, a
+# single-threaded static archive intended to be linked into downstream
+# emscripten projects (e.g. unifrac-binaries/duckdb-miint).
+#
+# Backend: Eigen header-only (fetched via scripts/fetch_eigen.sh).
+# No OpenMP, no pthread, no GPU, no dlopen, no CPU-arch dispatch.
+#
+# Expected toolchain (activated emsdk on PATH):
+# emcc, em++, emar
+#
+# Per-translation-unit rules for the CPU sources that are SHARED with
+# the native build are defined in src/Makefile via the `skbb_cpu_tu`
+# canned recipe (one definition emits both the .o and .wasm.o rules).
+# This file owns:
+# - the WASM compiler vars (used by skbb_cpu_tu)
+# - the WASM-only TU (Eigen linalg backend)
+# - the WASM_OBJS list and libskbb_wasm.a archive
+# - the WASM test infrastructure
+#
+# Invocation: this Makefile fragment is included from src/Makefile and
+# must be run with src/ as the working directory. All test source paths
+# (e.g. `tests/wasm/test_*.cpp`) are relative to src/, and `-I.` from
+# WASM_CXXFLAGS resolves `#include "tests/wasm/..."` against src/.
+
+WASM_REPO_ROOT := $(abspath $(dir $(lastword $(MAKEFILE_LIST)))/../../)
+WASM_EIGEN_INC := $(WASM_REPO_ROOT)/.wasm-cache/eigen/include
+
+# Compiler-independent target flags. The macros override the native Makefile's
+# x86 dispatch (we explicitly don't want x86-v3/v4 object files under WASM
+# regardless of what the host machine looks like) and select the Eigen
+# backend in linalg_backend_eigen.cpp.
+WASM_CXX := em++
+WASM_AR := emar
+WASM_CXXFLAGS := -std=c++17 -O3 -Wall -I. \
+ -I$(WASM_EIGEN_INC) \
+ -DSKBB_WASM=1 \
+ -DSKBB_BLAS_BACKEND_EIGEN=1 \
+ -DNOGPU=1 \
+ -fno-exceptions \
+ -Wno-unknown-pragmas
+
+# Object list for the WASM archive. Stems must match the .wasm.o targets
+# emitted by skbb_cpu_tu in src/Makefile, plus the Eigen-backed linalg
+# implementation defined below.
+WASM_OBJS := \
+ util_rand.wasm.o \
+ skbb_detect_acc.wasm.o \
+ skbb_accapi_cpu.wasm.o \
+ dist_permanova.wasm.o \
+ permanova_cpu.wasm.o \
+ ord_pcoa.wasm.o \
+ ord_linalg_backend_eigen.wasm.o \
+ skbb_extern_util.wasm.o \
+ skbb_extern_distance.wasm.o \
+ skbb_extern_ordination.wasm.o
+
+# WASM-only TU: the Eigen-backed linalg backend has no native equivalent
+# (native uses linalg_backend_lapacke.cpp), so this rule isn't generated
+# by the symmetric skbb_cpu_tu macro.
+ord_linalg_backend_eigen.wasm.o: ordination/linalg_backend_eigen.cpp ordination/linalg_backend.hpp
+ $(WASM_CXX) $(WASM_CXXFLAGS) -c $< -o $@
+
+libskbb_wasm.a: $(WASM_OBJS)
+ rm -f $@
+ $(WASM_AR) rcs $@ $(WASM_OBJS)
+
+# ---- WASM test binaries ----
+# Emscripten links a .wasm + .js pair; node runs the .js, which loads
+# the sibling .wasm. NODERAWFS=0 is the default — we don't touch disk.
+WASM_TEST_LDFLAGS := -sEXIT_RUNTIME=1 -sALLOW_MEMORY_GROWTH=1 \
+ -sENVIRONMENT=node -sNODERAWFS=0
+
+# The public headers refer to themselves as "scikit-bio-binaries/util.h".
+# Provide that prefix via a staged include directory.
+wasm_test_include_stage:
+ mkdir -p .wasm-test-include/scikit-bio-binaries
+ cp extern/util.h .wasm-test-include/scikit-bio-binaries/util.h
+ cp extern/distance.h .wasm-test-include/scikit-bio-binaries/distance.h
+ cp extern/ordination.h .wasm-test-include/scikit-bio-binaries/ordination.h
+
+test_smoke_wasm.js: libskbb_wasm.a tests/wasm/test_smoke.cpp wasm_test_include_stage
+ $(WASM_CXX) $(WASM_CXXFLAGS) -I.wasm-test-include \
+ tests/wasm/test_smoke.cpp libskbb_wasm.a \
+ $(WASM_TEST_LDFLAGS) -o $@
+
+test_permanova_wasm.js: libskbb_wasm.a tests/wasm/test_permanova_wasm.cpp \
+ tests/wasm/permanova_inputs.hpp \
+ tests/wasm/expected/permanova_expected.h \
+ wasm_test_include_stage
+ $(WASM_CXX) $(WASM_CXXFLAGS) -I.wasm-test-include \
+ tests/wasm/test_permanova_wasm.cpp libskbb_wasm.a \
+ $(WASM_TEST_LDFLAGS) -o $@
+
+test_center_wasm.js: libskbb_wasm.a tests/wasm/test_center_wasm.cpp \
+ wasm_test_include_stage
+ $(WASM_CXX) $(WASM_CXXFLAGS) -I.wasm-test-include \
+ tests/wasm/test_center_wasm.cpp libskbb_wasm.a \
+ $(WASM_TEST_LDFLAGS) -o $@
+
+test_pcoa_wasm.js: libskbb_wasm.a tests/wasm/test_pcoa_wasm.cpp \
+ tests/wasm/pcoa_inputs.hpp \
+ tests/wasm/expected/pcoa_expected.h \
+ wasm_test_include_stage
+ $(WASM_CXX) $(WASM_CXXFLAGS) -I.wasm-test-include \
+ tests/wasm/test_pcoa_wasm.cpp libskbb_wasm.a \
+ $(WASM_TEST_LDFLAGS) -o $@
+
+wasm_test: test_smoke_wasm.js test_permanova_wasm.js test_center_wasm.js test_pcoa_wasm.js
+ @echo "--- smoke ---"
+ node test_smoke_wasm.js
+ @echo "--- permanova ---"
+ node test_permanova_wasm.js
+ @echo "--- center ---"
+ node test_center_wasm.js
+ @echo "--- pcoa ---"
+ node test_pcoa_wasm.js
+
+# Native-build helpers that produce the expected-value headers consumed by
+# the WASM tests above. These are not committed to the repo (per reviewer
+# guidance); make wasm_test depends on them through the `expected/*.h`
+# prerequisites and they are regenerated whenever missing or out of date.
+# Both generators require the native build (libskbb_cpu.a) plus BLAS/LAPACKE
+# for the PCoA generator.
+tests/wasm/expected/permanova_expected.h: tests/wasm/generate_permanova_expected.cpp \
+ tests/wasm/permanova_inputs.hpp \
+ libskbb_cpu.a
+ $(CXX) $(CXXFLAGS) \
+ tests/wasm/generate_permanova_expected.cpp libskbb_cpu.a \
+ $(LDFLAGS) -o generate_permanova_expected.exe
+ mkdir -p tests/wasm/expected
+ OMP_NUM_THREADS=1 ./generate_permanova_expected.exe > $@
+
+tests/wasm/expected/pcoa_expected.h: tests/wasm/generate_pcoa_expected.cpp \
+ tests/wasm/pcoa_inputs.hpp \
+ libskbb_cpu.a
+ $(CXX) $(CXXFLAGS) \
+ tests/wasm/generate_pcoa_expected.cpp libskbb_cpu.a \
+ $(LDFLAGS) $(BLASLIB) -o generate_pcoa_expected.exe
+ mkdir -p tests/wasm/expected
+ OMP_NUM_THREADS=1 ./generate_pcoa_expected.exe > $@
+
+wasm_clean:
+ rm -f libskbb_wasm.a *.wasm.o *_wasm.js *_wasm.wasm
+ rm -f generate_permanova_expected.exe generate_pcoa_expected.exe
+ rm -rf .wasm-test-include tests/wasm/expected
+
+.PHONY: wasm_clean wasm_test wasm_test_include_stage
View it on GitLab: https://salsa.debian.org/med-team/scikit-bio-binaries/-/commit/ccba9dc03c036ebcb02e1e2ad48334f5cfedfccf
--
View it on GitLab: https://salsa.debian.org/med-team/scikit-bio-binaries/-/commit/ccba9dc03c036ebcb02e1e2ad48334f5cfedfccf
You're receiving this email because of your account on salsa.debian.org. Manage all notifications: https://salsa.debian.org/-/profile/notifications | Help: https://salsa.debian.org/help
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/debian-med-commit/attachments/20260907/c5921fea/attachment-0001.htm>
More information about the debian-med-commit
mailing list