diff --git a/.github/workflows/build_and_test.yml b/.github/workflows/build_and_test.yml index 6ce1c1100..04c2f9ae7 100644 --- a/.github/workflows/build_and_test.yml +++ b/.github/workflows/build_and_test.yml @@ -232,12 +232,12 @@ jobs: CCACHE_DIR: ${{ github.workspace }}/.ccache run: ccache -s || true - # Run the suite under BOTH thread settings. OMP_NUM_THREADS overrides the requested per-solve - # thread count (parallel::EffectiveThreadCount), so the serial and threaded code paths genuinely - # differ; running both on every platform catches threading bugs (e.g. clone-per-thread) that a - # single setting would hide. `--parallel` is ctest's *process* parallelism (independent of the - # per-solve OMP count), so each leg still runs the test executables concurrently. - - name: Test (threaded + OMP_NUM_THREADS=1) + # Run the suite under BOTH thread settings. BERTINI_NUM_THREADS overrides the requested + # per-solve thread count (parallel::EffectiveThreadCount), so the serial and threaded code paths + # genuinely differ; running both on every platform catches threading bugs (e.g. clone-per-thread) + # that a single setting would hide. `--parallel` is ctest's *process* parallelism (independent + # of the per-solve thread count), so each leg still runs the test executables concurrently. + - name: Test (threaded + BERTINI_NUM_THREADS=1) working-directory: build env: LD_LIBRARY_PATH: /tmp/boost-base/lib:/tmp/boost-base/lib64 @@ -245,24 +245,28 @@ jobs: run: | set -e jobs=$(nproc 2>/dev/null || sysctl -n hw.ncpu) - echo "::group::ctest with OMP_NUM_THREADS unset (per-solve threading on)" - env -u OMP_NUM_THREADS ctest --output-on-failure --timeout 600 --parallel "$jobs" + echo "::group::ctest with BERTINI_NUM_THREADS unset (per-solve threading on)" + env -u BERTINI_NUM_THREADS ctest --output-on-failure --timeout 600 --parallel "$jobs" echo "::endgroup::" - echo "::group::ctest with OMP_NUM_THREADS=1 (per-solve serial)" - OMP_NUM_THREADS=1 ctest --output-on-failure --timeout 600 --parallel "$jobs" + echo "::group::ctest with BERTINI_NUM_THREADS=1 (per-solve serial)" + BERTINI_NUM_THREADS=1 ctest --output-on-failure --timeout 600 --parallel "$jobs" echo "::endgroup::" # The published wheels ship WITHOUT MPI (we can't know the user's MPI), but this C++ build is # MPI-enabled (openmpi is installed above), so we verify the MPI runtime path actually works # on each OS we support -- and, crucially, that it gives the SAME answer as the serial solve. - # The blackbox CLI reads a file named `input` from the cwd and writes the solution count as - # the first line of `main_data`. mpirun -n 2 = one manager + one worker (the minimal pool). + # The blackbox CLI reads a file named `input` from the cwd and writes its results to + # `main_data`. mpirun -n 2 = one manager + one worker (the minimal pool). # # Run BOTH start-system families, because DistributeSystems broadcasts the start system to the # workers and a dropped serialized member is invisible until that exact start system is sent. # small.b2 (one affine group) infers total-degree; small_mhom.b2 (two groups) infers the # multihomogeneous start system -- whose serialize() was once incomplete, crashing every worker. - - name: MPI smoke test (CLI under mpirun) + # + # With a fixed randomseed, a path's result does not depend on the thread or rank that ran it + # (ADR-0071), so main_data must be BYTE-IDENTICAL across: serial, threaded without MPI, and + # MPI with one and with two threads per worker -- all set through BERTINI_NUM_THREADS. + - name: MPI and thread identity test (CLI) working-directory: build/core env: LD_LIBRARY_PATH: /tmp/boost-base/lib:/tmp/boost-base/lib64 @@ -270,14 +274,18 @@ jobs: run: | set -e for case in small small_mhom; do - cp "$GITHUB_WORKSPACE/benchmark/inputs/${case}.b2" input - ./bertini2 >/dev/null 2>&1 - serial=$(head -n1 main_data); rm -f main_data - mpirun -n 2 --bind-to none ./bertini2 >/dev/null 2>&1 - parallel=$(head -n1 main_data); rm -f main_data - echo "${case}: solution count serial=$serial mpi(-n 2)=$parallel" - test -n "$serial" - test "$serial" = "$parallel" # MPI must match serial, not merely run + { printf 'CONFIG\nrandomseed: 20261005;\n'; tail -n +2 "$GITHUB_WORKSPACE/benchmark/inputs/${case}.b2"; } > input + BERTINI_NUM_THREADS=1 ./bertini2 >/dev/null 2>&1; mv main_data serial + BERTINI_NUM_THREADS=3 ./bertini2 >/dev/null 2>&1; mv main_data threaded + BERTINI_NUM_THREADS=1 mpirun -n 2 --bind-to none ./bertini2 >/dev/null 2>&1; mv main_data mpi_1 + BERTINI_NUM_THREADS=2 mpirun -n 2 --bind-to none ./bertini2 >/dev/null 2>&1; mv main_data mpi_2 + echo "${case}: solution count $(head -n1 serial)" + test -n "$(head -n1 serial)" + for other in threaded mpi_1 mpi_2; do + cmp serial "$other" # identical, not merely the same count + echo " ${other}: identical to serial" + done + rm -f serial threaded mpi_1 mpi_2 done test_cpp_windows: @@ -360,17 +368,17 @@ jobs: ccache --show-stats # Run under both per-solve thread settings (see the unix C++ test job for the rationale). - - name: Test (threaded + OMP_NUM_THREADS=1) + - name: Test (threaded + BERTINI_NUM_THREADS=1) shell: pwsh working-directory: build run: | - Remove-Item Env:\OMP_NUM_THREADS -ErrorAction SilentlyContinue - Write-Host "::group::ctest with OMP_NUM_THREADS unset (per-solve threading on)" + Remove-Item Env:\BERTINI_NUM_THREADS -ErrorAction SilentlyContinue + Write-Host "::group::ctest with BERTINI_NUM_THREADS unset (per-solve threading on)" ctest --output-on-failure --timeout 600 --parallel if ($LASTEXITCODE -ne 0) { exit $LASTEXITCODE } Write-Host "::endgroup::" - $env:OMP_NUM_THREADS = '1' - Write-Host "::group::ctest with OMP_NUM_THREADS=1 (per-solve serial)" + $env:BERTINI_NUM_THREADS = '1' + Write-Host "::group::ctest with BERTINI_NUM_THREADS=1 (per-solve serial)" ctest --output-on-failure --timeout 600 --parallel if ($LASTEXITCODE -ne 0) { exit $LASTEXITCODE } Write-Host "::endgroup::" diff --git a/CHANGELOG.md b/CHANGELOG.md index ad173cdd0..94acf44dd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -339,6 +339,30 @@ the identity of the algorithm that answered it (#420). ### Changed +- **Seed 0 is a seed** (#475, ADR-0072). `bertini.set_random_seed(0)` and `SetGlobalSeed(0)` + meant "draw a seed from entropy", so a script that set seed 0 to make a run reproducible got + a different run every time. 0 is now an ordinary seed, as in numpy and Python's `random`. + Entropy is asked for by omitting the seed: `set_random_seed()` or `set_random_seed(None)` in + Python, `SetGlobalSeedFromEntropy()` in C++. `set_random_seed` now returns the seed in + effect, so the run can be reproduced from it. The classic input format keeps Bertini 1's + meaning: `randomseed: 0`, its default, still draws from entropy in the command-line program. + Seed 0 is therefore the one value that does not carry between a Python script and a classic + input file. +- **The thread count is `BERTINI_NUM_THREADS`, and it defaults to every available CPU.** + b2 used to read `OMP_NUM_THREADS`, although it uses no OpenMP: its threads are its own + `std::thread` pool. That variable also sets the thread count of numpy's OpenBLAS, so one + setting tied the two libraries together, and pinning b2 to one thread for a reproducible run + pinned numpy's linear algebra as well. b2 now reads `BERTINI_NUM_THREADS` and ignores + `OMP_NUM_THREADS`. The order is the variable if set, then the solver's `num_threads`, then + the number of CPUs the process may run on, which on Linux respects the affinity mask that + `taskset`, cpusets and a launcher's binding set. That last default is OpenMP's, so threaded + runs come without configuration. + **MPI ranks change the most.** A rank used one thread unless `OMP_NUM_THREADS` was set; it + now follows the same rule as a standalone solve, `num_threads` included. How many CPUs a rank + may run on depends on how the launcher binds it, so several ranks on one machine can together + start more threads than the machine has cores. Set `BERTINI_NUM_THREADS` whenever you run + more than one rank per machine. A script or job file that set `OMP_NUM_THREADS` for b2 + should set `BERTINI_NUM_THREADS` instead. - **An accuracy estimate names its coordinates** (ADR-0069). A solution's metadata has two accuracy estimates, both the distance between the endgame's last two approximations of the root. `accuracy_estimate` is renamed **`accuracy_estimate_internal_coords`**: it is in the @@ -460,6 +484,34 @@ the identity of the algorithm that answered it (#420). ### Fixed +- **A threaded solve gives the same bits as a serial one, and a path no longer depends on the + paths before it** (#378, ADR-0071). Each path's random draws were already deterministic, + but state that is not random flowed from one path into the next on a reused tracker or + endgame. Under a thread pool the scheduler chose each thread's predecessors, so a seeded + threaded solve differed from run to run in step counts, precision and the last digits of + its endpoints, and algorithms built on those endpoints differed with them. A serial solve + depended on path order. There were five sources: + - an adaptive track built its first step at the precision the previous track ended in; + - a drop to double left the step size and the time at the old multiple precision; + - the condition-number probe was rounded in place at each change of precision; + - the endgame's c/k probe was drawn once per endgame object, by whichever path first needed + it, or redrawn mid-path after an escalation; + - an endgame run in double left an earlier run's multiprecision times in place, so + `final_time_used` reported another path's last time. + Now a solve at any thread count, under MPI or not, is bit-identical path by path to a + serial solve with the same seed: every endpoint and every metadata field but the + wall-clock time. A tracker or endgame driven by hand gives a track the same result as a + fresh one would. Results differ from 4.0.0.dev1 in the last bits and, occasionally, in step + counts. +- **`max_precision_used` is the highest precision a path used, over all of its tracking** + (ADR-0071). The endgame tracks a path in many calls, and the solver's record started afresh + with every one, so a precision raised in an earlier endgame sub-track and lowered before the + last went unreported: seeded cyclic-5 had paths that reported 30 digits after tracking at 40. + The tracker's `FirstPrecisionRecorder` and `MinMaxPrecisionRecorder` now record each track + from its very start, the start-point refinement included, and keep recording every track + once attached. `FirstPrecisionRecorder` used to unsubscribe at its first increase and + report that track ever after, and a track whose initialization failed was reported with the + previous track's values. A new `PathPrecisionRecorder` records a whole path. - **A copied or unpickled system combines with its original** (ADR-0068). `copy.deepcopy`, `copy.copy` and `pickle` go through an archive, and loading from one built variables outside the variable factory: a copy was over fresh variables with the same names, and diff --git a/CLAUDE.md b/CLAUDE.md index 1b9c969df..f4a0ab5be 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -72,6 +72,18 @@ ctest --test-dir build/core ./build/core/test_settings ``` +**Threads.** b2 reads its thread count from `BERTINI_NUM_THREADS` alone: the variable if set, +else the solver's `num_threads`, else every CPU the process may run on (ADR-0070). It does *not* +read `OMP_NUM_THREADS`, which is numpy's OpenBLAS setting. CI runs the C++ suite with the variable +unset and again with `BERTINI_NUM_THREADS=1`; do the same locally (the bert2 MCP's test tools take +`threads: "both"`). Under `mpirun`, several ranks on one machine each default to all the CPUs they +may run on, so set the variable there. + +**Seeds.** To `set_random_seed` / `SetGlobalSeed`, 0 is an ordinary seed, as in numpy; +`set_random_seed()` (or `None`) draws one from entropy and returns it. A classic input file's +`randomseed: 0` keeps Bertini 1's meaning, "draw from entropy" (ADR-0072). A reproducible test sets +a seed explicitly. + ### Python Tests `pytest` is the single way to run the Python tests (the suites are plain pytest @@ -123,7 +135,7 @@ The project has three layers, built in order: 1. **`core/`** -- C++ shared library (`libbertini2`). Header-only-heavy design under `core/include/bertini2/`. Key subsystems: - `function_tree/` -- Expression tree (nodes, operators, symbols) for *building and representing* polynomial systems. Nodes no longer evaluate: node-level recursive evaluation was removed -- the **SLP (`straight_line_program`) is the sole evaluator** (compile a system once, evaluate the compiled program). `Function`/`Handle` are gone; `NamedExpression` is the sole root node. See ADR-0027 (SLP: immutable Program + per-thread Memory) and ADR-0028 (named-node taxonomy). - `system/` -- Polynomial system construction, start systems (`start/total_degree_linear_product.hpp`, `start/mhom.hpp`), patches, slices. **Variables are canonical by name** (`Variable::Make` interns), and every operation that combines systems (`Concatenate`, `MakeHomotopy`, `MakeMovingHomotopy`, `+=`) makes ONE check, `CheckVariableStructuresMatch`: same variables, same groups, same order. Do not compare counts, and do not give an operation its own comparison. A `System` or `Slice` loaded from an archive re-interns itself as the last step of the load, so copies and unpickled systems recombine with their originals; a new archivable class that holds nodes needs the same step. See ADR-0068. - - `trackers/` -- Path tracking (fixed-precision and adaptive-precision trackers, predictors, Newton correctors) + - `trackers/` -- Path tracking (fixed-precision and adaptive-precision trackers, predictors, Newton correctors). **A path's result depends on the systems, the settings, the seed and its index -- never on what a reused tracker or endgame ran before it**, so threaded, MPI and serial solves are bit-identical path by path (ADR-0071). A track or run sets what it uses from its own inputs, not from a leftover member; a random direction is stored as drawn and every working copy derived from it, never rounded in place; a per-path direction is drawn at run start after the per-path reseed, never lazily. A new member of a tracker or endgame needs a stated place where each track gives it its value. - `endgames/` -- Power series and Cauchy endgames for singular endpoint handling - `nag_algorithms/` -- Higher-level algorithms (zero-dim solve; numerical irreducible decomposition is *framework scaffolding* -- not yet implemented, its `Solve()` throws) - `records/` -- The structured output directory (record schema `b2rec/1`, spec at `docs/records/b2rec-1.md`): durable, self-describing run records with full provenance and resume-by-recall. `OutputDirectory` = append-only JSONL `history/` + content-addressed `definitions/` under kind folders (`systems/`, `configs/`, `givens/`, ...); `README.txt`/`INDEX.txt`/`results.json` are derived, rebuildable views, never truth. Solvers write through the emission seam in `nag_algorithms` (`RecordTo(...)` or the ambient `BERTINI_RECORDS_DIR`); `bertini.solve` is *ensure-answered* -- paths already recorded for an identical ask (system digest + settings digest + seed) are **recalled**, not recomputed. See ADR-0042..0047 and "Persistent digests" below. A records **session lasts the whole process**: `OutputDirectory::Shared()` keeps one instance -- one history file -- per directory per process however many solves attach, and a directory's files close as soon as no solve is using it (ADR-0066). Do not return it to weak ownership: a parameter sweep then claims a session file per solve and dies after 25 in one second. @@ -188,7 +200,7 @@ Rules that follow: ## CI/CD -- `.github/workflows/build_and_test.yml` -- Builds wheels on Ubuntu/macOS/Windows and runs tests. Triggered by pull requests and pushes to `develop`/`main`. Docs-only changes (`**/*.md`, `python/docs/**`, `Doxyfile`, `VERSION`) are in `paths-ignore` — but beware: for `pull_request` events GitHub evaluates the filter against the PR's **entire** diff, so a docs-only *push to a PR that also carries code* still re-runs CI (and, via the concurrency group, cancels the PR's in-flight run — batch docs commits with code pushes). Only pushes to `develop`/`main` and PRs whose whole diff is docs-only are skipped. MPI is verified per-platform: the CLI smoke test requires `mpirun -n 2` to produce the *same solution count* as serial, for both start-system families. **Wheel builds run CONCURRENTLY with the C++ test jobs, not behind them** — they are independent full compiles (prebuilt deps + cibuildwheel vs. boost-base-from-source + ctest; different ccache namespaces), so the only gate they share is the cheap `doc_lint`. A C++-test failure still fails the run at the workflow level (and blocks publish); it just no longer serializes the wheels behind ~7 min (unix) / ~22 min (windows) of tests. Do **not** re-add `test_cpp_*` to the wheel jobs' `needs`. +- `.github/workflows/build_and_test.yml` -- Builds wheels on Ubuntu/macOS/Windows and runs tests. Triggered by pull requests and pushes to `develop`/`main`. Docs-only changes (`**/*.md`, `python/docs/**`, `Doxyfile`, `VERSION`) are in `paths-ignore` — but beware: for `pull_request` events GitHub evaluates the filter against the PR's **entire** diff, so a docs-only *push to a PR that also carries code* still re-runs CI (and, via the concurrency group, cancels the PR's in-flight run — batch docs commits with code pushes). Only pushes to `develop`/`main` and PRs whose whole diff is docs-only are skipped. The C++ suite runs twice per platform, with `BERTINI_NUM_THREADS` unset (threaded) and `=1` (serial). MPI and threading are verified per-platform by the CLI identity test: with a fixed `randomseed`, `main_data` must be *byte-identical* across a serial run, a threaded run, and `mpirun -n 2` with one and with two threads per worker, for both start-system families. **Wheel builds run CONCURRENTLY with the C++ test jobs, not behind them** — they are independent full compiles (prebuilt deps + cibuildwheel vs. boost-base-from-source + ctest; different ccache namespaces), so the only gate they share is the cheap `doc_lint`. A C++-test failure still fails the run at the workflow level (and blocks publish); it just no longer serializes the wheels behind ~7 min (unix) / ~22 min (windows) of tests. Do **not** re-add `test_cpp_*` to the wheel jobs' `needs`. - `.github/workflows/doc_lint.yml` -- A cheap wheel-free doc-correctness gate that **the build matrix depends on** (it runs first; if it fails, nothing compiles). Two parallel jobs, one per language — **run both locally before pushing docs/source changes**, or CI will bounce the whole build: - **C++ (`bash tools/doclint.sh`)** — needs `doxygen` on PATH (`brew install doxygen`). Two passes: (1) *correctness* — `@param` names must match signatures, no doc blocks on removed signatures, no unresolved `\ref`/`\cite` (zero tolerance); (2) *undocumented ratchet* — the count in `tools/doc_undocumented_baseline.txt` may only **decrease** (currently `0`, so **every new public C++ entity — including each struct data member — needs a Doxygen comment**, e.g. `///< ...`). `bash tools/doclint.sh --update-baseline` locks in a reduction. - **Python (`python tools/py_doclint.py`)** — pure-stdlib `ast`, no build/wheel/pip. Same two-pass shape over `python/bertini/**/*.py`: (1) *correctness* — a numpydoc `Parameters` entry naming a parameter absent from the signature fails (zero tolerance); (2) *undocumented ratchet* — undocumented public entities (module / public class / public function / public method; dunders excluded) counted against `tools/py_doc_undocumented_baseline.txt`, may only **decrease**. `python tools/py_doclint.py --update-baseline` locks in a reduction. Only lints the pure-Python layer; native `_pybertini` docstrings (Boost.Python strings in `python_bindings/src/*_export.cpp`) are invisible to a wheel-free tool. diff --git a/benchmark/README.md b/benchmark/README.md index a9070c14f..bf8b95586 100644 --- a/benchmark/README.md +++ b/benchmark/README.md @@ -31,7 +31,7 @@ python benchmark/run_benchmark.py \ `--no-mpi` runs the solver directly (no `mpirun`) and sweeps the thread counts you list. The run at `threads=1` is the serial baseline; every other run's speedup is reported relative to it. -The thread count is passed to the solver through `OMP_NUM_THREADS`, which the script sets for +The thread count is passed to the solver through `BERTINI_NUM_THREADS`, which the script sets for you — don't set it yourself. Example output: @@ -79,7 +79,7 @@ python benchmark/run_benchmark.py \ --threads 1 2 4 ``` -This launches each combination under `mpirun -n ` with `OMP_NUM_THREADS=` per +This launches each combination under `mpirun -n ` with `BERTINI_NUM_THREADS=` per rank. `--bind-to none` is passed to `mpirun` automatically so threaded workers can spread across a node's cores without affinity conflicts. The total worker count in the table is `ranks × threads`. diff --git a/benchmark/comparison/README.md b/benchmark/comparison/README.md index 5fe6edc79..cc4033f0d 100644 --- a/benchmark/comparison/README.md +++ b/benchmark/comparison/README.md @@ -26,7 +26,7 @@ Two phases: 1. **Serial** — one process, one thread. The MPI-built binaries are simply run directly (a single MPI process), which is the serial baseline. -2. **MPI, single thread** — `mpirun -n N` with `OMP_NUM_THREADS=1`, sweeping `--ranks`. +2. **MPI, single thread** — `mpirun -n N` with `BERTINI_NUM_THREADS=1`, sweeping `--ranks`. ## Correctness/agreement comes for free — and is *not* part of the timing diff --git a/benchmark/comparison/run_comparison.py b/benchmark/comparison/run_comparison.py index f01dd220e..f49c8a3ad 100644 --- a/benchmark/comparison/run_comparison.py +++ b/benchmark/comparison/run_comparison.py @@ -7,7 +7,7 @@ solvers by construction. Two phases: * serial -- one process, one thread, every system through every solver; - * MPI single-thread -- sweep --ranks with OMP_NUM_THREADS=1, each solver under `mpirun -n N`. + * MPI single-thread -- sweep --ranks with BERTINI_NUM_THREADS=1, each solver under `mpirun -n N`. We get a correctness/agreement check for free: every solver should report the same solution count, equal to the system's known expected count. That check is done AFTER timing, from output files, and diff --git a/benchmark/comparison/solvers.py b/benchmark/comparison/solvers.py index 804e7db37..bea79099d 100644 --- a/benchmark/comparison/solvers.py +++ b/benchmark/comparison/solvers.py @@ -49,7 +49,7 @@ def _time_run(argv, input_text, threads, timeout): f.write(input_text) env = os.environ.copy() - env["OMP_NUM_THREADS"] = str(threads) # single-thread for this benchmark phase + env["BERTINI_NUM_THREADS"] = str(threads) # b2's path threads for this phase t0 = time.perf_counter() try: diff --git a/benchmark/run_benchmark.py b/benchmark/run_benchmark.py index 78e20d9e0..6e01cf715 100755 --- a/benchmark/run_benchmark.py +++ b/benchmark/run_benchmark.py @@ -16,7 +16,7 @@ Requirements: - mpirun must be on PATH - bertini2 must be compiled with MPI support (BERTINI2_HAVE_MPI) - - OMP_NUM_THREADS controls threads per MPI rank + - BERTINI_NUM_THREADS controls threads per MPI rank (the script sets it per run) """ import argparse @@ -40,7 +40,7 @@ def parse_args(): p.add_argument("--ranks", nargs="+", type=int, default=[1, 2, 4], help="MPI rank counts to sweep (default: 1 2 4)") p.add_argument("--threads", nargs="+", type=int, default=[1], - help="OMP thread counts to sweep (default: 1)") + help="Thread counts per rank to sweep, via BERTINI_NUM_THREADS (default: 1)") p.add_argument("--output", default="benchmark_results.csv", help="CSV output file (default: benchmark_results.csv)") p.add_argument("--timeout", type=float, default=600.0, @@ -68,7 +68,7 @@ def run_once(bertini2_path, input_file, ranks, threads, mpirun_cmd, mpirun_args, Run bertini2 with the given parallelism settings in a fresh temp directory. When use_mpi is True the solver is launched under `mpirun -n ` (MPI across ranks, - OMP_NUM_THREADS threads within each rank). When use_mpi is False the binary is run directly + BERTINI_NUM_THREADS threads within each rank). When use_mpi is False the binary is run directly -- a pure shared-memory thread sweep that needs no MPI at all; ranks is ignored (always 1). Returns (wall_time_s, solutions_found) or (float('nan'), -1) on failure. @@ -79,9 +79,9 @@ def run_once(bertini2_path, input_file, ranks, threads, mpirun_cmd, mpirun_args, shutil.copy2(input_file, dest_input) env = os.environ.copy() - # OMP_NUM_THREADS drives the worker-thread count in both modes: per MPI rank under mpirun, - # and the whole shared-memory solve when run directly. - env["OMP_NUM_THREADS"] = str(threads) + # BERTINI_NUM_THREADS drives the worker-thread count in both modes: per MPI rank under + # mpirun, and the whole shared-memory solve when run directly. + env["BERTINI_NUM_THREADS"] = str(threads) if use_mpi: cmd = [mpirun_cmd, *shlex.split(mpirun_args), "-n", str(ranks), diff --git a/benchmark/submit.slurm.sh b/benchmark/submit.slurm.sh index af1bf365f..e6a49145e 100755 --- a/benchmark/submit.slurm.sh +++ b/benchmark/submit.slurm.sh @@ -5,7 +5,7 @@ #SBATCH --job-name=bertini2_benchmark #SBATCH --nodes=4 # total number of nodes #SBATCH --ntasks-per-node=1 # one MPI rank per node -#SBATCH --cpus-per-task=8 # CPU cores per rank (= OMP_NUM_THREADS) +#SBATCH --cpus-per-task=8 # CPU cores per rank (the script sets BERTINI_NUM_THREADS per run) #SBATCH --time=02:00:00 #SBATCH --output=benchmark_%j.log #SBATCH --error=benchmark_%j.err diff --git a/core/include/bertini2/endgames/base_endgame.hpp b/core/include/bertini2/endgames/base_endgame.hpp index f5d9bcb99..c1a814a6c 100644 --- a/core/include/bertini2/endgames/base_endgame.hpp +++ b/core/include/bertini2/endgames/base_endgame.hpp @@ -146,22 +146,46 @@ class EndgameBase : BCT target_time_{}; ///< The final target time (default 0); set via SetTargetTime(). /** - The random direction the c/k estimate is projected along. Generated once per size+precision - and reused for the whole endgame, so consecutive estimates differ only because the samples - differ. A fresh random probe every call made the stabilization test noisy -- it could certify - the operating zone spuriously -- and churned mpfr allocations. See z_notes/20260629. + The random direction the c/k estimate is projected along, as drawn: in double precision, so + that it is exact in every lane and at every multiple precision. Drawn at the start of every + run (RefreshCOverKProbe) and fixed for the run, so consecutive estimates differ only because + the samples differ. A fresh random probe every call made the stabilization test noisy -- it + could certify the operating zone spuriously -- and churned mpfr allocations. See + z_notes/20260629. */ + mutable Vec c_over_k_probe_as_drawn_; + + /// The c/k probe in each numeric type the endgame computes in: exact conversions of + /// c_over_k_probe_as_drawn_, re-derived when their size or precision is wrong, never rounded. mutable TupOfVec c_over_k_probe_; + /** + \brief Draw the c/k probe direction for a run. + + Called at the start of every run, after whatever reseed the caller does for the path, so the + direction is a function of the random-number stream at that point and nothing else. + \param size The number of coordinates the probe must have. + */ + // Drawn per run, never lazily on first use: a lazy draw lands in whichever path first needs + // the probe, so the direction a path used, and the draws it consumed, would depend on the + // paths the same endgame ran before it, and on which thread ran them (#378, ADR-0071). + void RefreshCOverKProbe(unsigned size) const + { + c_over_k_probe_as_drawn_.resize(size); + for (unsigned ii = 0; ii < size; ++ii) + c_over_k_probe_as_drawn_(ii) = RandomUnit(); + std::apply([](auto&... working){ (working.resize(0), ...); }, c_over_k_probe_); + } /** \brief The fixed random direction the c/k estimate is projected along. - Regenerated only when the size changes; re-precisioned in place otherwise, so the direction - itself is preserved across a precision increase. + An exact conversion of the direction as drawn, at the requested precision. Converted, not + rounded in place: the drawn values are doubles, which every multiple precision holds exactly, + so the direction is the same in both lanes and survives any change of precision. \tparam ComplexT The complex number type. \param size The number of coordinates the probe must have. @@ -172,15 +196,21 @@ class EndgameBase : Vec const& GetCOverKProbe(unsigned size, unsigned prec) const { using bertini::Precision; + if (static_cast(c_over_k_probe_as_drawn_.size()) != size) + RefreshCOverKProbe(size); // samples of another size than this run began with + auto& probe = std::get >(c_over_k_probe_); - if (static_cast(probe.size()) != size) + bool stale = static_cast(probe.size()) != size; + if constexpr (!std::is_same::value) + stale = stale || Precision(probe) != prec; + if (stale) { probe.resize(size); for (unsigned ii = 0; ii < size; ++ii) - probe(ii) = RandomUnit(); + probe(ii) = ComplexT(c_over_k_probe_as_drawn_(ii)); + if constexpr (!std::is_same::value) + Precision(probe, prec); } - if (Precision(probe) != prec) - Precision(probe, prec); return probe; } @@ -313,6 +343,8 @@ class EndgameBase : // endgame concatenates every path into one undifferentiated heap. Initializing was // declared, and handled by GoryDetailLogger, but emitted by nothing. NotifyObservers(Initializing(this->AsFlavor())); + // a run's c/k direction is drawn here, at its start, never inherited from an earlier run + RefreshCOverKProbe(static_cast(start_point.size())); auto prec = Precision(start_point); BCT t = start_time_; Precision(t, prec); BCT t0 = target_time_; Precision(t0, prec); diff --git a/core/include/bertini2/endgames/cauchy.hpp b/core/include/bertini2/endgames/cauchy.hpp index e0b459d03..dec39843f 100644 --- a/core/include/bertini2/endgames/cauchy.hpp +++ b/core/include/bertini2/endgames/cauchy.hpp @@ -1410,6 +1410,12 @@ class CauchyEndgame : this->adaptive_numeric_type_active_ = true; struct Disarmer { bool& flag; ~Disarmer(){ flag = false; } } disarm{this->adaptive_numeric_type_active_}; + // Both lanes start empty: a phase clears only the lane it runs in, so a path whose endgame + // stays in double would otherwise leave an earlier path's multiprecision samples and times + // in place for anything that reads that lane first (#378). + ClearTimesAndSamples(); + ClearTimesAndSamples(); + // Start in the precision the tracker handed us at the endgame boundary: double for the easy // majority, already-mpfr for the few paths that escalated before the endgame. Migrate UP only. this->current_endgame_precision_ = std::max(DoublePrecision(), Precision(start_point)); @@ -1667,7 +1673,6 @@ class CauchyEndgame : this->template CrossSampsUp<>(pseg_samples_, newprec); this->template CrossTimesUp<>(cauchy_times_, newprec); this->template CrossSampsUp<>(cauchy_samples_, newprec); - this->template CrossVecUp<> (this->c_over_k_probe_, newprec); } else { @@ -1675,8 +1680,6 @@ class CauchyEndgame : tracking::adaptive::SetPrecision(std::get>(pseg_samples_), newprec); tracking::adaptive::SetPrecision(std::get>(cauchy_times_), newprec); tracking::adaptive::SetPrecision(std::get>(cauchy_samples_), newprec); - auto& pm = std::get>(this->c_over_k_probe_); - if (pm.size() > 0) Precision(pm, newprec); } if (this->final_approximation_.size() > 0) Precision(this->final_approximation_, newprec); diff --git a/core/include/bertini2/endgames/powerseries.hpp b/core/include/bertini2/endgames/powerseries.hpp index 08d8ccde6..caaac3c1c 100644 --- a/core/include/bertini2/endgames/powerseries.hpp +++ b/core/include/bertini2/endgames/powerseries.hpp @@ -911,6 +911,13 @@ class PowerSeriesEndgame : this->adaptive_numeric_type_active_ = true; struct Disarmer { bool& flag; ~Disarmer(){ flag = false; } } disarm{this->adaptive_numeric_type_active_}; + // Both lanes start empty. A segment clears only the lane it runs in, so a path whose + // endgame stays in double would otherwise leave an earlier path's multiprecision times in + // place -- and LatestTime reads that lane first whenever it is non-empty, reporting the + // earlier path's last time as this one's (#378). + ClearTimesAndSamples(); + ClearTimesAndSamples(); + this->current_endgame_precision_ = std::max(DoublePrecision(), Precision(start_point)); // ---- SETUP (restart-in-mpfr on escalation; bounded) ---- diff --git a/core/include/bertini2/nag_algorithms/common/config.hpp b/core/include/bertini2/nag_algorithms/common/config.hpp index b18371866..fab254dff 100644 --- a/core/include/bertini2/nag_algorithms/common/config.hpp +++ b/core/include/bertini2/nag_algorithms/common/config.hpp @@ -215,10 +215,10 @@ struct ZeroDimConfig unsigned initial_ambient_precision = DefaultPrecision(); ///< Initial ambient (working) precision for the solve. unsigned max_num_crossed_path_resolve_attempts = 2; ///< The maximum number of times to attempt to re-solve crossed paths at the endgame boundary. - /// Number of worker threads for a shared-memory (non-MPI) solve. 0 = auto - /// (std::thread::hardware_concurrency); 1 = serial (no thread pool). Overridden by the - /// OMP_NUM_THREADS environment variable when set. See parallel::EffectiveThreadCount. - /// Under MPI the per-rank thread count comes from OMP_NUM_THREADS, not this field. + /// Number of worker threads tracking paths, in a standalone solve and on each MPI rank. + /// 0 = auto (every CPU this process may run on, as OpenMP's default); 1 = serial (no + /// thread pool). Overridden by the BERTINI_NUM_THREADS environment variable when set. + /// See parallel::EffectiveThreadCount. unsigned num_threads = 0; /// Wall-clock budget for each path, in seconds; 0 (the default) means none. A path that has not @@ -254,14 +254,15 @@ struct MetaConfig }; /** -Global RNG seed for reproducible runs. random_seed == 0 (the default) draws from -std::random_device and reports the effective seed so the run can be reproduced. -Set via `randomseed: N;` in the classic Bertini input file or bertini.set_random_seed(N) -in Python. Must be applied before system construction (gamma, patch, TD-constants). +The RNG seed a classic Bertini input file asks for, via `randomseed: N;`. As in Bertini 1, +random_seed == 0 (the default) means "draw from entropy": the command-line program then draws +one and reports it, so the run can be reproduced. The library API differs: to +SetGlobalSeed and bertini.set_random_seed, 0 is an ordinary seed. Must be applied +before system construction (gamma, patch, TD-constants). */ struct RandomConfig { - unsigned long random_seed = 0; ///< The RNG seed; 0 draws from std::random_device and reports the effective seed. + unsigned long random_seed = 0; ///< The classic input's seed; 0 means draw one from entropy and report it. }; diff --git a/core/include/bertini2/nag_algorithms/zero_dim_solve.hpp b/core/include/bertini2/nag_algorithms/zero_dim_solve.hpp index 07a889f87..07dc53f84 100644 --- a/core/include/bertini2/nag_algorithms/zero_dim_solve.hpp +++ b/core/include/bertini2/nag_algorithms/zero_dim_solve.hpp @@ -757,7 +757,10 @@ run the endgame, classify the endpoints, report. See the forward-declare doc ab PreSolveChecks(); PreSolveSetup(); - const int n_threads = parallel::WorkerThreadCount(); + // threads per rank: the same rule as the standalone solve (BERTINI_NUM_THREADS, + // then num_threads, then every CPU this rank may run on) + const int n_threads = static_cast( + parallel::EffectiveThreadCount(this->template Get().num_threads)); // One round = dispatch a set of path indices, each executed as a WHOLE path // (pre-endgame + endgame) by a worker. Round 0 is every path; later rounds re-run @@ -866,7 +869,7 @@ run the endgame, classify the endpoints, report. See the forward-declare doc ab // Trackers/Endgames have no default ctor, so aggregate-init via new. std::unique_ptr s( new PathThreadState{ Clone(GetTracker().GetSystem()), Clone(TargetSystem()), - GetTracker(), GetEndgame(), {}, {} }); + GetTracker(), GetEndgame(), {} }); s->tracker.SetSystem(s->sys); s->endgame.SetTracker(s->tracker); SetThreadPrecision(this->template Get().initial_ambient_precision); @@ -1197,8 +1200,9 @@ run the endgame, classify the endpoints, report. See the forward-declare doc ab indices_to_run = RecallRecordedPaths(all_indices); } - // num_threads: 0 = auto (hardware_concurrency), 1 = serial, N = N threads; - // OMP_NUM_THREADS overrides. n_threads <= 1 takes the pool-free serial path. + // num_threads: 0 = auto (every CPU this process may run on), 1 = serial, N = N + // threads; BERTINI_NUM_THREADS overrides. n_threads <= 1 takes the pool-free + // serial path. const unsigned n_threads = parallel::EffectiveThreadCount(this->template Get().num_threads); @@ -1619,8 +1623,7 @@ run the endgame, classify the endpoints, report. See the forward-declare doc ab struct BeforeEGContext { TrackerType& tracker; - tracking::FirstPrecisionRecorder& first_prec_rec; - tracking::MinMaxPrecisionRecorder& min_max_prec; + tracking::PathPrecisionRecorder& path_prec; }; struct DuringEGContext { @@ -1630,17 +1633,29 @@ run the endgame, classify the endpoints, report. See the forward-declare doc ab SystemType const& target_sys; TrackerType& tracker; EndgameType& endgame; - tracking::FirstPrecisionRecorder& first_prec_rec; - tracking::MinMaxPrecisionRecorder& min_max_prec; + tracking::PathPrecisionRecorder& path_prec; }; BeforeEGContext MemberBeforeEGContext() { - return BeforeEGContext{ GetTracker(), first_prec_rec_, min_max_prec_ }; + return BeforeEGContext{ GetTracker(), path_prec_ }; } DuringEGContext MemberDuringEGContext() { - return DuringEGContext{ TargetSystem(), GetTracker(), GetEndgame(), first_prec_rec_, min_max_prec_ }; + return DuringEGContext{ TargetSystem(), GetTracker(), GetEndgame(), path_prec_ }; + } + + /// Fold the path's precision record so far into its metadata: called after each stage. + static void HarvestPathPrecision(SolutionMetaData& smd, + tracking::PathPrecisionRecorder const& rec) + { + if (rec.DidPrecisionIncrease() && !smd.precision_changed) + { + smd.precision_changed = true; + smd.time_of_first_prec_increase = rec.TimeOfIncrease(); + } + using std::max; + smd.max_precision_used = max(smd.max_precision_used, rec.MaxPrecision()); } /** @@ -1673,12 +1688,6 @@ run the endgame, classify the endpoints, report. See the forward-declare doc ab { ReseedThisThread(static_cast(soln_ind)); - if (tracking::TrackerTraits::IsAdaptivePrec) - { - ctx.tracker.AddObserver(ctx.first_prec_rec); - ctx.tracker.AddObserver(ctx.min_max_prec); - } - auto& smd = solution_final_metadata_[soln_ind]; smd.path_index = soln_ind; smd.solution_index = soln_ind; @@ -1754,19 +1763,7 @@ run the endgame, classify the endpoints, report. See the forward-declare doc ab smd.pre_endgame_success_code = tracking_success; if (tracking::TrackerTraits::IsAdaptivePrec) - { - if (ctx.first_prec_rec.DidPrecisionIncrease()) - { - smd.precision_changed = true; - smd.time_of_first_prec_increase = ctx.first_prec_rec.TimeOfIncrease(); - } - else - ctx.tracker.RemoveObserver(ctx.first_prec_rec); - ctx.tracker.RemoveObserver(ctx.min_max_prec); - using std::max; - smd.max_precision_used = - max(smd.max_precision_used, ctx.min_max_prec.MaxPrecision()); - } + HarvestPathPrecision(smd, ctx.path_prec); StampWhereThePathGotTo(smd, ctx.tracker, tracking_success, ctx.tracker.CurrentTime()); } @@ -1807,7 +1804,7 @@ run the endgame, classify the endpoints, report. See the forward-declare doc ab used by both the standalone (MPI-less) threaded solve and the threaded MPI worker. Each std::thread owns one: a homotopy copy (tracked by `tracker`), a target-system copy (residual evaluation mutates System precision state), a Tracker and Endgame, and its own - precision observers so observer-derived metadata matches serial runs. Build a + precision recorder so observer-derived metadata matches serial runs. Build a DuringEGContext from it to drive ExecuteOnePath. */ struct PathThreadState @@ -1816,12 +1813,11 @@ run the endgame, classify the endpoints, report. See the forward-declare doc ab System target_sys; // for function residuals / dehomogenization TrackerType tracker; EndgameType endgame; - tracking::FirstPrecisionRecorder first_prec_rec; - tracking::MinMaxPrecisionRecorder min_max_prec; + tracking::PathPrecisionRecorder path_prec; DuringEGContext Context() { - return DuringEGContext{ target_sys, tracker, endgame, first_prec_rec, min_max_prec }; + return DuringEGContext{ target_sys, tracker, endgame, path_prec }; } }; @@ -1840,11 +1836,11 @@ run the endgame, classify the endpoints, report. See the forward-declare doc ab { return [this]() -> std::unique_ptr { - // Trackers/Endgames have no default ctor, so aggregate-init via new; {},{} default - // the two precision recorders. + // Trackers/Endgames have no default ctor, so aggregate-init via new; {} defaults + // the path's precision recorder. std::unique_ptr s( new PathThreadState{ Clone(GetTracker().GetSystem()), Clone(TargetSystem()), - GetTracker(), GetEndgame(), {}, {} }); + GetTracker(), GetEndgame(), {} }); s->tracker.SetSystem(s->sys); // re-point the copy at its own System s->endgame.SetTracker(s->tracker); // ... and the endgame at that tracker SetThreadPrecision(this->template Get().initial_ambient_precision); @@ -1966,12 +1962,23 @@ run the endgame, classify the endpoints, report. See the forward-declare doc ab this->NotifyObservers(PathStarted(*this, static_cast(soln_ind), exec_tracker)); + // One precision record for the whole path, across the track to the boundary and + // every endgame sub-track: started here, read after each stage, detached at the end. + constexpr bool adaptive = tracking::TrackerTraits::IsAdaptivePrec; + if (adaptive) + { + ctx.path_prec.Reset(); + ctx.tracker.AddObserver(ctx.path_prec); + } + ctx.tracker.SetTrackingTolerance(midpath_retrack_tolerance_); - ExecuteBeforeEG(BeforeEGContext{ ctx.tracker, ctx.first_prec_rec, ctx.min_max_prec }, soln_ind, start_point); + ExecuteBeforeEG(BeforeEGContext{ ctx.tracker, ctx.path_prec }, soln_ind, start_point); ReadSolveTimeoutAsAnInterrupt(solution_final_metadata_[soln_ind].pre_endgame_success_code, path_start_clock); if (solution_final_metadata_[soln_ind].pre_endgame_success_code != SuccessCode::Success) { + if (adaptive) + ctx.tracker.RemoveObserver(ctx.path_prec); ctx.tracker.ClearMaxWallClockTime(); stamp_path_time(); this->NotifyObservers(PathComplete(*this, static_cast(soln_ind), exec_tracker, @@ -1981,6 +1988,8 @@ run the endgame, classify the endpoints, report. See the forward-declare doc ab ctx.tracker.SetTrackingTolerance(this->template Get().newton_during_endgame); ExecuteDuringEG(ctx, soln_ind); + if (adaptive) + ctx.tracker.RemoveObserver(ctx.path_prec); ReadSolveTimeoutAsAnInterrupt(solution_final_metadata_[soln_ind].endgame_success_code, path_start_clock); ctx.tracker.ClearMaxWallClockTime(); // the budget was this path's; do not carry it to the next @@ -2118,12 +2127,6 @@ run the endgame, classify the endpoints, report. See the forward-declare doc ab ReseedThisThread(static_cast(soln_ind) + static_cast(num_start_points_)); auto& smd = solution_final_metadata_[soln_ind]; - if (tracking::TrackerTraits::IsAdaptivePrec) - { - if (!smd.precision_changed) - ctx.tracker.AddObserver(ctx.first_prec_rec); - ctx.tracker.AddObserver(ctx.min_max_prec); - } const auto& bdry_point = solutions_at_endgame_boundary_[soln_ind].path_point; @@ -2155,28 +2158,12 @@ run the endgame, classify the endpoints, report. See the forward-declare doc ab eg_success == SuccessCode::Success ? ctx.endgame.LatestTime() : ctx.tracker.CurrentTime()); - // Harvest the AMP observers REGARDLESS of the endgame outcome: the precision a - // path used is a fact about the tracking that happened, and it is exactly the - // FAILED paths -- e.g. slow divergers escalating in the mp lane -- whose - // precision honesty matters most for diagnostics. (Previously the failure - // path returned before this harvest, so a path that escalated to 70 digits - // and then failed reported its pre-endgame precision, e.g. 16.) + // Harvest the path's precision record REGARDLESS of the endgame outcome: the + // precision a path used is a fact about the tracking that happened, and it is + // exactly the FAILED paths -- e.g. slow divergers escalating in the mp lane -- + // whose precision honesty matters most for diagnostics. if (tracking::TrackerTraits::IsAdaptivePrec) - { - if (!smd.precision_changed) - { - if (ctx.first_prec_rec.DidPrecisionIncrease()) - { - smd.precision_changed = true; - smd.time_of_first_prec_increase = ctx.first_prec_rec.TimeOfIncrease(); - } - ctx.tracker.RemoveObserver(ctx.first_prec_rec); - } - ctx.tracker.RemoveObserver(ctx.min_max_prec); - using std::max; - smd.max_precision_used = - max(smd.max_precision_used, ctx.min_max_prec.MaxPrecision()); - } + HarvestPathPrecision(smd, ctx.path_prec); // an unsuccessful endgame has no final approximation, so the final-point-dependent // metadata cannot be computed. @@ -3091,8 +3078,7 @@ run the endgame, classify the endpoints, report. See the forward-declare doc ab /// observers used during tracking // i feel like these should be factored out into some policy class which prescribes how they are used, so that the actions taken are customizable. - tracking::FirstPrecisionRecorder first_prec_rec_; ///< Records the time of the first precision increase on each path. - tracking::MinMaxPrecisionRecorder min_max_prec_; ///< Records the min/max precision used on each path. + tracking::PathPrecisionRecorder path_prec_; ///< Records the precision each path used, across all its tracks. /// function objects used during the algorithm diff --git a/core/include/bertini2/parallel/thread_pool.hpp b/core/include/bertini2/parallel/thread_pool.hpp index 5b5513ccb..12aa65ff8 100644 --- a/core/include/bertini2/parallel/thread_pool.hpp +++ b/core/include/bertini2/parallel/thread_pool.hpp @@ -59,6 +59,10 @@ of the thread. // "MPI-less threading" plan: threading must be available in builds (and wheels) compiled // without any MPI installation. +#if defined(__linux__) +#include // sched_getaffinity, CPU_COUNT (AvailableCpuCount) +#endif + #include #include // std::getenv, std::atoi #include @@ -74,22 +78,50 @@ namespace bertini { namespace parallel { +/** +\brief The number of CPUs this process may run on. + +On Linux, the size of the process's CPU affinity mask, so a process bound to a subset of +the machine (an MPI rank that `mpirun` bound, a job step confined to a cpuset) counts only +its own CPUs, as OpenMP's default thread count does. Elsewhere, and +if the mask cannot be read, std::thread::hardware_concurrency(). Always >= 1. +*/ +inline unsigned AvailableCpuCount() +{ +#if defined(__linux__) + cpu_set_t mask; + CPU_ZERO(&mask); + if (sched_getaffinity(0, sizeof(mask), &mask) == 0) + { + const int n = CPU_COUNT(&mask); + if (n >= 1) + return static_cast(n); + } +#endif + const unsigned hw = std::thread::hardware_concurrency(); + return hw >= 1 ? hw : 1u; +} + + /** \brief Resolve a configured thread count to the number of worker threads to actually spawn. Precedence (highest first): - 1. The OMP_NUM_THREADS environment variable, if set to an integer >= 1. HPC schedulers - (SLURM) set this from --cpus-per-task, so honoring it keeps the MPI+threads hybrid and - the standalone threaded solve behaving the same under a job allocation. - 2. The `configured` value, if >= 1 (e.g. a num_threads config knob set by the user). - 3. std::thread::hardware_concurrency() when `configured == 0` ("auto"). - -Always returns >= 1 (a return of 1 means "run serially, no pool"). hardware_concurrency() -can report 0 on exotic platforms; we clamp that to 1. + 1. The BERTINI_NUM_THREADS environment variable, if set to an integer >= 1. + 2. The `configured` value, if >= 1 (the `num_threads` config field). + 3. AvailableCpuCount() when `configured == 0` ("auto"): every CPU this process may run on, + the same default OpenMP uses, so a solve is threaded without any configuration. + +The same rule serves the standalone threaded solve and each MPI rank. OMP_NUM_THREADS has +no effect. + +Always returns >= 1 (a return of 1 means "run serially, no pool"). */ +// b2 reads its own variable, never OMP_NUM_THREADS: numpy's OpenBLAS takes its thread count from +// that one, so a shared name would tie the two libraries' thread counts together (ADR-0070). inline unsigned EffectiveThreadCount(unsigned configured = 0) { - if (const char* env = std::getenv("OMP_NUM_THREADS")) + if (const char* env = std::getenv("BERTINI_NUM_THREADS")) { int n = std::atoi(env); if (n >= 1) @@ -99,8 +131,7 @@ inline unsigned EffectiveThreadCount(unsigned configured = 0) if (configured >= 1) return configured; - unsigned hw = std::thread::hardware_concurrency(); - return hw >= 1 ? hw : 1u; + return AvailableCpuCount(); } diff --git a/core/include/bertini2/parallel/worker.hpp b/core/include/bertini2/parallel/worker.hpp index 75065d28a..7e2d45fe9 100644 --- a/core/include/bertini2/parallel/worker.hpp +++ b/core/include/bertini2/parallel/worker.hpp @@ -43,7 +43,6 @@ Each tracking thread owns its state (System copy + Tracker copy) built by state_ #include "bertini2/parallel/mpi_include.hpp" #include -#include // std::getenv #include #include #include @@ -172,25 +171,6 @@ void RunWorkerLoopThreaded( } -/** -\brief Read the thread count for a worker rank from OMP_NUM_THREADS. - -Returns the value of OMP_NUM_THREADS if set and >= 1, otherwise 1 (serial). -HPC schedulers (SLURM) set OMP_NUM_THREADS automatically from --cpus-per-task. -*/ -inline int WorkerThreadCount() -{ - const char* env = std::getenv("OMP_NUM_THREADS"); - if (env) - { - int n = std::atoi(env); - if (n >= 1) - return n; - } - return 1; -} - - } // namespace parallel } // namespace bertini diff --git a/core/include/bertini2/random.hpp b/core/include/bertini2/random.hpp index 813afa998..5bdf9634a 100644 --- a/core/include/bertini2/random.hpp +++ b/core/include/bertini2/random.hpp @@ -50,16 +50,35 @@ namespace bertini std::mt19937& ThreadEngine(); /** - Set the global RNG seed. seed == 0 draws from std::random_device and stores - the effective (non-zero) seed so it can be retrieved and reproduced later. - Must be called before system construction (gamma, patch, TD-constants) to make - those setup draws deterministic. + Set the global RNG seed. Every value is a seed, 0 included, as in numpy and Python's + random; SetGlobalSeedFromEntropy asks for an unpredictable one. Must be called before + system construction (gamma, patch, TD-constants) to make those setup draws deterministic. + + The classic input format differs: there `randomseed: 0`, its default, means "draw from + entropy", as in Bertini 1, and the command-line program reads it that way. */ + // 0 is a seed here and entropy in the classic format: each follows its users' convention (ADR-0072) void SetGlobalSeed(unsigned long seed); /** - Returns the effective global seed. If SetGlobalSeed has never been called, draws - from entropy on first call and caches the result. + Set the global RNG seed to one drawn from std::random_device, and return it, so the run + can be reproduced by setting that seed. + */ + unsigned long SetGlobalSeedFromEntropy(); + + /** + Apply a classic input file's `randomseed`, as Bertini 1 does: 0, the default, draws a seed + from entropy; any other value is the seed. Only the classic input format reads 0 this way; + to SetGlobalSeed it is an ordinary seed. + + \param randomseed The value of `randomseed` in the classic input file. + \return The seed now in effect. + */ + unsigned long ApplyClassicRandomSeed(unsigned long randomseed); + + /** + Returns the effective global seed. If no seed has been set, draws one from entropy on + the first call and keeps it. */ unsigned long GetGlobalSeed(); @@ -95,7 +114,8 @@ namespace bertini rekeys to each derived seed, so consecutive seedless solves form a deterministic seed chain from the initial master: same master, same sequence of solves, same seeds -- on every platform (the value is masked to 32 bits for LLP64 Windows). - Never returns 0 (SetGlobalSeed treats 0 as "draw from entropy"). + Never returns 0, so a recorded seed replays through a classic input file, where + `randomseed: 0` means "draw from entropy". \return A nonzero seed for this solve, ready for SetGlobalSeed and the ask record. */ diff --git a/core/include/bertini2/trackers/amp_tracker.hpp b/core/include/bertini2/trackers/amp_tracker.hpp index 35ffde545..089b9a8a8 100644 --- a/core/include/bertini2/trackers/amp_tracker.hpp +++ b/core/include/bertini2/trackers/amp_tracker.hpp @@ -537,12 +537,21 @@ namespace bertini{ endtime_ = end_time; endtime_.precision(initial_precision_); - current_stepsize_.precision(initial_precision_); + // The first step is built at THIS track's start precision, and the precision is set + // after the assignment, never before it. current_precision_ still holds wherever the + // previous track ended, and an assignment takes its source's precision + // (preserve_related_precision), so building it at current_precision_ rounded the + // exact initial_step_size at the previous track's precision. A track that starts in + // double then rounded that a second time, to a double, and the last bit of its first + // step -- and so its whole step sequence -- depended on which path the tracker had + // tracked before: thread scheduling made the threaded solve irreproducible, and path + // order made a serial solve depend on it (#378). if (reinitialize_stepsize_) { real_mp segment_length = abs(start_time-end_time)/Get().min_num_steps; - SetStepSize(min(real_mp(Get().initial_step_size, current_precision_),segment_length)); + SetStepSize(min(real_mp(Get().initial_step_size, initial_precision_),segment_length)); } + current_stepsize_.precision(initial_precision_); // populate the current space value with the start point, in appropriate precision if (initial_precision_==DoublePrecision()) @@ -1508,7 +1517,12 @@ namespace bertini{ for (unsigned ii=0; ii >(current_space_)(ii) = complex_dbl(source_point(ii)); - endtime_.precision(DoublePrecision()); // i question this one 2021-04-12 + // The same internals every other conversion re-precisions: the step size, the + // time, its increment, the end time, and the predictor's and corrector's work + // space. Re-precisioning only the end time here left the rest at whatever + // multiple precision the track came down from, so a track's double-precision + // steps depended on the precision history before them (#378). + AdjustInternalsPrecision(DoublePrecision()); } /** diff --git a/core/include/bertini2/trackers/explicit_predictors.hpp b/core/include/bertini2/trackers/explicit_predictors.hpp index dedc49b20..75875be6d 100644 --- a/core/include/bertini2/trackers/explicit_predictors.hpp +++ b/core/include/bertini2/trackers/explicit_predictors.hpp @@ -342,7 +342,8 @@ namespace bertini{ std::get< Vec >(step_temp_).resize(numTotalFunctions_); std::get< Vec >(step_temp_).resize(numTotalFunctions_); std::get< Vec >(rand_temp_) = RandomOfUnits(numVariables_); - std::get< Vec >(rand_temp_) = RandomOfUnits(numVariables_); + probe_as_drawn_ = RandomOfUnits(numVariables_); + std::get< Vec >(rand_temp_) = probe_as_drawn_; std::get< Vec >(solve_temp_).resize(numVariables_); std::get< Vec >(solve_temp_).resize(numVariables_); @@ -376,7 +377,9 @@ namespace bertini{ */ void SetConditionProbe(std::tuple< Vec, Vec > const& probe) { - rand_temp_ = probe; + std::get< Vec >(rand_temp_) = std::get< Vec >(probe); + probe_as_drawn_ = std::get< Vec >(probe); + std::get< Vec >(rand_temp_) = probe_as_drawn_; } @@ -402,7 +405,13 @@ namespace bertini{ Precision(std::get< Mat >(dh_dx_0_),new_precision); Precision(std::get< Mat >(dh_dx_temp_),new_precision); Precision(std::get< Vec >(step_temp_),new_precision); - Precision(std::get< Vec >(rand_temp_),new_precision); + // the probe is re-derived from the stored copy, never rounded in place: a drop + // in precision would otherwise destroy digits that a later rise cannot restore + { + auto& probe = std::get< Vec >(rand_temp_); + probe = probe_as_drawn_; + Precision(probe, new_precision); + } Precision(std::get< Vec >(solve_temp_),new_precision); Precision(std::get< Vec >(stage_pt_temp_),new_precision); Precision(std::get< Vec >(err_temp_),new_precision); @@ -960,6 +969,7 @@ namespace bertini{ mutable std::tuple< Vec, Vec > step_temp_; // reused scratch for FullStep stage accumulation mutable std::tuple< Vec, Vec > rand_temp_; // reused scratch: random RHS for norm_J_inverse + Vec probe_as_drawn_; ///< the probe direction at the precision it was drawn at; rand_temp_'s multiprecision half is this at the working precision mutable std::tuple< Vec, Vec > solve_temp_; // reused scratch: LU solve result mutable std::tuple< Vec, Vec > stage_pt_temp_; // reused scratch: RK stage point (current_space + delta_t*temp) mutable std::tuple< Vec, Vec > err_temp_; // reused scratch: embedded-RK error estimate vector diff --git a/core/include/bertini2/trackers/newton_corrector.hpp b/core/include/bertini2/trackers/newton_corrector.hpp index 028de82e6..be4287030 100644 --- a/core/include/bertini2/trackers/newton_corrector.hpp +++ b/core/include/bertini2/trackers/newton_corrector.hpp @@ -86,7 +86,13 @@ namespace bertini{ Precision(std::get< Vec >(f_temp_), new_precision); Precision(std::get< Vec >(step_temp_), new_precision); Precision(std::get< Mat >(J_temp_), new_precision); - Precision(std::get< Vec >(rand_temp_), new_precision); + // the probe is re-derived from the stored copy, never rounded in place: a drop + // in precision would otherwise destroy digits that a later rise cannot restore + { + auto& probe = std::get< Vec >(rand_temp_); + probe = probe_as_drawn_; + Precision(probe, new_precision); + } Precision(std::get< Vec >(solve_temp_), new_precision); std::get< linalg::PartialPivLU >(LU_).ChangePrecision(new_precision); @@ -137,17 +143,20 @@ namespace bertini{ void RefreshRandomDirection() { std::get< Vec >(rand_temp_) = RandomOfUnits(numVariables_); - std::get< Vec >(rand_temp_) = RandomOfUnits(numVariables_); + probe_as_drawn_ = RandomOfUnits(numVariables_); + std::get< Vec >(rand_temp_) = probe_as_drawn_; } /** - \brief The condition-number probe direction (both precisions), so the tracker can share - the SAME per-path probe with the predictor -- making the predictor's and corrector's - ||J^{-1}|| estimates use one consistent direction. \see RefreshRandomDirection. + \brief The condition-number probe direction as drawn (both precisions), so the tracker + can share the SAME per-path probe with the predictor -- making the predictor's and + corrector's ||J^{-1}|| estimates use one consistent direction. The multiprecision + half is the stored copy, at the precision it was drawn at, not the working copy the + corrector holds at its current precision. \see RefreshRandomDirection. */ - std::tuple< Vec, Vec > const& ConditionProbe() const + std::tuple< Vec, Vec > ConditionProbe() const { - return rand_temp_; + return { std::get< Vec >(rand_temp_), probe_as_drawn_ }; } @@ -304,6 +313,7 @@ namespace bertini{ std::tuple< Vec, Vec > step_temp_; std::tuple< Mat, Mat > J_temp_; std::tuple< Vec, Vec > rand_temp_; // reused scratch: random RHS for norm_J_inverse + Vec probe_as_drawn_; ///< the probe direction at the precision it was drawn at; rand_temp_'s multiprecision half is this at the working precision std::tuple< Vec, Vec > solve_temp_; // reused scratch: LU solve result std::tuple< linalg::PartialPivLU, linalg::PartialPivLU > LU_; diff --git a/core/include/bertini2/trackers/observers.hpp b/core/include/bertini2/trackers/observers.hpp index 4f09f1f7e..c3755ad05 100644 --- a/core/include/bertini2/trackers/observers.hpp +++ b/core/include/bertini2/trackers/observers.hpp @@ -40,37 +40,69 @@ namespace bertini { namespace tracking{ - /// \brief Observer that records the starting precision and the first precision increase during a track. + /** + \brief Observer that records the starting precision and the first precision increase of + each track: each call to TrackPath, from its start. + + A track begins with the tracker's Initializing event, which comes before the start point + is refined, so an increase that the refinement needs counts as the track's first. The + record starts afresh with each track, and the observer stays attached: attach it once and + it reports every track, including one whose initialization fails. + + For a whole path, which the endgame tracks in many calls, see PathPrecisionRecorder. + */ + // Reset on Initializing, never on TrackingStarted, and never unsubscribe: a reset at + // TrackingStarted comes after the refinement and misses its increases, and an observer + // that unsubscribes reports one track forever after (#378). template class FirstPrecisionRecorder : public TypedObserver< FirstPrecisionRecorder, TrackerT, + Initializing::EventEmitterType, + typename TrackerTraits::BaseComplexT>, TrackingStarted::EventEmitterType>, PrecisionChanged::EventEmitterType> > { BOOST_TYPE_INDEX_REGISTER_CLASS using EmitterT = typename TrackerTraits::EventEmitterType; + using BCT = typename TrackerTraits::BaseComplexT; public: - /// \brief On tracking start, record the starting precision. - ObserveResult OnEvent(TrackingStarted const& e) + /// \brief A track begins: the record starts afresh. + ObserveResult OnEvent(Initializing const&) { + have_start_ = false; precision_increased_ = false; - starting_precision_ = e.Get().CurrentPrecision(); + starting_precision_ = 0; + next_precision_ = 0; + time_of_first_increase_ = typename TrackerTraits::BaseComplexT{}; return ObserveResult::KeepObserving; } - /// \brief On a precision increase, record it and unsubscribe (only the first increase is wanted). + /// \brief Tracking starts after the start point is refined: the start precision, unless a refinement change already gave it. + ObserveResult OnEvent(TrackingStarted const& e) + { + if (!have_start_) + { + starting_precision_ = e.Get().CurrentPrecision(); + have_start_ = true; + } + return ObserveResult::KeepObserving; + } + + /// \brief A change of precision: the first one gives the start precision, and the first increase is recorded. ObserveResult OnEvent(PrecisionChanged const& e) { - auto next = e.Next(); - if (next > e.Previous()) + if (!have_start_) + { + starting_precision_ = e.Previous(); + have_start_ = true; + } + if (!precision_increased_ && e.Next() > e.Previous()) { precision_increased_ = true; - next_precision_ = next; + next_precision_ = e.Next(); time_of_first_increase_ = e.Get().CurrentTime(); - // done: ask to be dropped instead of mutating the list mid-dispatch. - return ObserveResult::Unsubscribe; } return ObserveResult::KeepObserving; } @@ -103,41 +135,65 @@ namespace bertini { private: - unsigned starting_precision_; - unsigned next_precision_; - bool precision_increased_; - typename TrackerTraits::BaseComplexT time_of_first_increase_; + bool have_start_ = false; ///< Whether this track's start precision is known yet. + unsigned starting_precision_ = 0; ///< The precision this track started at. + unsigned next_precision_ = 0; ///< The precision after this track's first increase. + bool precision_increased_ = false; ///< Whether precision increased during this track. + typename TrackerTraits::BaseComplexT time_of_first_increase_{}; ///< The time of this track's first increase. }; - /// \brief Observer that records the minimum and maximum precision reached during a track. + /** + \brief Observer that records the minimum and maximum precision of each track: each call + to TrackPath, from its start. + + A track begins with the tracker's Initializing event, before the start point is refined, + so the precision the track started at, and any the refinement moved through, count. + + For a whole path, which the endgame tracks in many calls, see PathPrecisionRecorder. + */ template class MinMaxPrecisionRecorder : public TypedObserver< MinMaxPrecisionRecorder, TrackerT, + Initializing::EventEmitterType, + typename TrackerTraits::BaseComplexT>, TrackingStarted::EventEmitterType>, PrecisionChanged::EventEmitterType> > { BOOST_TYPE_INDEX_REGISTER_CLASS using EmitterT = typename TrackerTraits::EventEmitterType; + using BCT = typename TrackerTraits::BaseComplexT; + + void Include(unsigned p) + { + if (p < min_precision_) + min_precision_ = p; + if (p > max_precision_) + max_precision_ = p; + } public: - /// \brief On tracking start, seed the min and max with the starting precision. + /// \brief A track begins: the record starts afresh. + ObserveResult OnEvent(Initializing const&) + { + min_precision_ = std::numeric_limits::max(); + max_precision_ = 0; + return ObserveResult::KeepObserving; + } + + /// \brief Tracking starts: include the precision it starts at. ObserveResult OnEvent(TrackingStarted const& e) { - min_precision_ = e.Get().CurrentPrecision(); - max_precision_ = e.Get().CurrentPrecision(); + Include(e.Get().CurrentPrecision()); return ObserveResult::KeepObserving; } - /// \brief On a precision change, update the running min and max. + /// \brief A change of precision: include both ends of it. ObserveResult OnEvent(PrecisionChanged const& e) { - auto next_precision = e.Next(); - if (next_precision < min_precision_) - min_precision_ = next_precision; - if (next_precision > max_precision_) - max_precision_ = next_precision; + Include(e.Previous()); + Include(e.Next()); return ObserveResult::KeepObserving; } @@ -170,6 +226,87 @@ namespace bertini { }; + /** + \brief Observer that records the precision a whole path used, across every call to + TrackPath that tracks it: the track to the endgame boundary and every endgame sub-track. + + Unlike the per-track recorders, it never starts afresh on its own: its owner calls Reset + when a path begins, and reads it when the path ends. The zero-dim solver uses it for each + path's precision_changed, time_of_first_prec_increase and max_precision_used. + */ + // A per-track recorder starts afresh at every endgame sub-track, so a precision raised in + // an earlier sub-track and lowered before the last would go unreported (#378). + template + class PathPrecisionRecorder + : public TypedObserver< PathPrecisionRecorder, TrackerT, + TrackingStarted::EventEmitterType>, + PrecisionChanged::EventEmitterType> > + { BOOST_TYPE_INDEX_REGISTER_CLASS + + using EmitterT = typename TrackerTraits::EventEmitterType; + + void Include(unsigned p) + { + if (p < min_precision_) + min_precision_ = p; + if (p > max_precision_) + max_precision_ = p; + } + + public: + + /// \brief Start the record of a new path. + void Reset() + { + precision_increased_ = false; + time_of_first_increase_ = typename TrackerTraits::BaseComplexT{}; + min_precision_ = std::numeric_limits::max(); + max_precision_ = 0; + } + + /// \brief A track of the path starts: include the precision it starts at. + ObserveResult OnEvent(TrackingStarted const& e) + { + Include(e.Get().CurrentPrecision()); + return ObserveResult::KeepObserving; + } + + /// \brief A change of precision: include both ends, and record the path's first increase. + ObserveResult OnEvent(PrecisionChanged const& e) + { + Include(e.Previous()); + Include(e.Next()); + if (!precision_increased_ && e.Next() > e.Previous()) + { + precision_increased_ = true; + time_of_first_increase_ = e.Get().CurrentTime(); + } + return ObserveResult::KeepObserving; + } + + /// \brief Whether precision increased anywhere on the path. + bool DidPrecisionIncrease() const { return precision_increased_; } + + /// \brief The time of the path's first precision increase. + typename TrackerTraits::BaseComplexT TimeOfIncrease() const { return time_of_first_increase_; } + + /// \brief The lowest precision the path used; the maximum unsigned value if it saw none. + unsigned MinPrecision() const { return min_precision_; } + + /// \brief The highest precision the path used; 0 if it saw none. + unsigned MaxPrecision() const { return max_precision_; } + + virtual ~PathPrecisionRecorder() = default; + + private: + + bool precision_increased_ = false; ///< Whether precision increased anywhere on the path. + typename TrackerTraits::BaseComplexT time_of_first_increase_{}; ///< The time of the path's first increase. + unsigned min_precision_ = std::numeric_limits::max(); ///< The lowest precision seen. + unsigned max_precision_ = 0; ///< The highest precision seen. + }; + + /// \brief Observer that collects the precision at every tracking event into a vector. template class PrecisionAccumulator : public Observer diff --git a/core/src/basics/random.cpp b/core/src/basics/random.cpp index 5bfab9a63..e614928b7 100644 --- a/core/src/basics/random.cpp +++ b/core/src/basics/random.cpp @@ -31,6 +31,7 @@ #include #include #include +#include namespace bertini { @@ -44,7 +45,17 @@ inline uint64_t splitmix64(uint64_t x) return x ^ (x >> 31); } +// The global seed, and whether one has been chosen. Every value is a seed, 0 included, so +// "not chosen yet" is its own flag; the mutex only guards the one-time entropy draw. std::atomic g_global_seed{0}; +std::atomic g_global_seed_chosen{false}; +std::mutex g_global_seed_mutex; + +unsigned long EntropySeed() +{ + std::random_device rd; + return static_cast(rd()); +} // Domain tags for stream derivation. Distinct domains (and distinct indices within a domain) // produce distinct engine states, so no two streams ever coincide and no path/worker stream can @@ -82,26 +93,38 @@ std::mt19937& ThreadEngine() { return g_thread_engine; } unsigned long GetGlobalSeed() { - unsigned long s = g_global_seed.load(std::memory_order_relaxed); - if (s == 0) { - std::random_device rd; - s = static_cast(rd()); - if (s == 0) s = 1; - unsigned long expected = 0; - if (!g_global_seed.compare_exchange_strong(expected, s, std::memory_order_relaxed)) - s = g_global_seed.load(std::memory_order_relaxed); + if (g_global_seed_chosen.load(std::memory_order_acquire)) + return g_global_seed.load(std::memory_order_relaxed); + + // never set: choose one from entropy, once, however many threads ask at the same time + std::lock_guard lock(g_global_seed_mutex); + if (!g_global_seed_chosen.load(std::memory_order_relaxed)) + { + g_global_seed.store(EntropySeed(), std::memory_order_relaxed); + g_global_seed_chosen.store(true, std::memory_order_release); } - return s; + return g_global_seed.load(std::memory_order_relaxed); +} + +unsigned long SetGlobalSeedFromEntropy() +{ + unsigned long const seed = EntropySeed(); + SetGlobalSeed(seed); + return seed; +} + +unsigned long ApplyClassicRandomSeed(unsigned long randomseed) +{ + if (randomseed == 0) + return SetGlobalSeedFromEntropy(); + SetGlobalSeed(randomseed); + return randomseed; } void SetGlobalSeed(unsigned long seed) { - if (seed == 0) { - std::random_device rd; - seed = static_cast(rd()); - if (seed == 0) seed = 1; - } g_global_seed.store(seed, std::memory_order_relaxed); + g_global_seed_chosen.store(true, std::memory_order_release); g_solve_ordinal.store(0, std::memory_order_relaxed); // (master, ordinal) restarts here // the setup stream: domain = setup, index = 0. Path/worker streams use other domains, so none // of them can ever reproduce this stream (the old ReseedThisThread(0) == SetGlobalSeed collision). @@ -126,7 +149,7 @@ unsigned long DerivedWorkerSeed(uint64_t worker_index) { uint64_t s = static_cast(GetGlobalSeed()); uint64_t h = splitmix64(s ^ kDomainWorker ^ splitmix64(worker_index)); - if (h == 0) h = 1; // SetGlobalSeed treats 0 as "draw from entropy"; avoid that + if (h == 0) h = 1; // classic input reads randomseed 0 as "draw from entropy"; never hand out 0 return static_cast(h); } @@ -145,7 +168,8 @@ unsigned long DeriveSolveSeed() uint64_t const h = splitmix64(master ^ kDomainSolve ^ splitmix64(ordinal)); unsigned long s = static_cast(h & 0xFFFFFFFFull); if (s == 0) s = static_cast(h >> 32); // deterministic nonzero fallback - if (s == 0) s = 1; // 0 means "entropy" to SetGlobalSeed + if (s == 0) s = 1; // so a recorded seed replays through a classic + // input file, where randomseed 0 means entropy return s; } diff --git a/core/src/blackbox/argc_argv.cpp b/core/src/blackbox/argc_argv.cpp index e4483922a..79a58b257 100644 --- a/core/src/blackbox/argc_argv.cpp +++ b/core/src/blackbox/argc_argv.cpp @@ -53,16 +53,16 @@ ParsedArgs ParseArgcArgv(int argc, char** argv) "\n" "Parallelism:\n" #ifdef BERTINI2_HAVE_MPI - " MPI ranks: mpirun --bind-to none -n N bertini2 [input_file]\n" - " Threads/rank: OMP_NUM_THREADS=T mpirun --bind-to none -n N bertini2 [input_file]\n" - " Total capacity = N ranks x T threads. OMP_NUM_THREADS defaults to 1.\n" - " Note: threading uses std::thread, not OpenMP. OMP_NUM_THREADS is reused\n" - " as the thread-count variable because HPC schedulers (e.g. SLURM) set it\n" - " automatically from --cpus-per-task, so no extra configuration is needed.\n" + " MPI ranks: mpirun -n N bertini2 [input_file]\n" + " Threads/rank: BERTINI_NUM_THREADS=T mpirun --bind-to none -n N bertini2 [input_file]\n" + " Total capacity = N ranks x T threads. Unset, each rank uses every CPU it\n" + " may run on, which depends on how mpirun binds it; with several ranks on\n" + " one machine, set BERTINI_NUM_THREADS so they do not start more threads\n" + " than the machine has cores.\n" #else " This build was compiled without MPI. Rebuild with MPI present for\n" " multi-rank parallelism (it is auto-detected at configure time).\n" - " Thread count: set OMP_NUM_THREADS (default: 1). Uses std::thread, not OpenMP.\n" + " Thread count: BERTINI_NUM_THREADS (default: every CPU this process may run on).\n" #endif ; std::exit(0); diff --git a/core/src/blackbox/main_mode_switch.cpp b/core/src/blackbox/main_mode_switch.cpp index a60bdd080..efed57ce7 100644 --- a/core/src/blackbox/main_mode_switch.cpp +++ b/core/src/blackbox/main_mode_switch.cpp @@ -53,17 +53,19 @@ int RunZeroDim(std::string const& config_str, std::string const& input_str) // Parse and apply the RNG seed before any setup draws (gamma, TD-constants, patch). auto rand_cfg = parsing::classic::FillConfigStruct(config_str); + // In the classic format randomseed 0, the default, means "draw from entropy", as in + // Bertini 1; to the library API 0 is an ordinary seed (ADR-0072). #ifdef BERTINI2_HAVE_MPI - // Manager sets the seed (possibly from entropy), then broadcasts the effective - // (non-zero) seed to workers so all ranks share identical per-path streams. + // Manager sets the seed (possibly from entropy), then broadcasts the effective seed to + // workers so all ranks share identical per-path streams. if (parallel::IsManager()) - SetGlobalSeed(rand_cfg.random_seed); + ApplyClassicRandomSeed(rand_cfg.random_seed); unsigned long effective_seed = GetGlobalSeed(); MPI_Bcast(&effective_seed, 1, MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD); if (!parallel::IsManager()) SetGlobalSeed(effective_seed); #else - SetGlobalSeed(rand_cfg.random_seed); + ApplyClassicRandomSeed(rand_cfg.random_seed); #endif std::cout << "bertini: random seed = " << GetGlobalSeed() << "\n"; diff --git a/core/test/classes/seeded_randomness_test.cpp b/core/test/classes/seeded_randomness_test.cpp index 1f3b6b302..ab0c8a84c 100644 --- a/core/test/classes/seeded_randomness_test.cpp +++ b/core/test/classes/seeded_randomness_test.cpp @@ -103,7 +103,7 @@ BOOST_AUTO_TEST_CASE(derived_solve_seeds_capture_the_stream_deterministically) SetGlobalSeed(42); auto const s1 = DeriveSolveSeed(); auto const s2 = DeriveSolveSeed(); - BOOST_CHECK(s1 != 0u); // 0 means "entropy" to SetGlobalSeed + BOOST_CHECK(s1 != 0u); // 0 means "entropy" to a classic input file BOOST_CHECK(s1 != s2); // consecutive solves get distinct seeds BOOST_CHECK(s1 <= 0xFFFFFFFFul); // 32-bit portable (Windows unsigned long) BOOST_CHECK(s2 <= 0xFFFFFFFFul); @@ -177,6 +177,53 @@ BOOST_AUTO_TEST_CASE(bounded_modulus_factories_are_seeded_and_real_when_asked) BOOST_CHECK(RandomRealBoundedModulus() != r1); // (not the seed-independent orthonormal footgun) } +// ---- seed 0 (#475, ADR-0072) ---- + +// 0 is a seed like any other, as in numpy and Python's random. It used to mean "draw one from +// entropy", so a script that set seed 0 to make a run reproducible got a different run each time. +BOOST_AUTO_TEST_CASE(seed_zero_is_an_ordinary_seed) +{ + SetGlobalSeed(0); + BOOST_CHECK_EQUAL(GetGlobalSeed(), 0ul); + auto const r = RandomRat(); + auto const m = RandomMp(30); + auto const digest = SeededHomotopy(0).ContentDigest().Hex(); + + SetGlobalSeed(0); + BOOST_CHECK_EQUAL(GetGlobalSeed(), 0ul); + BOOST_CHECK_EQUAL(RandomRat(), r); + BOOST_CHECK_EQUAL(RandomMp(30), m); + BOOST_CHECK_EQUAL(SeededHomotopy(0).ContentDigest().Hex(), digest); + + BOOST_CHECK(SeededHomotopy(1).ContentDigest().Hex() != digest); +} + +// Entropy is its own request, and it says which seed it chose, so the run can be reproduced. +BOOST_AUTO_TEST_CASE(a_seed_from_entropy_is_reported_and_reproduces_the_run) +{ + auto const chosen = SetGlobalSeedFromEntropy(); + BOOST_CHECK_EQUAL(GetGlobalSeed(), chosen); + auto const r = RandomRat(); + + SetGlobalSeed(chosen); + BOOST_CHECK_EQUAL(RandomRat(), r); +} + +// The classic input format keeps Bertini 1's meaning: randomseed 0, its default, draws a seed +// from entropy; any other value is the seed. +BOOST_AUTO_TEST_CASE(classic_randomseed_zero_draws_from_entropy) +{ + BOOST_CHECK_EQUAL(ApplyClassicRandomSeed(5), 5ul); + BOOST_CHECK_EQUAL(GetGlobalSeed(), 5ul); + + auto const first = ApplyClassicRandomSeed(0); + BOOST_CHECK_EQUAL(GetGlobalSeed(), first); + auto const second = ApplyClassicRandomSeed(0); + BOOST_CHECK_EQUAL(GetGlobalSeed(), second); + // two draws from entropy agree with probability about 2^-32 + BOOST_CHECK(first != second); +} + // ---- the acceptance test: same seed => digest-identical homotopies ---- BOOST_AUTO_TEST_CASE(same_seed_builds_digest_identical_homotopies) diff --git a/core/test/endgames/generic_cauchy_test.hpp b/core/test/endgames/generic_cauchy_test.hpp index ab5829523..751371219 100644 --- a/core/test/endgames/generic_cauchy_test.hpp +++ b/core/test/endgames/generic_cauchy_test.hpp @@ -1601,6 +1601,91 @@ BOOST_AUTO_TEST_CASE(cauchy_multiple_variables) +/** +An endgame run does not depend on the run before it (#378). + +An endgame is reused path after path -- by the zero-dim solver on each thread, and by anyone +driving one by hand. Its c/k probe direction used to be drawn the first time the endgame needed +it and kept, so a run's direction came from whichever run went first; and a run in double after +one in multiple precision read the earlier run's times. The same run, after another under a +different random stream, must match a fresh endgame's run bit for bit, given the same reseed. +*/ +BOOST_AUTO_TEST_CASE(cauchy_run_does_not_depend_on_the_run_before_it) +{ + DefaultPrecision(ambient_precision); + bertini::SetGlobalSeed(20261005); // ReseedThisThread derives from it; entropy would vary the runs + + bertini::System sys; + Var x = Variable::Make("x"), t = Variable::Make("t"), y = Variable::Make("y"); + sys.AddVariableGroup(VariableGroup{x, y}); + sys.AddPathVariable(t); + sys.AddFunction((pow(x-1,3))*(1-t) + (pow(x,3) + 1)*t); + sys.AddFunction((pow(y-1,2))*(1-t) + (pow(y,2) + 1)*t); + + // A tracker draws its condition-number probe when it is built, so every tracker here is built + // from the same point of the random stream: they differ only in what they have tracked. + auto make_tracker = [&sys]() + { + bertini::ReseedThisThread(7); + TrackerType tracker(sys); + tracker.Setup(TestedPredictor, 1e-5, 1e5, + bertini::tracking::SteppingConfig(), bertini::tracking::NewtonConfig()); + tracker.PrecisionSetup(PrecisionConfig(sys)); + return tracker; + }; + + BCT const boundary_time = ComplexFromString(".1"); + Vec boundary_point(2); + boundary_point << ComplexFromString("5.000000000000001e-01", "9.084258952712920e-17"), + ComplexFromString("9.000000000000001e-01", "4.358898943540673e-01"); + + auto used_tracker = make_tracker(); + TestedEGType used(used_tracker); + used.SetBoundaryTime(boundary_time); + bertini::ReseedThisThread(1); + BOOST_REQUIRE(used.Run(boundary_point) == SuccessCode::Success); + bertini::ReseedThisThread(2); + auto const used_code = used.Run(boundary_point); + + auto fresh_tracker = make_tracker(); + TestedEGType fresh(fresh_tracker); + fresh.SetBoundaryTime(boundary_time); + bertini::ReseedThisThread(2); + auto const fresh_code = fresh.Run(boundary_point); + + BOOST_REQUIRE(fresh_code == SuccessCode::Success); + BOOST_CHECK(used_code == fresh_code); + BOOST_CHECK_EQUAL(used.CycleNumber(), fresh.CycleNumber()); + BOOST_CHECK(used.LatestTime() == fresh.LatestTime()); + BOOST_CHECK_EQUAL(used.ApproximateError(), fresh.ApproximateError()); + auto const& a = used.template FinalApproximation(); + auto const& b = fresh.template FinalApproximation(); + BOOST_REQUIRE_EQUAL(a.size(), b.size()); + for (Eigen::Index ii = 0; ii < a.size(); ++ii) + { + BOOST_CHECK_EQUAL(Precision(a(ii)), Precision(b(ii))); + BOOST_CHECK(a(ii) == b(ii)); + } + + // The c/k estimate reads the samples only through the probe direction, so on the same samples + // the two endgames agree exactly when their directions do. The run above need not show a + // difference in direction: on this path the estimate tips no decision. The successive + // differences, (1, 0) and (0.9, 0.1), are not parallel, so the estimate depends on the + // direction; and their ratio along any unit-modulus direction lies in [0.8, 1], so it is never + // clipped to 1. + bertini::SampCont samples; + for (auto const& [first, second] : {std::pair{"0", "0"}, + {"1", "0"}, {"1.9", "0.1"}}) + { + Vec s(2); + s << ComplexFromString(first), ComplexFromString(second); + samples.push_back(s); + } + BOOST_CHECK(used.ComputeCOverK(samples) == fresh.ComputeCOverK(samples)); +} + + + /** Test to see if we can compute Cauchy samples around a nonzero target time. */ diff --git a/core/test/endgames/generic_pseg_test.hpp b/core/test/endgames/generic_pseg_test.hpp index 576568f16..e5c4246c4 100644 --- a/core/test/endgames/generic_pseg_test.hpp +++ b/core/test/endgames/generic_pseg_test.hpp @@ -1124,6 +1124,86 @@ BOOST_AUTO_TEST_CASE(full_run_multiple_variables) +/** +A power-series run does not depend on the run before it (#378). + +The same run, after another under a different random stream, must match a fresh endgame's run +bit for bit, given the same reseed. See cauchy_run_does_not_depend_on_the_run_before_it: the +c/k probe direction and the sample times are shared state of every endgame flavor, and a run +in double after one in multiple precision read the earlier run's times. +*/ +BOOST_AUTO_TEST_CASE(power_series_run_does_not_depend_on_the_run_before_it) +{ + DefaultPrecision(ambient_precision); + bertini::SetGlobalSeed(20261005); // ReseedThisThread derives from it; entropy would vary the runs + + bertini::System sys; + Var x = Variable::Make("x"), t = Variable::Make("t"), y = Variable::Make("y"); + sys.AddVariableGroup(VariableGroup{x, y}); + sys.AddPathVariable(t); + sys.AddFunction((pow(x-1,3))*(1-t) + (pow(x,3) + 1)*t); + sys.AddFunction((pow(y-1,2))*(1-t) + (pow(y,2) + 1)*t); + + // A tracker draws its condition-number probe when it is built, so every tracker here is built + // from the same point of the random stream: they differ only in what they have tracked. + auto make_tracker = [&sys]() + { + bertini::ReseedThisThread(7); + TrackerType tracker(sys); + tracker.Setup(TestedPredictor, 1e-6, 1e5, + bertini::tracking::SteppingConfig(), bertini::tracking::NewtonConfig()); + tracker.PrecisionSetup(PrecisionConfig(sys)); + return tracker; + }; + + BCT const boundary_time = ComplexFromString(".1"); + Vec boundary_point(2); + boundary_point << ComplexFromString("5.000000000000001e-01", "9.084258952712920e-17"), + ComplexFromString("9.000000000000001e-01", "4.358898943540673e-01"); + + auto used_tracker = make_tracker(); + TestedEGType used(used_tracker); + used.SetBoundaryTime(boundary_time); + bertini::ReseedThisThread(1); + BOOST_REQUIRE(used.Run(boundary_point) == SuccessCode::Success); + bertini::ReseedThisThread(2); + auto const used_code = used.Run(boundary_point); + + auto fresh_tracker = make_tracker(); + TestedEGType fresh(fresh_tracker); + fresh.SetBoundaryTime(boundary_time); + bertini::ReseedThisThread(2); + auto const fresh_code = fresh.Run(boundary_point); + + BOOST_REQUIRE(fresh_code == SuccessCode::Success); + BOOST_CHECK(used_code == fresh_code); + BOOST_CHECK_EQUAL(used.CycleNumber(), fresh.CycleNumber()); + BOOST_CHECK(used.LatestTime() == fresh.LatestTime()); + BOOST_CHECK_EQUAL(used.ApproximateError(), fresh.ApproximateError()); + auto const& a = used.template FinalApproximation(); + auto const& b = fresh.template FinalApproximation(); + BOOST_REQUIRE_EQUAL(a.size(), b.size()); + for (Eigen::Index ii = 0; ii < a.size(); ++ii) + { + BOOST_CHECK_EQUAL(Precision(a(ii)), Precision(b(ii))); + BOOST_CHECK(a(ii) == b(ii)); + } + + // The c/k estimate reads the samples only through the probe direction; see the same check in + // cauchy_run_does_not_depend_on_the_run_before_it for why these samples. + bertini::SampCont samples; + for (auto const& [first, second] : {std::pair{"0", "0"}, + {"1", "0"}, {"1.9", "0.1"}}) + { + Vec s(2); + s << ComplexFromString(first), ComplexFromString(second); + samples.push_back(s); + } + BOOST_CHECK(used.ComputeCOverK(samples) == fresh.ComputeCOverK(samples)); +} + + + /** The function that runs the power series endgame is called PSEG. PSEG takes an endgame_time value and an endgame_space value that is diff --git a/core/test/nag_algorithms/threaded_solve.cpp b/core/test/nag_algorithms/threaded_solve.cpp index 31287a146..ba9e7719b 100644 --- a/core/test/nag_algorithms/threaded_solve.cpp +++ b/core/test/nag_algorithms/threaded_solve.cpp @@ -28,15 +28,21 @@ as a serial one -- threading is a performance feature, never a correctness chang on every platform (no MPI required), which is what gives us cross-OS threading coverage. */ +#include #include +#include +#include #include #include #include +#include +#include #include #include #include "bertini2/parallel/thread_pool.hpp" +#include "bertini2/random.hpp" #include "bertini2/nag_algorithms/zero_dim_solve.hpp" #include "bertini2/endgames.hpp" #include "bertini2/system/start_systems.hpp" @@ -121,16 +127,72 @@ BOOST_AUTO_TEST_CASE(thread_pool_state_is_per_thread) BOOST_CHECK_GE(total, static_cast(N)); } -BOOST_AUTO_TEST_CASE(effective_thread_count_is_sane) +namespace { +// Set (value != nullptr) or unset an environment variable for the life of the object, and put +// back whatever was there. POSIX setenv/unsetenv, or _putenv_s on Windows (where an empty value +// unsets). +class ScopedEnv { - using bertini::parallel::EffectiveThreadCount; - // An explicit positive request is honored (when OMP_NUM_THREADS is unset in the test env). - if (std::getenv("OMP_NUM_THREADS") == nullptr) + std::string name_; + bool had_ = false; + std::string old_; + + static void Put(std::string const& name, char const* value) + { +#ifdef _WIN32 + _putenv_s(name.c_str(), value ? value : ""); +#else + if (value) + setenv(name.c_str(), value, 1); + else + unsetenv(name.c_str()); +#endif + } + +public: + ScopedEnv(std::string name, char const* value) : name_(std::move(name)) { - BOOST_CHECK_EQUAL(EffectiveThreadCount(1), 1u); - BOOST_CHECK_EQUAL(EffectiveThreadCount(3), 3u); - BOOST_CHECK_GE(EffectiveThreadCount(0), 1u); // auto -> hardware_concurrency, clamped >= 1 + if (char const* prev = std::getenv(name_.c_str())) + { + had_ = true; + old_ = prev; + } + Put(name_, value); } + ~ScopedEnv() { Put(name_, had_ ? old_.c_str() : nullptr); } + ScopedEnv(ScopedEnv const&) = delete; + ScopedEnv& operator=(ScopedEnv const&) = delete; +}; +} // namespace + +BOOST_AUTO_TEST_CASE(effective_thread_count_is_sane) +{ + using bertini::parallel::EffectiveThreadCount; + ScopedEnv no_override("BERTINI_NUM_THREADS", nullptr); + // an explicit positive request is honored; auto is every available CPU, at least one + BOOST_CHECK_EQUAL(EffectiveThreadCount(1), 1u); + BOOST_CHECK_EQUAL(EffectiveThreadCount(3), 3u); + BOOST_CHECK_EQUAL(EffectiveThreadCount(0), bertini::parallel::AvailableCpuCount()); + BOOST_CHECK_GE(bertini::parallel::AvailableCpuCount(), 1u); +} + +BOOST_AUTO_TEST_CASE(bertini_num_threads_overrides_the_configured_count) +{ + using bertini::parallel::EffectiveThreadCount; + ScopedEnv two("BERTINI_NUM_THREADS", "2"); + BOOST_CHECK_EQUAL(EffectiveThreadCount(3), 2u); + BOOST_CHECK_EQUAL(EffectiveThreadCount(0), 2u); +} + +// b2 reads its own variable. OMP_NUM_THREADS also sets numpy's OpenBLAS thread count, and until +// 4.0 b2 read it too, coupling two unrelated settings; it must not reach b2 any more. +BOOST_AUTO_TEST_CASE(omp_num_threads_does_not_reach_b2) +{ + using bertini::parallel::EffectiveThreadCount; + ScopedEnv no_override("BERTINI_NUM_THREADS", nullptr); + ScopedEnv omp("OMP_NUM_THREADS", "1"); + BOOST_CHECK_EQUAL(EffectiveThreadCount(3), 3u); + BOOST_CHECK_EQUAL(EffectiveThreadCount(0), bertini::parallel::AvailableCpuCount()); } BOOST_AUTO_TEST_SUITE_END() // thread_pool @@ -320,16 +382,16 @@ BOOST_AUTO_TEST_CASE(event_tracker_is_member_tracker_when_serial) BOOST_AUTO_TEST_CASE(event_tracker_is_a_clone_when_threaded) { using namespace bertini; - // OMP_NUM_THREADS overrides the requested thread count (parallel::EffectiveThreadCount), so under - // OMP_NUM_THREADS=1 a "threaded" solve actually runs serially on the member tracker, and the - // clone-per-thread expectation below does not hold. CI runs the test suite both with and without - // OMP_NUM_THREADS=1; the threaded path is exercised in the non-serial leg, so skip here when the - // environment forces serial. + // BERTINI_NUM_THREADS overrides the requested thread count (parallel::EffectiveThreadCount), so + // under BERTINI_NUM_THREADS=1 a "threaded" solve actually runs serially on the member tracker, + // and the clone-per-thread expectation below does not hold. CI runs the test suite both with + // and without BERTINI_NUM_THREADS=1; the threaded path is exercised in the non-serial leg, so + // skip here when the environment forces serial. if (parallel::EffectiveThreadCount(4) <= 1) { - BOOST_TEST_MESSAGE("event_tracker_is_a_clone_when_threaded: skipped -- OMP_NUM_THREADS forces " - "a serial run (threaded clone behavior is covered by the OMP_NUM_THREADS!=1 " - "CI leg)."); + BOOST_TEST_MESSAGE("event_tracker_is_a_clone_when_threaded: skipped -- BERTINI_NUM_THREADS " + "forces a serial run (threaded clone behavior is covered by the " + "BERTINI_NUM_THREADS!=1 CI leg)."); return; } auto sys = TwoCubics(); @@ -349,4 +411,212 @@ BOOST_AUTO_TEST_CASE(event_tracker_is_a_clone_when_threaded) BOOST_CHECK(cap.trackers_seen.size() <= 4u); // at most one clone per worker thread } + +namespace { + +// cyclic-5: 120 total-degree paths, 70 finite roots, the rest diverging. Enough paths, and +// enough precision changes, that every thread count hands each thread a different sequence of +// predecessor paths. +bertini::System Cyclic5() +{ + using namespace bertini; + std::vector> z; + for (int i = 0; i < 5; ++i) + z.push_back(Variable::Make("z" + std::to_string(i))); + System sys; + using Nd = std::shared_ptr; + for (int k = 1; k < 5; ++k) + { + Nd total; + for (int i = 0; i < 5; ++i) + { + Nd term = z[i]; + for (int j = 1; j < k; ++j) + term = term * z[(i + j) % 5]; + total = total ? total + term : term; + } + sys.AddFunction(total); + } + sys.AddFunction(z[0] * z[1] * z[2] * z[3] * z[4] - 1); + sys.AddVariableGroup(VariableGroup(z.begin(), z.end())); + return sys; +} + +// Every path's endpoint and its whole metadata record, from a seeded adaptive-precision solve. +template +struct PathByPath +{ + using ZD = bertini::algorithm::ZeroDimSolver; + std::vector> endpoints; + std::vector> metadata; + + explicit PathByPath(unsigned num_threads) + { + bertini::SetGlobalSeed(20261005); + auto sys = Cyclic5(); + ZD zd(sys); + zd.DefaultSetup(); + auto cfg = zd.template Get(); + cfg.num_threads = num_threads; + zd.Set(cfg); + zd.Solve(); + endpoints = zd.SolutionsInternalCoords(); + metadata = zd.SolutionMetadata(); + } +}; + +template +void CheckBitIdenticalAcrossThreadCounts() +{ + PathByPath const serial(1); + BOOST_REQUIRE_EQUAL(serial.endpoints.size(), 120u); + + for (unsigned nt : {2u, 3u, 8u}) + { + PathByPath const threaded(nt); + BOOST_REQUIRE_EQUAL(threaded.endpoints.size(), serial.endpoints.size()); + for (std::size_t i = 0; i < serial.endpoints.size(); ++i) + { + BOOST_TEST_CONTEXT("path " << i << " at " << nt << " threads") + { + auto const& a = serial.endpoints[i]; + auto const& b = threaded.endpoints[i]; + BOOST_REQUIRE_EQUAL(a.size(), b.size()); + for (Eigen::Index j = 0; j < a.size(); ++j) + { + BOOST_CHECK_EQUAL(a(j).precision(), b(j).precision()); + BOOST_CHECK(a(j) == b(j)); + } + auto md_a = serial.metadata[i]; // operator== is not const + BOOST_CHECK_MESSAGE(md_a == threaded.metadata[i], + "metadata differs:\n" << serial.metadata[i] + << "\nversus\n" << threaded.metadata[i]); + } + } + } +} + +} // namespace + + +/** +A threaded solve is bit-identical to the serial solve, path by path (#378). + +A path's result is a function of the systems, the settings, the seed and its index -- not of +the thread that ran it, nor of the paths that thread ran before it. Endpoints are compared +digit for digit and precision for precision, and every metadata field (step counts, precision +record, condition number, residuals, accuracy estimates) exactly; only the wall-clock time is +excluded. Matching within a tolerance, as the tests above do, cannot see this: before the fix +the endpoints of a threaded solve agreed with the serial ones to many digits and still differed. + +Under BERTINI_NUM_THREADS=1 every one of these solves is serial and the test proves nothing; +the CI leg without it is the one that counts. +*/ +BOOST_AUTO_TEST_CASE(threaded_solve_is_bit_identical_to_serial_power_series) +{ + using TrackerA = bertini::tracking::AMPTracker; + CheckBitIdenticalAcrossThreadCounts::PSEG>(); +} + +BOOST_AUTO_TEST_CASE(threaded_solve_is_bit_identical_to_serial_cauchy) +{ + using TrackerA = bertini::tracking::AMPTracker; + CheckBitIdenticalAcrossThreadCounts::Cauchy>(); +} + + +namespace { + +// The highest precision the watched tracker reports, in any event: the start of a track, and +// both ends of a change of precision. +struct TrackerPrecisionWatch : public bertini::Observer +{ + unsigned highest = 0; + + bertini::ObserveResult Observe(bertini::AnyEvent const& e) override + { + using namespace bertini::tracking; + if (auto s = dynamic_cast const*>(&e)) + highest = std::max(highest, s->Get().CurrentPrecision()); + else if (auto c = dynamic_cast const*>(&e)) + highest = std::max({highest, c->Previous(), c->Next()}); + return bertini::ObserveResult::KeepObserving; + } +}; + +// Starts the watch afresh at each path, and keeps its reading when the path completes. +struct PerPathPrecision : public bertini::Observer +{ + TrackerPrecisionWatch& watch; + std::map highest_by_path; + + explicit PerPathPrecision(TrackerPrecisionWatch& w) : watch(w) {} + + bertini::ObserveResult Observe(bertini::AnyEvent const& e) override + { + using namespace bertini::algorithm; + if (dynamic_cast const*>(&e)) + watch.highest = 0; + else if (auto c = dynamic_cast const*>(&e)) + highest_by_path[c->PathIndex()] = watch.highest; + return bertini::ObserveResult::KeepObserving; + } +}; + +template +void CheckMaxPrecisionUsedCoversEveryTrack() +{ + using namespace bertini; + // serial, so every path runs on the solver's own tracker, where the watch is attached + thread_pool::ScopedEnv serial("BERTINI_NUM_THREADS", "1"); + + SetGlobalSeed(20261005); + auto sys = Cyclic5(); + algorithm::ZeroDimSolver zd(sys); + zd.DefaultSetup(); + + TrackerPrecisionWatch watch; + PerPathPrecision per_path(watch); + zd.GetTracker().AddObserver(watch); + zd.AddObserver(per_path); + zd.Solve(); + + auto const& md = zd.SolutionMetadata(); + BOOST_REQUIRE_EQUAL(per_path.highest_by_path.size(), md.size()); + bool some_path_raised_precision = false; + for (std::size_t i = 0; i < md.size(); ++i) + { + BOOST_TEST_CONTEXT("path " << i) + { + BOOST_CHECK_EQUAL(md[i].max_precision_used, per_path.highest_by_path.at(i)); + some_path_raised_precision = some_path_raised_precision || md[i].max_precision_used > 16; + } + } + BOOST_CHECK(some_path_raised_precision); // the premise: the check is not about doubles alone +} + +} // namespace + + +/** +A path's max_precision_used is the highest precision any of its tracks used (#378 follow-up). + +The endgame tracks a path in many calls to TrackPath. The solver's record of a path's precision +used to start afresh with every call, so a precision raised in an earlier endgame sub-track and +lowered before the last went unreported. Checked against every precision the tracker reports +while tracking each path, on a serial solve so that the solver's own tracker runs every path. +*/ +BOOST_AUTO_TEST_CASE(max_precision_used_covers_every_track_of_a_path_power_series) +{ + using TrackerA = bertini::tracking::AMPTracker; + CheckMaxPrecisionUsedCoversEveryTrack::PSEG>(); +} + +BOOST_AUTO_TEST_CASE(max_precision_used_covers_every_track_of_a_path_cauchy) +{ + using TrackerA = bertini::tracking::AMPTracker; + CheckMaxPrecisionUsedCoversEveryTrack::Cauchy>(); +} + BOOST_AUTO_TEST_SUITE_END() diff --git a/core/test/tracking_basics/amp_tracker_test.cpp b/core/test/tracking_basics/amp_tracker_test.cpp index 61b1c1334..0e4ce5702 100644 --- a/core/test/tracking_basics/amp_tracker_test.cpp +++ b/core/test/tracking_basics/amp_tracker_test.cpp @@ -25,8 +25,10 @@ #include #include +#include "bertini2/random.hpp" #include "bertini2/system/start_systems.hpp" #include "bertini2/trackers/tracker.hpp" +#include "bertini2/trackers/observers.hpp" #include "bertini2/records/solver_recording.hpp" // CanonicalName, for readable test messages @@ -1158,6 +1160,306 @@ BOOST_AUTO_TEST_CASE(set_start_precision_overrides_to_higher_precision) BOOST_CHECK(abs(y_end(0) - mpfr(0)) < 1e-5); } +namespace { + +// Tracked from all ones at t = 1. Cube root: y in (1-t)(y^3 - 2) + t(y^3 - 1), to the cube +// root of 2 at t = 0, a smooth path that drops to double on its own. Square root: (x, y) in +// x^2 - t, y - x, toward the double root at t = 0; the Jacobian's condition grows like 1/x, so +// a track to t = 1e-20 at tracking tolerance 1e-12 ends in multiple precision (30 digits), +// where the condition estimate uses the multiprecision probe direction. (One unknown would +// not do: a 1x1 Jacobian has condition number 1. And at a loose tolerance the corrector +// accepts points well off the path near the root, so the track never needs the digits.) +struct TestHomotopy +{ + enum class Kind { CubeRoot, SquareRoot }; + + Var x = Variable::Make("x"); + Var y = Variable::Make("y"); + Var t = Variable::Make("t"); + System sys; + double tracking_tolerance; + mpfr t_end; + + explicit TestHomotopy(Kind kind = Kind::CubeRoot) + : tracking_tolerance(kind == Kind::CubeRoot ? 1e-5 : 1e-12), + t_end(kind == Kind::CubeRoot ? mpfr(0) : mpfr("1e-20")) + { + if (kind == Kind::CubeRoot) + { + sys.AddFunction((1 - t) * (pow(y, 3) - 2) + t * (pow(y, 3) - 1)); + sys.AddVariableGroup(VariableGroup{y}); + } + else + { + sys.AddFunction(pow(x, 2) - t); + sys.AddFunction(y - x); + sys.AddVariableGroup(VariableGroup{x, y}); + } + sys.AddPathVariable(t); + } + + // Every tracker made here draws the same condition-number probe direction: the same seed, + // at the same precision. + bertini::tracking::AMPTracker MakeTracker() const + { + using namespace bertini::tracking; + DefaultPrecision(30); + bertini::ReseedThisThread(20261005); + AMPTracker tracker(sys); + tracker.Setup(Predictor::RKF45, tracking_tolerance, 1e5, SteppingConfig(), NewtonConfig()); + tracker.PrecisionSetup(AMPConfigFrom(sys)); + return tracker; + } +}; + +// Everything a track leaves behind that a caller can read: the code, the step count, the +// precision, the step size, the latest condition-number estimate (which is computed from the +// probe direction) and the endpoint, compared bit for bit and precision for precision. +struct TrackOutcome +{ + bertini::SuccessCode code; + unsigned steps; + unsigned precision; + real_mp stepsize; + double condition_number; + Vec endpoint; +}; + +TrackOutcome Track(bertini::tracking::AMPTracker const& tracker, unsigned start_digits, + mpfr const& t_end) +{ + DefaultPrecision(start_digits); + Vec start(tracker.NumVariables()); + for (Eigen::Index ii = 0; ii < start.size(); ++ii) + start(ii) = mpfr(1); + Vec endpoint; + auto code = tracker.TrackPath(endpoint, mpfr(1), t_end, start); + return { code, tracker.NumTotalStepsTaken(), tracker.CurrentPrecision(), + tracker.CurrentStepsize(), static_cast(tracker.LatestConditionNumber()), + endpoint }; +} + +void CheckIdentical(TrackOutcome const& a, TrackOutcome const& b) +{ + BOOST_CHECK(a.code == b.code); + BOOST_CHECK_EQUAL(a.steps, b.steps); + BOOST_CHECK_EQUAL(a.precision, b.precision); + BOOST_CHECK_EQUAL(a.condition_number, b.condition_number); + BOOST_CHECK_EQUAL(a.stepsize.precision(), b.stepsize.precision()); + BOOST_CHECK(a.stepsize == b.stepsize); + BOOST_REQUIRE_EQUAL(a.endpoint.size(), b.endpoint.size()); + for (Eigen::Index ii = 0; ii < a.endpoint.size(); ++ii) + { + BOOST_CHECK_EQUAL(a.endpoint(ii).precision(), b.endpoint(ii).precision()); + BOOST_CHECK(a.endpoint(ii) == b.endpoint(ii)); + } +} + +} // namespace + + +/** +A track's result does not depend on what the tracker tracked before it (#378). + +A tracker is reused path after path -- by the zero-dim solver on each thread, by every +endgame, and by anyone driving one by hand. Until 4.0 the first step of a track was built at +the precision the PREVIOUS track ended in, and kept that precision, so the same path tracked +by the same tracker gave different bits depending on its predecessor; a threaded solve, where +the scheduler picks each thread's predecessors, was irreproducible. Every pairing of where the +previous track ended and where this one starts: above, at and below it, double and not. +*/ +BOOST_AUTO_TEST_CASE(a_track_does_not_depend_on_the_track_before_it) +{ + using Kind = TestHomotopy::Kind; + // (digits the previous track starts at, digits this track starts at) + std::vector> const pairings{ + {60, 16}, {60, 30}, {16, 30}, {30, 16}, {30, 60}, {40, 30}, + // below the 30 digits the probe direction was drawn at: a drop that rounds it + {20, 30}, {20, 60}, {20, 16}}; + + // the cube root ends in double; the square root, stopped just short of its double root, + // ends in multiple precision + for (auto const kind : {Kind::CubeRoot, Kind::SquareRoot}) + { + TestHomotopy h(kind); + mpfr const& t_end = h.t_end; + + for (auto const& [before_digits, digits] : pairings) + { + BOOST_TEST_CONTEXT((kind == Kind::CubeRoot ? "cube root" : "square root") + << ", previous track starts at " << before_digits + << " digits, this one at " << digits) + { + auto used = h.MakeTracker(); + // preserved, the previous track ends where it started -- elsewhere than this one + used.PrecisionPreservation(true); + auto const before = Track(used, before_digits, mpfr(0.5)); + used.PrecisionPreservation(false); + BOOST_REQUIRE(before.code == bertini::SuccessCode::Success); + BOOST_REQUIRE_NE(before.precision, digits); + + auto const after_another = Track(used, digits, t_end); + + auto fresh = h.MakeTracker(); + auto const first_ever = Track(fresh, digits, t_end); + + BOOST_CHECK(first_ever.code == bertini::SuccessCode::Success); + if (kind == Kind::SquareRoot) // the premise of this half + BOOST_CHECK_GT(first_ever.precision, 16u); + CheckIdentical(after_another, first_ever); + } + } + } +} + + +/** +The step size is held at the tracker's working precision. + +An assignment takes its source's precision (preserve_related_precision), so a step size +assigned from a number of another precision carries that precision into every product it +enters. Checked after tracks that end in double and in multiple precision, each after a +predecessor at another precision. +*/ +BOOST_AUTO_TEST_CASE(the_step_size_is_at_the_working_precision_after_a_track) +{ + using Kind = TestHomotopy::Kind; + for (auto const kind : {Kind::CubeRoot, Kind::SquareRoot}) + { + TestHomotopy h(kind); + mpfr const& t_end = h.t_end; + auto tracker = h.MakeTracker(); + for (unsigned digits : {60u, 30u, 16u, 40u, 16u, 20u}) + { + BOOST_TEST_CONTEXT((kind == Kind::CubeRoot ? "cube root" : "square root") + << ", starting at " << digits << " digits") + { + auto const outcome = Track(tracker, digits, t_end); + BOOST_REQUIRE(outcome.code == bertini::SuccessCode::Success); + BOOST_CHECK_EQUAL(outcome.stepsize.precision(), outcome.precision); + } + } + } +} + + +namespace { + +// x^2 + (1-t)x, y^2 + (1-t)y. At t = 1 the start point (0, 0) is a double root: its refinement +// raises the precision until it gives up with SingularStartPoint, so the track fails during its +// initialization. At t = 1/2 the point (-1/2, -1/2) is a regular root, tracked to t = 1/4. +struct SingularAndRegularStarts +{ + Var x = Variable::Make("x"); + Var y = Variable::Make("y"); + Var t = Variable::Make("t"); + System sys; + + SingularAndRegularStarts() + { + sys.AddFunction(pow(x, 2) + (1 - t) * x); + sys.AddFunction(pow(y, 2) + (1 - t) * y); + sys.AddPathVariable(t); + sys.AddVariableGroup(VariableGroup{x, y}); + } + + bertini::tracking::AMPTracker MakeTracker() const + { + using namespace bertini::tracking; + AMPTracker tracker(sys); + tracker.Setup(Predictor::RKF45, 1e-5, 1e5, SteppingConfig(), NewtonConfig()); + tracker.PrecisionSetup(AMPConfigFrom(sys)); + return tracker; + } + + static bertini::SuccessCode TrackSingular(bertini::tracking::AMPTracker const& tracker) + { + DefaultPrecision(30); + Vec start(2); + start << mpfr(0), mpfr(0); + Vec end; + return tracker.TrackPath(end, mpfr(1), mpfr(0), start); + } + + static bertini::SuccessCode TrackRegular(bertini::tracking::AMPTracker const& tracker) + { + DefaultPrecision(16); + Vec start(2); + start << mpfr("-0.5"), mpfr("-0.5"); + Vec end; + return tracker.TrackPath(end, mpfr("0.5"), mpfr("0.25"), start); + } +}; + +} // namespace + + +/** +A precision recorder attached once reports every track. + +Each recorder starts afresh when a track begins and stays attached, so after a track whose +start-point refinement raised the precision, the next track, which raises nothing, says so. +*/ +BOOST_AUTO_TEST_CASE(a_precision_recorder_attached_once_reports_every_track) +{ + using namespace bertini::tracking; + SingularAndRegularStarts h; + auto tracker = h.MakeTracker(); + FirstPrecisionRecorder first; + MinMaxPrecisionRecorder min_max; + tracker.AddObserver(first); + tracker.AddObserver(min_max); + + BOOST_REQUIRE(SingularAndRegularStarts::TrackSingular(tracker) == bertini::SuccessCode::SingularStartPoint); + BOOST_CHECK(first.DidPrecisionIncrease()); + BOOST_CHECK_EQUAL(first.StartPrecision(), 30u); + BOOST_CHECK_EQUAL(min_max.MinPrecision(), 30u); + BOOST_CHECK_GT(min_max.MaxPrecision(), 30u); + + BOOST_REQUIRE(SingularAndRegularStarts::TrackRegular(tracker) == bertini::SuccessCode::Success); + BOOST_CHECK(!first.DidPrecisionIncrease()); + BOOST_CHECK_EQUAL(first.StartPrecision(), 16u); + BOOST_CHECK_EQUAL(min_max.MinPrecision(), 16u); + BOOST_CHECK_EQUAL(min_max.MaxPrecision(), 16u); +} + + +/** +A track whose initialization fails is recorded as itself, not as the track before it. + +The start-point refinement comes before TrackingStarted, and a failed refinement ends the track +before TrackingStarted is ever sent. The record starts at the track's Initializing event, so a +recorder that saw an earlier track reports the failed one exactly as a fresh recorder does. +*/ +BOOST_AUTO_TEST_CASE(a_track_that_fails_to_initialize_is_recorded_as_itself) +{ + using namespace bertini::tracking; + SingularAndRegularStarts h; + + auto used = h.MakeTracker(); + FirstPrecisionRecorder used_first; + MinMaxPrecisionRecorder used_min_max; + used.AddObserver(used_first); + used.AddObserver(used_min_max); + BOOST_REQUIRE(SingularAndRegularStarts::TrackRegular(used) == bertini::SuccessCode::Success); + BOOST_REQUIRE(SingularAndRegularStarts::TrackSingular(used) == bertini::SuccessCode::SingularStartPoint); + + auto fresh = h.MakeTracker(); + FirstPrecisionRecorder fresh_first; + MinMaxPrecisionRecorder fresh_min_max; + fresh.AddObserver(fresh_first); + fresh.AddObserver(fresh_min_max); + BOOST_REQUIRE(SingularAndRegularStarts::TrackSingular(fresh) == bertini::SuccessCode::SingularStartPoint); + + BOOST_CHECK_EQUAL(used_first.DidPrecisionIncrease(), fresh_first.DidPrecisionIncrease()); + BOOST_CHECK_EQUAL(used_first.StartPrecision(), fresh_first.StartPrecision()); + BOOST_CHECK_EQUAL(used_first.NextPrecision(), fresh_first.NextPrecision()); + BOOST_CHECK(used_first.TimeOfIncrease() == fresh_first.TimeOfIncrease()); + BOOST_CHECK_EQUAL(used_min_max.MinPrecision(), fresh_min_max.MinPrecision()); + BOOST_CHECK_EQUAL(used_min_max.MaxPrecision(), fresh_min_max.MaxPrecision()); +} + + BOOST_AUTO_TEST_CASE(set_start_precision_clear_restores_default_behavior) { DefaultPrecision(30); diff --git a/core/test/tracking_basics/newton_correct_test.cpp b/core/test/tracking_basics/newton_correct_test.cpp index 9a5d3d958..e136622ce 100644 --- a/core/test/tracking_basics/newton_correct_test.cpp +++ b/core/test/tracking_basics/newton_correct_test.cpp @@ -712,4 +712,43 @@ BOOST_AUTO_TEST_CASE(circle_line_one_corrector_step_double) BOOST_CHECK(success_code==bertini::SuccessCode::FailedToConverge); } + +/** +The condition-number probe direction survives a drop in precision. + +The probe is drawn once per path and held fixed for the whole track, while the working +precision rises and falls under it. It used to be rounded in place at each change, so a drop +to fewer digits destroyed digits that the next rise could not restore, and the direction a +path used depended on the precisions it had passed through -- for a tracker reused path after +path, on the paths before it. The direction is stored as drawn, and every working copy is +taken from that. +*/ +BOOST_AUTO_TEST_CASE(the_probe_direction_survives_a_drop_in_precision) +{ + DefaultPrecision(50); + bertini::System sys; + Var x = Variable::Make("x"), y = Variable::Make("y"), t = Variable::Make("t"); + sys.AddVariableGroup(VariableGroup{x, y}); + sys.AddPathVariable(t); + sys.AddFunction(t * (pow(x, 2) - 1) + (1 - t) * (pow(x, 2) + pow(y, 2) - 4)); + sys.AddFunction(t * (y - 1) + (1 - t) * (2 * x + 5 * y)); + + NewtonCorrector corrector(sys); + corrector.ChangePrecision(50); + corrector.RefreshRandomDirection(); + Vec const drawn = std::get>(corrector.ConditionProbe()); + + for (unsigned digits : {20u, 60u, 30u, 50u}) + { + corrector.ChangePrecision(digits); + Vec const now = std::get>(corrector.ConditionProbe()); + BOOST_REQUIRE_EQUAL(now.size(), drawn.size()); + for (Eigen::Index ii = 0; ii < drawn.size(); ++ii) + { + BOOST_CHECK_EQUAL(now(ii).precision(), drawn(ii).precision()); + BOOST_CHECK(now(ii) == drawn(ii)); + } + } +} + BOOST_AUTO_TEST_SUITE_END() diff --git a/docs/adr/0070-b2-reads-its-own-thread-count-variable.md b/docs/adr/0070-b2-reads-its-own-thread-count-variable.md new file mode 100644 index 000000000..1f94f23ea --- /dev/null +++ b/docs/adr/0070-b2-reads-its-own-thread-count-variable.md @@ -0,0 +1,51 @@ +# ADR-0070: b2 reads its own thread-count variable, and defaults to every available CPU + +**Status:** Accepted +**Date:** 2026-10-05 + +## Context + +b2 tracks paths on a pool of `std::thread` workers. It uses no OpenMP, but until 4.0 it took +its thread count from `OMP_NUM_THREADS`, because that variable is familiar and because SLURM +sets it from `--cpus-per-task`. The default differed by entry point: a standalone solve used +every CPU, and an MPI rank used one thread. + +`OMP_NUM_THREADS` is not b2's alone. numpy's OpenBLAS reads it when `OPENBLAS_NUM_THREADS` is +unset, and so does any OpenMP library in the same process. Two things followed. + +- Pinning b2 to one thread for a reproducible run pinned numpy's linear algebra as well, and + pinning numpy pinned b2. An experiment that wanted one and not the other could not say so. +- A user who set the variable for one library changed the other without knowing it. + +## Decision + +b2 reads **`BERTINI_NUM_THREADS`** and ignores `OMP_NUM_THREADS`. One function, +`parallel::EffectiveThreadCount`, decides the count for a standalone solve and for each MPI +rank, in this order: + +1. `BERTINI_NUM_THREADS`, when set to a positive integer; +2. the solver's `ZeroDimConfig::num_threads`, when positive; +3. `parallel::AvailableCpuCount()`: the CPUs the process may run on. On Linux that is the + affinity mask (`sched_getaffinity`), so `taskset`, cpusets and a launcher's binding are + respected; elsewhere it is `std::thread::hardware_concurrency()`. + +The default is the one GCC's OpenMP runtime uses, so a user gets threaded runs without +configuring anything. + +## Consequences + +- **Do not read `OMP_NUM_THREADS` again**, not even as a fallback. A fallback brings back the + coupling with numpy for everyone who sets it for numpy's sake. +- **Do not give MPI ranks a different default.** One rule for both entry points is the point; + a rank that should be serial is told so by the variable or the config. +- Under MPI, several ranks on one machine can each see most of its cores, and together start + more threads than there are cores. The CLI help and the "Solving at scale" tutorial tell + users to set `BERTINI_NUM_THREADS` when they run more than one rank per machine. That is the + same burden OpenMP places on hybrid MPI programs. +- A job script that set `OMP_NUM_THREADS` for b2 must change. The 4.0.0 CHANGELOG says so. +- `AvailableCpuCount` counts CPUs, not CPU quota. A container limited by a CFS quota + (`docker --cpus`) without a cpuset still reports every CPU it can be scheduled on. +- Tools that need a serial, reproducible run (`tools/refresh_doc_artifacts.py`) set + `BERTINI_NUM_THREADS=1` for b2 and `OMP_NUM_THREADS=1` for numpy, each for its own library. +- Pinned by `effective_thread_count_is_sane`, `bertini_num_threads_overrides_the_configured_count` + and `omp_num_threads_does_not_reach_b2` in `core/test/nag_algorithms/threaded_solve.cpp`. diff --git a/docs/adr/0071-a-path-does-not-depend-on-the-paths-before-it.md b/docs/adr/0071-a-path-does-not-depend-on-the-paths-before-it.md new file mode 100644 index 000000000..24428f6fc --- /dev/null +++ b/docs/adr/0071-a-path-does-not-depend-on-the-paths-before-it.md @@ -0,0 +1,92 @@ +# ADR-0071: A path does not depend on the paths tracked before it + +**Status:** Accepted +**Date:** 2026-10-05 + +## Context + +A tracker and an endgame are reused path after path: by the zero-dim solver, on its own +objects in a serial solve and on one clone per thread in a threaded one, and by anyone driving +them by hand. Each path's random draws were already deterministic, because the solver reseeds +the thread's generator from the global seed and the path index before each stage +(`ReseedThisThread`, ADR-0044). Results were still irreproducible under threads (#378), because +state that is not random flowed from one path into the next: + +- The adaptive tracker built a track's first step at `current_precision_`, which still held the + precision the previous track ended in, and set the step's precision before the assignment + that replaced it. +- The downward conversion to double re-precisioned only the end time. The step size, the time + and its increment kept the multiple precision the track came down from. +- The predictor and corrector rounded the condition-number probe in place at every change of + precision, so a drop destroyed digits that a later rise could not restore. +- The endgame's c/k probe was drawn the first time an endgame needed it and kept for the + object's life. An escalation then moved it out of the double lane, so the next run in double + drew a new one mid-path. +- An adaptive endgame cleared only the lane a run used, so a run in double left an earlier run's + multiprecision times in place, and `LatestTime` read them. + +With one thread, a path's predecessors are fixed by the path order, so a serial solve was +reproducible and still depended on that order. With a pool, the scheduler picks each thread's +predecessors. + +## Decision + +A path's result is a function of the systems, the settings, the seed and the path's index, and +of nothing a previous path left behind. The rule holds whatever runs the paths: one thread or +many, one process or MPI ranks. It is also met by state hygiene in the tracker and endgame, +not by building fresh objects for each path. A fresh object would hide the problem from the +solver and leave it in place for anyone reusing a tracker or an endgame by hand. + +- Whatever a track or run sets up, it sets from its own inputs: the start precision, the start + point, the stored configuration. Never from the value a previous track left in a member. +- A random direction is stored as drawn. Every working copy is derived from that copy at the + precision it is needed in, never rounded in place. The endgame's c/k probe is drawn in double, + so its value is exact in every lane and at every precision. +- A random direction that should be fresh per path is drawn at the start of the run, after the + per-path reseed, never lazily mid-path. +- Containers kept per lane are cleared in every lane when a run begins. +- **A record has a stated span, and it starts afresh exactly where that span starts.** The + tracker's precision recorders, `FirstPrecisionRecorder` and `MinMaxPrecisionRecorder`, record + one track. They start afresh at the track's `Initializing` event, before the start point is + refined, and they never unsubscribe, so a recorder attached once reports every track, failed + ones included. Before, they started afresh at `TrackingStarted`, which a track that fails to + initialize never sends and which comes after the refinement; and `FirstPrecisionRecorder` + unsubscribed at its first increase. The solver's per-path metadata (`precision_changed`, + `time_of_first_prec_increase`, `max_precision_used`) comes from `PathPrecisionRecorder`, which + spans every track of a path and starts afresh only when the solver resets it at the path's + start. Before, the solver read a per-track recorder, so a precision raised in an earlier + endgame sub-track and lowered before the last went unreported. + +## Consequences + +- **Do not store a value that the next path will read without resetting it.** When adding a + member to a tracker or an endgame, decide where each track or run gives it its value. +- **Do not round a stored random direction in place.** Derive the working copy from the stored + copy. +- **Do not draw random numbers lazily on first use.** The draw then lands in whichever path + happens to be first. +- **Do not replace this with a fresh tracker per path.** That would pass the solver's tests and + leave the trap for hand users. +- **Do not let an observer unsubscribe to save work, or start afresh on `TrackingStarted`.** The + first leaves a recorder reporting one track forever after; the second misses the start-point + refinement and leaves a failed track's record stale. +- **Do not read a per-track recorder for a per-path quantity.** The endgame tracks a path in many + calls. +- Pinned, for the recorders, by `a_precision_recorder_attached_once_reports_every_track` and + `a_track_that_fails_to_initialize_is_recorded_as_itself` (`amp_tracker_test.cpp`), and by + `max_precision_used_covers_every_track_of_a_path_*` (`threaded_solve.cpp`); each fails against + the old behavior. +- Results changed from 4.0.0.dev1 in the last bits and, occasionally, in step counts. The old + results depended on path order, so no earlier result is the reference. +- Pinned by tests at every level, each of which fails against the code without its fix: + - `a_track_does_not_depend_on_the_track_before_it` and + `the_step_size_is_at_the_working_precision_after_a_track` (`amp_tracker_test.cpp`); + - `the_probe_direction_survives_a_drop_in_precision` (`newton_correct_test.cpp`); + - `cauchy_run_does_not_depend_on_the_run_before_it` and + `power_series_run_does_not_depend_on_the_run_before_it`, which run for every tracker type; + - `threaded_solve_is_bit_identical_to_serial_*` (`threaded_solve.cpp`), which compares + endpoints and every metadata field across 1, 2, 3 and 8 threads. + Python counterparts are in `python/test/zero_dim/reproducibility_test.py`, and the + distributed-equals-serial test is in `python/test/parallel/test_mpi_zerodim.py`. +- The guarantee is bit-identity for one build on one machine. Different compilers, libraries + or instruction sets may round differently. diff --git a/docs/adr/0072-seed-zero-is-a-seed.md b/docs/adr/0072-seed-zero-is-a-seed.md new file mode 100644 index 000000000..7516da9e2 --- /dev/null +++ b/docs/adr/0072-seed-zero-is-a-seed.md @@ -0,0 +1,42 @@ +# ADR-0072: Seed 0 is a seed; the classic input format keeps Bertini 1's meaning + +**Status:** Accepted +**Date:** 2026-10-05 + +## Context + +`SetGlobalSeed(0)`, and so `bertini.set_random_seed(0)`, meant "draw a seed from entropy". A +user who set seed 0 to make a run reproducible got a different run every time (#475). 0 is the +first value people try, and in numpy (`default_rng(0)`, `np.random.seed(0)`) and in Python's +`random` (`random.seed(0)`) it is an ordinary seed. Those libraries ask for entropy with `None` +or no argument. + +Bertini 1's input file means something else. There `randomseed: 0` is the default, and it means +"choose a seed". b2's command-line program reads Bertini 1 input files, and a file that says 0, +or says nothing, expects a fresh seed. + +## Decision + +- **In the library API, every value is a seed, 0 included.** `SetGlobalSeed(0)` sets seed 0. + In Python, `set_random_seed(0)` sets seed 0. +- **Entropy is its own request.** `SetGlobalSeedFromEntropy()` in C++, and `set_random_seed()` + or `set_random_seed(None)` in Python, draw a seed and return it, so the run can be reproduced + from that value. `GetGlobalSeed()` draws one on first use if none was set, and keeps it. +- **The classic input format keeps Bertini 1's meaning.** `ApplyClassicRandomSeed` reads + `randomseed: 0` as "draw from entropy" and any other value as the seed. The command-line + program applies a classic input's seed through it. + +## Consequences + +- Seed 0 does not carry between the two interfaces. `set_random_seed(0)` is reproducible, while + `randomseed: 0` in a classic input file is not. This is the one value on which they disagree, + and the docstrings of both say so. +- **Do not make 0 mean entropy in the library API again**, and do not make the classic format's 0 + an ordinary seed: each interface follows the convention its users bring to it. +- "Not set yet" is a flag of its own in `random.cpp`, not the value 0. +- The derived seeds (`DeriveSolveSeed`, `DerivedWorkerSeed`) still never return 0. That way a + seed recorded from a library solve replays through a classic input file too. +- Pinned by `seed_zero_is_an_ordinary_seed`, + `a_seed_from_entropy_is_reported_and_reproduces_the_run` and + `classic_randomseed_zero_draws_from_entropy` (`core/test/classes/seeded_randomness_test.cpp`), + and in Python by `python/test/random/seeding_test.py`, which includes the reproducer from #475. diff --git a/docs/adr/README.md b/docs/adr/README.md index 8d5cd598b..350c6ef2f 100644 --- a/docs/adr/README.md +++ b/docs/adr/README.md @@ -86,3 +86,6 @@ Each ADR follows the template: | [0067](0067-the-endgames-differ-in-root-estimation-only.md) | The endgames differ in root estimation ONLY: the security ceiling (two consecutive endpoint approximations above `max_norm`, at security level 0) is asked BEFORE acceptance in the Cauchy endgame as it is in the power series one, and an approximation above the ceiling is never accepted -- Cauchy had returned `Success` for finite roots above the ceiling that power series truncated; Cauchy's count still runs over operating-zone rounds only, since an out-of-zone mean is not an endpoint approximation (that condition keeps cyclic-6 at 156); don't move acceptance back above security, don't drop the `!above_ceiling` clause | Core / endgames | | [0068](0068-combining-systems-compares-the-variables.md) | Combining systems compares THE VARIABLES, and a loaded system holds the live ones: `System` and `Slice` re-intern as the last step of loading from an archive (deep copies, pickles and the MPI broadcast came back over fresh variables of the same names and could not be recombined with their original); `Concatenate`, `MakeHomotopy`, `MakeMovingHomotopy` and `+=` make ONE check, variable by variable, group by group, in order -- the builders had counted, and blended `(x, y)` with `(y, x)` by position; a system that declares no variables takes the other's; in Python a node compares by identity, not by wrapper; don't go back to counting, don't remove the load hooks, don't build a node outside `Make` | Core / system, function tree; Python bindings | | [0069](0069-an-accuracy-estimate-names-its-coordinates.md) | An accuracy estimate NAMES ITS COORDINATES: `accuracy_estimate_internal_coords` (what `final_tolerance` is compared with; reads as digits) and `accuracy_estimate_user_coords` (the absolute error in the user's units); the bare `accuracy_estimate` held the internal one and read as the user's, so it is RETIRED -- reading it raises a `RuntimeError` naming both, not an `AttributeError`, which `getattr` with a default would swallow; retired and not handed to the other estimate, because the same name returning a different number breaks every reader silently; the records reader accepts the earlier key; don't re-add the bare name as an alias | Core / nag_algorithms, records; Python bindings | +| [0070](0070-b2-reads-its-own-thread-count-variable.md) | b2 reads ITS OWN thread-count variable, `BERTINI_NUM_THREADS`, and ignores `OMP_NUM_THREADS`, which numpy's OpenBLAS also reads (pinning one pinned the other); one rule for a standalone solve and each MPI rank: the variable, then `num_threads`, then every CPU the process may run on (the affinity mask on Linux), OpenMP's default; MPI ranks no longer default to one thread; don't read `OMP_NUM_THREADS` again, even as a fallback, and don't give ranks a different default | Core / parallel, nag_algorithms; CLI; docs | +| [0071](0071-a-path-does-not-depend-on-the-paths-before-it.md) | A path DOES NOT DEPEND ON THE PATHS BEFORE IT: its result is a function of the systems, settings, seed and index, met by state hygiene in the tracker and endgame rather than a fresh object per path (which would leave the trap for hand users); a track sets what it uses from its own inputs, never a leftover member (the first step was built at the previous track's precision); random directions are stored as drawn and every working copy derived, never rounded in place; per-path directions are drawn at run start after the reseed, never lazily; per-lane containers cleared in every lane; threaded, MPI and serial solves are bit-identical; a record starts afresh exactly where its span starts: the per-track precision recorders at `Initializing`, never unsubscribing, and the solver's per-path metadata from `PathPrecisionRecorder`, reset only at path start | Core / trackers, endgames, nag_algorithms | +| [0072](0072-seed-zero-is-a-seed.md) | SEED 0 IS A SEED in the library API (`SetGlobalSeed(0)`, `set_random_seed(0)`), as in numpy and Python's random; entropy is its own request (`SetGlobalSeedFromEntropy()`, `set_random_seed()` / `None`) and returns the seed it chose; the classic input format keeps Bertini 1's meaning, `randomseed: 0` = draw from entropy, via `ApplyClassicRandomSeed`; so seed 0 alone does not carry between the two interfaces; "not set" is its own flag; derived seeds still avoid 0 so they replay through classic input; don't make the API's 0 mean entropy again | Core / random, blackbox; Python bindings | diff --git a/python/bertini/nag_algorithm/__init__.py b/python/bertini/nag_algorithm/__init__.py index 84a419c3a..e77a73b9a 100644 --- a/python/bertini/nag_algorithm/__init__.py +++ b/python/bertini/nag_algorithm/__init__.py @@ -1231,7 +1231,7 @@ def parameter_sweep(make_system, generic_parameters, target_parameters, comm : mpi4py communicator, optional If given, ``target_parameters`` is split across the ranks (each rank solves the generic once, then tracks its slice -- with the path tracking inside each solve threaded across the - rank's cores via ``OMP_NUM_THREADS``). The collected results are gathered and **every rank + rank's cores, set with ``BERTINI_NUM_THREADS``). The collected results are gathered and **every rank returns the full list in the original order**. This is the two-level model: MPI across parameter points, threads across paths -- one extra argument, no manual scatter/gather. mptype, endgame diff --git a/python/docs/source/tutorials/parallelism/parallel_parameter_homotopy/index.rst b/python/docs/source/tutorials/parallelism/parallel_parameter_homotopy/index.rst index 02b029c3e..194735058 100644 --- a/python/docs/source/tutorials/parallelism/parallel_parameter_homotopy/index.rst +++ b/python/docs/source/tutorials/parallelism/parallel_parameter_homotopy/index.rst @@ -135,7 +135,7 @@ Run it serially, or across ranks with threads inside each (the ``--bind-to none` ``--map-by`` notes from :doc:`/tutorials/parallelism/solving_at_scale/index` apply -- and always run a *file*, never a heredoc):: python parallel_parameter_homotopy.py --grid 28 --save schlogl.svg # serial + threads - OMP_NUM_THREADS=3 mpirun -n 4 --bind-to none \ + BERTINI_NUM_THREADS=3 mpirun -n 4 --bind-to none \ python parallel_parameter_homotopy.py --grid 28 --save schlogl.svg # 4 ranks x 3 threads Both produce the **identical** map -- parallelism changes the wall-clock, never the answer: diff --git a/python/docs/source/tutorials/parallelism/solving_at_scale/index.rst b/python/docs/source/tutorials/parallelism/solving_at_scale/index.rst index cddf51e59..6d83ff287 100644 --- a/python/docs/source/tutorials/parallelism/solving_at_scale/index.rst +++ b/python/docs/source/tutorials/parallelism/solving_at_scale/index.rst @@ -13,7 +13,7 @@ the parallel payoff. This tutorial shows the two knobs: * **MPI ranks** -- a *manager-worker* pool, possibly spread across many machines. Pass a communicator to :func:`solve` and rank ``0`` hands paths out to the rest. -* **threads per rank** -- each worker can track several paths at once. Set ``OMP_NUM_THREADS``. +* **threads per rank** -- each worker can track several paths at once. Set ``BERTINI_NUM_THREADS``. Ideally, the two knobs scale *multiplicatively*: ``R`` worker ranks of ``T`` threads each track ``R × T`` paths at a time. We demonstrate on two driving problems -- a **cyclic-n** system (hundreds of cheap @@ -45,8 +45,8 @@ communicator. the manager says "done". Because the manager does not track, a distributed run needs **at least two ranks** (one manager + -one worker). The serial baseline is just :func:`solve` with no communicator. The minimal program -is:: +one worker). A single process is just :func:`solve` with no communicator; it still tracks paths on +every CPU it may run on, and ``BERTINI_NUM_THREADS=1`` makes it serial. The minimal program is:: from mpi4py import MPI import bertini as pb @@ -57,7 +57,7 @@ is:: if comm.Get_size() > 1: solver.solve(communicator=comm) # rank 0 dispatches; the others track else: - solver.solve() # serial: a lone manager would have no workers + solver.solve() # one process: a lone manager would have no workers if pb.parallel.is_manager(): # only rank 0 has the collected solutions print(len(solver.all_solutions()), "solutions") @@ -68,11 +68,13 @@ The :mod:`bertini.parallel` helpers -- :func:`~bertini.parallel.rank`, holding the solution list, so guard your reporting and your correctness checks with ``is_manager()``. -Launch it with ``mpirun``:: +Launch it with ``mpirun``. The examples in this section set ``BERTINI_NUM_THREADS=1`` so that +each worker tracks one path at a time and the rank count alone sets the parallelism; threads are +the subject of a later section:: - python my_solve.py # serial baseline (1 process) - mpirun -n 2 python my_solve.py # 1 manager + 1 worker (the smallest parallel run) - mpirun -n 9 python my_solve.py # 1 manager + 8 workers + BERTINI_NUM_THREADS=1 python my_solve.py # serial baseline (1 process) + BERTINI_NUM_THREADS=1 mpirun -n 2 python my_solve.py # 1 manager + 1 worker + BERTINI_NUM_THREADS=1 mpirun -n 9 python my_solve.py # 1 manager + 8 workers .. warning:: @@ -86,7 +88,7 @@ Launch it with ``mpirun``:: ("not enough slots"); add ``--map-by :OVERSUBSCRIBE`` to let the manager share a core (it barely uses one). For example, to run 8 workers on an 8-core box:: - mpirun -n 9 --map-by :OVERSUBSCRIBE python my_solve.py + BERTINI_NUM_THREADS=1 mpirun -n 9 --map-by :OVERSUBSCRIBE python my_solve.py Scaling a cyclic system ======================= @@ -129,6 +131,7 @@ Run the ladder (serial, then 2, 4, 8 workers): .. code-block:: console + $ export BERTINI_NUM_THREADS=1 # one path per worker at a time $ python python/examples/solve_cyclic.py --n 6 # serial baseline $ mpirun -n 3 python python/examples/solve_cyclic.py --n 6 # 2 workers $ mpirun -n 5 python python/examples/solve_cyclic.py --n 6 # 4 workers @@ -176,6 +179,7 @@ Run the same ladder on a sizable matrix: .. code-block:: console + $ export BERTINI_NUM_THREADS=1 # one path per worker at a time $ python python/examples/solve_eigenvalues.py --size 24 # serial baseline $ mpirun -n 3 python python/examples/solve_eigenvalues.py --size 24 # 2 workers $ mpirun -n 5 python python/examples/solve_eigenvalues.py --size 24 # 4 workers @@ -203,14 +207,16 @@ Using both threads and processes ====================================== Each worker rank can itself track several paths concurrently, one per thread, set with the -``OMP_NUM_THREADS`` environment variable (default 1). This stacks with MPI: ``R`` worker ranks of -``T`` threads each give ``R × T`` paths in flight. +``BERTINI_NUM_THREADS`` environment variable. Unset, a rank uses one thread per CPU it may run on, +the same default OpenMP uses. This stacks with MPI: ``R`` worker ranks of ``T`` threads each give +``R × T`` paths in flight. -The launch needs one extra flag, ``--bind-to none``. By default ``mpirun`` pins each rank to a -single core, which would strangle that rank's threads onto one core; ``--bind-to none`` lets a -rank's threads spread across the machine:: +Set ``BERTINI_NUM_THREADS`` whenever you run more than one rank on a machine. How many CPUs a rank +may run on depends on how ``mpirun`` binds it, and ranks that each see most of the machine would +together start several threads per core. Launch with ``--bind-to none`` as well: it lets a +rank's threads spread across the machine instead of sharing the cores the rank was bound to:: - OMP_NUM_THREADS=4 mpirun -n 3 --bind-to none python python/examples/solve_cyclic.py --n 6 + BERTINI_NUM_THREADS=4 mpirun -n 3 --bind-to none python python/examples/solve_cyclic.py --n 6 # 4 threads x (3 ranks = 1 manager + 2 workers) = 8 paths at a time Why have two knobs for the same cores? MPI ranks can live on diff --git a/python/docs/source/tutorials/performance_benchmarking/bertini1_vs_bertini2_timing/b1_vs_b2_timing.py b/python/docs/source/tutorials/performance_benchmarking/bertini1_vs_bertini2_timing/b1_vs_b2_timing.py index cc6bf053d..eb5fc14de 100644 --- a/python/docs/source/tutorials/performance_benchmarking/bertini1_vs_bertini2_timing/b1_vs_b2_timing.py +++ b/python/docs/source/tutorials/performance_benchmarking/bertini1_vs_bertini2_timing/b1_vs_b2_timing.py @@ -109,7 +109,7 @@ def time_solver(exe, input_text, repeats=3, timeout=600): d = tempfile.mkdtemp() try: open(os.path.join(d, "input"), "w").write(input_text) - env = dict(os.environ, OMP_NUM_THREADS="1") + env = dict(os.environ, BERTINI_NUM_THREADS="1") t0 = time.perf_counter() r = subprocess.run([os.path.abspath(exe), "input"], cwd=d, env=env, stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL, timeout=timeout) diff --git a/python/docs/source/tutorials/performance_benchmarking/bertini1_vs_bertini2_timing/index.rst b/python/docs/source/tutorials/performance_benchmarking/bertini1_vs_bertini2_timing/index.rst index 1200a2ad3..f89c18e98 100644 --- a/python/docs/source/tutorials/performance_benchmarking/bertini1_vs_bertini2_timing/index.rst +++ b/python/docs/source/tutorials/performance_benchmarking/bertini1_vs_bertini2_timing/index.rst @@ -14,7 +14,7 @@ build each system once in pybertini and emit it as a Bertini-1 classic input wit :py:meth:`System.to_classic_input`, which writes the equations *and* the tracking knobs (precision mode, predictor, tolerances, step cadence). Both solvers read that one file, each builds its own random start system, and solves. We time only the solver subprocess (never parsing/IO) and keep the -fastest of a few repeats, serial (``OMP_NUM_THREADS=1``): +fastest of a few repeats, serial (``BERTINI_NUM_THREADS=1``): .. literalinclude:: b1_vs_b2_timing.py :language: python diff --git a/python/examples/parallel_parameter_homotopy.py b/python/examples/parallel_parameter_homotopy.py index f0f16149c..c137f56c2 100644 --- a/python/examples/parallel_parameter_homotopy.py +++ b/python/examples/parallel_parameter_homotopy.py @@ -7,7 +7,7 @@ *outer* loop; the points are independent, so no manager-worker dispatch is needed -- just split the grid and gather the answers). * **threads across paths** -- inside each point's solve, the path tracking runs on all the cores - the rank is given (``OMP_NUM_THREADS``). + the rank may run on, or on ``BERTINI_NUM_THREADS`` of them. The model is the **Schlogl** reaction network, whose steady state is a cubic in the species concentration ``X``:: @@ -22,9 +22,10 @@ Run it (per the MPI notes in the "Solving at scale" tutorial -- always run a *file*, never a heredoc; ``--bind-to none`` lets each rank's threads spread across cores):: - python parallel_parameter_homotopy.py # serial baseline - mpirun -n 4 python parallel_parameter_homotopy.py # 4 ranks, 1 thread each - OMP_NUM_THREADS=3 mpirun -n 4 --bind-to none \ + BERTINI_NUM_THREADS=1 python parallel_parameter_homotopy.py # serial baseline + BERTINI_NUM_THREADS=1 mpirun -n 4 \ + python parallel_parameter_homotopy.py # 4 ranks, 1 thread each + BERTINI_NUM_THREADS=3 mpirun -n 4 --bind-to none \ python parallel_parameter_homotopy.py # 4 ranks x 3 threads The result (a counts grid, and optionally a heatmap PNG) is assembled on rank 0. diff --git a/python/examples/solve_cyclic.py b/python/examples/solve_cyclic.py index 950cc06c0..dceb668e7 100644 --- a/python/examples/solve_cyclic.py +++ b/python/examples/solve_cyclic.py @@ -4,14 +4,14 @@ The cyclic-n polynomials are a classic benchmark for polynomial system solving. Total-degree homotopy tracks one path per Bezout count (the product of the degrees), and the work is pleasantly parallel: each path is independent. Bertini hands the paths out from rank 0 (the -manager) to the worker ranks, and -- if OMP_NUM_THREADS > 1 -- each worker runs several tracking -threads. +manager) to the worker ranks, and each worker runs several tracking threads: by default one per +CPU it may run on, or ``BERTINI_NUM_THREADS`` of them. Run it:: - python solve_cyclic.py --n 6 # serial baseline + BERTINI_NUM_THREADS=1 python solve_cyclic.py --n 6 # serial baseline mpirun -n 5 python solve_cyclic.py --n 6 # 1 manager + 4 workers - OMP_NUM_THREADS=4 mpirun -n 3 --bind-to none python solve_cyclic.py --n 6 # 2 workers x 4 threads + BERTINI_NUM_THREADS=4 mpirun -n 3 --bind-to none python solve_cyclic.py --n 6 # 2 workers x 4 threads Only rank 0 prints; it checks the distinct finite-solution count -- taken from the solver's own solution metadata -- against the known cyclic-n value. @@ -80,7 +80,7 @@ def main(): report = solver.report() print('cyclic-{}: ranks={} threads/rank={} paths tracked={} finite solutions={} wall={:.1f}s'.format( - args.n, comm.Get_size(), os.environ.get('OMP_NUM_THREADS', '1'), + args.n, comm.Get_size(), os.environ.get('BERTINI_NUM_THREADS', 'auto'), report.num_paths_tracked, report.num_finite_solutions, elapsed)) if not report.all_paths_resolved: # a failed path OR an unresolved crossing -- show what and why diff --git a/python/examples/solve_eigenvalues.py b/python/examples/solve_eigenvalues.py index d7755e3c9..ed40abdd2 100644 --- a/python/examples/solve_eigenvalues.py +++ b/python/examples/solve_eigenvalues.py @@ -10,9 +10,9 @@ Run it:: - python solve_eigenvalues.py --size 24 # serial baseline + BERTINI_NUM_THREADS=1 python solve_eigenvalues.py --size 24 # serial baseline mpirun -n 5 python solve_eigenvalues.py --size 24 # 1 manager + 4 workers - OMP_NUM_THREADS=4 mpirun -n 3 --bind-to none python solve_eigenvalues.py --size 24 + BERTINI_NUM_THREADS=4 mpirun -n 3 --bind-to none python solve_eigenvalues.py --size 24 Only rank 0 prints; it checks the recovered eigenvalues against ``numpy.linalg.eigvals`` -- correctness, not just a path count. @@ -108,7 +108,7 @@ def main(): distinct.append(g) print('eigen-{}: ranks={} threads/rank={} paths={} recovered={} wall={:.1f}s'.format( - args.size, comm.Get_size(), os.environ.get('OMP_NUM_THREADS', '1'), + args.size, comm.Get_size(), os.environ.get('BERTINI_NUM_THREADS', 'auto'), args.size, len(distinct), elapsed)) # Correctness: every numpy eigenvalue is matched by a homotopy-recovered one (closest-match, diff --git a/python/test/parallel/test_mpi_zerodim.py b/python/test/parallel/test_mpi_zerodim.py index a815534c8..8c6f029ef 100644 --- a/python/test/parallel/test_mpi_zerodim.py +++ b/python/test/parallel/test_mpi_zerodim.py @@ -174,6 +174,34 @@ def _eigen_system(A): return sys +def _cyclic5_path_by_path(comm=None): + pb.random.set_random_seed(20261005) + solver = ZeroDimSolver(_cyclic_system(5)) + if comm is not None: + solver.solve(communicator=comm) + else: + solver.solve() + if not pb.parallel.is_manager(): + return None + return [([repr(c) for c in point], + (int(md.endgame_success_code), int(md.num_successful_steps), int(md.num_failed_steps), + int(md.max_precision_used), repr(md.final_time_used), repr(md.condition_number))) + for point, md in zip(solver.all_solutions(), solver.solution_metadata())] + + +def test_distributed_solve_is_bit_identical_to_serial(): + """A path's result does not depend on which rank ran it, nor which thread, nor what ran + before it there (#378): the distributed solve matches a serial one bit for bit, path by path. + Under plain pytest (one rank) it compares two serial solves.""" + comm = MPI.COMM_WORLD + serial = _cyclic5_path_by_path() + distributed = _cyclic5_path_by_path(comm if comm.Get_size() > 1 else None) + if pb.parallel.is_manager(): + assert len(serial) == len(distributed) == 120 + for index, (a, b) in enumerate(zip(serial, distributed)): + assert a == b, 'path {} differs between the serial and the distributed solve'.format(index) + + def test_distributed_mhom_eigenvalues_match_known_count(): pb.random.set_random_seed(3) # small symmetric integer matrix: real, generically distinct spectrum, n paths under mhom. diff --git a/python/test/random/seeding_test.py b/python/test/random/seeding_test.py index 8f1248043..98dfffd49 100644 --- a/python/test/random/seeding_test.py +++ b/python/test/random/seeding_test.py @@ -64,6 +64,36 @@ def test_different_seed_changes_multiprecision_draws(): assert first != other +def test_seed_zero_is_an_ordinary_seed(): + """0 is a seed like any other, as in numpy and Python's random (#475, ADR-0072).""" + assert pb.random.set_random_seed(0) == 0 + assert pb.random.get_random_seed() == 0 + first = _draw_sequence() + + pb.random.set_random_seed(0) + assert first == _draw_sequence() + + pb.random.set_random_seed(1) + assert first != _draw_sequence() + + +@pytest.mark.parametrize("call", [lambda: pb.random.set_random_seed(), + lambda: pb.random.set_random_seed(None)], + ids=["no argument", "None"]) +def test_entropy_is_asked_for_by_omitting_the_seed_and_the_seed_is_returned(call): + chosen = call() + assert pb.random.get_random_seed() == chosen + first = _draw_sequence() + + pb.random.set_random_seed(chosen) # the returned seed reproduces the run + assert first == _draw_sequence() + + +def test_two_entropy_seeds_differ(): + # equal with probability about 2^-32 + assert pb.random.set_random_seed() != pb.random.set_random_seed() + + def test_unit_draws_also_honor_the_seed(): pb.random.set_random_seed(99) first = [repr(pb.random.complex_unit()) for _ in range(4)] @@ -120,6 +150,31 @@ def test_same_seed_reproduces_patched_solution_coordinates(): assert np.linalg.norm(va - vb) < 1e-10 +@pytest.mark.parametrize("num_threads", [1, 4]) +def test_seed_zero_reproduces_a_solve_path_by_path(num_threads, monkeypatch): + """The reproducer from #475: re-solving after set_random_seed(0) changed which path diverged + and how many steps each took, from trial to trial, single-threaded too.""" + monkeypatch.delenv('BERTINI_NUM_THREADS', raising=False) + x, y = pb.Variable('x'), pb.Variable('y') + sys = pb.System() + sys.add_variable_group(pb.VariableGroup([x, y])) + sys.add_function(x * y - 1) + sys.add_function(x**2 - y) + + trials = [] + for _ in range(3): + pb.random.set_random_seed(0) + solver = ZeroDimSolver(sys) + cfg = solver.get_config(pb.nag_algorithm.ZeroDimConfig) + cfg.num_threads = num_threads + solver.set_config(cfg) + solver.solve() + trials.append([([repr(c) for c in point], int(md.num_successful_steps), + int(md.num_failed_steps), bool(md.is_finite)) + for point, md in zip(solver.all_solutions(), solver.solution_metadata())]) + assert trials[0] == trials[1] == trials[2] + + def test_user_solutions_are_seed_independent(): """The patch differs by seed, but the actual solutions do not.""" for seed in (31415, 27182): diff --git a/python/test/zero_dim/reproducibility_test.py b/python/test/zero_dim/reproducibility_test.py new file mode 100644 index 000000000..5cb6dd390 --- /dev/null +++ b/python/test/zero_dim/reproducibility_test.py @@ -0,0 +1,149 @@ +"""A path's result does not depend on what ran before it, nor on how many threads ran (#378). + +Trackers and endgames are reused path after path: by the solver on each thread, and by anyone +driving one by hand. These tests ask for bit-identical results -- every coordinate compared +through its full-precision repr, and every metadata field that is not a wall-clock time -- where +the older tests asked only for the same solution set within a tolerance, which this bug passed. + +The C++ suite is the gate (core/test/nag_algorithms/threaded_solve.cpp, +core/test/tracking_basics/amp_tracker_test.cpp, the generic endgame tests); these confirm the +same holds through the Python interface. +""" + +import numpy as np +import pytest + +import bertini as pb +import bertini.multiprec as mp +from bertini import ZeroDimSolver +from bertini.tracking import AMPTracker, amp_config_from, Predictor, SteppingConfig, NewtonConfig + + +def _cyclic(n): + xs = pb.variables('x', n) + sys = pb.System() + sys.add_variable_group(pb.VariableGroup(list(xs))) + for k in range(1, n): + total = None + for i in range(n): + term = None + for j in range(k): + v = xs[(i + j) % n] + term = v if term is None else term * v + total = term if total is None else total + term + sys.add_function(total) + prod = None + for v in xs: + prod = v if prod is None else prod * v + sys.add_function(prod - 1) + return sys + + +# every metadata field a path's computation determines; path_time_seconds is a wall-clock time +_FIELDS = ('endgame_success_code', 'pre_endgame_success_code', 'num_successful_steps', + 'num_failed_steps', 'max_precision_used', 'precision_changed', + 'time_of_first_prec_increase', 'final_time_used', 'condition_number', + 'accuracy_estimate_internal_coords', 'accuracy_estimate_user_coords', 'cycle_num', + 'function_residual', 'newton_residual', 'multiplicity', 'is_finite', 'is_real', + 'is_singular') + + +def _path_by_path(endgame, num_threads): + pb.random.set_random_seed(20261005) + solver = ZeroDimSolver(_cyclic(5), endgame=endgame) + cfg = solver.get_config(pb.nag_algorithm.ZeroDimConfig) + cfg.num_threads = num_threads + solver.set_config(cfg) + solver.solve() + out = [] + for point, md in zip(solver.all_solutions(), solver.solution_metadata()): + out.append(([repr(c) for c in point], + {f: repr(getattr(md, f)) for f in _FIELDS})) + return out + + +@pytest.mark.parametrize("endgame", ["powerseries", "cauchy"]) +def test_threaded_solve_is_bit_identical_to_serial(endgame, monkeypatch): + # the variable would override num_threads and make every solve here serial + monkeypatch.delenv('BERTINI_NUM_THREADS', raising=False) + serial = _path_by_path(endgame, 1) + assert len(serial) == 120 + for num_threads in (2, 4): + threaded = _path_by_path(endgame, num_threads) + assert len(threaded) == len(serial) + for index, (a, b) in enumerate(zip(serial, threaded)): + assert a[0] == b[0], 'path {} endpoint differs at {} threads'.format(index, num_threads) + for field in _FIELDS: + assert a[1][field] == b[1][field], \ + 'path {} {} differs at {} threads: {} vs {}'.format( + index, field, num_threads, a[1][field], b[1][field]) + + +@pytest.mark.parametrize("endgame", ["powerseries", "cauchy"]) +def test_max_precision_used_covers_every_track_of_a_path(endgame): + """A path's max_precision_used is at least every precision its tracker stepped at, across the + track to the endgame boundary and every endgame sub-track. (On this solve it once reported 30 + for a path whose endgame had run at 40: the record started afresh with every sub-track.)""" + pb.random.set_random_seed(20261005) + solver = ZeroDimSolver(_cyclic(5), endgame=endgame) + collector = pb.nag_algorithm.SolutionPathCollector() + solver.add_observer(collector) + solver.solve() + + metadata = solver.solution_metadata() + assert len(collector.series) == len(metadata) == 120 + for path in collector.series: + stepped = path.diagnostics()[:, 2] # the precision column, at every accepted step + if len(stepped) == 0: + continue + reported = int(metadata[path.path_index].max_precision_used) + assert reported >= int(stepped.max()), \ + 'path {}: max_precision_used {} below the {} it stepped at'.format( + path.path_index, reported, int(stepped.max())) + + +def _double_root_tracker(): + """x^2 - t, y - x: toward a double root at t = 0, tracked at tolerance 1e-12 so that a track + to t = 1e-20 ends in multiple precision. The seed is set before building, so every tracker + made here draws the same condition-number probe.""" + x, y, t = pb.Variable('x'), pb.Variable('y'), pb.Variable('t') + s = pb.System() + s.add_function(x**2 - t) + s.add_function(y - x) + s.add_variable_group(pb.VariableGroup([x, y])) + s.add_path_variable(t) + mp.default_precision(30) + pb.random.set_random_seed(20261005) + tracker = AMPTracker(s) + tracker.setup(Predictor.RKF45, 1e-12, 1e5, SteppingConfig(), NewtonConfig()) + tracker.precision_setup(amp_config_from(s)) + return s, tracker + + +def _track(tracker, digits, t_end): + mp.default_precision(digits) + start = np.array([mp.complex_mp(1), mp.complex_mp(1)]) + end = np.array([mp.complex_mp(0)] * 2) + code = tracker.track_path(end, mp.complex_mp(1), mp.complex_mp(t_end), start) + return (code, tracker.num_total_steps_taken(), tracker.current_precision(), + repr(tracker.current_stepsize()), repr(tracker.latest_condition_number()), + [repr(c) for c in end]) + + +@pytest.mark.parametrize("before_digits, digits", [(60, 16), (16, 30), (20, 30), (30, 60)]) +def test_a_track_does_not_depend_on_the_track_before_it(before_digits, digits): + """A hand-driven tracker, reused: the same track after another must match a fresh tracker's.""" + _, used = _double_root_tracker() + used.precision_preservation(True) # the previous track ends where it started + before = _track(used, before_digits, '0.5') + used.precision_preservation(False) + assert before[0] == pb.SuccessCode.Success + assert before[2] != digits + + after_another = _track(used, digits, '1e-20') + _, fresh = _double_root_tracker() + first_ever = _track(fresh, digits, '1e-20') + + assert first_ever[0] == pb.SuccessCode.Success + assert first_ever[2] > 16 # ends in multiple precision, where the mp probe is used + assert after_another == first_ever diff --git a/python/test/zero_dim/threaded_solve_test.py b/python/test/zero_dim/threaded_solve_test.py index 7894f96cd..243979440 100644 --- a/python/test/zero_dim/threaded_solve_test.py +++ b/python/test/zero_dim/threaded_solve_test.py @@ -124,6 +124,58 @@ def Observe(self, e): assert not any(saw_member) # never the member tracker -> a clone +def _threads_that_ran_paths(num_threads): + """Solve with `num_threads` configured; return the set of OS threads that started a path. + + A path's PathStarted event is delivered on the thread that runs the path: the calling thread in + a serial solve, a pool worker in a threaded one.""" + import threading + from bertini._pybertini import nag_algorithms as nag + + solver = ZeroDimSolver(_two_cubics(), mptype='amp') + cfg = solver.get_config(pb.nag_algorithm.ZeroDimConfig) + cfg.num_threads = num_threads + solver.set_config(cfg) + threads = set() + + class Probe(nag.observers.CustomObserver): + def Observe(self, e): + if isinstance(e, nag.observers.PathStarted): + threads.add(threading.get_ident()) + + solver.add_observer(Probe()) + solver.solve() + return threads + + +def test_bertini_num_threads_one_overrides_a_threaded_config(monkeypatch): + """BERTINI_NUM_THREADS=1 wins over num_threads=4: every path runs on the calling thread.""" + import threading + monkeypatch.setenv('BERTINI_NUM_THREADS', '1') + assert _threads_that_ran_paths(4) == {threading.get_ident()} + + +def test_bertini_num_threads_overrides_a_serial_config(monkeypatch): + """BERTINI_NUM_THREADS=3 wins over num_threads=1: the paths run on at most three pool threads, + none of them the calling thread.""" + import threading + monkeypatch.setenv('BERTINI_NUM_THREADS', '3') + threads = _threads_that_ran_paths(1) + assert threading.get_ident() not in threads + assert 1 <= len(threads) <= 3 + + +def test_omp_num_threads_does_not_reach_bertini(monkeypatch): + """OMP_NUM_THREADS is numpy's (OpenBLAS's) setting, not b2's (ADR-0070): with it at 1 and + BERTINI_NUM_THREADS unset, a solve configured for 4 threads still runs on pool threads.""" + import threading + monkeypatch.delenv('BERTINI_NUM_THREADS', raising=False) + monkeypatch.setenv('OMP_NUM_THREADS', '1') + threads = _threads_that_ran_paths(4) + assert threading.get_ident() not in threads + assert 1 <= len(threads) <= 4 + + def test_solution_path_collector_under_threads(): """The two-level meta-observer collects exactly one series per path under threading, because it now attaches to event.tracker() (the clone that actually runs the path).""" diff --git a/python_bindings/src/random_export.cpp b/python_bindings/src/random_export.cpp index 6735bf1f0..d1060c729 100644 --- a/python_bindings/src/random_export.cpp +++ b/python_bindings/src/random_export.cpp @@ -75,13 +75,23 @@ void ExportRandom(){ (arg("rows"), arg("cols"), arg("real")=false), "A rows x cols random conjugate-orthonormal matrix (its rows orthonormal -- QR-factored from a square matrix of units, then truncated; perfectly conditioned), at the current default precision. real=True yields a real orthogonal matrix (entries with zero imaginary part). Returns a complex_mp matrix."); - def("set_random_seed", &bertini::SetGlobalSeed, boost::python::arg("seed") = 0ul, - "Set the global RNG seed (0 = draw from entropy). Call before constructing any " - "homotopy or solver to get reproducible results. The effective seed (which may " - "differ from 0 when entropy is used) is retrievable via get_random_seed()."); + def("set_random_seed", + +[](object seed) -> unsigned long { + if (seed.is_none()) + return bertini::SetGlobalSeedFromEntropy(); + unsigned long const s = extract(seed); + bertini::SetGlobalSeed(s); + return s; + }, + boost::python::arg("seed") = object(), + "Set the global RNG seed, and return it. Every non-negative integer is a seed, 0 " + "included, as in numpy and Python's random; set_random_seed() or set_random_seed(None) " + "draws one from entropy, and the returned value reproduces the run. Call before " + "constructing any homotopy or solver to get reproducible results. (In a classic Bertini " + "input file, randomseed: 0 instead means draw from entropy, as in Bertini 1.)"); def("get_random_seed", &bertini::GetGlobalSeed, - "Return the effective global RNG seed. If set_random_seed has not been called, " - "draws from entropy on first call and caches the result."); + "Return the global RNG seed in effect. If none has been set, draws one from entropy " + "on the first call and keeps it."); def("derive_solve_seed", &bertini::DeriveSolveSeed, "Capture the session stream's current position as an effective per-solve seed: " "one draw off the current stream (advancing it), nonzero, 32-bit portable. " diff --git a/python_bindings/src/zero_dim_configs_export.cpp b/python_bindings/src/zero_dim_configs_export.cpp index eef638944..aba50ba95 100644 --- a/python_bindings/src/zero_dim_configs_export.cpp +++ b/python_bindings/src/zero_dim_configs_export.cpp @@ -125,9 +125,9 @@ namespace bertini{ "How many times to re-track crossed paths (with tightened settings) at the endgame " "boundary before giving up. 0 = detect and report only, do not re-track. Default 2.") .def_readwrite("num_threads", &ZeroDimConfig::num_threads, - "Worker threads for a shared-memory (non-MPI) solve. 0 = auto " - "(all available cores), 1 = serial (no thread pool), N = N threads. The " - "OMP_NUM_THREADS environment variable overrides this. Threading needs no MPI and " + "Worker threads tracking paths, in a standalone solve and on each MPI rank. 0 = auto " + "(every CPU this process may run on), 1 = serial (no thread pool), N = N threads. The " + "BERTINI_NUM_THREADS environment variable overrides this. Threading needs no MPI and " "no free-threaded Python: the heavy tracking runs in C++ with the GIL released.") .def_readwrite("max_path_wall_clock_duration", &ZeroDimConfig::max_path_wall_clock_duration, "Wall-clock budget for EACH path, in seconds; 0 (the default) means none. A path that " diff --git a/tools/refresh_doc_artifacts.py b/tools/refresh_doc_artifacts.py index 7372fd519..88fe961d7 100644 --- a/tools/refresh_doc_artifacts.py +++ b/tools/refresh_doc_artifacts.py @@ -149,7 +149,10 @@ def _env_with_determinism(): rc = tempfile.NamedTemporaryFile("w", suffix="matplotlibrc", delete=False) rc.write(f"svg.hashsalt: {SVG_HASHSALT}\n") rc.close() - env = dict(os.environ, MATPLOTLIBRC=rc.name, MPLBACKEND="Agg", OMP_NUM_THREADS="1") + # serial b2 (BERTINI_NUM_THREADS) and serial numpy (OMP_NUM_THREADS, which OpenBLAS reads): + # a figure must not depend on thread scheduling in either library + env = dict(os.environ, MATPLOTLIBRC=rc.name, MPLBACKEND="Agg", BERTINI_NUM_THREADS="1", + OMP_NUM_THREADS="1") return env, rc.name diff --git a/tools/run_with_watchdog.sh b/tools/run_with_watchdog.sh index 60713eb48..675d69abb 100755 --- a/tools/run_with_watchdog.sh +++ b/tools/run_with_watchdog.sh @@ -13,7 +13,7 @@ # tools/run_with_watchdog.sh [args...] # Examples: # tools/run_with_watchdog.sh 30 ./build/core/test_nag_algorithms -# OMP_NUM_THREADS=1 PYTHONPATH=python tools/run_with_watchdog.sh 60 python3 my_solve.py +# BERTINI_NUM_THREADS=1 PYTHONPATH=python tools/run_with_watchdog.sh 60 python3 my_solve.py # # Exit code: the command's exit code if it finished in time; 124 if it was killed by the watchdog # (the conventional timeout(1) code). diff --git a/tools/update_scaling_timings.py b/tools/update_scaling_timings.py index 9d1f5c224..e3617a002 100644 --- a/tools/update_scaling_timings.py +++ b/tools/update_scaling_timings.py @@ -74,7 +74,7 @@ def launch(script, args, nprocs, omp, bind_none, timeout): Aborts the whole tool if the run fails or does not print its ``OK:`` correctness line -- we never want to publish a number from a solve that did not verify. """ - env = dict(os.environ, OMP_NUM_THREADS=str(omp)) + env = dict(os.environ, BERTINI_NUM_THREADS=str(omp)) if nprocs <= 1: cmd = [sys.executable, script, *args] else: @@ -87,7 +87,7 @@ def launch(script, args, nprocs, omp, bind_none, timeout): cmd += ["--bind-to", "none"] cmd += [sys.executable, script, *args] - shown = (f"OMP_NUM_THREADS={omp} " if omp != 1 else "") + " ".join(cmd) + shown = f"BERTINI_NUM_THREADS={omp} " + " ".join(cmd) print(f" $ {shown}", flush=True) proc = subprocess.run(cmd, env=env, cwd=REPO, capture_output=True, text=True, timeout=timeout) out = proc.stdout + proc.stderr