From 3b88ab28efafd7a65710cced599f71a77da7728c Mon Sep 17 00:00:00 2001 From: xihaian251 Date: Thu, 22 Oct 2026 22:30:24 +0800 Subject: [PATCH 1/6] Make BSP-OT candidate plans depend only on the seed and their index ot.bsp.compute_bspot_bijection documents `seed` as a reproducibility knob, but the BSP-OT engines are `thread_local std::mt19937` objects constructed once per thread from `global_seed` (BSP-OT_header_only.h:549-553), so the `global_seed = seed` assignment in BSPOT_wrap reaches the engine of the first call made on that thread only: a repeated call continued the previous stream and returned different plans, different intermediary stacks and a different cost (Issue #897). The same lines also cost solution quality, because the candidate plans are generated by `#pragma omp parallel for` and each worker built its engine from the same `global_seed`: iterations at equal positions in their chunks consumed identical subsequences and produced identical matchings, so `n_plans=64` yielded ceil(n_plans / n_threads) distinct candidates on this machine and `MergePlansNoPar` walked over the copies. Re-seed the engine at the start of each iteration from a SplitMix64 mix of (user_seed, plan index) and pass the seed down as a parameter, so that candidate plan i is a function of the arguments alone. The mixing must not be additive: with `user_seed + task`, plan t of seed s is plan t-1 of seed s+1, so changing the seed would leave nearly the whole candidate set unchanged. The matching object is now constructed per iteration, which removes the only mutable state shared between the iterations of this loop, and the now dead `global_seed = seed` write is dropped; those two thread-safety edits do not change the result for a given seed. --- ot/bsp/bsp_wrapper.cpp | 47 ++++++++++++++++++++++++++++-------------- 1 file changed, 31 insertions(+), 16 deletions(-) diff --git a/ot/bsp/bsp_wrapper.cpp b/ot/bsp/bsp_wrapper.cpp index 3867d6305..982aa5a9d 100644 --- a/ot/bsp/bsp_wrapper.cpp +++ b/ot/bsp/bsp_wrapper.cpp @@ -2,6 +2,18 @@ #include "BSP-OT_header_only.h" +// Seed of the mt19937 that generates candidate plan `task` of a solve started with +// `user_seed`. The mixing must not be additive: with `user_seed + task`, plan `t` of +// seed `s` is plan `t-1` of seed `s+1`, so changing the seed would leave almost the +// whole candidate set untouched. SplitMix64 finalizer gives unrelated streams. +static std::uint32_t derive_seed(std::uint64_t user_seed, std::uint64_t task) { + std::uint64_t z = user_seed + 0x9E3779B97F4A7C15ull * (task + 1ull); + z = (z ^ (z >> 30)) * 0xBF58476D1CE4E5B9ull; + z = (z ^ (z >> 27)) * 0x94D049BB133111EBull; + return static_cast((z ^ (z >> 31)) & 0xFFFFFFFFull); +} + + template BSPOT::Points UnLinearize(double* data,int n,int d) { return Eigen::Map>(data, d, n).template cast(); @@ -31,14 +43,19 @@ std::function makeCost(const BSPOT::Points& A,const template -std::vector computeBSPMatchings_dim(const BSPOT::Points& A,const BSPOT::Points& B,int nb_plans, bool gaussian = true){ +std::vector computeBSPMatchings_dim(const BSPOT::Points& A,const BSPOT::Points& B,int nb_plans, bool gaussian = true, std::uint64_t seed = 0){ using namespace BSPOT; std::vector plans(nb_plans); - BijectiveBSPMatching BSP(A,B); #pragma omp parallel for for (int i = 0; i < nb_plans; i++) { + // Candidate plan i must depend on (seed, i) alone: the thread_local engine is + // re-seeded here, and the matching object is constructed per iteration so that + // this loop shares no mutable state between iterations. Neither the OpenMP team + // size nor the plans generated by earlier calls can then reach the result. + rng().seed(derive_seed(seed, static_cast(i))); + BijectiveBSPMatching BSP(A,B); if (gaussian){ plans[i] = BSP.computeGaussianMatching(); } @@ -57,10 +74,10 @@ BSPOT::BijectiveMatching MergeBijections(const std::vector -double BSPOT_wrap_dim(int n, int d, double *X, double *Y, uint64_t nb_plans,std::vector& plans, BSPOT::BijectiveMatching& plan,int lp_power,const BSPOT::BijectiveMatching& T0,bool gaussian) { +double BSPOT_wrap_dim(int n, int d, double *X, double *Y, uint64_t nb_plans,std::vector& plans, BSPOT::BijectiveMatching& plan,int lp_power,const BSPOT::BijectiveMatching& T0,bool gaussian,uint64_t seed) { auto A = UnLinearize(X,n,d); auto B = UnLinearize(Y,n,d); - plans = computeBSPMatchings_dim(A,B,nb_plans,gaussian); + plans = computeBSPMatchings_dim(A,B,nb_plans,gaussian,seed); auto cost_func = makeCost(A,B,lp_power); plan = MergeBijections(plans,cost_func,T0); @@ -70,8 +87,6 @@ double BSPOT_wrap_dim(int n, int d, double *X, double *Y, uint64_t nb_plans,std: double BSPOT_wrap(int n, int d, double *X, double *Y, uint64_t nb_plans, int *plans_ptr, int *final_plan_ptr,int lp_power,int* initial_plan,bool gaussian,int seed) { using namespace BSPOT; - global_seed = seed; - BijectiveMatching T0; if (initial_plan){ @@ -87,16 +102,16 @@ double BSPOT_wrap(int n, int d, double *X, double *Y, uint64_t nb_plans, int *pl switch (d) { - case 2: {cost = BSPOT_wrap_dim<2>(n,d,X,Y,nb_plans,plans,plan,lp_power,T0, gaussian);break;} - case 3: {cost = BSPOT_wrap_dim<3>(n,d,X,Y,nb_plans,plans,plan,lp_power,T0, gaussian);break;} - case 4: {cost = BSPOT_wrap_dim<4>(n,d,X,Y,nb_plans,plans,plan,lp_power,T0, gaussian);break;} - case 5: {cost = BSPOT_wrap_dim<5>(n,d,X,Y,nb_plans,plans,plan,lp_power,T0, gaussian);break;} - case 6: {cost = BSPOT_wrap_dim<6>(n,d,X,Y,nb_plans,plans,plan,lp_power,T0, gaussian);break;} - case 7: {cost = BSPOT_wrap_dim<7>(n,d,X,Y,nb_plans,plans,plan,lp_power,T0, gaussian);break;} - case 8: {cost = BSPOT_wrap_dim<8>(n,d,X,Y,nb_plans,plans,plan,lp_power,T0, gaussian);break;} - case 9: {cost = BSPOT_wrap_dim<9>(n,d,X,Y,nb_plans,plans,plan,lp_power,T0, gaussian);break;} - case 10: {cost = BSPOT_wrap_dim<10>(n,d,X,Y,nb_plans,plans,plan,lp_power,T0, gaussian);break;} - default: {cost = BSPOT_wrap_dim<-1>(n,d,X,Y,nb_plans,plans,plan,lp_power,T0, gaussian);break;} + case 2: {cost = BSPOT_wrap_dim<2>(n,d,X,Y,nb_plans,plans,plan,lp_power,T0, gaussian, static_cast(seed));break;} + case 3: {cost = BSPOT_wrap_dim<3>(n,d,X,Y,nb_plans,plans,plan,lp_power,T0, gaussian, static_cast(seed));break;} + case 4: {cost = BSPOT_wrap_dim<4>(n,d,X,Y,nb_plans,plans,plan,lp_power,T0, gaussian, static_cast(seed));break;} + case 5: {cost = BSPOT_wrap_dim<5>(n,d,X,Y,nb_plans,plans,plan,lp_power,T0, gaussian, static_cast(seed));break;} + case 6: {cost = BSPOT_wrap_dim<6>(n,d,X,Y,nb_plans,plans,plan,lp_power,T0, gaussian, static_cast(seed));break;} + case 7: {cost = BSPOT_wrap_dim<7>(n,d,X,Y,nb_plans,plans,plan,lp_power,T0, gaussian, static_cast(seed));break;} + case 8: {cost = BSPOT_wrap_dim<8>(n,d,X,Y,nb_plans,plans,plan,lp_power,T0, gaussian, static_cast(seed));break;} + case 9: {cost = BSPOT_wrap_dim<9>(n,d,X,Y,nb_plans,plans,plan,lp_power,T0, gaussian, static_cast(seed));break;} + case 10: {cost = BSPOT_wrap_dim<10>(n,d,X,Y,nb_plans,plans,plan,lp_power,T0, gaussian, static_cast(seed));break;} + default: {cost = BSPOT_wrap_dim<-1>(n,d,X,Y,nb_plans,plans,plan,lp_power,T0, gaussian, static_cast(seed));break;} } std::copy(plan.getPlan().begin(), plan.getPlan().end(), final_plan_ptr); From 7fc28104472f0c929fbc7d76fc1a40dabbe0f143 Mon Sep 17 00:00:00 2001 From: xihaian251 Date: Thu, 22 Oct 2026 22:30:25 +0800 Subject: [PATCH 2/6] Add tests for BSP-OT seed reproducibility Repeated calls, interleaved history, a fresh-interpreter reference, invariance to OMP_NUM_THREADS, and the property that candidate plan i is determined by (seed, i), which is checked by comparing the first rows of the candidate stacks returned for n_plans=4 and n_plans=37. Five of the six fail on current master. The candidate matchings are deliberately not asserted to be pairwise distinct: independent streams can legitimately return the same bijection -- for a similarity transform every BSP candidate is the identity -- so distinctness is a property of the instance and of the standard library's distributions rather than of the seeding, and #897 reports it as a quality metric. No historical cost or permutation is hardcoded, so the tests do not freeze the new random sequence. The children report the extension hash they loaded and the parent asserts it matches, so a child importing another build cannot silently pass. --- test/test_bspot.py | 176 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 176 insertions(+) diff --git a/test/test_bspot.py b/test/test_bspot.py index 862923e12..de2fc8782 100644 --- a/test/test_bspot.py +++ b/test/test_bspot.py @@ -4,6 +4,13 @@ # # License: MIT License +import hashlib +import json +import os +import subprocess +import sys +import warnings + import ot import ot.bsp import numpy as np @@ -135,3 +142,172 @@ def test_bsp_ot_torch_backend(): # cost should be zero np.testing.assert_allclose(cost_new, 0, atol=1e-5) + + +def _two_clouds(n=100, data_seed=0): + rng = np.random.RandomState(data_seed) + return rng.randn(n, 2), rng.randn(n, 2) + + +def _bsp_fingerprint(A, B, n_plans, seed): + cost, perm, perms = ot.bsp.compute_bspot_bijection(A, B, n_plans, seed=seed) + return float(np.asarray(cost)), np.asarray(perm), np.asarray(perms) + + +def test_bsp_seed_repeated_calls_are_reproducible(): + # ``seed`` is documented as "Random seed for reproducibility", so two calls with + # the same arguments must return the same matchings, not just a similar cost. + A, B = _two_clouds() + cost0, perm0, perms0 = _bsp_fingerprint(A, B, 32, 0) + cost1, perm1, perms1 = _bsp_fingerprint(A, B, 32, 0) + + np.testing.assert_array_equal(perm0, perm1) + np.testing.assert_array_equal(perms0, perms1) + np.testing.assert_allclose(cost0, cost1, rtol=0, atol=0) + + +def test_bsp_seed_is_independent_of_call_history(): + # A preceding call with another seed must not disturb the seeded stream. + A, B = _two_clouds() + _bsp_fingerprint(A, B, 16, 7) + cost0, perm0, perms0 = _bsp_fingerprint(A, B, 32, 0) + _bsp_fingerprint(A, B, 16, 11) + cost1, perm1, perms1 = _bsp_fingerprint(A, B, 32, 0) + + np.testing.assert_array_equal(perm0, perm1) + np.testing.assert_array_equal(perms0, perms1) + np.testing.assert_allclose(cost0, cost1, rtol=0, atol=0) + + +_FRESH_PROCESS_SNIPPET = """ +import hashlib, json, sys +import numpy as np +import ot +from ot.bsp import bsp_wrap + +seed, n_plans = int(sys.argv[1]), int(sys.argv[2]) +rng = np.random.RandomState(0) +A, B = rng.randn(100, 2), rng.randn(100, 2) +out = [] +for _ in range(2): + cost, perm, perms = ot.bsp.compute_bspot_bijection(A, B, n_plans, seed=seed) + perm, perms = np.asarray(perm), np.asarray(perms) + out.append({ + "cost": float(np.asarray(cost)), + "perm": hashlib.sha256(perm.tobytes()).hexdigest(), + "perms": hashlib.sha256(perms.tobytes()).hexdigest(), + }) +print(json.dumps({"results": out, "ot_file": ot.__file__, + "ext_sha256": hashlib.sha256(open(bsp_wrap.__file__, "rb").read()).hexdigest()})) +""" + + +def _ext_sha256(): + from ot.bsp import bsp_wrap + + return hashlib.sha256(open(bsp_wrap.__file__, "rb").read()).hexdigest() + + +def _bsp_in_fresh_process(seed, n_plans, omp_threads): + """Run the fingerprint in a fresh interpreter. + + Returns the reported results, or None together with the raw process output when the + child did not get as far as printing them. + """ + root = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + env = dict(os.environ) + env["PYTHONPATH"] = root + os.pathsep + env.get("PYTHONPATH", "") + if omp_threads is not None: + env["OMP_NUM_THREADS"] = str(omp_threads) + proc = subprocess.run( + [sys.executable, "-c", _FRESH_PROCESS_SNIPPET, str(seed), str(n_plans)], + cwd=root, + env=env, + text=True, + capture_output=True, + check=False, + ) + if proc.returncode != 0: + return None, proc.returncode, proc.stderr + report = json.loads(proc.stdout) + # a child that imported a different POT build would prove nothing + assert report["ext_sha256"] == _ext_sha256(), "child loaded another build" + return report["results"], proc.returncode, proc.stderr + + +def _run_in_fresh_process(seed=0, n_plans=32, omp_threads=None, attempts=3): + # A child interpreter can die while importing POT on a host that is out of memory; + # that shows up as a plain exit code 1 with a Python traceback and is retried. + # Anything else -- a native crash such as 0xC0000005 or a signal -- is reported at + # once, since a real defect is deterministic and retrying it would only hide it. + failure = None + for _ in range(attempts): + results, returncode, stderr = _bsp_in_fresh_process(seed, n_plans, omp_threads) + if results is not None: + return results + failure = (returncode, stderr) + if returncode != 1: + break + raise AssertionError( + f"BSP subprocess failed with rc={failure[0]}\n{failure[1]}" + ) + + +def test_bsp_seed_reproducible_in_fresh_process(): + # Reproducibility must not depend on whether the solver ran earlier in the + # session, so the reference is a first call in a brand new interpreter. + ref = _run_in_fresh_process() + assert ref[0] == ref[1], "two first calls in one fresh interpreter already disagree" + + A, B = _two_clouds() + _bsp_fingerprint(A, B, 16, 7) # make the current session non-pristine + cost, perm, perms = _bsp_fingerprint(A, B, 32, 0) + assert cost == ref[0]["cost"] + assert hashlib.sha256(np.asarray(perm).tobytes()).hexdigest() == ref[0]["perm"] + assert hashlib.sha256(np.asarray(perms).tobytes()).hexdigest() == ref[0]["perms"] + + +def test_bsp_result_independent_of_openmp_thread_count(): + # The candidate plans are generated by an OpenMP loop. How the loop happens to be + # split over threads is not part of the solver, so the answer must not change with + # the number of available threads (CI images and laptops differ here). Staying + # parallel is a performance property, checked by the benchmarks, not by this test. + ref = _run_in_fresh_process(omp_threads=1) + for threads in (2, 4): + got = _run_in_fresh_process(omp_threads=threads) + assert got == ref, f"OMP_NUM_THREADS={threads} gives a different result" + + +def test_bsp_candidate_plans_are_determined_by_seed_and_index(): + # The property the seeding gives: candidate plan i depends on + # (A, B, p, gaussian_slicing, seed, i) alone, so asking for 4 plans and asking for 37 + # must agree on the first four rows, whatever the OpenMP chunking made of them and + # whatever was computed before. Whether the plans are pairwise distinct is a property + # of the instance and of the standard library's distributions, not of the seeding -- + # for a similarity transform every BSP candidate legitimately is the identity -- so it + # is not asserted here. + A, B = _two_clouds() + _, _, few = _bsp_fingerprint(A, B, 4, 0) + _bsp_fingerprint(A, B, 16, 5) # unrelated call in between + _, _, many = _bsp_fingerprint(A, B, 37, 0) + np.testing.assert_array_equal(few, many[:4]) + _, _, few_again = _bsp_fingerprint(A, B, 4, 0) + np.testing.assert_array_equal(few, few_again) + + +def test_bsp_permutation_structure_and_cost(): + # Structural invariants that must hold whatever random sequence is drawn. + # compute_bspot_bijection warns when the cost it recomputes from ``perm`` disagrees + # with the cost the solver reports for it, so warnings are turned into errors here: + # that compares the solver's own bookkeeping against the permutation it returned. + A, B = _two_clouds(n=80, data_seed=3) + M = ot.dist(A, B, metric="sqeuclidean") + with warnings.catch_warnings(): + warnings.simplefilter("error") + for seed in (0, 1, 2): + cost, perm, perms = _bsp_fingerprint(A, B, 16, seed) + perm = np.asarray(perm) + np.testing.assert_array_equal(np.sort(perm), np.arange(80)) + np.testing.assert_allclose(cost, np.mean(M[np.arange(80), perm]), rtol=1e-9) + for row in np.asarray(perms): + np.testing.assert_array_equal(np.sort(row), np.arange(80)) From 3173c824e6db0936dc651ca76223ae08f3495391 Mon Sep 17 00:00:00 2001 From: xihaian251 Date: Thu, 22 Oct 2026 22:30:25 +0800 Subject: [PATCH 3/6] Note the BSP-OT seeding fix in RELEASES.md Required by .github/workflows/check_release.yml, which fails any pull request that does not touch RELEASES.md. --- RELEASES.md | 1 + 1 file changed, 1 insertion(+) diff --git a/RELEASES.md b/RELEASES.md index d01ff927f..5aa4c8265 100644 --- a/RELEASES.md +++ b/RELEASES.md @@ -18,6 +18,7 @@ #### Closed issues +- Fix the `seed` of `ot.bsp.compute_bspot_bijection` (and `random_state` of `ot.solve_sample(..., method="bsp")`) not being reproducible: the BSP-OT candidate matchings were generated from a `thread_local` `std::mt19937` that is seeded once per OpenMP thread, so a repeated call continued the previous stream and returned different plans, costs and merged bijections, and every thread of a team started from the same state, which collapsed the `n_plans` candidates down to `ceil(n_plans / n_threads)` distinct matchings. Each candidate plan now draws its own stream, derived from the seed and the plan index, so the result is reproducible and independent of the call history and of the number of OpenMP threads (Issue #897) - Fix duplicate function and class descriptions in generated module API pages while preserving top-level aliases and example backreferences (Issue #527) - Fix internal arrays created on the torch default device instead of the input device in `ot.partial.partial_wasserstein`, `ot.wasserstein_circle` (and `ot.sliced_wasserstein_sphere`), the high-dimensional Bures-Wasserstein functions in `ot.gaussian`, `ot.gmm`, `ot.lowrank_sinkhorn` with `init="deterministic"`, the semi-relaxed (F)GW barycenters and `ot.utils.projection_sparse_simplex`, which made them fail on GPU inputs. Add a `meta_default_device` test fixture to catch such allocations on CPU (PR #883, Issue #852) - Allow `NumpyBackend.seed` to adopt an existing `np.random.RandomState` instance and remove NumPy-specific random sampling paths in sliced utilities (PR #849, Issue #848) From 577bc95bebe832aab4eab0ba92024db09cdd823e Mon Sep 17 00:00:00 2001 From: xihaian251 Date: Wed, 28 Oct 2026 11:53:30 +0800 Subject: [PATCH 4/6] test: handle noisy subprocess stdout in BSP reproducibility tests --- test/test_bspot.py | 79 +++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 74 insertions(+), 5 deletions(-) diff --git a/test/test_bspot.py b/test/test_bspot.py index de2fc8782..8c92ce61e 100644 --- a/test/test_bspot.py +++ b/test/test_bspot.py @@ -179,13 +179,17 @@ def test_bsp_seed_is_independent_of_call_history(): np.testing.assert_allclose(cost0, cost1, rtol=0, atol=0) +# Third-party imports write to stdout on import (e.g. the KeOps warning printed on +# machines without CUDA), so the child prefixes its record with a unique marker. +_RESULT_MARKER = "__POT_BSP_JSON_V1__=" + _FRESH_PROCESS_SNIPPET = """ import hashlib, json, sys import numpy as np import ot from ot.bsp import bsp_wrap -seed, n_plans = int(sys.argv[1]), int(sys.argv[2]) +seed, n_plans, marker = int(sys.argv[1]), int(sys.argv[2]), sys.argv[3] rng = np.random.RandomState(0) A, B = rng.randn(100, 2), rng.randn(100, 2) out = [] @@ -197,8 +201,9 @@ def test_bsp_seed_is_independent_of_call_history(): "perm": hashlib.sha256(perm.tobytes()).hexdigest(), "perms": hashlib.sha256(perms.tobytes()).hexdigest(), }) -print(json.dumps({"results": out, "ot_file": ot.__file__, - "ext_sha256": hashlib.sha256(open(bsp_wrap.__file__, "rb").read()).hexdigest()})) +print(marker + json.dumps({"results": out, "ot_file": ot.__file__, + "ext_sha256": hashlib.sha256(open(bsp_wrap.__file__, "rb").read()).hexdigest()}), + flush=True) """ @@ -208,6 +213,37 @@ def _ext_sha256(): return hashlib.sha256(open(bsp_wrap.__file__, "rb").read()).hexdigest() +def _parse_fresh_process_report(stdout, stderr, returncode): + """Return the JSON record printed by the child after ``_RESULT_MARKER``. + + The marker is searched for anywhere in stdout: the child shares stdout with + third-party imports, so its record is neither first nor necessarily last. + """ + context = ( + f"returncode={returncode}\n" + f"--- child stdout ---\n{stdout!r}\n" + f"--- child stderr ---\n{stderr!r}" + ) + occurrences = stdout.count(_RESULT_MARKER) + if occurrences == 0: + raise AssertionError( + f"child stdout contains no {_RESULT_MARKER!r} record\n{context}" + ) + if occurrences > 1: + raise AssertionError( + f"child stdout contains {occurrences} {_RESULT_MARKER!r} records\n{context}" + ) + start = stdout.index(_RESULT_MARKER) + len(_RESULT_MARKER) + end = stdout.find("\n", start) + payload = stdout[start:] if end == -1 else stdout[start:end] + try: + return json.loads(payload) + except json.JSONDecodeError as exc: + raise AssertionError( + f"invalid JSON in {_RESULT_MARKER!r} record: {payload!r}\n{context}" + ) from exc + + def _bsp_in_fresh_process(seed, n_plans, omp_threads): """Run the fingerprint in a fresh interpreter. @@ -220,7 +256,14 @@ def _bsp_in_fresh_process(seed, n_plans, omp_threads): if omp_threads is not None: env["OMP_NUM_THREADS"] = str(omp_threads) proc = subprocess.run( - [sys.executable, "-c", _FRESH_PROCESS_SNIPPET, str(seed), str(n_plans)], + [ + sys.executable, + "-c", + _FRESH_PROCESS_SNIPPET, + str(seed), + str(n_plans), + _RESULT_MARKER, + ], cwd=root, env=env, text=True, @@ -229,7 +272,7 @@ def _bsp_in_fresh_process(seed, n_plans, omp_threads): ) if proc.returncode != 0: return None, proc.returncode, proc.stderr - report = json.loads(proc.stdout) + report = _parse_fresh_process_report(proc.stdout, proc.stderr, proc.returncode) # a child that imported a different POT build would prove nothing assert report["ext_sha256"] == _ext_sha256(), "child loaded another build" return report["results"], proc.returncode, proc.stderr @@ -253,6 +296,32 @@ def _run_in_fresh_process(seed=0, n_plans=32, omp_threads=None, attempts=3): ) +def test_fresh_process_report_parsing_ignores_foreign_stdout(): + report = {"results": [{"cost": 0.5}], "ot_file": "ot.py", "ext_sha256": "abc"} + stdout = ( + "[KeOps] Warning : CUDA libraries not found or could not be loaded; " + "Switching to CPU only." + + _RESULT_MARKER + + json.dumps(report) + + "\ntrailing log line\n" + ) + assert _parse_fresh_process_report(stdout, "", 0) == report + + +def test_fresh_process_report_parsing_diagnoses_protocol_errors(): + record = _RESULT_MARKER + '{"results": []}' + + with pytest.raises(AssertionError, match="contains no") as missing: + _parse_fresh_process_report("only foreign output\n", "boom", 0) + assert "boom" in str(missing.value) + + with pytest.raises(AssertionError, match="contains 2"): + _parse_fresh_process_report(record + "\n" + record + "\n", "", 0) + + with pytest.raises(AssertionError, match="invalid JSON"): + _parse_fresh_process_report(_RESULT_MARKER + "{not json}\n", "", 0) + + def test_bsp_seed_reproducible_in_fresh_process(): # Reproducibility must not depend on whether the solver ran earlier in the # session, so the reference is a first call in a brand new interpreter. From fd4706b30e3c50900ad91ce6ebd2a371e350c657 Mon Sep 17 00:00:00 2001 From: xihaian251 Date: Tue, 3 Nov 2026 11:21:57 +0800 Subject: [PATCH 5/6] Document the BSP-OT seed contract and trim the review surface State in the seed docstring that candidate plans are reproducible for a given seed within a fixed build, and that a loop refining a solution with initial_perm must change the seed each round, since reusing it repeats the same candidates and stalls refinement. Shorten the RELEASES.md entry to match its neighbours and credit PR #898, and drop repeated explanations from the seeding tests. --- RELEASES.md | 2 +- ot/bsp/bsp_ot.py | 6 +++++- test/test_bspot.py | 31 ++++++++++++++----------------- 3 files changed, 20 insertions(+), 19 deletions(-) diff --git a/RELEASES.md b/RELEASES.md index 2ae028ab3..b4e356f79 100644 --- a/RELEASES.md +++ b/RELEASES.md @@ -19,7 +19,7 @@ #### Closed issues - Fix overflow to NaN in `ot.stochastic.coordinate_grad_semi_dual` by shifting before exponentiating, as `ot.stochastic.c_transform_entropic` already does. The factor cancels in the normalisation, so results are unchanged where the old code did not overflow (PR #866, Issue #264) -- Fix the `seed` of `ot.bsp.compute_bspot_bijection` (and `random_state` of `ot.solve_sample(..., method="bsp")`) not being reproducible: the BSP-OT candidate matchings were generated from a `thread_local` `std::mt19937` that is seeded once per OpenMP thread, so a repeated call continued the previous stream and returned different plans, costs and merged bijections, and every thread of a team started from the same state, which collapsed the `n_plans` candidates down to `ceil(n_plans / n_threads)` distinct matchings. Each candidate plan now draws its own stream, derived from the seed and the plan index, so the result is reproducible and independent of the call history and of the number of OpenMP threads (Issue #897) +- Fix `ot.bsp.compute_bspot_bijection` (and `random_state` of `ot.solve_sample(..., method="bsp")`) not being reproducible for a given `seed`: its candidate plans came from a `thread_local` `std::mt19937` seeded once per OpenMP thread, so the result depended on the call history and on the number of threads. Each plan now draws its own stream, derived from the seed and the plan index (PR #898, Issue #897) - Fix duplicate function and class descriptions in generated module API pages while preserving top-level aliases and example backreferences (Issue #527) - Fix internal arrays created on the torch default device instead of the input device in `ot.partial.partial_wasserstein`, `ot.wasserstein_circle` (and `ot.sliced_wasserstein_sphere`), the high-dimensional Bures-Wasserstein functions in `ot.gaussian`, `ot.gmm`, `ot.lowrank_sinkhorn` with `init="deterministic"`, the semi-relaxed (F)GW barycenters and `ot.utils.projection_sparse_simplex`, which made them fail on GPU inputs. Add a `meta_default_device` test fixture to catch such allocations on CPU (PR #883, Issue #852) - Allow `NumpyBackend.seed` to adopt an existing `np.random.RandomState` instance and remove NumPy-specific random sampling paths in sliced utilities (PR #849, Issue #848) diff --git a/ot/bsp/bsp_ot.py b/ot/bsp/bsp_ot.py index d7d13a7b5..cca23dfef 100644 --- a/ot/bsp/bsp_ot.py +++ b/ot/bsp/bsp_ot.py @@ -42,7 +42,11 @@ def compute_bspot_bijection( Comes with a cubic complexity with dimension. If 'auto', then the heuristic is used for dimensions smaller than 64, and disabled for larger dimensions. seed : int, optional - Random seed for reproducibility (default 0). + Random seed for reproducibility (default 0). Candidate plans are + reproducible for a given seed within a fixed build, so a loop that + refines a solution with ``initial_perm`` should use a different seed + at each round: reusing the same seed regenerates the same candidate + plans and can stall the refinement. Returns ------- diff --git a/test/test_bspot.py b/test/test_bspot.py index 8c92ce61e..e366d18a8 100644 --- a/test/test_bspot.py +++ b/test/test_bspot.py @@ -155,8 +155,8 @@ def _bsp_fingerprint(A, B, n_plans, seed): def test_bsp_seed_repeated_calls_are_reproducible(): - # ``seed`` is documented as "Random seed for reproducibility", so two calls with - # the same arguments must return the same matchings, not just a similar cost. + # Two calls with the same arguments must return the same matchings, not + # just a similar cost. A, B = _two_clouds() cost0, perm0, perms0 = _bsp_fingerprint(A, B, 32, 0) cost1, perm1, perms1 = _bsp_fingerprint(A, B, 32, 0) @@ -216,8 +216,8 @@ def _ext_sha256(): def _parse_fresh_process_report(stdout, stderr, returncode): """Return the JSON record printed by the child after ``_RESULT_MARKER``. - The marker is searched for anywhere in stdout: the child shares stdout with - third-party imports, so its record is neither first nor necessarily last. + The marker is searched for anywhere in stdout, since the child's record is + neither first nor necessarily last. """ context = ( f"returncode={returncode}\n" @@ -337,10 +337,9 @@ def test_bsp_seed_reproducible_in_fresh_process(): def test_bsp_result_independent_of_openmp_thread_count(): - # The candidate plans are generated by an OpenMP loop. How the loop happens to be - # split over threads is not part of the solver, so the answer must not change with - # the number of available threads (CI images and laptops differ here). Staying - # parallel is a performance property, checked by the benchmarks, not by this test. + # How the candidate-plan loop happens to be split over threads is not part of the + # solver, so the answer must not change with the number of available threads. + # Staying parallel is a performance property, not asserted here. ref = _run_in_fresh_process(omp_threads=1) for threads in (2, 4): got = _run_in_fresh_process(omp_threads=threads) @@ -350,11 +349,10 @@ def test_bsp_result_independent_of_openmp_thread_count(): def test_bsp_candidate_plans_are_determined_by_seed_and_index(): # The property the seeding gives: candidate plan i depends on # (A, B, p, gaussian_slicing, seed, i) alone, so asking for 4 plans and asking for 37 - # must agree on the first four rows, whatever the OpenMP chunking made of them and - # whatever was computed before. Whether the plans are pairwise distinct is a property - # of the instance and of the standard library's distributions, not of the seeding -- - # for a similarity transform every BSP candidate legitimately is the identity -- so it - # is not asserted here. + # must agree on the first four rows. Whether the plans are pairwise distinct is a + # property of the instance and of the standard library's distributions, not of the + # seeding -- for a similarity transform every candidate legitimately is the identity + # -- so it is not asserted here. A, B = _two_clouds() _, _, few = _bsp_fingerprint(A, B, 4, 0) _bsp_fingerprint(A, B, 16, 5) # unrelated call in between @@ -365,10 +363,9 @@ def test_bsp_candidate_plans_are_determined_by_seed_and_index(): def test_bsp_permutation_structure_and_cost(): - # Structural invariants that must hold whatever random sequence is drawn. - # compute_bspot_bijection warns when the cost it recomputes from ``perm`` disagrees - # with the cost the solver reports for it, so warnings are turned into errors here: - # that compares the solver's own bookkeeping against the permutation it returned. + # Invariants that must hold whatever random sequence is drawn. The solver warns when + # the cost it recomputes from ``perm`` disagrees with the cost it reports for it, so + # warnings are turned into errors here. A, B = _two_clouds(n=80, data_seed=3) M = ot.dist(A, B, metric="sqeuclidean") with warnings.catch_warnings(): From df935ee9339cb5bc9ecadee54f88a718fd86a3ea Mon Sep 17 00:00:00 2001 From: B Date: Tue, 3 Nov 2026 17:30:45 +0800 Subject: [PATCH 6/6] docs: clarify BSP-OT multithread performance regression --- RELEASES.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/RELEASES.md b/RELEASES.md index b4e356f79..9f3d4a13e 100644 --- a/RELEASES.md +++ b/RELEASES.md @@ -19,7 +19,7 @@ #### Closed issues - Fix overflow to NaN in `ot.stochastic.coordinate_grad_semi_dual` by shifting before exponentiating, as `ot.stochastic.c_transform_entropic` already does. The factor cancels in the normalisation, so results are unchanged where the old code did not overflow (PR #866, Issue #264) -- Fix `ot.bsp.compute_bspot_bijection` (and `random_state` of `ot.solve_sample(..., method="bsp")`) not being reproducible for a given `seed`: its candidate plans came from a `thread_local` `std::mt19937` seeded once per OpenMP thread, so the result depended on the call history and on the number of threads. Each plan now draws its own stream, derived from the seed and the plan index (PR #898, Issue #897) +- Fix `ot.bsp.compute_bspot_bijection` (and `random_state` of `ot.solve_sample(..., method="bsp")`) not being reproducible for a given `seed`: its candidate plans came from a `thread_local` `std::mt19937` seeded once per OpenMP thread, so the result depended on the call history and on the number of threads. Every thread started from the same state, so the previous version could underperform with multiple threads. Each plan now draws its own stream, derived from the seed and the plan index (PR #898, Issue #897) - Fix duplicate function and class descriptions in generated module API pages while preserving top-level aliases and example backreferences (Issue #527) - Fix internal arrays created on the torch default device instead of the input device in `ot.partial.partial_wasserstein`, `ot.wasserstein_circle` (and `ot.sliced_wasserstein_sphere`), the high-dimensional Bures-Wasserstein functions in `ot.gaussian`, `ot.gmm`, `ot.lowrank_sinkhorn` with `init="deterministic"`, the semi-relaxed (F)GW barycenters and `ot.utils.projection_sparse_simplex`, which made them fail on GPU inputs. Add a `meta_default_device` test fixture to catch such allocations on CPU (PR #883, Issue #852) - Allow `NumpyBackend.seed` to adopt an existing `np.random.RandomState` instance and remove NumPy-specific random sampling paths in sliced utilities (PR #849, Issue #848)