From 2cb3a49ab40c777add34b87b7526ce5df0524264 Mon Sep 17 00:00:00 2001
From: sou-cheng-choi <sou.cheng.choi@stanfordalumni.org>
Date: Mon, 31 Aug 2026 10:58:19 +0800
Subject: [PATCH 1/3] tandardize test filenames and improve formatter
 diagnostics

---
 CONTRIBUTING.md                               |    9 +-
 docs/api/discrete_distributions.md            |    5 +-
 docs/mpmc-compatibility.md                    |    7 +-
 docs/tests.md                                 |   42 +
 makefile                                      |   15 +
 qmcpy/accumulate_data/__init__.py             |    1 -
 scripts/check_test_style.py                   |  131 ++
 scripts/flatten_qmcpy_imports.py              |   27 +-
 scripts/remove_trailing_whitespace.py         |   14 +-
 scripts/unwrap_markdown.py                    |   22 +-
 test/README.md                                |   42 +
 test/test_check_links.py                      |  178 ---
 test/test_check_removed_urls.py               |  134 --
 test/test_copulas.py                          | 1362 -----------------
 ...stribs.py => test_dd_discrete_distribs.py} |    0
 test/test_dd_dummy_sampler.py                 |  107 ++
 test/test_dummy_sampler.py                    |  111 --
 ...test_integrate.py => test_ee_integrate.py} |    0
 test/{test_keister.py => test_ee_keister.py}  |    0
 ...st_pi_problem.py => test_ee_pi_problem.py} |    0
 test/test_fast_transform_fallbacks.py         |   46 -
 test/test_financial_option_quick.py           |   94 --
 test/test_flatten_qmcpy_imports.py            |  330 ----
 test/test_ft_fast_transform_fallbacks.py      |   52 +
 test/test_ig_financial_option_quick.py        |  102 ++
 ...st_integrands.py => test_ig_integrands.py} |    0
 test/{test_option.py => test_ig_option.py}    |    0
 ...test_option_ml.py => test_ig_option_ml.py} |    0
 test/test_install_mpmc_pyg.py                 |   85 -
 test/{test_kernels.py => test_kn_kernels.py}  |    0
 test/test_mpmc_optional_imports.py            |  100 --
 test/test_plot_and_stop.py                    |  155 --
 test/test_product_measure.py                  |  277 ----
 ...ate_data.py => test_sc_accumulate_data.py} |    0
 ...ubbayes_vec.py => test_sc_cubbayes_vec.py} |    0
 ...iteria.py => test_sc_stopping_criteria.py} |    0
 test/test_scipy_wrapper_custom.py             |  302 ----
 test/test_sr_check_links.py                   |  198 +++
 test/test_sr_check_removed_urls.py            |  165 ++
 test/test_sr_flatten_qmcpy_imports.py         |  351 +++++
 test/test_sr_install_mpmc_pyg.py              |   92 ++
 test/test_sr_mpmc_optional_imports.py         |  114 ++
 test/test_sr_unwrap_markdown.py               |   99 ++
 test/test_tm_copulas.py                       | 1271 +++++++++++++++
 test/test_tm_product_measure.py               |  271 ++++
 test/test_tm_scipy_wrapper_custom.py          |  292 ++++
 ...e_measures.py => test_tm_true_measures.py} |    0
 test/test_unwrap_markdown.py                  |   89 --
 test/test_ut_plot_and_stop.py                 |  167 ++
 test/{test_util.py => test_ut_util.py}        |    0
 50 files changed, 3556 insertions(+), 3303 deletions(-)
 delete mode 100644 qmcpy/accumulate_data/__init__.py
 create mode 100755 scripts/check_test_style.py
 delete mode 100644 test/test_check_links.py
 delete mode 100644 test/test_check_removed_urls.py
 delete mode 100644 test/test_copulas.py
 rename test/{test_discrete_distribs.py => test_dd_discrete_distribs.py} (100%)
 create mode 100644 test/test_dd_dummy_sampler.py
 delete mode 100644 test/test_dummy_sampler.py
 rename test/{test_integrate.py => test_ee_integrate.py} (100%)
 rename test/{test_keister.py => test_ee_keister.py} (100%)
 rename test/{test_pi_problem.py => test_ee_pi_problem.py} (100%)
 delete mode 100644 test/test_fast_transform_fallbacks.py
 delete mode 100644 test/test_financial_option_quick.py
 delete mode 100644 test/test_flatten_qmcpy_imports.py
 create mode 100644 test/test_ft_fast_transform_fallbacks.py
 create mode 100644 test/test_ig_financial_option_quick.py
 rename test/{test_integrands.py => test_ig_integrands.py} (100%)
 rename test/{test_option.py => test_ig_option.py} (100%)
 rename test/{test_option_ml.py => test_ig_option_ml.py} (100%)
 delete mode 100644 test/test_install_mpmc_pyg.py
 rename test/{test_kernels.py => test_kn_kernels.py} (100%)
 delete mode 100644 test/test_mpmc_optional_imports.py
 delete mode 100644 test/test_plot_and_stop.py
 delete mode 100644 test/test_product_measure.py
 rename test/{test_accumulate_data.py => test_sc_accumulate_data.py} (100%)
 rename test/{test_cubbayes_vec.py => test_sc_cubbayes_vec.py} (100%)
 rename test/{test_stopping_criteria.py => test_sc_stopping_criteria.py} (100%)
 delete mode 100644 test/test_scipy_wrapper_custom.py
 create mode 100644 test/test_sr_check_links.py
 create mode 100644 test/test_sr_check_removed_urls.py
 create mode 100644 test/test_sr_flatten_qmcpy_imports.py
 create mode 100644 test/test_sr_install_mpmc_pyg.py
 create mode 100644 test/test_sr_mpmc_optional_imports.py
 create mode 100644 test/test_sr_unwrap_markdown.py
 create mode 100644 test/test_tm_copulas.py
 create mode 100644 test/test_tm_product_measure.py
 create mode 100644 test/test_tm_scipy_wrapper_custom.py
 rename test/{test_true_measures.py => test_tm_true_measures.py} (100%)
 delete mode 100644 test/test_unwrap_markdown.py
 create mode 100644 test/test_ut_plot_and_stop.py
 rename test/{test_util.py => test_ut_util.py} (100%)

diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index a4dc7a59a..63feb062a 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -60,10 +60,7 @@ While `dev` contains the most complete set of install dependencies, a number of
 pip install -e ".[dev]"
 ~~~
 
-The `dev` extra includes QMCPy's PyPI-hosted MPMC dependencies. MPMC additionally
-requires a platform-specific `pyg_lib` wheel that is not available from PyPI.
-After installing `dev`, let the QMCPy installer select the wheel page matching
-the installed PyTorch build:
+The `dev` extra includes QMCPy's PyPI-hosted MPMC dependencies. MPMC additionally requires a platform-specific `pyg_lib` wheel that is not available from PyPI. After installing `dev`, let the QMCPy installer select the wheel page matching the installed PyTorch build:
 
 ~~~bash
 qmcpy-install-mpmc
@@ -145,6 +142,10 @@ make tests
 
 Please see the targets in the makefile for more granular control over tests.
 
+### Test file layout
+
+Unit tests live flat in `test/`, named `test_<area>_<topic>.py` where `<area>` is a short code for the `qmcpy` subpackage under test (`tm` true_measure, `dd` discrete_distribution, `sc` stopping_criterion, `ig` integrand, ...) or a cross-cutting bucket (`ee`, `sr`). So `pytest test/ -k test_tm_` runs every true-measure test. New files should also be written as a `unittest.TestCase` subclass rather than bare `def test_*` functions. `make check_test_style` lists any file that breaks either convention (informational; also runs inside `make format`; `STRICT=--strict` makes it fail). The full area table is in [`test/README.md`](test/README.md#test-file-organization).
+
 ## Documentation
 
 ### Ensure `pyreverse` Is On Your PATH
diff --git a/docs/api/discrete_distributions.md b/docs/api/discrete_distributions.md
index daa782cb6..49b60fbb4 100644
--- a/docs/api/discrete_distributions.md
+++ b/docs/api/discrete_distributions.md
@@ -70,10 +70,7 @@ python -m pip install "qmcpy[mpmc]"
 qmcpy-install-mpmc
 ```
 
-The second command selects the `pyg_lib` wheel page matching the installed
-PyTorch and accelerator builds. For GPU support or platform-specific wheels,
-see the [PyTorch installation guide](https://pytorch.org/get-started/locally/)
-and the [PyTorch Geometric installation guide](https://pytorch-geometric.readthedocs.io/en/latest/install/installation.html).
+The second command selects the `pyg_lib` wheel page matching the installed PyTorch and accelerator builds. For GPU support or platform-specific wheels, see the [PyTorch installation guide](https://pytorch.org/get-started/locally/) and the [PyTorch Geometric installation guide](https://pytorch-geometric.readthedocs.io/en/latest/install/installation.html).
 
 ::: qmcpy.discrete_distribution.mpmc.mpmc.MPMC
 
diff --git a/docs/mpmc-compatibility.md b/docs/mpmc-compatibility.md
index 92d6f1526..4e9bbf424 100644
--- a/docs/mpmc-compatibility.md
+++ b/docs/mpmc-compatibility.md
@@ -36,17 +36,14 @@ This gives one place to enforce modern MPMC compatibility without forcing the en
 
 ## Local Developer Commands
 
-Install the usual test and MPMC extras first, then add the platform-specific
-PyG runtime with QMCPy's installed helper command:
+Install the usual test and MPMC extras first, then add the platform-specific PyG runtime with QMCPy's installed helper command:
 
 ```bash
 python -m pip install -e ".[test,test_torch,test_gpytorch,test_botorch,mpmc]"
 qmcpy-install-mpmc
 ```
 
-The `mpmc` extra contains dependencies available from PyPI. The helper handles
-`pyg_lib` separately because its wheel page depends on the installed PyTorch
-version and accelerator build, which standard project metadata cannot select.
+The `mpmc` extra contains dependencies available from PyPI. The helper handles `pyg_lib` separately because its wheel page depends on the installed PyTorch version and accelerator build, which standard project metadata cannot select.
 
 Then run the MPMC-specific checks:
 
diff --git a/docs/tests.md b/docs/tests.md
index 1120146dc..73fb226be 100644
--- a/docs/tests.md
+++ b/docs/tests.md
@@ -18,6 +18,48 @@ This document describes the available test targets in the Makefile for QMCSoftwa
 | `make delcoverage` | Reset coverage tracking | Instant | Start fresh coverage analysis |
 
 
+## Test File Organization
+
+Unit tests live flat in `test/` (no subpackage subfolders). Every file is named:
+
+```
+test_<area>_<topic>.py
+```
+
+`<area>` is a short code for the `qmcpy` subpackage under test, or a cross-cutting bucket:
+
+| area | scope |
+|------|-------|
+| `dd` | `qmcpy/discrete_distribution` |
+| `ft` | `qmcpy/fast_transform` |
+| `ig` | `qmcpy/integrand` |
+| `kn` | `qmcpy/kernel` |
+| `sc` | `qmcpy/stopping_criterion` |
+| `tm` | `qmcpy/true_measure` |
+| `ut` | `qmcpy/util` |
+| `ee` | end-to-end / cross-cutting pipeline (`integrate()`, worked problems such as Keister and pi) |
+| `sr` | `scripts/` tooling, packaging, and docs checks |
+
+This keeps related tests adjacent when the directory is sorted, and lets you run one area at a time:
+
+```bash
+python -m pytest test/ -k test_tm_      # every true_measure test
+make unittests PYTEST_EXTRA_ARGS="-k test_sc_"
+```
+
+Notebook tests are separate: they live in `test/booktests/` as `tb_*.py` and are generated from `demos/` (see `test/booktests/README.md`).
+
+### Conventions checked by `make check_test_style`
+
+1. **Area prefix** — the filename must start with a recognized `test_<area>_` prefix from the table above.
+2. **Object class** — write a test file as one or more `unittest.TestCase` subclasses rather than bare `def test_*` pytest functions. A class groups related assertions under a name (so `pytest -k TestCubMCG` selects them and a failure report names the group), shares construction through `setUp` / `setUpClass` / `self.addCleanup`, and runs identically under `pytest`, `python -m unittest`, and the coverage and booktest runners without depending on pytest fixtures. Most of the suite already follows this; a few legacy files still use bare functions and new files should not.
+
+`make check_test_style` lists any violation and is informational (exit 0). It also runs as part of `make format`. To make it fail instead — for a pre-commit hook or CI gate — pass `--strict`:
+
+```bash
+STRICT=--strict make check_test_style
+```
+
 ## Detailed Descriptions
 
 ## Scope
diff --git a/makefile b/makefile
index 57e4f66e6..924c0e9b5 100644
--- a/makefile
+++ b/makefile
@@ -40,6 +40,16 @@ clean_local_only_files:
 clean_coverage:
 	rm -fr artifacts/coverage/ .coverage* test/booktests/.coverage*
 
+TEST_STYLE_PATH ?= test
+# Check test/test_*.py against two suite conventions: (1) written as a
+# unittest.TestCase subclass ("object class"), not bare pytest functions;
+# (2) named test_<area>_*.py where <area> is the qmcpy subpackage under test
+# (dd ft ig kn sc tm ut) or a cross-cutting bucket (ee sr).
+# Informational by default; pass --strict to make it fail
+# (e.g. STRICT=--strict make check_test_style).
+check_test_style:
+	@$(PYTHON) scripts/check_test_style.py $(TEST_STYLE_PATH) $(STRICT)
+
 ##########################################################
 # Doctests
 ##########################################################
@@ -406,8 +416,13 @@ MARKDOWN_UNWRAP_PATH ?= $(FORMAT_PATH)
 
 format:
 	$(MAKE) flatten_qmcpy_imports
+	@echo ""
 	$(MAKE) markdown-unwrap MARKDOWN_UNWRAP_PATH="$(MARKDOWN_UNWRAP_PATH)"
+	@echo ""
 	$(MAKE) rm_trailing_whitespace FORMAT_PATH="$(FORMAT_PATH)"
+	@echo ""
+	$(MAKE) check_test_style
+	@echo ""
 
 flatten_qmcpy_imports:
 	$(PYTHON) scripts/flatten_qmcpy_imports.py
diff --git a/qmcpy/accumulate_data/__init__.py b/qmcpy/accumulate_data/__init__.py
deleted file mode 100644
index 6f913ef81..000000000
--- a/qmcpy/accumulate_data/__init__.py
+++ /dev/null
@@ -1 +0,0 @@
-"""Accumulate data module."""
diff --git a/scripts/check_test_style.py b/scripts/check_test_style.py
new file mode 100755
index 000000000..9fdda9fa0
--- /dev/null
+++ b/scripts/check_test_style.py
@@ -0,0 +1,131 @@
+#!/usr/bin/env python3
+"""Check ``test/test_*.py`` files against two suite conventions.
+
+1. **Object class.** A test file should be written as a ``unittest.TestCase``
+   subclass, not as bare ``def test_*`` pytest functions.  The numeric-correctness
+   backbone (``test_tm_true_measures.py``, ``test_sc_stopping_criteria.py``,
+   ``test_dd_discrete_distribs.py``, ...) already follows this; newer per-measure
+   and tooling files do not.  The split is listed so it stays visible in review.
+
+2. **Area prefix.** A test file should be named ``test_<area>_<rest>.py`` where
+   ``<area>`` marks the ``qmcpy`` subpackage under test (or a cross-cutting
+   bucket).  Recognized areas:
+
+       dd    discrete_distribution    tm    true_measure
+       ft    fast_transform           ut    util
+       ig    integrand                ee    end-to-end / cross-cutting pipeline
+       kn    kernel                   sr    scripts/ tooling, packaging, docs checks
+       sc    stopping_criterion
+
+Usage:
+    python scripts/check_test_style.py [TEST_DIR] [--strict] [--quiet]
+
+TEST_DIR defaults to ``test``.  With ``--strict`` the exit code is non-zero when
+any file violates either convention (so it can gate CI); otherwise it is always
+0 and the output is informational.
+"""
+import ast
+import re
+import sys
+from pathlib import Path
+
+AREA_PREFIXES = {
+    "dd": "discrete_distribution",
+    "ft": "fast_transform",
+    "ig": "integrand",
+    "kn": "kernel",
+    "sc": "stopping_criterion",
+    "tm": "true_measure",
+    "ut": "util",
+    "ee": "end-to-end / cross-cutting pipeline",
+    "sr": "scripts/ tooling, packaging, docs checks",
+}
+AREA_RE = re.compile(r"^test_(?:" + "|".join(sorted(AREA_PREFIXES)) + r")_.+\.py$")
+
+
+def _area_ok(path):
+    """True if the filename starts with a recognized ``test_<area>_`` prefix."""
+    return bool(AREA_RE.match(path.name))
+
+
+def _subclasses_testcase(node):
+    """True if a ClassDef lists ``TestCase`` / ``unittest.TestCase`` as a base."""
+    for base in node.bases:
+        if isinstance(base, ast.Attribute) and base.attr == "TestCase":
+            return True
+        if isinstance(base, ast.Name) and base.id == "TestCase":
+            return True
+    return False
+
+
+def classify(path):
+    """Return (has_testcase_class, has_test_callables)."""
+    tree = ast.parse(path.read_text(encoding="utf-8"), filename=str(path))
+    nodes = list(ast.walk(tree))
+    has_class = any(
+        isinstance(n, ast.ClassDef) and _subclasses_testcase(n) for n in nodes
+    )
+    has_tests = any(
+        isinstance(n, (ast.FunctionDef, ast.AsyncFunctionDef))
+        and n.name.startswith("test_")
+        for n in nodes
+    )
+    return has_class, has_tests
+
+
+def main(argv):
+    strict = "--strict" in argv
+    quiet = "--quiet" in argv
+    positional = [a for a in argv if not a.startswith("-")]
+    test_dir = Path(positional[0]) if positional else Path("test")
+
+    files = sorted(test_dir.glob("test_*.py"))
+    if not files:
+        print(f"no test_*.py files under {test_dir}/", file=sys.stderr)
+        return 1
+
+    class_based, function_based, no_tests = [], [], []
+    for f in files:
+        has_class, has_tests = classify(f)
+        if has_class:
+            class_based.append(f)
+        elif has_tests:
+            function_based.append(f)
+        else:
+            no_tests.append(f)
+
+    misnamed = [f for f in files if not _area_ok(f)]
+
+    if not quiet:
+        print(f"{len(class_based)}/{len(files)} file(s) use a unittest.TestCase class")
+    if function_based:
+        print(
+            f"{len(function_based)} file(s) use bare pytest functions "
+            f"(no unittest.TestCase class):"
+        )
+        for f in function_based:
+            print(f"  {f.as_posix()}")
+    elif not quiet:
+        print("  no bare-function test files found")
+    if no_tests and not quiet:
+        print(f"{len(no_tests)} file(s) define no test_* callables:")
+        for f in no_tests:
+            print(f"  {f.as_posix()}")
+
+    if not quiet:
+        print(
+            f"{len(files) - len(misnamed)}/{len(files)} file(s) use a "
+            f"test_<area>_ prefix ({', '.join(sorted(AREA_PREFIXES))})"
+        )
+    if misnamed:
+        print(f"  {len(misnamed)} file(s) have no recognized test_<area>_ prefix:")
+        for f in misnamed:
+            print(f"  {f.as_posix()}")
+    elif not quiet:
+        print("  no misnamed test files found")
+
+    return 1 if (strict and (function_based or misnamed)) else 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main(sys.argv[1:]))
diff --git a/scripts/flatten_qmcpy_imports.py b/scripts/flatten_qmcpy_imports.py
index 033952768..f48a5cbd1 100644
--- a/scripts/flatten_qmcpy_imports.py
+++ b/scripts/flatten_qmcpy_imports.py
@@ -907,7 +907,7 @@ def main(argv: list[str] | None = None) -> int:
             file=sys.stderr,
         )
 
-    changed_files = 0
+    changed = []  # list of (display_path, import_count)
     changed_imports = 0
     for path in targets:
         original = path.read_bytes()
@@ -919,29 +919,26 @@ def main(argv: list[str] | None = None) -> int:
         if not count:
             continue
 
-        changed_files += 1
+        changed.append((_display_path(path, repository_root), count))
         changed_imports += count
         if not args.check:
             path.write_bytes(updated)
-        action = "Would update" if args.check else "Updated"
-        import_label = "import" if count == 1 else "imports"
-        print(
-            f"{action}: {_display_path(path, repository_root)} "
-            f"({count} {import_label})"
-        )
 
-    if changed_files:
-        action = "need updates" if args.check else "updated"
+    action = "would update" if args.check else "updated"
+    if changed:
+        file_label = "file" if len(changed) == 1 else "files"
         import_label = "import" if changed_imports == 1 else "imports"
-        file_label = "file" if changed_files == 1 else "files"
         print(
-            f"{changed_imports} {import_label} in "
-            f"{changed_files} {file_label} {action}."
+            f"qmcpy imports {action}: {len(changed)} {file_label}, "
+            f"{changed_imports} {import_label}:"
         )
+        for display_path, count in sorted(changed):
+            per = "import" if count == 1 else "imports"
+            print(f"  {display_path} ({count} {per})")
     else:
-        print("All eligible QMCPy imports already use the top-level package.")
+        print(f"  qmcpy imports {action}: 0 files")
 
-    return int(args.check and changed_files > 0)
+    return int(args.check and bool(changed))
 
 
 if __name__ == "__main__":
diff --git a/scripts/remove_trailing_whitespace.py b/scripts/remove_trailing_whitespace.py
index e8def8d17..bbdc65f58 100644
--- a/scripts/remove_trailing_whitespace.py
+++ b/scripts/remove_trailing_whitespace.py
@@ -125,11 +125,17 @@ def main() -> int:
     parser.add_argument("paths", nargs="+", help="tracked files or directories to process")
     args = parser.parse_args()
 
-    changed = [
-        path for path in iter_source_files(args.paths) if remove_trailing_whitespace(path, args.check)
-    ]
+    changed = sorted(
+        path for path in iter_source_files(args.paths)
+        if remove_trailing_whitespace(path, args.check)
+    )
     action = "would update" if args.check else "updated"
-    print(f"trailing whitespace {action}: {len(changed)} file(s)")
+    if changed:
+        print(f"trailing whitespace {action}: {len(changed)} file(s):")
+        for path in changed:
+            print(f"  {path}")
+    else:
+        print(f"  trailing whitespace {action}: 0 file(s)")
     return int(args.check and bool(changed))
 
 
diff --git a/scripts/unwrap_markdown.py b/scripts/unwrap_markdown.py
index 7841d8cc7..b02004c8f 100755
--- a/scripts/unwrap_markdown.py
+++ b/scripts/unwrap_markdown.py
@@ -280,23 +280,31 @@ def main() -> int:
         print("error: no .md or .ipynb files found", file=sys.stderr)
         return 2
 
-    changed_files = 0
+    changed_paths = []
     changed_cells = 0
     for path in targets:
         suffix = path.suffix.lower()
         if suffix == ".md":
-            changed = process_markdown_file(path, args.check)
-            changed_files += int(changed)
+            if process_markdown_file(path, args.check):
+                changed_paths.append(path)
         elif suffix == ".ipynb":
             changed, cell_count = process_notebook(path, args.check)
-            changed_files += int(changed)
+            if changed:
+                changed_paths.append(path)
             changed_cells += cell_count
 
     mode = "would update" if args.check else "updated"
-    print(
-        f"markdown unwrap {mode}: {changed_files} file(s), {changed_cells} markdown cell(s)",
+    summary = (
+        f"markdown unwrap {mode}: {len(changed_paths)} file(s), "
+        f"{changed_cells} markdown cell(s)"
     )
-    return 1 if args.check and changed_files else 0
+    if changed_paths:
+        print(summary + ":")
+        for path in sorted(changed_paths):
+            print(f"  {path}")
+    else:
+        print("  " + summary)
+    return 1 if args.check and changed_paths else 0
 
 
 if __name__ == "__main__":
diff --git a/test/README.md b/test/README.md
index 1120146dc..73fb226be 100644
--- a/test/README.md
+++ b/test/README.md
@@ -18,6 +18,48 @@ This document describes the available test targets in the Makefile for QMCSoftwa
 | `make delcoverage` | Reset coverage tracking | Instant | Start fresh coverage analysis |
 
 
+## Test File Organization
+
+Unit tests live flat in `test/` (no subpackage subfolders). Every file is named:
+
+```
+test_<area>_<topic>.py
+```
+
+`<area>` is a short code for the `qmcpy` subpackage under test, or a cross-cutting bucket:
+
+| area | scope |
+|------|-------|
+| `dd` | `qmcpy/discrete_distribution` |
+| `ft` | `qmcpy/fast_transform` |
+| `ig` | `qmcpy/integrand` |
+| `kn` | `qmcpy/kernel` |
+| `sc` | `qmcpy/stopping_criterion` |
+| `tm` | `qmcpy/true_measure` |
+| `ut` | `qmcpy/util` |
+| `ee` | end-to-end / cross-cutting pipeline (`integrate()`, worked problems such as Keister and pi) |
+| `sr` | `scripts/` tooling, packaging, and docs checks |
+
+This keeps related tests adjacent when the directory is sorted, and lets you run one area at a time:
+
+```bash
+python -m pytest test/ -k test_tm_      # every true_measure test
+make unittests PYTEST_EXTRA_ARGS="-k test_sc_"
+```
+
+Notebook tests are separate: they live in `test/booktests/` as `tb_*.py` and are generated from `demos/` (see `test/booktests/README.md`).
+
+### Conventions checked by `make check_test_style`
+
+1. **Area prefix** — the filename must start with a recognized `test_<area>_` prefix from the table above.
+2. **Object class** — write a test file as one or more `unittest.TestCase` subclasses rather than bare `def test_*` pytest functions. A class groups related assertions under a name (so `pytest -k TestCubMCG` selects them and a failure report names the group), shares construction through `setUp` / `setUpClass` / `self.addCleanup`, and runs identically under `pytest`, `python -m unittest`, and the coverage and booktest runners without depending on pytest fixtures. Most of the suite already follows this; a few legacy files still use bare functions and new files should not.
+
+`make check_test_style` lists any violation and is informational (exit 0). It also runs as part of `make format`. To make it fail instead — for a pre-commit hook or CI gate — pass `--strict`:
+
+```bash
+STRICT=--strict make check_test_style
+```
+
 ## Detailed Descriptions
 
 ## Scope
diff --git a/test/test_check_links.py b/test/test_check_links.py
deleted file mode 100644
index 6a2b3dc52..000000000
--- a/test/test_check_links.py
+++ /dev/null
@@ -1,178 +0,0 @@
-import ssl
-import sys
-import urllib.error
-from unittest.mock import patch
-
-from scripts import check_links
-
-
-def _http_error(url, code):
-    return urllib.error.HTTPError(url, code, "test response", {}, None)
-
-
-def test_head_success_is_reachable():
-    with patch.object(check_links.urllib.request, "urlopen", return_value=object()) as urlopen:
-        assert check_links._check_one("https://example.test", timeout=1) is None
-
-    assert urlopen.call_count == 1
-    assert urlopen.call_args.args[0].get_method() == "HEAD"
-
-
-def test_get_success_after_head_failure_is_reachable():
-    url = "https://example.test"
-    with patch.object(
-        check_links.urllib.request,
-        "urlopen",
-        side_effect=[_http_error(url, 405), object()],
-    ) as urlopen:
-        assert check_links._check_one(url, timeout=1) is None
-
-    assert urlopen.call_count == 2
-    assert urlopen.call_args_list[1].args[0].get_method() == "GET"
-
-
-def test_not_found_and_gone_gets_are_broken():
-    for code in (404, 410):
-        url = f"https://example.test/{code}"
-        with patch.object(
-            check_links.urllib.request,
-            "urlopen",
-            side_effect=[_http_error(url, code), _http_error(url, code)],
-        ):
-            assert check_links._check_one(url, timeout=1) == (
-                "broken",
-                f"{url} -- HTTP {code}",
-            )
-
-
-def test_bot_block_and_rate_limit_are_warnings():
-    for code in (403, 429):
-        url = f"https://example.test/{code}"
-        with patch.object(
-            check_links.urllib.request,
-            "urlopen",
-            side_effect=[_http_error(url, code), _http_error(url, code)],
-        ):
-            severity, message = check_links._check_one(url, timeout=1)
-
-        assert severity == "warning"
-        assert f"HTTP {code}" in message
-
-
-def test_tls_and_timeout_failures_are_warnings():
-    failures = (
-        ssl.SSLCertVerificationError("certificate verify failed"),
-        TimeoutError("timed out"),
-    )
-    for failure in failures:
-        with patch.object(
-            check_links.urllib.request,
-            "urlopen",
-            side_effect=[failure, failure],
-        ):
-            severity, message = check_links._check_one(
-                "https://example.test", timeout=1
-            )
-
-        assert severity == "warning"
-        assert str(failure) in message
-
-
-def test_external_results_are_separated_and_duplicate_urls_checked_once(tmp_path):
-    (tmp_path / "page.html").write_text(
-        '<a href="https://example.test/missing">missing</a>'
-        '<a href="https://example.test/missing">duplicate</a>'
-        '<a href="https://example.test/blocked">blocked</a>',
-        encoding="utf-8",
-    )
-
-    def result_for(url, _timeout):
-        if url.endswith("/missing"):
-            return "broken", f"{url} -- HTTP 404"
-        return "warning", f"{url} -- HTTP 403"
-
-    with patch.object(check_links, "_check_one", side_effect=result_for) as check_one:
-        broken, warnings = check_links.check_external(tmp_path, workers=1)
-
-    assert check_one.call_count == 2
-    assert broken == [
-        "https://example.test/missing -- HTTP 404 (seen on page.html)"
-    ]
-    assert warnings == [
-        "https://example.test/blocked -- HTTP 403 (seen on page.html)"
-    ]
-
-
-def test_internal_links_strip_site_url_deployment_path(tmp_path):
-    target = tmp_path / "target"
-    target.mkdir()
-    (target / "index.html").write_text(
-        '<h2 id="section">Target</h2>', encoding="utf-8"
-    )
-    (tmp_path / "index.html").write_text(
-        '<a href="/QMCSoftware/target/">root-relative</a>'
-        '<a href="https://qmcsoftware.github.io/QMCSoftware/target/#section">absolute</a>',
-        encoding="utf-8",
-    )
-
-    assert (
-        check_links.check_internal(
-            tmp_path, site_url="https://qmcsoftware.github.io/QMCSoftware/"
-        )
-        == []
-    )
-
-
-def test_external_check_skips_same_site_urls(tmp_path):
-    (tmp_path / "page.html").write_text(
-        '<a href="https://qmcsoftware.github.io/QMCSoftware/target/">same</a>'
-        '<a href="https://example.test/target/">external</a>',
-        encoding="utf-8",
-    )
-
-    with patch.object(check_links, "_check_one", return_value=None) as check_one:
-        broken, warnings = check_links.check_external(
-            tmp_path,
-            workers=1,
-            site_url="https://qmcsoftware.github.io/QMCSoftware/",
-        )
-
-    assert broken == []
-    assert warnings == []
-    assert check_one.call_count == 1
-    assert check_one.call_args.args[0] == "https://example.test/target/"
-
-
-def test_external_warnings_do_not_make_main_fail(tmp_path, monkeypatch, capsys):
-    monkeypatch.setattr(sys, "argv", ["check_links.py", str(tmp_path), "--external"])
-    monkeypatch.setattr(
-        check_links, "check_internal", lambda _site_dir, site_url=None: []
-    )
-    monkeypatch.setattr(
-        check_links,
-        "check_external",
-        lambda _site_dir, site_url=None: (
-            [],
-            ["https://example.test -- HTTP 403"],
-        ),
-    )
-
-    assert check_links.main() == 0
-    assert "0 broken link(s), 1 warning(s)" in capsys.readouterr().out
-
-
-def test_confirmed_external_breakage_makes_main_fail(tmp_path, monkeypatch):
-    monkeypatch.setattr(sys, "argv", ["check_links.py", str(tmp_path), "--external"])
-    monkeypatch.setattr(
-        check_links, "check_internal", lambda _site_dir, site_url=None: []
-    )
-    monkeypatch.setattr(
-        check_links,
-        "check_external",
-        lambda _site_dir, site_url=None: (
-            ["https://example.test -- HTTP 404"],
-            [],
-        ),
-    )
-
-    assert check_links.main() == 1
diff --git a/test/test_check_removed_urls.py b/test/test_check_removed_urls.py
deleted file mode 100644
index 9374e34ad..000000000
--- a/test/test_check_removed_urls.py
+++ /dev/null
@@ -1,134 +0,0 @@
-import sys
-import urllib.error
-from unittest.mock import patch
-
-from scripts import check_removed_urls as cru
-
-SITE = "https://qmcsoftware.github.io/QMCSoftware/"
-
-
-def _sitemap(*paths):
-    locs = "".join(f"<loc>{SITE}{path}</loc>" for path in paths)
-    return f'<?xml version="1.0" encoding="UTF-8"?><urlset>{locs}</urlset>'
-
-
-def _config(redirect_maps=None):
-    plugins = ["material/search", {"mkdocs-jupyter": {"execute": False}}]
-    if redirect_maps is not None:
-        plugins.append({"redirects": {"redirect_maps": redirect_maps}})
-    return {"site_url": SITE, "plugins": plugins}
-
-
-def _run(tmp_path, monkeypatch, sitemap_paths, redirect_maps=None, extra_argv=()):
-    """Run main() offline against a temp sitemap and a temp docs/ tree."""
-    docs = tmp_path / "docs"
-    docs.mkdir(parents=True)
-    (docs / "README.md").write_text("home", encoding="utf-8")
-    (docs / "good_practices.md").write_text("page", encoding="utf-8")
-    sitemap = tmp_path / "sitemap.xml"
-    sitemap.write_text(_sitemap(*sitemap_paths), encoding="utf-8")
-
-    monkeypatch.setattr(cru, "read_config", lambda *a, **k: _config(redirect_maps))
-    monkeypatch.setattr(sys, "argv", [
-        "check_removed_urls.py", "--sitemap", str(sitemap), "--docs-dir", str(docs),
-        *extra_argv,
-    ])
-    return cru.main()
-
-
-def test_url_path_and_source_round_trip(tmp_path):
-    for source, url_path in [("blogs/scipywrapper/index.md", "blogs/scipywrapper/"),
-                             ("good_practices.md", "good_practices/"),
-                             ("demos/quickstart.ipynb", "demos/quickstart/"),
-                             ("index.md", ""), ("README.md", "")]:
-        assert cru.url_path_for_source(source) == url_path
-
-    for source in ("README.md", "good_practices.md", "demos/quickstart.ipynb",
-                   "api/index.md"):
-        path = tmp_path / source
-        path.parent.mkdir(parents=True, exist_ok=True)
-        path.write_text("page", encoding="utf-8")
-        assert cru.source_exists(cru.url_path_for_source(source), tmp_path)
-    assert not cru.source_exists("blogs/scipywrapper/", tmp_path)
-
-
-def test_redirect_maps_reads_the_plugin_and_tolerates_its_absence():
-    entry = {"blogs/x/index.md": "https://qmcsoftware.org/blogs/x/"}
-    assert cru.redirect_maps(_config(entry)) == entry
-    assert cru.redirect_maps(_config()) == {}
-    assert cru.redirect_maps({}) == {}
-
-
-def test_published_paths_separates_foreign_urls():
-    sitemap = _sitemap("", "good_practices/").replace(
-        "</urlset>", "<loc>https://example.test/other/</loc></urlset>")
-
-    assert cru.published_paths(sitemap, SITE) == (
-        ["", "good_practices/"], ["https://example.test/other/"])
-
-
-def test_http_status_falls_back_to_get_when_head_is_unsupported():
-    url = "https://example.test"
-    error = urllib.error.HTTPError(url, 405, "test response", {}, None)
-    response = type("Response", (), {"status": 200, "__enter__": lambda s: s,
-                                     "__exit__": lambda s, *a: False})()
-    with patch.object(cru.urllib.request, "urlopen",
-                      side_effect=[error, response]) as urlopen:
-        assert cru.http_status(url, timeout=1) == "200"
-
-    assert urlopen.call_count == 2
-    assert urlopen.call_args_list[1].args[0].get_method() == "GET"
-
-
-def test_removed_page_without_redirect_is_flagged(tmp_path, monkeypatch, capsys):
-    code = _run(tmp_path, monkeypatch, ["", "good_practices/", "blogs/scipywrapper/"])
-    out = capsys.readouterr().out
-
-    assert code == 1
-    assert "1 removed with no redirect" in out
-    assert f"[ORPHAN]   {SITE}blogs/scipywrapper/" in out
-    assert "blogs/scipywrapper/index.md: <new URL or page>" in out
-
-
-def test_removed_page_covered_by_a_redirect_passes(tmp_path, monkeypatch, capsys):
-    code = _run(
-        tmp_path, monkeypatch, ["", "good_practices/", "blogs/scipywrapper/"],
-        redirect_maps={
-            "blogs/scipywrapper/index.md": "https://qmcsoftware.org/blogs/scipywrapper/"},
-    )
-    out = capsys.readouterr().out
-
-    assert code == 0
-    assert "0 removed with no redirect" in out
-    assert "[redirect]" in out and "[ORPHAN]" not in out
-
-
-def test_intact_site_passes(tmp_path, monkeypatch, capsys):
-    assert _run(tmp_path, monkeypatch, ["", "good_practices/"]) == 0
-    assert "2 still have a page source" in capsys.readouterr().out
-
-
-def test_verify_redirects_follows_the_target_status(tmp_path, monkeypatch, capsys):
-    redirects = {"blogs/x/index.md": "https://qmcsoftware.org/blogs/x/"}
-    for status, expected_code in [("200", 0), ("404", 1)]:
-        monkeypatch.setattr(cru, "http_status", lambda *a, **k: status)
-        code = _run(tmp_path / status, monkeypatch, ["", "blogs/x/"],
-                    redirect_maps=redirects, extra_argv=("--verify-redirects",))
-        out = capsys.readouterr().out
-
-        assert code == expected_code
-        assert status in out
-        # The URL itself is covered, so a failure is the target, not an orphan.
-        assert "[ORPHAN]" not in out
-
-
-def test_unreachable_sitemap_fails_unless_offline_is_allowed(tmp_path, monkeypatch, capsys):
-    monkeypatch.setattr(cru, "read_config", lambda *a, **k: _config())
-    argv = ["check_removed_urls.py", "--sitemap", str(tmp_path / "absent.xml")]
-
-    monkeypatch.setattr(sys, "argv", argv)
-    assert cru.main() == 1
-
-    monkeypatch.setattr(sys, "argv", argv + ["--allow-offline"])
-    assert cru.main() == 0
-    assert "skipping the check" in capsys.readouterr().out
diff --git a/test/test_copulas.py b/test/test_copulas.py
deleted file mode 100644
index 2f4bd06ba..000000000
--- a/test/test_copulas.py
+++ /dev/null
@@ -1,1362 +0,0 @@
-import warnings
-
-import numpy as np
-import pytest
-import scipy.stats as stats
-
-from qmcpy import (
-    AbstractCopula,
-    ClaytonCopula,
-    DigitalNetB2,
-    FrankCopula,
-    GaussianCopula,
-    GumbelCopula,
-    StudentTCopula,
-)
-
-from qmcpy.true_measure.copula import (
-    AbstractCopula as ModuleAbstractCopula,
-    _apply_marginal_ppfs,
-    _build_marginal_range,
-    _clip_unit_interval,
-    _marginal_cdfs_and_logpdf,
-    _validate_correlation_matrix,
-    _validate_dimension,
-    _validate_marginals,
-)
-
-from qmcpy.util import DimensionError, MethodImplementationError, ParameterError
-
-
-class PPFOnlyMarginal:
-    def ppf(self, u):
-        return np.asarray(u, dtype=float)
-
-
-class NonCallablePPFMarginal:
-    ppf = 1.0
-
-
-class UnitPDFMarginal:
-    def ppf(self, u):
-        return np.asarray(u, dtype=float)
-
-    def cdf(self, x):
-        return np.asarray(x, dtype=float)
-
-    def pdf(self, x):
-        return np.ones_like(np.asarray(x, dtype=float))
-
-
-class CDFOnlyMarginal(PPFOnlyMarginal):
-    def cdf(self, x):
-        return np.asarray(x, dtype=float)
-
-
-class BadIntervalMarginal(PPFOnlyMarginal):
-    def interval(self, confidence):
-        raise ValueError("interval unavailable")
-
-
-class BadRangeMarginal:
-    def ppf(self, u):
-        raise ValueError("ppf unavailable")
-
-
-def _equicorrelation(d, rho):
-    corr = np.full((d, d), rho, dtype=float)
-    np.fill_diagonal(corr, 1.0)
-    return corr
-
-
-def _make_copula(copula_cls, dimension=2, marginals=None, correlation=None, seed=7):
-    if marginals is None:
-        marginals = [stats.norm()] * dimension
-    if correlation is None:
-        correlation = np.eye(dimension)
-
-    kwargs = {}
-    if copula_cls is StudentTCopula:
-        kwargs["df"] = 4
-    if copula_cls is ClaytonCopula:
-        kwargs["theta"] = 2.0
-    if copula_cls is FrankCopula:
-        kwargs["theta"] = 5.0
-    if copula_cls is GumbelCopula:
-        kwargs["theta"] = 2.0
-
-    common = {
-        "sampler": DigitalNetB2(dimension, seed=seed),
-        "marginals": marginals,
-        **kwargs,
-    }
-    if copula_cls in [ClaytonCopula, FrankCopula, GumbelCopula]:
-        return copula_cls(**common)
-    return copula_cls(correlation=correlation, **common)
-
-
-# Base AbstractCopula and helper tests
-
-
-def test_abstract_copula_is_importable_from_public_module_path():
-    assert ModuleAbstractCopula is AbstractCopula
-
-
-def test_public_api_imports_and_normal_usage():
-    for copula_cls in [
-        GaussianCopula,
-        StudentTCopula,
-        ClaytonCopula,
-        FrankCopula,
-        GumbelCopula,
-    ]:
-        assert issubclass(copula_cls, AbstractCopula)
-
-        tm = _make_copula(copula_cls)
-        x = tm(8)
-        x_gen = tm.gen_samples(8)
-        v = tm.gen_copula_samples(8)
-
-        assert x.shape == (8, 2)
-        assert x_gen.shape == (8, 2)
-        assert v.shape == (8, 2)
-        assert np.all(np.isfinite(x))
-        assert np.all(np.isfinite(x_gen))
-        assert np.all((0 <= v) & (v <= 1))
-
-
-def test_abstract_copula_rejects_unimplemented_transform():
-    tm = AbstractCopula(
-        DigitalNetB2(2, seed=101),
-        marginals=[stats.uniform(), stats.uniform()],
-    )
-
-    with pytest.raises(MethodImplementationError):
-        tm.copula_transform(np.full((3, 2), 0.5))
-
-
-def test_abstract_copula_rejects_invalid_sampler():
-    with pytest.raises(ParameterError, match="sampler"):
-        AbstractCopula(object(), marginals=[stats.uniform()])
-
-
-def test_validate_marginals_error_branches():
-    with pytest.raises(ParameterError, match="marginals"):
-        _validate_marginals(None)
-
-    with pytest.raises(ParameterError, match="at least one"):
-        _validate_marginals([])
-
-    with pytest.raises(ParameterError, match="ppf"):
-        _validate_marginals([NonCallablePPFMarginal()])
-
-
-def test_validate_dimension_error_branches():
-    with pytest.raises(DimensionError, match="integer dimension"):
-        _validate_dimension(object(), [stats.uniform()])
-
-    with pytest.raises(DimensionError, match="marginals"):
-        _validate_dimension(3, [stats.uniform(), stats.uniform()])
-
-
-def test_apply_marginal_ppfs_clips_endpoints_and_checks_dimension():
-    transformed = _apply_marginal_ppfs(
-        np.array([[0.0, 1.0], [1.0, 0.0]]),
-        [stats.norm(), stats.norm()],
-    )
-
-    assert transformed.shape == (2, 2)
-    assert np.all(np.isfinite(transformed))
-
-    with pytest.raises(DimensionError, match="marginals"):
-        _apply_marginal_ppfs(np.full((2, 3), 0.5), [stats.uniform(), stats.uniform()])
-
-
-def test_marginal_range_falls_back_when_interval_or_ppf_fails():
-    ranges = _build_marginal_range([BadIntervalMarginal(), BadRangeMarginal()])
-
-    assert ranges.shape == (2, 2)
-    assert np.all(np.isfinite(ranges[0]))
-    np.testing.assert_allclose(ranges[1], [-np.inf, np.inf])
-
-
-def test_marginal_cdfs_and_logpdf_pdf_branch_and_errors():
-    x = np.array([[0.25, 0.75], [0.4, 0.6]])
-    u, log_density = _marginal_cdfs_and_logpdf(
-        x,
-        [UnitPDFMarginal(), UnitPDFMarginal()],
-    )
-
-    np.testing.assert_allclose(u, x)
-    np.testing.assert_allclose(log_density, np.zeros(2))
-
-    with pytest.raises(ParameterError, match="cdf"):
-        _marginal_cdfs_and_logpdf(x, [PPFOnlyMarginal(), UnitPDFMarginal()])
-
-    with pytest.raises(ParameterError, match="pdf"):
-        _marginal_cdfs_and_logpdf(x, [CDFOnlyMarginal(), UnitPDFMarginal()])
-
-
-def test_validate_correlation_matrix_rejects_nonfinite_values():
-    with pytest.raises(ValueError, match="finite"):
-        _validate_correlation_matrix([[1.0, np.nan], [np.nan, 1.0]], 2)
-
-
-def test_clip_unit_interval_uses_machine_epsilon():
-    clipped = _clip_unit_interval(np.array([0.0, 0.5, 1.0]))
-    eps = np.finfo(float).eps
-
-    np.testing.assert_allclose(clipped, [eps, 0.5, 1.0 - eps])
-
-
-@pytest.mark.parametrize(
-    "copula_cls",
-    [GaussianCopula, StudentTCopula, ClaytonCopula, GumbelCopula, FrankCopula],
-)
-def test_copula_transform_outputs_dependent_uniforms_in_unit_cube(copula_cls):
-    tm = _make_copula(copula_cls, dimension=3)
-    u = np.array(
-        [
-            [0.1, 0.3, 0.7],
-            [0.5, 0.5, 0.5],
-            [0.9, 0.8, 0.2],
-        ]
-    )
-
-    v = tm.copula_transform(u)
-
-    assert v.shape == u.shape
-    assert np.all(np.isfinite(v))
-    assert np.all((0.0 <= v) & (v <= 1.0))
-
-
-@pytest.mark.parametrize(
-    "copula_cls,dimension",
-    [
-        (GaussianCopula, 3),
-        (StudentTCopula, 3),
-        (ClaytonCopula, 3),
-        (FrankCopula, 3),
-        (GumbelCopula, 3),
-    ],
-)
-def test_copula_sample_shapes_are_preserved(copula_cls, dimension):
-    tm = _make_copula(copula_cls, dimension=dimension, seed=9)
-
-    one = tm(1)
-    many = tm(8)
-    batched_transform = tm._transform(np.full((2, 3, dimension), 0.5))
-
-    assert one.shape == (1, dimension)
-    assert many.shape == (8, dimension)
-    assert batched_transform.shape == (2, 3, dimension)
-    assert np.all(np.isfinite(one))
-    assert np.all(np.isfinite(many))
-    assert np.all(np.isfinite(batched_transform))
-
-
-# Elliptical copulas
-
-
-def test_output_shape_with_nonnormal_marginals():
-    tm = GaussianCopula(
-        sampler=DigitalNetB2(2, seed=7),
-        marginals=[stats.beta(a=2, b=5), stats.gamma(a=3, scale=2)],
-        correlation=[[1.0, 0.4], [0.4, 1.0]],
-    )
-
-    x = tm(16)
-
-    assert x.shape == (16, 2)
-
-
-def test_finite_output_for_normal_marginals():
-    tm = GaussianCopula(
-        sampler=DigitalNetB2(2, seed=11),
-        marginals=[stats.norm(), stats.norm(loc=1.0, scale=2.0)],
-        correlation=[[1.0, -0.3], [-0.3, 1.0]],
-    )
-
-    x = tm(128)
-
-    assert np.all(np.isfinite(x))
-
-
-def test_return_weights_shape_when_marginal_densities_available():
-    tm = GaussianCopula(
-        sampler=DigitalNetB2(2, seed=12),
-        marginals=[stats.norm(), stats.gamma(a=2.0)],
-        correlation=[[1.0, 0.25], [0.25, 1.0]],
-    )
-
-    x, weights = tm(32, return_weights=True)
-
-    assert x.shape == (32, 2)
-    assert weights.shape == (32,)
-    assert np.all(np.isfinite(weights))
-    assert np.all(weights > 0.0)
-
-
-def test_identity_correlation_matches_independent_marginal_transforms():
-    marginals = [stats.norm(loc=-1.0, scale=2.0), stats.gamma(a=2.0, scale=3.0)]
-    tm = GaussianCopula(
-        sampler=DigitalNetB2(2, seed=13),
-        marginals=marginals,
-        correlation=np.eye(2),
-    )
-    u = np.array([[0.2, 0.7], [0.4, 0.8], [0.9, 0.1]])
-
-    x = tm._transform(u)
-    expected = np.column_stack(
-        [marginals[j].ppf(u[:, j]) for j in range(len(marginals))]
-    )
-
-    np.testing.assert_allclose(x, expected, rtol=1e-12, atol=1e-12)
-
-
-def test_positive_correlation_produces_positive_dependence():
-    rho = 0.75
-    tm = GaussianCopula(
-        sampler=DigitalNetB2(2, seed=17),
-        marginals=[stats.norm(), stats.norm()],
-        correlation=[[1.0, rho], [rho, 1.0]],
-    )
-
-    x = tm(4096)
-    empirical_corr = np.corrcoef(x.T)[0, 1]
-
-    assert empirical_corr > 0.5
-    assert abs(empirical_corr - rho) < 0.2
-
-
-@pytest.mark.parametrize("copula_cls", [GaussianCopula, StudentTCopula])
-@pytest.mark.parametrize("dimension", [1, 3, 5])
-def test_elliptical_copulas_support_general_dimensions(copula_cls, dimension):
-    correlation = _equicorrelation(dimension, 0.25)
-    tm = _make_copula(
-        copula_cls,
-        dimension=dimension,
-        marginals=[stats.norm()] * dimension,
-        correlation=correlation,
-        seed=19,
-    )
-
-    x = tm(16)
-    one = tm(1)
-
-    assert x.shape == (16, dimension)
-    assert one.shape == (1, dimension)
-    assert np.all(np.isfinite(x))
-    assert np.all(np.isfinite(one))
-
-
-@pytest.mark.parametrize("copula_cls", [GaussianCopula, StudentTCopula])
-def test_elliptical_copulas_handle_valid_near_singular_correlation(copula_cls):
-    dimension = 5
-    tm = _make_copula(
-        copula_cls,
-        dimension=dimension,
-        marginals=[stats.norm()] * dimension,
-        correlation=_equicorrelation(dimension, 0.999),
-        seed=20,
-    )
-
-    x = tm(32)
-
-    assert x.shape == (32, dimension)
-    assert np.all(np.isfinite(x))
-
-
-@pytest.mark.parametrize("copula_cls", [GaussianCopula, StudentTCopula])
-def test_elliptical_copulas_reject_singular_correlation(copula_cls):
-    with pytest.raises(ValueError, match="positive definite"):
-        _make_copula(
-            copula_cls,
-            dimension=3,
-            marginals=[stats.norm(), stats.norm(), stats.norm()],
-            correlation=np.ones((3, 3)),
-            seed=22,
-        )
-
-
-@pytest.mark.parametrize(
-    "copula_cls",
-    [GaussianCopula, StudentTCopula, ClaytonCopula, FrankCopula, GumbelCopula],
-)
-def test_distribution_dimension_matches_number_of_marginals(copula_cls):
-    tm = _make_copula(
-        copula_cls,
-        dimension=5,
-        marginals=[
-            stats.norm(),
-            stats.beta(a=2, b=5),
-            stats.gamma(a=3),
-            stats.expon(),
-            stats.lognorm(s=0.5),
-        ],
-        correlation=np.eye(5),
-    )
-
-    x = tm(32)
-
-    assert x.shape == (32, 5)
-    assert np.all(np.isfinite(x))
-
-
-@pytest.mark.parametrize("copula_cls", [GaussianCopula, StudentTCopula])
-def test_invalid_dimension_mismatches_raise(copula_cls):
-    with pytest.raises(DimensionError, match="marginals"):
-        _make_copula(
-            copula_cls,
-            dimension=2,
-            marginals=[stats.norm(), stats.norm(), stats.norm()],
-            correlation=np.eye(2),
-        )
-
-    with pytest.raises(ValueError, match="shape"):
-        _make_copula(
-            copula_cls,
-            dimension=2,
-            marginals=[stats.norm(), stats.norm()],
-            correlation=np.eye(3),
-        )
-
-    with pytest.raises(ValueError, match="square"):
-        _make_copula(
-            copula_cls,
-            dimension=2,
-            marginals=[stats.norm(), stats.norm()],
-            correlation=[[1.0, 0.2, 0.3], [0.2, 1.0, 0.4]],
-        )
-
-
-@pytest.mark.parametrize("copula_cls", [ClaytonCopula, FrankCopula, GumbelCopula])
-def test_archimedean_dimension_mismatch_raises_dimension_error(copula_cls):
-    with pytest.raises(DimensionError, match="marginals"):
-        _make_copula(
-            copula_cls,
-            dimension=2,
-            marginals=[stats.norm(), stats.norm(), stats.norm()],
-        )
-
-
-@pytest.mark.parametrize(
-    "copula_cls",
-    [GaussianCopula, StudentTCopula],
-)
-@pytest.mark.parametrize(
-    "correlation",
-    [
-        [[1.0, 0.2], [0.3, 1.0]],
-        [[1.0, 0.2], [0.2, 0.9]],
-        [[1.0, 1.2], [1.2, 1.0]],
-    ],
-)
-def test_invalid_correlation_matrices_raise_value_error(copula_cls, correlation):
-    with pytest.raises(ValueError):
-        _make_copula(
-            copula_cls,
-            dimension=2,
-            marginals=[stats.norm(), stats.norm()],
-            correlation=correlation,
-        )
-
-
-def test_marginal_length_mismatch_raises_dimension_error():
-    with pytest.raises(DimensionError, match="marginals"):
-        GaussianCopula(
-            sampler=DigitalNetB2(2, seed=21),
-            marginals=[stats.norm()],
-            correlation=np.eye(2),
-        )
-
-
-def test_marginal_without_ppf_raises_clear_error():
-    class NoPPF:
-        pass
-
-    with pytest.raises(ParameterError, match="ppf"):
-        GaussianCopula(
-            sampler=DigitalNetB2(1, seed=23),
-            marginals=[NoPPF()],
-            correlation=[[1.0]],
-        )
-
-
-@pytest.mark.parametrize(
-    "copula_cls",
-    [GaussianCopula, StudentTCopula, ClaytonCopula, FrankCopula, GumbelCopula],
-)
-def test_common_scipy_frozen_marginals_work(copula_cls):
-    tm = _make_copula(
-        copula_cls,
-        dimension=5,
-        marginals=[
-            stats.norm(),
-            stats.beta(a=2, b=5),
-            stats.gamma(a=3),
-            stats.expon(),
-            stats.lognorm(s=0.5),
-        ],
-        correlation=np.eye(5),
-        seed=47,
-    )
-
-    x = tm(128)
-
-    assert x.shape == (128, 5)
-    assert np.all(np.isfinite(x))
-
-
-@pytest.mark.parametrize(
-    "copula_cls",
-    [GaussianCopula, StudentTCopula, ClaytonCopula, FrankCopula, GumbelCopula],
-)
-def test_endpoint_uniforms_are_clipped_to_finite_outputs(copula_cls):
-    tm = _make_copula(
-        copula_cls,
-        dimension=5,
-        marginals=[
-            stats.norm(),
-            stats.beta(a=2, b=5),
-            stats.gamma(a=3),
-            stats.expon(),
-            stats.lognorm(s=0.5),
-        ],
-        correlation=np.eye(5),
-        seed=53,
-    )
-    u = np.array(
-        [
-            [0.0, 1.0, 0.0, 1.0, 0.5],
-            [1.0, 0.0, 1.0, 0.0, 0.5],
-        ]
-    )
-
-    x = tm._transform(u)
-
-    assert x.shape == (2, 5)
-    assert np.all(np.isfinite(x))
-
-
-def test_student_t_copula_output_shape_and_finite_values():
-    tm = StudentTCopula(
-        sampler=DigitalNetB2(2, seed=29),
-        marginals=[stats.norm(), stats.gamma(a=3.0, scale=2.0)],
-        correlation=[[1.0, 0.5], [0.5, 1.0]],
-        df=4,
-    )
-
-    x = tm(128)
-
-    assert x.shape == (128, 2)
-    assert np.all(np.isfinite(x))
-
-
-def test_student_t_copula_positive_correlation_produces_positive_dependence():
-    tm = StudentTCopula(
-        sampler=DigitalNetB2(2, seed=31),
-        marginals=[stats.norm(), stats.norm()],
-        correlation=[[1.0, 0.7], [0.7, 1.0]],
-        df=5,
-    )
-
-    x = tm(4096)
-    empirical_corr = np.corrcoef(x.T)[0, 1]
-
-    assert empirical_corr > 0.45
-
-
-def test_student_t_copula_has_stronger_joint_tail_than_gaussian_copula():
-    rho = 0.7
-    df = 4
-    n = 2**12
-    marginals = [stats.norm(), stats.norm()]
-    correlation = [[1.0, rho], [rho, 1.0]]
-
-    gaussian = GaussianCopula(
-        sampler=DigitalNetB2(2, seed=101),
-        marginals=marginals,
-        correlation=correlation,
-    )
-    student_t = StudentTCopula(
-        sampler=DigitalNetB2(2, seed=101),
-        marginals=marginals,
-        correlation=correlation,
-        df=df,
-    )
-
-    x_gaussian = gaussian(n)
-    x_student_t = student_t(n)
-    threshold = stats.norm.ppf(0.99)
-
-    def joint_tail_rate(x):
-        tail_0 = x[:, 0] > threshold
-        return np.mean(x[tail_0, 1] > threshold)
-
-    gaussian_tail = joint_tail_rate(x_gaussian)
-    student_t_tail = joint_tail_rate(x_student_t)
-
-    assert student_t_tail > gaussian_tail + 0.08
-
-
-def test_student_t_copula_return_weights_shape_when_density_available():
-    tm = StudentTCopula(
-        sampler=DigitalNetB2(2, seed=37),
-        marginals=[stats.norm(), stats.gamma(a=2.0)],
-        correlation=[[1.0, 0.3], [0.3, 1.0]],
-        df=6,
-    )
-
-    x, weights = tm(32, return_weights=True)
-
-    assert x.shape == (32, 2)
-    assert weights.shape == (32,)
-    assert np.all(np.isfinite(weights))
-    assert np.all(weights > 0.0)
-
-
-@pytest.mark.parametrize("df", [1.0, 100.0])
-def test_student_t_copula_boundary_df_values_are_finite(df):
-    dimension = 3
-    tm = StudentTCopula(
-        sampler=DigitalNetB2(dimension, seed=39),
-        marginals=[stats.norm()] * dimension,
-        correlation=_equicorrelation(dimension, 0.4),
-        df=df,
-    )
-
-    x = tm(128)
-
-    assert x.shape == (128, dimension)
-    assert np.all(np.isfinite(x))
-
-
-def test_student_t_copula_large_df_is_close_to_gaussian_copula():
-    rho = 0.6
-    correlation = [[1.0, rho], [rho, 1.0]]
-    marginals = [stats.norm(), stats.norm()]
-    gaussian = GaussianCopula(
-        sampler=DigitalNetB2(2, seed=40),
-        marginals=marginals,
-        correlation=correlation,
-    )
-    student_t = StudentTCopula(
-        sampler=DigitalNetB2(2, seed=40),
-        marginals=marginals,
-        correlation=correlation,
-        df=100,
-    )
-
-    x_gaussian = gaussian(4096)
-    x_student_t = student_t(4096)
-    corr_gaussian = np.corrcoef(x_gaussian.T)[0, 1]
-    corr_student_t = np.corrcoef(x_student_t.T)[0, 1]
-
-    assert abs(corr_student_t - corr_gaussian) < 0.02
-
-
-@pytest.mark.parametrize("df", [0, -1, np.inf, "not-a-number"])
-def test_student_t_copula_invalid_df_raises_parameter_error(df):
-    with pytest.raises(ParameterError, match="df"):
-        StudentTCopula(
-            sampler=DigitalNetB2(2, seed=41),
-            marginals=[stats.norm(), stats.norm()],
-            correlation=np.eye(2),
-            df=df,
-        )
-
-
-def test_student_t_copula_marginal_without_ppf_raises_clear_error():
-    class NoPPF:
-        pass
-
-    with pytest.raises(ParameterError, match="ppf"):
-        StudentTCopula(
-            sampler=DigitalNetB2(1, seed=43),
-            marginals=[NoPPF()],
-            correlation=[[1.0]],
-            df=4,
-        )
-
-
-# Archimedean copulas
-
-
-def test_clayton_copula_output_shape_and_finite_values():
-    tm = ClaytonCopula(
-        sampler=DigitalNetB2(2, seed=57),
-        marginals=[stats.norm(), stats.gamma(a=3.0, scale=2.0)],
-        theta=2.0,
-    )
-
-    x = tm(128)
-
-    assert x.shape == (128, 2)
-    assert np.all(np.isfinite(x))
-
-
-def test_clayton_copula_return_weights_shape_when_density_available():
-    tm = ClaytonCopula(
-        sampler=DigitalNetB2(3, seed=59),
-        marginals=[stats.norm(), stats.gamma(a=2.0), stats.expon()],
-        theta=1.5,
-    )
-
-    x, weights = tm(32, return_weights=True)
-
-    assert x.shape == (32, 3)
-    assert weights.shape == (32,)
-    assert np.all(np.isfinite(weights))
-    assert np.all(weights > 0.0)
-
-
-@pytest.mark.parametrize("theta", [0, -1, np.inf, "not-a-number"])
-def test_clayton_copula_invalid_theta_raises_parameter_error(theta):
-    with pytest.raises(ParameterError, match="theta"):
-        ClaytonCopula(
-            sampler=DigitalNetB2(2, seed=61),
-            marginals=[stats.norm(), stats.norm()],
-            theta=theta,
-        )
-
-
-@pytest.mark.parametrize("dimension", [2, 3, 5])
-def test_clayton_copula_supports_general_dimension(dimension):
-    tm = ClaytonCopula(
-        sampler=DigitalNetB2(dimension, seed=63),
-        marginals=[stats.norm()] * dimension,
-        theta=2.0,
-    )
-
-    x = tm(128)
-
-    assert x.shape == (128, dimension)
-    assert np.all(np.isfinite(x))
-
-
-def test_clayton_copula_marginal_without_ppf_raises_clear_error():
-    class NoPPF:
-        pass
-
-    with pytest.raises(ParameterError, match="ppf"):
-        ClaytonCopula(
-            sampler=DigitalNetB2(2, seed=67),
-            marginals=[stats.norm(), NoPPF()],
-            theta=2.0,
-        )
-
-
-@pytest.mark.parametrize(
-    "marginals",
-    [
-        [stats.norm(), stats.beta(a=2, b=5)],
-        [stats.gamma(a=3), stats.expon()],
-        [stats.lognorm(s=0.5), stats.norm()],
-    ],
-)
-def test_clayton_copula_common_scipy_frozen_marginals_work(marginals):
-    tm = ClaytonCopula(
-        sampler=DigitalNetB2(2, seed=69),
-        marginals=marginals,
-        theta=2.0,
-    )
-
-    x = tm(128)
-
-    assert x.shape == (128, 2)
-    assert np.all(np.isfinite(x))
-
-
-def test_clayton_copula_endpoint_uniforms_are_clipped_to_finite_outputs():
-    tm = ClaytonCopula(
-        sampler=DigitalNetB2(2, seed=70),
-        marginals=[stats.norm(), stats.lognorm(s=0.5)],
-        theta=2.0,
-    )
-    u = np.array([[0.0, 1.0], [1.0, 0.0]])
-
-    x = tm._transform(u)
-
-    assert x.shape == (2, 2)
-    assert np.all(np.isfinite(x))
-
-
-@pytest.mark.parametrize("dimension", [2, 3, 5])
-def test_clayton_copula_tiny_theta_is_near_independent(dimension):
-    marginals = [stats.uniform()] * dimension
-    tm = ClaytonCopula(
-        sampler=DigitalNetB2(dimension, seed=70),
-        marginals=marginals,
-        theta=1e-8,
-    )
-    u = np.array(
-        [
-            [0.2, 0.7, 0.4, 0.6, 0.8],
-            [0.4, 0.8, 0.9, 0.3, 0.2],
-            [0.9, 0.1, 0.3, 0.7, 0.5],
-        ]
-    )[:, :dimension]
-
-    x = tm._transform(u)
-
-    assert x.shape == (3, dimension)
-    assert np.all(np.isfinite(x))
-    np.testing.assert_allclose(x, u, atol=5e-6)
-
-
-@pytest.mark.parametrize("dimension", [2, 3, 5])
-@pytest.mark.parametrize("theta", [20.0, 50.0])
-def test_clayton_copula_large_theta_is_finite(theta, dimension):
-    tm = ClaytonCopula(
-        sampler=DigitalNetB2(dimension, seed=70),
-        marginals=[stats.norm()] * dimension,
-        theta=theta,
-    )
-
-    x = tm(128)
-
-    assert x.shape == (128, dimension)
-    assert np.all(np.isfinite(x))
-
-
-def test_clayton_copula_positive_dependence_behavior():
-    tm = ClaytonCopula(
-        sampler=DigitalNetB2(2, seed=71),
-        marginals=[stats.uniform(), stats.uniform()],
-        theta=2.0,
-    )
-
-    x = tm(4096)
-    empirical_corr = np.corrcoef(x.T)[0, 1]
-
-    assert empirical_corr > 0.45
-
-
-def test_clayton_copula_has_stronger_lower_tail_than_gaussian_copula():
-    theta = 2.0
-    n = 2**12
-    marginals = [stats.uniform(), stats.uniform()]
-    # Clayton Kendall tau is theta/(theta+2); convert to Gaussian rho.
-    rho = np.sin(np.pi * (theta / (theta + 2.0)) / 2.0)
-
-    clayton = ClaytonCopula(
-        sampler=DigitalNetB2(2, seed=73),
-        marginals=marginals,
-        theta=theta,
-    )
-    gaussian = GaussianCopula(
-        sampler=DigitalNetB2(2, seed=73),
-        marginals=marginals,
-        correlation=[[1.0, rho], [rho, 1.0]],
-    )
-
-    x_clayton = clayton(n)
-    x_gaussian = gaussian(n)
-    threshold = 0.05
-
-    def lower_tail_rate(x):
-        tail_0 = x[:, 0] < threshold
-        return np.mean(x[tail_0, 1] < threshold)
-
-    clayton_tail = lower_tail_rate(x_clayton)
-    gaussian_tail = lower_tail_rate(x_gaussian)
-
-    assert clayton_tail > gaussian_tail + 0.2
-
-
-def test_frank_copula_output_shape_for_two_dimensions():
-    tm = FrankCopula(
-        sampler=DigitalNetB2(2, seed=75),
-        marginals=[stats.norm(), stats.gamma(a=3.0, scale=2.0)],
-        theta=5.0,
-    )
-
-    x = tm(128)
-
-    assert x.shape == (128, 2)
-    assert np.all(np.isfinite(x))
-
-
-@pytest.mark.parametrize("dimension", [3, 5])
-def test_frank_copula_positive_theta_supports_higher_dimensions(dimension):
-    tm = FrankCopula(
-        sampler=DigitalNetB2(dimension, seed=76),
-        marginals=[stats.norm()] * dimension,
-        theta=5.0,
-    )
-
-    x = tm(128)
-
-    assert x.shape == (128, dimension)
-    assert np.all(np.isfinite(x))
-
-
-def test_frank_copula_return_weights_shape_when_density_available():
-    tm = FrankCopula(
-        sampler=DigitalNetB2(3, seed=77),
-        marginals=[stats.norm(), stats.gamma(a=2.0), stats.expon()],
-        theta=4.0,
-    )
-
-    x, weights = tm(32, return_weights=True)
-
-    assert x.shape == (32, 3)
-    assert weights.shape == (32,)
-    assert np.all(np.isfinite(weights))
-    assert np.all(weights > 0.0)
-
-
-@pytest.mark.parametrize("theta", [0, np.inf, -np.inf, "not-a-number"])
-def test_frank_copula_invalid_theta_raises_parameter_error(theta):
-    with pytest.raises(ParameterError, match="theta"):
-        FrankCopula(
-            sampler=DigitalNetB2(2, seed=78),
-            marginals=[stats.norm(), stats.norm()],
-            theta=theta,
-        )
-
-
-def test_frank_copula_negative_theta_rejected_above_two_dimensions():
-    with pytest.raises(ParameterError, match="d=2"):
-        FrankCopula(
-            sampler=DigitalNetB2(3, seed=79),
-            marginals=[stats.norm(), stats.norm(), stats.norm()],
-            theta=-2.0,
-        )
-
-
-def test_frank_copula_dimension_mismatch_raises_dimension_error():
-    with pytest.raises(DimensionError, match="marginals"):
-        FrankCopula(
-            sampler=DigitalNetB2(2, seed=80),
-            marginals=[stats.norm(), stats.norm(), stats.norm()],
-            theta=5.0,
-        )
-
-
-def test_frank_copula_marginal_without_ppf_raises_clear_error():
-    class NoPPF:
-        pass
-
-    with pytest.raises(ParameterError, match="ppf"):
-        FrankCopula(
-            sampler=DigitalNetB2(2, seed=82),
-            marginals=[stats.norm(), NoPPF()],
-            theta=5.0,
-        )
-
-
-def test_frank_copula_positive_dependence_behavior():
-    tm = FrankCopula(
-        sampler=DigitalNetB2(2, seed=84),
-        marginals=[stats.uniform(), stats.uniform()],
-        theta=6.0,
-    )
-
-    x = tm(4096)
-    empirical_corr = np.corrcoef(x.T)[0, 1]
-
-    assert empirical_corr > 0.45
-
-
-@pytest.mark.parametrize(
-    "theta,dimension",
-    [
-        (1e-8, 3),
-        (-1e-8, 2),
-    ],
-)
-def test_frank_copula_tiny_theta_is_close_to_independence(theta, dimension):
-    marginals = [stats.uniform()] * dimension
-    tm = FrankCopula(
-        sampler=DigitalNetB2(dimension, seed=86),
-        marginals=marginals,
-        theta=theta,
-    )
-    u = np.array(
-        [
-            [0.2, 0.7, 0.4, 0.6, 0.8],
-            [0.4, 0.8, 0.9, 0.3, 0.2],
-            [0.9, 0.1, 0.3, 0.7, 0.5],
-        ]
-    )[:, :dimension]
-
-    x = tm._transform(u)
-
-    assert x.shape == (3, dimension)
-    assert np.all(np.isfinite(x))
-    np.testing.assert_allclose(x, u, atol=5e-6)
-
-
-@pytest.mark.parametrize(
-    "theta,dimension",
-    [
-        (50.0, 5),
-        (-50.0, 2),
-    ],
-)
-def test_frank_copula_large_theta_is_finite(theta, dimension):
-    tm = FrankCopula(
-        sampler=DigitalNetB2(dimension, seed=87),
-        marginals=[stats.norm()] * dimension,
-        theta=theta,
-    )
-
-    x = tm(128)
-
-    assert x.shape == (128, dimension)
-    assert np.all(np.isfinite(x))
-
-
-def test_frank_copula_negative_theta_produces_negative_dependence_in_2d():
-    tm = FrankCopula(
-        sampler=DigitalNetB2(2, seed=88),
-        marginals=[stats.uniform(), stats.uniform()],
-        theta=-6.0,
-    )
-
-    x = tm(4096)
-    empirical_corr = np.corrcoef(x.T)[0, 1]
-
-    assert empirical_corr < -0.35
-
-
-def test_gumbel_copula_output_shape_and_finite_values():
-    tm = GumbelCopula(
-        sampler=DigitalNetB2(2, seed=79),
-        marginals=[stats.norm(), stats.gamma(a=3.0, scale=2.0)],
-        theta=2.0,
-    )
-
-    x = tm(128)
-
-    assert x.shape == (128, 2)
-    assert np.all(np.isfinite(x))
-
-
-def test_gumbel_copula_return_weights_shape_when_density_available():
-    tm = GumbelCopula(
-        sampler=DigitalNetB2(3, seed=81),
-        marginals=[stats.norm(), stats.gamma(a=2.0), stats.expon()],
-        theta=1.5,
-    )
-
-    x, weights = tm(32, return_weights=True)
-
-    assert x.shape == (32, 3)
-    assert weights.shape == (32,)
-    assert np.all(np.isfinite(weights))
-    assert np.all(weights > 0.0)
-
-
-@pytest.mark.parametrize("theta", [0, 0.5, -1, np.inf, "not-a-number"])
-def test_gumbel_copula_invalid_theta_raises_parameter_error(theta):
-    with pytest.raises(ParameterError, match="theta"):
-        GumbelCopula(
-            sampler=DigitalNetB2(2, seed=83),
-            marginals=[stats.norm(), stats.norm()],
-            theta=theta,
-        )
-
-
-def test_gumbel_copula_theta_one_is_independent_marginal_transform():
-    marginals = [stats.norm(loc=-1.0, scale=2.0), stats.gamma(a=2.0, scale=3.0)]
-    tm = GumbelCopula(
-        sampler=DigitalNetB2(2, seed=85),
-        marginals=marginals,
-        theta=1.0,
-    )
-    u = np.array([[0.2, 0.7], [0.4, 0.8], [0.9, 0.1]])
-
-    x = tm._transform(u)
-    expected = np.column_stack(
-        [marginals[j].ppf(u[:, j]) for j in range(len(marginals))]
-    )
-
-    np.testing.assert_allclose(x, expected, rtol=1e-12, atol=1e-12)
-
-
-@pytest.mark.parametrize("dimension", [2, 3, 5])
-def test_gumbel_copula_theta_close_to_one_is_near_independent(dimension):
-    marginals = [stats.uniform()] * dimension
-    tm = GumbelCopula(
-        sampler=DigitalNetB2(dimension, seed=85),
-        marginals=marginals,
-        theta=1.000001,
-    )
-    u = np.array(
-        [
-            [0.2, 0.7, 0.4, 0.6, 0.8],
-            [0.4, 0.8, 0.9, 0.3, 0.2],
-            [0.9, 0.1, 0.3, 0.7, 0.5],
-        ]
-    )[:, :dimension]
-
-    x = tm._transform(u)
-
-    assert x.shape == (3, dimension)
-    assert np.all(np.isfinite(x))
-    np.testing.assert_allclose(x, u, atol=5e-5)
-
-
-@pytest.mark.parametrize("dimension", [2, 3, 5])
-@pytest.mark.parametrize("theta", [20.0, 50.0])
-def test_gumbel_copula_large_theta_is_finite(theta, dimension):
-    tm = GumbelCopula(
-        sampler=DigitalNetB2(dimension, seed=86),
-        marginals=[stats.norm()] * dimension,
-        theta=theta,
-    )
-
-    x = tm(128)
-
-    assert x.shape == (128, dimension)
-    assert np.all(np.isfinite(x))
-
-
-@pytest.mark.parametrize("dimension", [2, 3, 5])
-def test_gumbel_copula_supports_general_dimension(dimension):
-    tm = GumbelCopula(
-        sampler=DigitalNetB2(dimension, seed=87),
-        marginals=[stats.norm()] * dimension,
-        theta=2.0,
-    )
-
-    x = tm(128)
-
-    assert x.shape == (128, dimension)
-    assert np.all(np.isfinite(x))
-
-
-def test_gumbel_copula_marginal_without_ppf_raises_clear_error():
-    class NoPPF:
-        pass
-
-    with pytest.raises(ParameterError, match="ppf"):
-        GumbelCopula(
-            sampler=DigitalNetB2(2, seed=89),
-            marginals=[stats.norm(), NoPPF()],
-            theta=2.0,
-        )
-
-
-@pytest.mark.parametrize(
-    "marginals",
-    [
-        [stats.norm(), stats.beta(a=2, b=5)],
-        [stats.gamma(a=3), stats.expon()],
-        [stats.lognorm(s=0.5), stats.norm()],
-    ],
-)
-def test_gumbel_copula_common_scipy_frozen_marginals_work(marginals):
-    tm = GumbelCopula(
-        sampler=DigitalNetB2(2, seed=91),
-        marginals=marginals,
-        theta=2.0,
-    )
-
-    x = tm(128)
-
-    assert x.shape == (128, 2)
-    assert np.all(np.isfinite(x))
-
-
-def test_gumbel_copula_endpoint_uniforms_are_clipped_to_finite_outputs():
-    tm = GumbelCopula(
-        sampler=DigitalNetB2(2, seed=93),
-        marginals=[stats.norm(), stats.lognorm(s=0.5)],
-        theta=2.0,
-    )
-    u = np.array([[0.0, 1.0], [1.0, 0.0]])
-
-    x = tm._transform(u)
-
-    assert x.shape == (2, 2)
-    assert np.all(np.isfinite(x))
-
-
-def test_gumbel_copula_positive_dependence_behavior():
-    tm = GumbelCopula(
-        sampler=DigitalNetB2(2, seed=95),
-        marginals=[stats.uniform(), stats.uniform()],
-        theta=2.0,
-    )
-
-    x = tm(4096)
-    empirical_corr = np.corrcoef(x.T)[0, 1]
-
-    assert empirical_corr > 0.45
-
-
-def test_gumbel_copula_has_stronger_upper_tail_than_gaussian_copula():
-    theta = 2.0
-    n = 2**12
-    marginals = [stats.uniform(), stats.uniform()]
-    # Gumbel Kendall tau is 1 - 1/theta; convert to Gaussian rho.
-    rho = np.sin(np.pi * (1.0 - 1.0 / theta) / 2.0)
-
-    gumbel = GumbelCopula(
-        sampler=DigitalNetB2(2, seed=97),
-        marginals=marginals,
-        theta=theta,
-    )
-    gaussian = GaussianCopula(
-        sampler=DigitalNetB2(2, seed=97),
-        marginals=marginals,
-        correlation=[[1.0, rho], [rho, 1.0]],
-    )
-
-    x_gumbel = gumbel(n)
-    x_gaussian = gaussian(n)
-    threshold = 0.95
-
-    def upper_tail_rate(x):
-        tail_0 = x[:, 0] > threshold
-        return np.mean(x[tail_0, 1] > threshold)
-
-    gumbel_tail = upper_tail_rate(x_gumbel)
-    gaussian_tail = upper_tail_rate(x_gaussian)
-
-    assert gumbel_tail > gaussian_tail + 0.15
-
-
-# Weights, fallback behavior, spawn, and edge cases
-
-
-@pytest.mark.parametrize(
-    "copula_cls",
-    [GaussianCopula, StudentTCopula, ClaytonCopula, GumbelCopula, FrankCopula],
-)
-def test_copula_weight_fallback_warns_once_when_density_methods_are_missing(
-    copula_cls,
-):
-    tm = _make_copula(
-        copula_cls,
-        dimension=2,
-        marginals=[PPFOnlyMarginal(), PPFOnlyMarginal()],
-    )
-    x = np.full((4, 2), 0.5)
-    expected_message = getattr(
-        tm,
-        "_missing_weight_warning_message",
-        f"{copula_cls.__name__} marginals must implement 'cdf' and "
-        "'pdf' or 'logpdf' to compute density weights. "
-        "Weights will be treated as 1.",
-    )
-
-    assert "_unit_weight_with_warning" not in copula_cls.__dict__
-    assert (
-        tm._unit_weight_with_warning.__func__
-        is AbstractCopula._unit_weight_with_warning
-    )
-
-    with pytest.warns(UserWarning) as warning_info:
-        weights = tm._weight(x)
-
-    with warnings.catch_warnings(record=True) as caught:
-        warnings.simplefilter("always")
-        second_weights = tm._weight(x)
-
-    np.testing.assert_allclose(weights, np.ones(4))
-    np.testing.assert_allclose(second_weights, np.ones(4))
-    assert str(warning_info[0].message) == expected_message
-    assert caught == []
-
-
-def test_student_t_weight_falls_back_when_multivariate_t_is_unavailable():
-    tm = StudentTCopula(
-        DigitalNetB2(2, seed=115),
-        marginals=[stats.norm(), stats.norm()],
-        correlation=np.eye(2),
-        df=4,
-    )
-    tm._mvt_scipy = None
-
-    with pytest.warns(UserWarning, match="Weights will be treated as 1"):
-        weights = tm._weight(np.full((3, 2), 0.25))
-
-    np.testing.assert_allclose(weights, np.ones(3))
-
-
-def test_gaussian_weight_uses_pdf_branch_when_logpdf_is_unavailable():
-    tm = GaussianCopula(
-        DigitalNetB2(2, seed=117),
-        marginals=[UnitPDFMarginal(), UnitPDFMarginal()],
-        correlation=[[1.0, 0.4], [0.4, 1.0]],
-    )
-
-    weights = tm._weight(np.array([[0.25, 0.5], [0.75, 0.5]]))
-
-    assert weights.shape == (2,)
-    assert np.all(np.isfinite(weights))
-    assert np.all(weights > 0.0)
-
-
-def test_gumbel_theta_one_weight_is_independent_marginal_density():
-    tm = GumbelCopula(
-        DigitalNetB2(2, seed=119),
-        marginals=[stats.gamma(a=2.0), stats.expon()],
-        theta=1.0,
-    )
-    x = np.array([[1.0, 0.5], [2.0, 1.5]])
-    expected = stats.gamma(a=2.0).pdf(x[:, 0]) * stats.expon().pdf(x[:, 1])
-
-    weights = tm._weight(x)
-
-    np.testing.assert_allclose(weights, expected)
-
-
-def test_gen_copula_samples_composed_transform_branch():
-    inner = GaussianCopula(
-        DigitalNetB2(2, seed=121),
-        marginals=[stats.uniform(), stats.uniform()],
-        correlation=[[1.0, 0.3], [0.3, 1.0]],
-    )
-    outer = ClaytonCopula(inner, marginals=[stats.uniform(), stats.uniform()], theta=1.5)
-
-    v = outer.gen_copula_samples(n_min=4, n_max=8)
-
-    assert v.shape == (4, 2)
-    assert np.all(np.isfinite(v))
-    assert np.all((0.0 <= v) & (v <= 1.0))
-
-
-@pytest.mark.parametrize(
-    "copula_cls",
-    [GaussianCopula, StudentTCopula, ClaytonCopula, GumbelCopula, FrankCopula],
-)
-def test_copula_spawn_same_dimension_and_reject_different_dimension(copula_cls):
-    tm = _make_copula(copula_cls, dimension=2)
-
-    spawned = tm.spawn(s=1, dimensions=[2])
-    assert len(spawned) == 1
-    assert isinstance(spawned[0], copula_cls)
-    assert spawned[0](4).shape == (4, 2)
-
-    with pytest.raises(DimensionError):
-        tm._spawn(DigitalNetB2(3, seed=123), 3)
-
-
-def test_frank_one_dimensional_weight_covers_zero_order_eulerian_term():
-    tm = FrankCopula(
-        DigitalNetB2(1, seed=125),
-        marginals=[UnitPDFMarginal()],
-        theta=3.0,
-    )
-
-    weights = tm._weight(np.array([[0.25], [0.75]]))
-
-    assert weights.shape == (2,)
-    assert np.all(np.isfinite(weights))
-    assert np.all(weights > 0.0)
-
-
-def test_frank_rejects_large_negative_theta_when_exponential_overflows():
-    with np.errstate(over="ignore"):
-        with pytest.raises(ParameterError, match="too close to 0 or too large"):
-            FrankCopula(
-                DigitalNetB2(2, seed=127),
-                marginals=[stats.uniform(), stats.uniform()],
-                theta=-1000.0,
-            )
diff --git a/test/test_discrete_distribs.py b/test/test_dd_discrete_distribs.py
similarity index 100%
rename from test/test_discrete_distribs.py
rename to test/test_dd_discrete_distribs.py
diff --git a/test/test_dd_dummy_sampler.py b/test/test_dd_dummy_sampler.py
new file mode 100644
index 000000000..50d5005b8
--- /dev/null
+++ b/test/test_dd_dummy_sampler.py
@@ -0,0 +1,107 @@
+import unittest
+
+import numpy as np
+
+from qmcpy import DummySampler
+from qmcpy.util import ParameterError
+
+
+PLACEHOLDER_ERROR = "construction placeholder"
+
+
+class TestDummySampler(unittest.TestCase):
+
+    def test_dummy_sampler_constructs_dimension_one(self):
+        sampler = DummySampler(1)
+
+        self.assertEqual(sampler.d, 1)
+        self.assertEqual(sampler.replications, 1)
+        self.assertTrue(sampler.no_replications)
+        self.assertEqual(sampler.mimics, "StdUniform")
+        self.assertEqual(sampler.parameters, [])
+
+    def test_dummy_sampler_constructs_larger_dimensions(self):
+        sampler = DummySampler(3, seed=7)
+
+        self.assertEqual(sampler.d, 3)
+        self.assertEqual(sampler.replications, 1)
+        self.assertTrue(sampler.no_replications)
+        self.assertTrue(np.array_equal(sampler.dvec, np.arange(3)))
+
+    def test_dummy_sampler_constructs_larger_dimension_with_replications(self):
+        sampler = DummySampler(4, replications=3, seed=7)
+
+        self.assertEqual(sampler.d, 4)
+        self.assertEqual(sampler.replications, 3)
+        self.assertFalse(sampler.no_replications)
+        self.assertTrue(np.array_equal(sampler.dvec, np.arange(4)))
+
+    def test_dummy_sampler_direct_sampling_raises_placeholder_error(self):
+        sampler = DummySampler(2)
+
+        with self.assertRaisesRegex(ParameterError, PLACEHOLDER_ERROR):
+            sampler(8)
+
+    def test_dummy_sampler_replicated_direct_sampling_raises_placeholder_error(self):
+        sampler = DummySampler(2, replications=3)
+
+        with self.assertRaisesRegex(ParameterError, PLACEHOLDER_ERROR):
+            sampler(8)
+
+    def test_dummy_sampler_supported_calling_conventions_raise_placeholder_error(self):
+        sampler = DummySampler(2)
+
+        with self.assertRaisesRegex(ParameterError, PLACEHOLDER_ERROR):
+            sampler(n=4)
+        with self.assertRaisesRegex(ParameterError, PLACEHOLDER_ERROR):
+            sampler(n_min=2, n_max=6)
+        with self.assertRaisesRegex(ParameterError, PLACEHOLDER_ERROR):
+            sampler(n=2, n_min=6)
+
+    def test_dummy_sampler_nonzero_n_min_raises_placeholder_error(self):
+        sampler = DummySampler(2)
+
+        with self.assertRaisesRegex(ParameterError, PLACEHOLDER_ERROR):
+            sampler(n_min=5, n_max=9)
+
+    def test_dummy_sampler_rejects_return_binary(self):
+        sampler = DummySampler(2)
+
+        with self.assertRaisesRegex(ParameterError, PLACEHOLDER_ERROR):
+            sampler(4, return_binary=True)
+
+    def test_dummy_sampler_internal_gen_samples_raises_placeholder_error(self):
+        sampler = DummySampler(2)
+
+        with self.assertRaisesRegex(ParameterError, PLACEHOLDER_ERROR):
+            sampler._gen_samples(n_min=5, n_max=9, return_binary=False, warn=True)
+
+    def test_dummy_sampler_spawn_preserves_relevant_fields(self):
+        sampler = DummySampler(2, replications=3, seed=11)
+
+        spawned = sampler.spawn(s=2, dimensions=[1, 5])
+
+        self.assertEqual([spawn.d for spawn in spawned], [1, 5])
+        self.assertEqual([spawn.replications for spawn in spawned], [3, 3])
+        self.assertTrue(all(isinstance(spawn, DummySampler) for spawn in spawned))
+
+    def test_dummy_sampler_spawn_without_explicit_replications(self):
+        sampler = DummySampler(2, seed=11)
+
+        spawned = sampler.spawn(s=1, dimensions=4)[0]
+
+        self.assertEqual(spawned.d, 4)
+        self.assertEqual(spawned.replications, 1)
+        self.assertTrue(spawned.no_replications)
+
+    def test_dummy_sampler_limits_are_enforced(self):
+        with self.assertRaisesRegex(ParameterError, "dimension greater than dimension limit"):
+            DummySampler(10_002)
+
+        sampler = DummySampler(1)
+        with self.assertRaisesRegex(ParameterError, "n_limit"):
+            sampler(n_min=0, n_max=2**32 + 1)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/test_dummy_sampler.py b/test/test_dummy_sampler.py
deleted file mode 100644
index 24bec84eb..000000000
--- a/test/test_dummy_sampler.py
+++ /dev/null
@@ -1,111 +0,0 @@
-import numpy as np
-import pytest
-
-from qmcpy import DummySampler
-from qmcpy.util import ParameterError
-
-
-PLACEHOLDER_ERROR = "construction placeholder"
-
-
-def test_dummy_sampler_constructs_dimension_one():
-    sampler = DummySampler(1)
-
-    assert sampler.d == 1
-    assert sampler.replications == 1
-    assert sampler.no_replications
-    assert sampler.mimics == "StdUniform"
-    assert sampler.parameters == []
-
-
-def test_dummy_sampler_constructs_larger_dimensions():
-    sampler = DummySampler(3, seed=7)
-
-    assert sampler.d == 3
-    assert sampler.replications == 1
-    assert sampler.no_replications
-    assert np.array_equal(sampler.dvec, np.arange(3))
-
-
-def test_dummy_sampler_constructs_larger_dimension_with_replications():
-    sampler = DummySampler(4, replications=3, seed=7)
-
-    assert sampler.d == 4
-    assert sampler.replications == 3
-    assert not sampler.no_replications
-    assert np.array_equal(sampler.dvec, np.arange(4))
-
-
-def test_dummy_sampler_direct_sampling_raises_placeholder_error():
-    sampler = DummySampler(2)
-
-    with pytest.raises(ParameterError, match=PLACEHOLDER_ERROR):
-        sampler(8)
-
-
-def test_dummy_sampler_replicated_direct_sampling_raises_placeholder_error():
-    sampler = DummySampler(2, replications=3)
-
-    with pytest.raises(ParameterError, match=PLACEHOLDER_ERROR):
-        sampler(8)
-
-
-def test_dummy_sampler_supported_calling_conventions_raise_placeholder_error():
-    sampler = DummySampler(2)
-
-    with pytest.raises(ParameterError, match=PLACEHOLDER_ERROR):
-        sampler(n=4)
-    with pytest.raises(ParameterError, match=PLACEHOLDER_ERROR):
-        sampler(n_min=2, n_max=6)
-    with pytest.raises(ParameterError, match=PLACEHOLDER_ERROR):
-        sampler(n=2, n_min=6)
-
-
-def test_dummy_sampler_nonzero_n_min_raises_placeholder_error():
-    sampler = DummySampler(2)
-
-    with pytest.raises(ParameterError, match=PLACEHOLDER_ERROR):
-        sampler(n_min=5, n_max=9)
-
-
-def test_dummy_sampler_rejects_return_binary():
-    sampler = DummySampler(2)
-
-    with pytest.raises(ParameterError, match=PLACEHOLDER_ERROR):
-        sampler(4, return_binary=True)
-
-
-def test_dummy_sampler_internal_gen_samples_raises_placeholder_error():
-    sampler = DummySampler(2)
-
-    with pytest.raises(ParameterError, match=PLACEHOLDER_ERROR):
-        sampler._gen_samples(n_min=5, n_max=9, return_binary=False, warn=True)
-
-
-def test_dummy_sampler_spawn_preserves_relevant_fields():
-    sampler = DummySampler(2, replications=3, seed=11)
-
-    spawned = sampler.spawn(s=2, dimensions=[1, 5])
-
-    assert [spawn.d for spawn in spawned] == [1, 5]
-    assert [spawn.replications for spawn in spawned] == [3, 3]
-    assert all(isinstance(spawn, DummySampler) for spawn in spawned)
-
-
-def test_dummy_sampler_spawn_without_explicit_replications():
-    sampler = DummySampler(2, seed=11)
-
-    spawned = sampler.spawn(s=1, dimensions=4)[0]
-
-    assert spawned.d == 4
-    assert spawned.replications == 1
-    assert spawned.no_replications
-
-
-def test_dummy_sampler_limits_are_enforced():
-    with pytest.raises(ParameterError, match="dimension greater than dimension limit"):
-        DummySampler(10_002)
-
-    sampler = DummySampler(1)
-    with pytest.raises(ParameterError, match="n_limit"):
-        sampler(n_min=0, n_max=2**32 + 1)
diff --git a/test/test_integrate.py b/test/test_ee_integrate.py
similarity index 100%
rename from test/test_integrate.py
rename to test/test_ee_integrate.py
diff --git a/test/test_keister.py b/test/test_ee_keister.py
similarity index 100%
rename from test/test_keister.py
rename to test/test_ee_keister.py
diff --git a/test/test_pi_problem.py b/test/test_ee_pi_problem.py
similarity index 100%
rename from test/test_pi_problem.py
rename to test/test_ee_pi_problem.py
diff --git a/test/test_fast_transform_fallbacks.py b/test/test_fast_transform_fallbacks.py
deleted file mode 100644
index 2dcbd2b7a..000000000
--- a/test/test_fast_transform_fallbacks.py
+++ /dev/null
@@ -1,46 +0,0 @@
-import numpy as np
-import pytest
-
-from qmcpy import (
-    fftbr,
-    fftbr_torch,
-    fwht,
-    fwht_torch,
-    ifftbr,
-    ifftbr_torch,
-    omega_fftbr,
-    omega_fftbr_torch,
-    omega_fwht,
-    omega_fwht_torch,
-)
-
-
-def test_non_torch_transforms_basic():
-    rng = np.random.default_rng(11)
-    x = rng.random(8) + 1j * rng.random(8)
-    y = fftbr(x)
-    assert y.shape == x.shape
-    xr = ifftbr(y)
-    assert xr.shape == x.shape
-
-    a = rng.random(8)
-    b = fwht(a)
-    assert b.shape == a.shape
-
-    omega = omega_fftbr(3)
-    assert omega.shape[0] == 2**3
-    omega2 = omega_fwht(3)
-    assert omega2.shape[0] == 2**3
-
-
-def test_torch_fallbacks_raise():
-    with pytest.raises(Exception):
-        fftbr_torch()
-    with pytest.raises(Exception):
-        ifftbr_torch()
-    with pytest.raises(Exception):
-        fwht_torch()
-    with pytest.raises(Exception):
-        omega_fftbr_torch()
-    with pytest.raises(Exception):
-        omega_fwht_torch()
diff --git a/test/test_financial_option_quick.py b/test/test_financial_option_quick.py
deleted file mode 100644
index bd9f3022b..000000000
--- a/test/test_financial_option_quick.py
+++ /dev/null
@@ -1,94 +0,0 @@
-import numpy as np
-
-from qmcpy import FinancialOption
-import qmcpy
-
-
-class SmallSampler(qmcpy.AbstractDiscreteDistribution):
-    def __init__(self, d=3):
-        super().__init__(
-            dimension=d, replications=1, seed=123, d_limit=100, n_limit=1024
-        )
-
-    def _gen_samples(self, n_min, n_max, return_binary=False, warn=True):
-        n = n_max - n_min
-        # return shape (replications, n, d)
-        arr = np.tile(np.linspace(0.1, 1.0, n)[:, None], (1, self.d))
-        return arr.reshape(self.replications, n, self.d)
-
-
-def test_financial_option_payoffs_and_exact():
-    sampler = SmallSampler(d=3)
-    fo = FinancialOption(
-        sampler,
-        option="EUROPEAN",
-        call_put="CALL",
-        volatility=0.5,
-        start_price=30,
-        strike_price=25,
-        interest_rate=0.01,
-        t_final=1,
-    )
-    gbm = np.array([[30.0, 28.0, 35.0]])
-    c = fo.payoff_european_call(gbm)
-    p = fo.payoff_european_put(gbm)
-    assert c.shape == (1,)
-    assert p.shape == (1,)
-
-    # Asian arithmetic trapezoidal
-    fo_asian = FinancialOption(
-        sampler,
-        option="ASIAN",
-        asian_mean="ARITHMETIC",
-        asian_mean_quadrature_rule="TRAPEZOIDAL",
-    )
-    gbm2 = np.array([[30.0, 32.0, 34.0]])
-    a_call = fo_asian.payoff_asian_arithmetic_trap_call(gbm2)
-    assert a_call.shape == (1,)
-
-    # geometric right call
-    fo_geo = FinancialOption(
-        sampler,
-        option="ASIAN",
-        asian_mean="GEOMETRIC",
-        asian_mean_quadrature_rule="RIGHT",
-    )
-    g_call = fo_geo.payoff_asian_geometric_right_call(np.array([[30.0, 30.0, 30.0]]))
-    assert g_call.shape == (1,)
-
-    # barrier options: up and down behaviors
-    fo_barrier_up = FinancialOption(
-        sampler, option="BARRIER", barrier_in_out="IN", barrier_price=25, start_price=20
-    )
-    gbm_up = np.array([[20.0, 26.0, 27.0]])
-    v = fo_barrier_up.payoff_barrier_in_up_call(gbm_up)
-    assert v.shape == (1,)
-
-    fo_barrier_out = FinancialOption(
-        sampler,
-        option="BARRIER",
-        barrier_in_out="OUT",
-        barrier_price=40,
-        start_price=30,
-    )
-    gbm_out = np.array([[30.0, 32.0, 33.0]])
-    v2 = fo_barrier_out.payoff_barrier_out_up_call(gbm_out)
-    assert v2.shape == (1,)
-
-    # lookback
-    fo_lb = FinancialOption(sampler, option="LOOKBACK")
-    lb = fo_lb.payoff_lookback_call(np.array([[10.0, 9.0, 12.0]]))
-    assert lb.shape == (1,)
-
-    # digital
-    fo_dig = FinancialOption(sampler, option="DIGITAL", digital_payout=5)
-    dig = fo_dig.payoff_digital_call(np.array([[10.0, 11.0, 12.0]]))
-    assert dig.shape == (1,)
-
-    # exact value for European should return a float
-    val = fo.get_exact_value()
-    assert np.isscalar(val)
-
-    # exact value for Asian geometric right
-    val2 = fo_geo.get_exact_value()
-    assert np.isscalar(val2)
diff --git a/test/test_flatten_qmcpy_imports.py b/test/test_flatten_qmcpy_imports.py
deleted file mode 100644
index c7fc36fdb..000000000
--- a/test/test_flatten_qmcpy_imports.py
+++ /dev/null
@@ -1,330 +0,0 @@
-import json
-from pathlib import Path
-
-from scripts.flatten_qmcpy_imports import (
-    _load_qmcpy_public_names,
-    flatten_imports,
-    main,
-)
-
-
-def _nested_import(module, imported):
-    return f"from {'qmcpy.' + module} import {imported}"
-
-
-def test_flatten_imports_basic():
-    source = (
-        _nested_import("integrand", "Keister")
-        + "\n"
-        + _nested_import("discrete_distribution.lattice", "Lattice as LD")
-        + "\nfrom qmcpy import DigitalNetB2\nimport qmcpy.util\n"
-    ).encode()
-
-    updated, count = flatten_imports(
-        source, frozenset({"DigitalNetB2", "Keister", "Lattice"})
-    )
-
-    assert count == 3
-    assert updated == (
-        b"from qmcpy import DigitalNetB2, Keister, Lattice as LD\n"
-        b"import qmcpy.util\n"
-    )
-
-
-def test_flatten_preserves_private():
-    source = (
-        _nested_import("_internal._helpers", "PublicHelper")
-        + "\n"
-        + _nested_import(
-            "true_measure.uniform_triangle",
-            "UniformTriangle, _UniformTriangleAdapter",
-        )
-        + "\n"
-        + _nested_import(
-            "true_measure.copula",
-            "(\n    AbstractCopula,\n    _validate_dimension,\n)",
-        )
-        + "\n"
-        + _nested_import("integrand", "Keister")
-        + "\n"
-    ).encode()
-
-    updated, count = flatten_imports(source, frozenset({"Keister"}))
-
-    assert count == 1
-    assert updated == source.replace(
-        _nested_import("integrand", "Keister").encode(),
-        b"from qmcpy import Keister",
-    )
-
-
-def test_private_module_splits_groups():
-    source = (
-        b"from qmcpy import Zeta\n"
-        b"from qmcpy._internal._helpers import PublicHelper\n"
-        b"from qmcpy import Alpha\n"
-    )
-
-    updated, count = flatten_imports(source)
-
-    assert count == 0
-    assert updated == source
-
-
-def test_flatten_preserves_util_imports():
-    source = (
-        b"from qmcpy.util import ParameterError\n"
-        b"from qmcpy.util.transforms import tf_exp\n"
-    )
-
-    updated, count = flatten_imports(source, frozenset({"ParameterError", "tf_exp"}))
-
-    assert (updated, count) == (source, 0)
-
-
-def test_flatten_keeps_nonpublic_names():
-    source = b"from qmcpy.stopping_criterion.pf_gp_ci import PFGPCIData\n"
-
-    updated, count = flatten_imports(source, frozenset({"PFGPCI"}))
-
-    assert (updated, count) == (source, 0)
-
-
-def test_flatten_no_public_api_noop():
-    source = (_nested_import("integrand", "Keister") + "\n").encode()
-
-    updated, count = flatten_imports(source)
-
-    assert (updated, count) == (source, 0)
-
-
-def test_flatten_preserve_str_literals():
-    source = b'text = """\nfrom qmcpy.integrand import Keister\n"""\n'
-
-    updated, count = flatten_imports(source, frozenset({"Keister"}))
-
-    assert count == 0
-    assert updated == source
-
-
-def test_python_string_protection_applies_to_every_rewrite_stage():
-    string_body = (
-        b'text = """\n'
-        b"from qmcpy.integrand import Keister\n"
-        b"from qmcpy import Zeta,Beta\n"
-        b"from qmcpy import Alpha\n"
-        b"from qmcpy import *\n"
-        b"from qmcpy import *\n"
-        b'"""\n'
-    )
-    source = string_body + b"from qmcpy.integrand import Keister\n"
-
-    updated, count = flatten_imports(source, frozenset({"Keister"}))
-
-    assert count == 1
-    assert updated == string_body + b"from qmcpy import Keister\n"
-
-
-def test_python_tokenize_failure_is_fail_closed():
-    source = b'"""unterminated\nfrom qmcpy.integrand import Keister\n'
-
-    assert flatten_imports(source, frozenset({"Keister"})) == (source, 0)
-
-
-def test_flatten_skip_star_expansion():
-    source = (
-        b"from qmcpy import *\n\n"
-        b"def f(Lattice):\n"
-        b"    return Lattice\n\n"
-        b"y = Keister(dimension=2)\n"
-        b"x = Lattice(dimension=2)\n"
-    )
-
-    updated, count = flatten_imports(source, frozenset({"Keister", "Lattice"}))
-
-    assert count == 0
-    assert updated == source
-
-
-def test_notebook_star_dedup():
-    notebook = {
-        "cells": [
-            {
-                "cell_type": "code",
-                "source": [
-                    _nested_import("integrand", "*") + "\n",
-                    _nested_import("true_measure", "*"),
-                ],
-            }
-        ]
-    }
-    source = json.dumps(notebook, indent=1).encode()
-
-    updated, count = flatten_imports(source, frozenset({"Keister"}))
-
-    assert count == 3
-    assert json.loads(updated)["cells"][0]["source"] == ["from qmcpy import *"]
-
-
-def test_named_imports_merge_sort():
-    source = (
-        b"from qmcpy import Zeta,Beta\n"
-        b"from qmcpy import Alpha\n"
-        b"\n"
-        b"from qmcpy import Gamma\n"
-    )
-
-    updated, count = flatten_imports(source)
-
-    assert count == 1
-    assert updated == (
-        b"from qmcpy import Alpha, Beta, Zeta\n"
-        b"\n"
-        b"from qmcpy import Gamma\n"
-    )
-    assert flatten_imports(updated) == (updated, 0)
-
-
-def test_merge_paren_and_single_line():
-    source = b"""from qmcpy import (
-    KernelDigShiftInvar,
-    KernelDigShiftInvarAdaptiveAlpha,
-    KernelDigShiftInvarCombined,
-    KernelShiftInvar,
-    KernelShiftInvarCombined,
-)
-from qmcpy import tf_exp_eps, tf_exp_eps_inv
-"""
-
-    updated, count = flatten_imports(source)
-
-    assert count == 1
-    assert updated == b"""from qmcpy import (
-    KernelDigShiftInvar,
-    KernelDigShiftInvarAdaptiveAlpha,
-    KernelDigShiftInvarCombined,
-    KernelShiftInvar,
-    KernelShiftInvarCombined,
-    tf_exp_eps,
-    tf_exp_eps_inv,
-)
-"""
-    assert flatten_imports(updated) == (updated, 0)
-
-
-def test_merge_same_scope_only():
-    source = (
-        b"if enabled:\n"
-        b"    from qmcpy import Zeta\n"
-        b"    from qmcpy import Alpha as First\n"
-        b"else:\n"
-        b"    from qmcpy import Beta\n"
-        b"from qmcpy import _Private\n"
-        b"from qmcpy import Gamma  # keep this comment\n"
-    )
-
-    updated, count = flatten_imports(source)
-
-    assert count == 1
-    assert updated == (
-        b"if enabled:\n"
-        b"    from qmcpy import Alpha as First, Zeta\n"
-        b"else:\n"
-        b"    from qmcpy import Beta\n"
-        b"from qmcpy import _Private\n"
-        b"from qmcpy import Gamma  # keep this comment\n"
-    )
-
-
-def test_notebook_named_merge():
-    notebook = {
-        "cells": [
-            {
-                "cell_type": "code",
-                "source": [
-                    "from qmcpy import Zeta\n",
-                    "from qmcpy import Alpha,Beta\n",
-                    "print(Alpha)\n",
-                ],
-            }
-        ]
-    }
-    source = json.dumps(notebook, indent=1).encode()
-
-    updated, count = flatten_imports(source)
-
-    assert count == 1
-    assert json.loads(updated)["cells"][0]["source"] == [
-        "from qmcpy import Alpha, Beta, Zeta\n",
-        "print(Alpha)\n",
-    ]
-    assert flatten_imports(updated) == (updated, 0)
-
-
-def test_notebook_flattens_nested_imports_only_in_code_cells():
-    nested_import = _nested_import("integrand", "Keister") + "\n"
-    metadata_import = _nested_import("true_measure", "Gaussian") + "\n"
-    string_literal = f'text = "{nested_import.rstrip()}"\n'
-    multiline_string = ['text = """\n', nested_import, '"""\n']
-    notebook = {
-        "metadata": {"source": [metadata_import]},
-        "cells": [
-            {"cell_type": "markdown", "source": [nested_import]},
-            {"cell_type": "code", "source": [nested_import]},
-            {"cell_type": "code", "source": [string_literal]},
-            {"cell_type": "code", "source": multiline_string},
-        ]
-    }
-    source = json.dumps(notebook, indent=1).encode()
-
-    updated, count = flatten_imports(source, frozenset({"Keister"}))
-
-    cells = json.loads(updated)["cells"]
-    assert count == 1
-    assert json.loads(updated)["metadata"]["source"] == [metadata_import]
-    assert cells[0]["source"] == [nested_import]
-    assert cells[1]["source"] == ["from qmcpy import Keister\n"]
-    assert cells[2]["source"] == [string_literal]
-    assert cells[3]["source"] == multiline_string
-    assert flatten_imports(updated, frozenset({"Keister"})) == (updated, 0)
-
-
-def test_markdown_import_examples_are_flattened(tmp_path):
-    path = tmp_path / "example.md"
-    path.write_bytes(
-        b'Example with unmatched prose delimiter: """\n\n'
-        b"```python\n"
-        b"from qmcpy.integrand import Keister\n"
-        b"```\n"
-    )
-
-    assert main([str(path)]) == 0
-    assert b"from qmcpy import Keister" in path.read_bytes()
-
-
-def test_check_mode_no_write(tmp_path):
-    path = tmp_path / "example.py"
-    original = (_nested_import("true_measure", "Gaussian") + "\n").encode()
-    path.write_bytes(original)
-
-    assert main(["--check", str(path)]) == 1
-    assert path.read_bytes() == original
-
-    assert main([str(path)]) == 0
-    assert path.read_bytes() == b"from qmcpy import Gaussian\n"
-    assert main(["--check", str(path)]) == 0
-
-
-def test_public_names_optional_free_stable():
-    repository_root = Path(__file__).resolve().parent.parent
-    names = _load_qmcpy_public_names(repository_root)
-
-    assert names is not None
-    assert "Gaussian" in names
-    assert "Keister" in names
-    # Optional dependencies are blocked in the probe context, so fallback
-    # exports are part of the deterministic name set.
-    assert "PFGPCI" in names
-    # Helpers that are deliberately not part of the top-level API.
-    assert "PFGPCIData" not in names
-    assert "TriangularDistribution" not in names
\ No newline at end of file
diff --git a/test/test_ft_fast_transform_fallbacks.py b/test/test_ft_fast_transform_fallbacks.py
new file mode 100644
index 000000000..d60266030
--- /dev/null
+++ b/test/test_ft_fast_transform_fallbacks.py
@@ -0,0 +1,52 @@
+import unittest
+
+import numpy as np
+
+from qmcpy import (
+    fftbr,
+    fftbr_torch,
+    fwht,
+    fwht_torch,
+    ifftbr,
+    ifftbr_torch,
+    omega_fftbr,
+    omega_fftbr_torch,
+    omega_fwht,
+    omega_fwht_torch,
+)
+
+
+class TestFastTransformFallbacks(unittest.TestCase):
+
+    def test_non_torch_transforms_basic(self):
+        rng = np.random.default_rng(11)
+        x = rng.random(8) + 1j * rng.random(8)
+        y = fftbr(x)
+        self.assertEqual(y.shape, x.shape)
+        xr = ifftbr(y)
+        self.assertEqual(xr.shape, x.shape)
+
+        a = rng.random(8)
+        b = fwht(a)
+        self.assertEqual(b.shape, a.shape)
+
+        omega = omega_fftbr(3)
+        self.assertEqual(omega.shape[0], 2**3)
+        omega2 = omega_fwht(3)
+        self.assertEqual(omega2.shape[0], 2**3)
+
+    def test_torch_fallbacks_raise(self):
+        with self.assertRaises(Exception):
+            fftbr_torch()
+        with self.assertRaises(Exception):
+            ifftbr_torch()
+        with self.assertRaises(Exception):
+            fwht_torch()
+        with self.assertRaises(Exception):
+            omega_fftbr_torch()
+        with self.assertRaises(Exception):
+            omega_fwht_torch()
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/test_ig_financial_option_quick.py b/test/test_ig_financial_option_quick.py
new file mode 100644
index 000000000..921701712
--- /dev/null
+++ b/test/test_ig_financial_option_quick.py
@@ -0,0 +1,102 @@
+import unittest
+
+import numpy as np
+
+from qmcpy import FinancialOption
+import qmcpy
+
+
+class SmallSampler(qmcpy.AbstractDiscreteDistribution):
+    def __init__(self, d=3):
+        super().__init__(
+            dimension=d, replications=1, seed=123, d_limit=100, n_limit=1024
+        )
+
+    def _gen_samples(self, n_min, n_max, return_binary=False, warn=True):
+        n = n_max - n_min
+        # return shape (replications, n, d)
+        arr = np.tile(np.linspace(0.1, 1.0, n)[:, None], (1, self.d))
+        return arr.reshape(self.replications, n, self.d)
+
+
+class TestFinancialOptionPayoffs(unittest.TestCase):
+
+    def test_financial_option_payoffs_and_exact(self):
+        sampler = SmallSampler(d=3)
+        fo = FinancialOption(
+            sampler,
+            option="EUROPEAN",
+            call_put="CALL",
+            volatility=0.5,
+            start_price=30,
+            strike_price=25,
+            interest_rate=0.01,
+            t_final=1,
+        )
+        gbm = np.array([[30.0, 28.0, 35.0]])
+        c = fo.payoff_european_call(gbm)
+        p = fo.payoff_european_put(gbm)
+        self.assertEqual(c.shape, (1,))
+        self.assertEqual(p.shape, (1,))
+
+        # Asian arithmetic trapezoidal
+        fo_asian = FinancialOption(
+            sampler,
+            option="ASIAN",
+            asian_mean="ARITHMETIC",
+            asian_mean_quadrature_rule="TRAPEZOIDAL",
+        )
+        gbm2 = np.array([[30.0, 32.0, 34.0]])
+        a_call = fo_asian.payoff_asian_arithmetic_trap_call(gbm2)
+        self.assertEqual(a_call.shape, (1,))
+
+        # geometric right call
+        fo_geo = FinancialOption(
+            sampler,
+            option="ASIAN",
+            asian_mean="GEOMETRIC",
+            asian_mean_quadrature_rule="RIGHT",
+        )
+        g_call = fo_geo.payoff_asian_geometric_right_call(np.array([[30.0, 30.0, 30.0]]))
+        self.assertEqual(g_call.shape, (1,))
+
+        # barrier options: up and down behaviors
+        fo_barrier_up = FinancialOption(
+            sampler, option="BARRIER", barrier_in_out="IN", barrier_price=25, start_price=20
+        )
+        gbm_up = np.array([[20.0, 26.0, 27.0]])
+        v = fo_barrier_up.payoff_barrier_in_up_call(gbm_up)
+        self.assertEqual(v.shape, (1,))
+
+        fo_barrier_out = FinancialOption(
+            sampler,
+            option="BARRIER",
+            barrier_in_out="OUT",
+            barrier_price=40,
+            start_price=30,
+        )
+        gbm_out = np.array([[30.0, 32.0, 33.0]])
+        v2 = fo_barrier_out.payoff_barrier_out_up_call(gbm_out)
+        self.assertEqual(v2.shape, (1,))
+
+        # lookback
+        fo_lb = FinancialOption(sampler, option="LOOKBACK")
+        lb = fo_lb.payoff_lookback_call(np.array([[10.0, 9.0, 12.0]]))
+        self.assertEqual(lb.shape, (1,))
+
+        # digital
+        fo_dig = FinancialOption(sampler, option="DIGITAL", digital_payout=5)
+        dig = fo_dig.payoff_digital_call(np.array([[10.0, 11.0, 12.0]]))
+        self.assertEqual(dig.shape, (1,))
+
+        # exact value for European should return a float
+        val = fo.get_exact_value()
+        self.assertTrue(np.isscalar(val))
+
+        # exact value for Asian geometric right
+        val2 = fo_geo.get_exact_value()
+        self.assertTrue(np.isscalar(val2))
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/test_integrands.py b/test/test_ig_integrands.py
similarity index 100%
rename from test/test_integrands.py
rename to test/test_ig_integrands.py
diff --git a/test/test_option.py b/test/test_ig_option.py
similarity index 100%
rename from test/test_option.py
rename to test/test_ig_option.py
diff --git a/test/test_option_ml.py b/test/test_ig_option_ml.py
similarity index 100%
rename from test/test_option_ml.py
rename to test/test_ig_option_ml.py
diff --git a/test/test_install_mpmc_pyg.py b/test/test_install_mpmc_pyg.py
deleted file mode 100644
index 41b1b4e0b..000000000
--- a/test/test_install_mpmc_pyg.py
+++ /dev/null
@@ -1,85 +0,0 @@
-"""Tests for the platform-specific MPMC dependency installer."""
-
-import subprocess
-from types import SimpleNamespace
-
-import pytest
-
-from qmcpy.util import install_mpmc_pyg
-
-
-def _torch(version="2.12.1+cpu", cuda=None, hip=None):
-    return SimpleNamespace(
-        __version__=version,
-        version=SimpleNamespace(cuda=cuda, hip=hip),
-    )
-
-
-def test_torch_versions_include_baseline_fallback():
-    """Wheel lookup tries an exact patch release, then its minor baseline."""
-    assert install_mpmc_pyg.torch_versions("2.12.1+cpu") == ["2.12.1", "2.12.0"]
-    assert install_mpmc_pyg.torch_versions("2.12.0") == ["2.12.0"]
-
-    with pytest.raises(RuntimeError, match="Unable to parse torch version"):
-        install_mpmc_pyg.torch_versions("development")
-
-
-@pytest.mark.parametrize(
-    ("torch_module", "expected"),
-    [
-        (_torch(), "cpu"),
-        (_torch(cuda="12.6"), "cu126"),
-        (_torch(cuda="13.0.1"), "cu130"),
-    ],
-)
-def test_accelerator_tag(torch_module, expected):
-    """PyTorch build metadata maps to the expected PyG wheel tag."""
-    assert install_mpmc_pyg.accelerator_tag(torch_module) == expected
-
-
-def test_accelerator_tag_rejects_rocm():
-    """The installer directs unsupported ROCm users to upstream guidance."""
-    with pytest.raises(RuntimeError, match="does not currently support ROCm"):
-        install_mpmc_pyg.accelerator_tag(_torch(hip="6.3"))
-
-
-def test_main_retries_with_torch_minor_baseline(monkeypatch):
-    """A missing exact wheel page falls back to the minor baseline page."""
-    calls = []
-
-    def fake_run(*args):
-        calls.append(args)
-        if args[-1].endswith("torch-2.12.1+cpu.html"):
-            raise subprocess.CalledProcessError(1, args)
-
-    monkeypatch.setattr(install_mpmc_pyg, "run", fake_run)
-
-    install_mpmc_pyg.main(_torch())
-
-    assert calls[0][-1] == "torch-geometric>=2.6.1"
-    assert calls[1][-1] == "https://data.pyg.org/whl/torch-2.12.1+cpu.html"
-    assert calls[2][-1] == "https://data.pyg.org/whl/torch-2.12.0+cpu.html"
-    assert "--only-binary" in calls[1]
-
-
-def test_main_explains_that_torch_must_be_installed(monkeypatch):
-    """Running the helper before installing the extra gives a useful error."""
-    def missing_torch(_name):
-        raise ModuleNotFoundError("No module named 'torch'", name="torch")
-
-    monkeypatch.setattr(install_mpmc_pyg.importlib, "import_module", missing_torch)
-
-    with pytest.raises(RuntimeError, match=r"install 'qmcpy\[mpmc\]'"):
-        install_mpmc_pyg.main()
-
-
-def test_main_reports_missing_wheel(monkeypatch):
-    """Exhausting candidate wheel pages reports the build that failed."""
-    def fail_pyg_lib(*args):
-        if "pyg_lib>=0.6.0" in args:
-            raise subprocess.CalledProcessError(1, args)
-
-    monkeypatch.setattr(install_mpmc_pyg, "run", fail_pyg_lib)
-
-    with pytest.raises(RuntimeError, match=r"torch 2\.12\.1\+cpu \(cpu\)"):
-        install_mpmc_pyg.main(_torch())
diff --git a/test/test_kernels.py b/test/test_kn_kernels.py
similarity index 100%
rename from test/test_kernels.py
rename to test/test_kn_kernels.py
diff --git a/test/test_mpmc_optional_imports.py b/test/test_mpmc_optional_imports.py
deleted file mode 100644
index 33b71b841..000000000
--- a/test/test_mpmc_optional_imports.py
+++ /dev/null
@@ -1,100 +0,0 @@
-import ast
-import builtins
-from pathlib import Path
-
-import pytest
-
-
-def _execute_optional_import(blocked_import):
-    repository_root = Path(__file__).resolve().parent.parent
-    init_path = repository_root / "qmcpy" / "__init__.py"
-    init_tree = ast.parse(init_path.read_text())
-    optional_import = next(
-        node
-        for node in init_tree.body
-        if isinstance(node, ast.Try)
-        and any(
-            isinstance(statement, ast.ImportFrom)
-            and statement.module == "discrete_distribution.mpmc"
-            for statement in node.body
-        )
-    )
-
-    import qmcpy
-
-    real_import = builtins.__import__
-
-    def guarded_import(name, globals=None, locals=None, fromlist=(), level=0):
-        missing_module = blocked_import(name, fromlist, level)
-        if missing_module is not None:
-            raise ModuleNotFoundError(
-                "blocked optional dependency",
-                name=missing_module,
-            )
-        return real_import(name, globals, locals, fromlist, level)
-
-    test_builtins = vars(builtins).copy()
-    test_builtins["__import__"] = guarded_import
-    namespace = {"__builtins__": test_builtins, "__package__": "qmcpy"}
-    module = ast.Module(body=[optional_import], type_ignores=[])
-    exec(compile(module, str(init_path), "exec"), namespace)
-    return namespace
-
-
-def test_mpmc_utils_remain_available_without_pyg():
-    pytest.importorskip("torch")
-
-    def block_pyg_models(name, fromlist, level):
-        if level == 1 and name == "discrete_distribution.mpmc.models":
-            return "torch_geometric"
-        return None
-
-    namespace = _execute_optional_import(block_pyg_models)
-
-    import qmcpy
-
-    assert namespace["mpmc_utils"] is qmcpy.mpmc_utils
-    assert namespace["mpmc_utils"].__name__ == (
-        "qmcpy.discrete_distribution.mpmc.utils"
-    )
-    assert "utils" not in namespace
-
-    with pytest.raises(ModuleNotFoundError, match="MPMC_net.*torch_geometric") as error:
-        namespace["MPMC_net"]()
-    assert error.value.name == "torch_geometric"
-
-
-def test_mpmc_placeholders_report_missing_torch():
-    def block_torch_utils(name, fromlist, level):
-        if (
-            level == 1
-            and name == "discrete_distribution.mpmc"
-            and "utils" in fromlist
-        ):
-            return "torch"
-        return None
-
-    namespace = _execute_optional_import(block_torch_utils)
-
-    with pytest.raises(ModuleNotFoundError, match="mpmc_utils.*torch") as error:
-        namespace["mpmc_utils"].L2star
-    assert error.value.name == "torch"
-
-    with pytest.raises(ModuleNotFoundError, match="MPMC_net.*torch") as error:
-        namespace["MPMC_net"]()
-    assert error.value.name == "torch"
-
-
-def test_mpmc_placeholder_missing_torch_scatter():
-    pytest.importorskip("torch")
-
-    def block_torch_scatter(name, fromlist, level):
-        if level == 1 and name == "discrete_distribution.mpmc.models":
-            return "torch_scatter"
-        return None
-
-    namespace = _execute_optional_import(block_torch_scatter)
-
-    with pytest.raises(ModuleNotFoundError, match="MPMC_net.*torch_scatter") as error:
-        namespace["MPMC_net"]()
-    assert error.value.name == "torch_scatter"
diff --git a/test/test_plot_and_stop.py b/test/test_plot_and_stop.py
deleted file mode 100644
index dd564cd9d..000000000
--- a/test/test_plot_and_stop.py
+++ /dev/null
@@ -1,155 +0,0 @@
-import sys
-import types
-import numpy as np
-import builtins
-import pytest
-
-import qmcpy
-from qmcpy import plot_proj
-from qmcpy.util import stop_notebook
-
-
-class FakeAxes:
-    def __init__(self):
-        self.removed = False
-        self.calls = []
-
-    def remove(self):
-        self.removed = True
-
-    def set_xlim(self, *a, **k):
-        self.calls.append(("set_xlim", a))
-
-    def set_ylim(self, *a, **k):
-        self.calls.append(("set_ylim", a))
-
-    def set_xticks(self, *a, **k):
-        self.calls.append(("set_xticks", a))
-
-    def set_yticks(self, *a, **k):
-        self.calls.append(("set_yticks", a))
-
-    def set_aspect(self, *a, **k):
-        self.calls.append(("set_aspect", a))
-
-    def grid(self, *a, **k):
-        self.calls.append(("grid", a))
-
-    def tick_params(self, *a, **k):
-        self.calls.append(("tick_params", a))
-
-    def set_xlabel(self, *a, **k):
-        self.calls.append(("set_xlabel", a))
-
-    def set_ylabel(self, *a, **k):
-        self.calls.append(("set_ylabel", a))
-
-    def scatter(self, *a, **k):
-        self.calls.append(("scatter", a))
-
-
-class FakeFig:
-    def __init__(self):
-        self.tl = False
-
-    def tight_layout(self, *a, **k):
-        self.tl = True
-
-
-def make_fake_matplotlib(nrows, ncols):
-    plt = types.ModuleType("matplotlib.pyplot")
-    plt.style = types.SimpleNamespace()
-    plt.style.use = lambda *a, **k: None
-    plt.rcParams = {
-        "font.family": "sans-serif",
-        "axes.prop_cycle": types.SimpleNamespace(
-            by_key=lambda: {"color": ["k", "b", "r"]}
-        ),
-    }
-
-    def subplots(nrows=1, ncols=1, figsize=None, squeeze=False):
-        fig = FakeFig()
-        ax = np.empty((nrows, ncols), dtype=object)
-        for i in range(nrows):
-            for j in range(ncols):
-                ax[i, j] = FakeAxes()
-        return fig, ax
-
-    plt.subplots = subplots
-    plt.suptitle = lambda *a, **k: None
-    return plt
-
-
-class DummySampler(qmcpy.AbstractDiscreteDistribution):
-    def __init__(self, d=2):
-        super().__init__(dimension=d, replications=1, seed=1, d_limit=10, n_limit=100)
-
-    def _gen_samples(self, n_min, n_max, return_binary=False, warn=True):
-        n = n_max - n_min
-        return np.tile(np.arange(n)[:, None] / max(1, n - 1), (1, 1, self.d)).reshape(
-            self.replications, n, self.d
-        )
-
-    def __repr__(self):
-        return "DummySampler"
-
-
-def test_plot_proj_with_fake_matplotlib_and_sampler(monkeypatch):
-    # Inject fake matplotlib.pyplot
-    fake_plt = make_fake_matplotlib(1, 1)
-    # Create a proper matplotlib package module with colors submodule
-    fake_matplotlib = types.ModuleType("matplotlib")
-    fake_matplotlib.pyplot = fake_plt
-    fake_matplotlib.colors = types.SimpleNamespace()
-    monkeypatch.setitem(sys.modules, "matplotlib.pyplot", fake_plt)
-    monkeypatch.setitem(sys.modules, "matplotlib", fake_matplotlib)
-
-    sampler = DummySampler(d=3)
-    fig, ax = plot_proj(
-        sampler,
-        n=4,
-        d_horizontal=1,
-        d_vertical=2,
-        math_ind=True,
-        marker_size=1,
-        figfac=1,
-    )
-    assert isinstance(fig, FakeFig)
-    assert isinstance(ax, np.ndarray)
-    # At least one axes should have scatter calls or be removed
-    found = False
-    for a in ax.flatten():
-        if getattr(a, "removed", False) or any(c[0] == "scatter" for c in a.calls):
-            found = True
-            break
-    assert found
-
-
-def test_plot_proj_with_callable_sampler(monkeypatch):
-    # sampler not instance of AbstractDiscreteDistribution -> uses t_i labels
-    fake_plt = make_fake_matplotlib(1, 1)
-    fake_matplotlib = types.ModuleType("matplotlib")
-    fake_matplotlib.pyplot = fake_plt
-    fake_matplotlib.colors = types.SimpleNamespace()
-    monkeypatch.setitem(sys.modules, "matplotlib.pyplot", fake_plt)
-    monkeypatch.setitem(sys.modules, "matplotlib", fake_matplotlib)
-
-    def sampler_callable(n):
-        return np.zeros((n, 1))
-
-    fig, ax = plot_proj(
-        sampler_callable, n=3, d_horizontal=0, d_vertical=0, math_ind=False
-    )
-    assert isinstance(fig, FakeFig)
-
-
-def test_stop_notebook_yes_and_no(monkeypatch):
-    # When input is 'yes' nothing should happen
-    monkeypatch.setattr(builtins, "input", lambda prompt="": "yes")
-    # Should not raise
-    stop_notebook("prompt")
-
-    # When input is not 'yes' should exit
-    monkeypatch.setattr(builtins, "input", lambda prompt="": "no")
-    with pytest.raises(SystemExit):
-        stop_notebook("prompt")
diff --git a/test/test_product_measure.py b/test/test_product_measure.py
deleted file mode 100644
index c91b313bd..000000000
--- a/test/test_product_measure.py
+++ /dev/null
@@ -1,277 +0,0 @@
-import numpy as np
-import pytest
-import scipy.stats as stats
-
-from qmcpy import (
-    AcceptanceRejection,
-    DigitalNetB2,
-    DummySampler,
-    Gaussian,
-    GaussianCopula,
-    ProductMeasure,
-    SciPyWrapper,
-    Uniform,
-    ZeroInflatedExpUniform,
-)
-from qmcpy.util import DimensionError, ParameterError
-
-
-def test_product_measure_zero_inflated_with_scipy_uniform_shape():
-    n = 32
-    marginals = [
-        ZeroInflatedExpUniform(DummySampler(1), p_zero=0.4, lam=1.5),
-        SciPyWrapper(DummySampler(1), stats.uniform(loc=2.0, scale=3.0)),
-    ]
-    tm = ProductMeasure(sampler=DigitalNetB2(2, seed=23), marginals=marginals)
-
-    x = tm(n)
-
-    assert x.shape == (n, 2)
-    assert np.any(x[:, 0] == 0.0)
-    assert np.all((2.0 <= x[:, 1]) & (x[:, 1] <= 5.0))
-
-
-def test_product_measure_replication_shape():
-    n = 16
-    r = 3
-    marginals = [
-        ZeroInflatedExpUniform(DummySampler(1), p_zero=0.4, lam=1.5),
-        Uniform(DummySampler(1), lower_bound=2.0, upper_bound=5.0),
-    ]
-    tm = ProductMeasure(
-        sampler=DigitalNetB2(2, seed=23, replications=r),
-        marginals=marginals,
-    )
-
-    x = tm(n)
-
-    assert x.shape == (r, n, 2)
-
-
-def test_product_measure_marginals_with_different_dimensions():
-    n = 32
-    marginals = [
-        Gaussian(
-            DummySampler(2),
-            mean=[1.0, -1.0],
-            covariance=[[2.0, 0.25], [0.25, 1.0]],
-        ),
-        ZeroInflatedExpUniform(DummySampler(1), p_zero=0.4, lam=1.5),
-    ]
-    tm = ProductMeasure(sampler=DigitalNetB2(3, seed=31), marginals=marginals)
-
-    x = tm(n)
-
-    assert tm.d == 3
-    assert np.array_equal(tm.marginal_dimensions, np.array([2, 1]))
-    assert x.shape == (n, 3)
-    assert np.all(np.isfinite(x[:, :2]))
-    assert np.all(x[:, 2] >= 0.0)
-
-
-def test_product_measure_block_split_range_and_weight_product():
-    n = 16
-    marginals = [
-        Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0),
-        Uniform(
-            DummySampler(2),
-            lower_bound=[20.0, 30.0],
-            upper_bound=[24.0, 36.0],
-        ),
-    ]
-    tm = ProductMeasure(sampler=DigitalNetB2(3, seed=41), marginals=marginals)
-
-    u = tm.discrete_distrib.gen_samples(n)
-    x = tm._transform(u)
-    x_call, jac = tm(n, return_weights=True)
-    expected = np.concatenate(
-        [
-            marginals[0]._jacobian_transform_r(u[..., :1], return_weights=False),
-            marginals[1]._jacobian_transform_r(u[..., 1:], return_weights=False),
-        ],
-        axis=-1,
-    )
-    expected_range = np.array([[10.0, 12.0], [20.0, 24.0], [30.0, 36.0]])
-
-    assert x.shape == (n, 3)
-    assert np.allclose(tm.range, expected_range)
-    assert np.allclose(x, expected)
-    assert np.all((10.0 <= x[:, 0]) & (x[:, 0] <= 12.0))
-    assert np.all((20.0 <= x[:, 1]) & (x[:, 1] <= 24.0))
-    assert np.all((30.0 <= x[:, 2]) & (x[:, 2] <= 36.0))
-    assert np.allclose(tm._weight(x), 1.0 / (2.0 * 4.0 * 6.0))
-    assert x_call.shape == (n, 3)
-    assert np.allclose(jac, 2.0 * 4.0 * 6.0)
-
-
-def test_product_measure_invalid_inputs():
-    with pytest.raises(ParameterError, match="nonempty list of marginals"):
-        ProductMeasure(sampler=DigitalNetB2(1, seed=7), marginals=[])
-
-    with pytest.raises(ParameterError, match="marginal"):
-        ProductMeasure(sampler=DigitalNetB2(1, seed=7), marginals=[object()])
-
-    with pytest.raises(ParameterError, match="AbstractDiscreteDistribution"):
-        ProductMeasure(sampler=object(), marginals=[Uniform(DummySampler(1))])
-
-    marginals = [Uniform(DummySampler(1))]
-    with pytest.raises(DimensionError, match="sum of marginal dimensions"):
-        ProductMeasure(sampler=DigitalNetB2(2, seed=7), marginals=marginals)
-
-
-def test_product_measure_rejects_non_dimension_preserving_marginal():
-    marginal = AcceptanceRejection(
-        DigitalNetB2(2, seed=7),
-        lambda x: np.ones(len(x)),
-        1.0,
-        1.0,
-    )
-
-    with pytest.raises(DimensionError, match="dimension-preserving"):
-        ProductMeasure(DigitalNetB2(2, seed=11), [marginal])
-
-
-def test_product_measure_spawn_preserves_marginal_blocks_and_replaces_outer_sampler():
-    marginals = [
-        Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0),
-        Uniform(
-            DummySampler(2),
-            lower_bound=[20.0, 30.0],
-            upper_bound=[24.0, 36.0],
-        ),
-    ]
-    tm = ProductMeasure(sampler=DigitalNetB2(3, seed=41), marginals=marginals)
-
-    spawn = tm.spawn(s=1)[0]
-
-    assert isinstance(spawn, ProductMeasure)
-    assert spawn.d == 3
-    assert spawn.marginals == tm.marginals
-    assert spawn.discrete_distrib is not tm.discrete_distrib
-    assert np.array_equal(spawn.marginal_dimensions, np.array([1, 2]))
-
-    with pytest.raises(DimensionError):
-        tm.spawn(s=1, dimensions=4)
-
-
-def test_product_measure_does_not_use_marginal_dummy_sampler_values():
-    marginals = [
-        Uniform(DummySampler(1), lower_bound=0.0, upper_bound=2.0),
-        Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0),
-    ]
-
-    with pytest.raises(ParameterError, match="construction placeholder"):
-        marginals[0].discrete_distrib(4)
-
-    tm = ProductMeasure(sampler=DigitalNetB2(2, seed=19), marginals=marginals)
-    x = tm(8)
-
-    assert x.shape == (8, 2)
-    assert np.all((0.0 <= x[:, 0]) & (x[:, 0] <= 2.0))
-    assert np.all((10.0 <= x[:, 1]) & (x[:, 1] <= 12.0))
-
-
-def test_product_measure_same_outer_seed_matches_different_outer_seed_changes():
-    marginals = [
-        Uniform(DummySampler(1), lower_bound=0.0, upper_bound=2.0),
-        Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0),
-    ]
-
-    first = ProductMeasure(sampler=DigitalNetB2(2, seed=101), marginals=marginals)(16)
-    same_outer = ProductMeasure(sampler=DigitalNetB2(2, seed=101), marginals=marginals)(16)
-    different_outer = ProductMeasure(sampler=DigitalNetB2(2, seed=102), marginals=marginals)(16)
-
-    assert np.array_equal(first, same_outer)
-    assert not np.array_equal(first, different_outer)
-
-
-def test_product_measure_replication_means_close_to_uniform_targets():
-    n = 1024
-    r = 4
-    marginals = [
-        Uniform(DummySampler(1), lower_bound=0.0, upper_bound=2.0),
-        Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0),
-    ]
-    tm = ProductMeasure(
-        sampler=DigitalNetB2(2, seed=101, replications=r),
-        marginals=marginals,
-    )
-
-    x = tm(n)
-    replication_means = x.mean(axis=1)
-
-    assert x.shape == (r, n, 2)
-    assert np.allclose(replication_means[:, 0], 1.0, atol=0.03)
-    assert np.allclose(replication_means[:, 1], 11.0, atol=0.03)
-
-
-def test_product_measure_with_scipywrapper_beta_marginal():
-    n = 64
-    marginals = [
-        Uniform(DummySampler(1), lower_bound=-1.0, upper_bound=1.0),
-        SciPyWrapper(DummySampler(1), stats.beta(a=2.0, b=5.0)),
-    ]
-    tm = ProductMeasure(sampler=DigitalNetB2(2, seed=71), marginals=marginals)
-
-    x = tm(n)
-
-    assert x.shape == (n, 2)
-    assert np.all((-1.0 <= x[:, 0]) & (x[:, 0] <= 1.0))
-    assert np.all((0.0 <= x[:, 1]) & (x[:, 1] <= 1.0))
-
-
-def test_product_measure_matches_equivalent_scipywrapper():
-    n = 128
-    seed = 55
-    scipy_marginals = [stats.norm(loc=0.0, scale=1.0), stats.gamma(a=2.0, scale=1.0)]
-    product_marginals = [
-        SciPyWrapper(DummySampler(1), scipy_marginals[0]),
-        SciPyWrapper(DummySampler(1), scipy_marginals[1]),
-    ]
-
-    product_samples = ProductMeasure(
-        sampler=DigitalNetB2(2, seed=seed),
-        marginals=product_marginals,
-    )(n)
-    scipy_samples = SciPyWrapper(DigitalNetB2(2, seed=seed), scipy_marginals)(n)
-
-    assert np.array_equal(product_samples, scipy_samples)
-
-
-def test_product_measure_with_gaussian_copula_marginal():
-    n = 64
-    copula = GaussianCopula(
-        DummySampler(2),
-        marginals=[stats.beta(a=2.0, b=5.0), stats.gamma(a=3.0, scale=1.0)],
-        correlation=[[1.0, 0.5], [0.5, 1.0]],
-    )
-    marginals = [copula, Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0)]
-    tm = ProductMeasure(sampler=DigitalNetB2(3, seed=81), marginals=marginals)
-
-    x = tm(n)
-
-    assert x.shape == (n, 3)
-    assert np.all((0.0 <= x[:, 0]) & (x[:, 0] <= 1.0))
-    assert np.all(x[:, 1] >= 0.0)
-    assert np.all((10.0 <= x[:, 2]) & (x[:, 2] <= 12.0))
-
-
-def test_product_measure_recursive_transform_sampling_supported_but_weights_restricted():
-    recursive_marginal = Uniform(
-        Uniform(DummySampler(1), lower_bound=0.0, upper_bound=1.0),
-        lower_bound=2.0,
-        upper_bound=4.0,
-    )
-    direct_marginal = Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0)
-    tm = ProductMeasure(
-        sampler=DigitalNetB2(2, seed=91),
-        marginals=[recursive_marginal, direct_marginal],
-    )
-
-    x = tm(16)
-
-    assert x.shape == (16, 2)
-    assert np.all((2.0 <= x[:, 0]) & (x[:, 0] <= 4.0))
-    assert np.all((10.0 <= x[:, 1]) & (x[:, 1] <= 12.0))
-    with pytest.raises(ParameterError, match="direct marginal"):
-        tm(16, return_weights=True)
diff --git a/test/test_accumulate_data.py b/test/test_sc_accumulate_data.py
similarity index 100%
rename from test/test_accumulate_data.py
rename to test/test_sc_accumulate_data.py
diff --git a/test/test_cubbayes_vec.py b/test/test_sc_cubbayes_vec.py
similarity index 100%
rename from test/test_cubbayes_vec.py
rename to test/test_sc_cubbayes_vec.py
diff --git a/test/test_stopping_criteria.py b/test/test_sc_stopping_criteria.py
similarity index 100%
rename from test/test_stopping_criteria.py
rename to test/test_sc_stopping_criteria.py
diff --git a/test/test_scipy_wrapper_custom.py b/test/test_scipy_wrapper_custom.py
deleted file mode 100644
index dd713e934..000000000
--- a/test/test_scipy_wrapper_custom.py
+++ /dev/null
@@ -1,302 +0,0 @@
-import warnings
-
-import pytest
-import numpy as np
-import scipy.stats as stats
-
-from qmcpy import DigitalNetB2, SciPyWrapper, StudentT, ZeroInflatedExpUniform
-
-from qmcpy.true_measure.triangular import TriangularDistribution
-from qmcpy.util import DimensionError, ParameterError
-
-
-MISSING_PDF_WARNING = "no 'pdf' or 'logpdf'"
-
-
-def _missing_pdf_warnings(caught):
-    return [
-        warning
-        for warning in caught
-        if issubclass(warning.category, UserWarning)
-        and MISSING_PDF_WARNING in str(warning.message)
-    ]
-
-
-def test_mvn_dependence_correlation_and_moment():
-    """
-    Check that passing a SciPy multivariate normal through SciPyWrapper
-    preserves correlation and the mixed moment E[X1 X2].
-    """
-    sampler = DigitalNetB2(2, seed=5)
-    rho_target = 0.7
-    cov = [[1.0, rho_target], [rho_target, 1.0]]
-    mvn = stats.multivariate_normal(mean=[0.0, 0.0], cov=cov)
-    tm_mvn = SciPyWrapper(sampler, scipy_distribs=mvn)
-
-    n = 4096
-    x = tm_mvn(n)
-
-    rho_hat = np.corrcoef(x.T)[0, 1]
-    est_moment = np.mean(x[:, 0] * x[:, 1])
-
-    assert np.isfinite(rho_hat)
-    assert np.isfinite(est_moment)
-
-    assert abs(rho_hat - rho_target) < 0.05
-    assert abs(est_moment - rho_target) < 0.05
-
-
-def test_triangular_custom_marginal_range_and_shape():
-    """
-    Make sure our custom triangular marginal behaves sensibly:
-    samples stay in the right interval and the empirical mean is close
-    to the analytic mean.
-    """
-    tri = TriangularDistribution(c=0.3, loc=-1.0, scale=2.0)
-    tm = SciPyWrapper(DigitalNetB2(1, seed=11), scipy_distribs=tri)
-
-    n = 4096
-    x = tm(n).ravel()
-
-    assert x.min() >= -1.1
-    assert x.max() <= 1.1
-
-    a = -1.0
-    b = 1.0
-    m = -1.0 + 0.3 * 2.0
-    true_mean = (a + b + m) / 3.0
-    emp_mean = x.mean()
-    assert abs(emp_mean - true_mean) < 0.05
-
-
-def test_zero_inflated_zero_rate():
-    """
-    Check that the zero-inflated exponential distribution preserves the
-    specified probability mass at X = 0.
-    """
-    p_zero = 0.4
-    sampler = DigitalNetB2(1, seed=17)
-    tm = ZeroInflatedExpUniform(sampler, p_zero=p_zero, lam=1.5)
-
-    n = 4096
-    samples = tm(n)
-    x = samples.ravel()
-    zero_rate = np.mean(x == 0.0)
-
-    assert samples.shape == (n, 1)
-    assert abs(zero_rate - p_zero) < 0.05
-
-
-def test_zero_inflated_replications_shape():
-    tm = ZeroInflatedExpUniform(
-        DigitalNetB2(1, seed=17, replications=2),
-        p_zero=0.4,
-        lam=1.5,
-    )
-
-    x = tm(8)
-
-    assert x.shape == (2, 8, 1)
-    assert np.all(x >= 0.0)
-
-
-@pytest.mark.parametrize("p_zero", [0.0, 1.0, -0.1, 1.1])
-def test_zero_inflated_rejects_invalid_p_zero(p_zero):
-    with pytest.raises(ParameterError, match="p_zero must be in"):
-        ZeroInflatedExpUniform(
-            DigitalNetB2(1, seed=17),
-            p_zero=p_zero,
-            lam=1.5,
-        )
-
-
-@pytest.mark.parametrize("lam", [0.0, -1.0])
-def test_zero_inflated_rejects_nonpositive_lam(lam):
-    with pytest.raises(ParameterError, match="lam must be positive"):
-        ZeroInflatedExpUniform(
-            DigitalNetB2(1, seed=17),
-            p_zero=0.4,
-            lam=lam,
-        )
-
-
-def test_zero_inflated_requires_one_dimensional_sampler():
-    with pytest.raises(
-        DimensionError,
-        match="requires a one-dimensional sampler",
-    ):
-        ZeroInflatedExpUniform(
-            DigitalNetB2(2, seed=17),
-            p_zero=0.4,
-            lam=1.5,
-        )
-
-
-def test_zero_inflated_inverse_transform_exact_values():
-    tm = ZeroInflatedExpUniform(
-        DigitalNetB2(1, seed=17),
-        p_zero=0.4,
-        lam=2.0,
-    )
-    u = np.array([[0.0], [0.2], [0.4], [0.7], [0.9]])
-
-    x = tm._transform(u)
-
-    assert x.shape == (5, 1)
-    assert np.array_equal(x[:3], np.zeros((3, 1)))
-    assert np.all(x[3:] > 0.0)
-
-    u_positive = u[3:, 0]
-    u_rescaled = (u_positive - 0.4) / 0.6
-    expected = -np.log1p(-u_rescaled) / 2.0
-    assert np.allclose(x[3:, 0], expected)
-
-
-def test_zero_inflated_inverse_transform_all_zero_branch():
-    tm = ZeroInflatedExpUniform(
-        DigitalNetB2(1, seed=17),
-        p_zero=0.4,
-        lam=2.0,
-    )
-    u = np.array([[0.0], [0.1], [0.4]])
-
-    x = tm._transform(u)
-
-    assert x.shape == (3, 1)
-    assert np.array_equal(x, np.zeros((3, 1)))
-
-
-def test_zero_inflated_inverse_transform_clips_one():
-    tm = ZeroInflatedExpUniform(
-        DigitalNetB2(1, seed=17),
-        p_zero=0.4,
-        lam=2.0,
-    )
-    u = np.array([[1.0]])
-
-    x = tm._transform(u)
-
-    assert x.shape == (1, 1)
-    assert np.isfinite(x).all()
-    assert x[0, 0] > 0.0
-
-
-def test_zero_inflated_construction_does_not_warn_about_missing_pdf():
-    with warnings.catch_warnings(record=True) as caught:
-        warnings.simplefilter("always")
-        tm = ZeroInflatedExpUniform(
-            DigitalNetB2(1, seed=17),
-            p_zero=0.4,
-            lam=1.5,
-        )
-
-    assert tm.d == 1
-    assert _missing_pdf_warnings(caught) == []
-
-
-def test_zero_inflated_sampling_does_not_warn_about_missing_pdf():
-    tm = ZeroInflatedExpUniform(DigitalNetB2(1, seed=17), p_zero=0.4, lam=1.5)
-
-    with warnings.catch_warnings(record=True) as caught:
-        warnings.simplefilter("always")
-        x = tm(8)
-
-    assert x.shape == (8, 1)
-    assert _missing_pdf_warnings(caught) == []
-
-
-def test_zero_inflated_return_weights_warns_once_for_missing_pdf():
-    tm = ZeroInflatedExpUniform(DigitalNetB2(1, seed=17), p_zero=0.4, lam=1.5)
-
-    with pytest.warns(UserWarning, match=MISSING_PDF_WARNING):
-        x, jac = tm(8, return_weights=True)
-
-    assert x.shape == (8, 1)
-    assert jac.shape == (8,)
-    assert np.allclose(jac, 1.0)
-
-    with warnings.catch_warnings(record=True) as caught:
-        warnings.simplefilter("always")
-        x_second, jac_second = tm(8, return_weights=True)
-
-    assert x_second.shape == (8, 1)
-    assert np.allclose(jac_second, 1.0)
-    assert _missing_pdf_warnings(caught) == []
-
-
-def test_zero_inflated_y_split_warns_and_uses_one_dimensional_interface():
-    with pytest.warns(DeprecationWarning, match="y_split"):
-        tm = ZeroInflatedExpUniform(
-            DigitalNetB2(1, seed=17),
-            p_zero=0.4,
-            lam=1.5,
-            y_split=0.5,
-        )
-
-    x = tm(4)
-
-    assert x.shape == (4, 1)
-    assert np.all(x >= 0.0)
-
-
-def test_zero_inflated_y_split_preserves_deprecated_two_dimensional_usage():
-    with pytest.warns(DeprecationWarning, match="2D zero-inflated"):
-        tm = ZeroInflatedExpUniform(
-            DigitalNetB2(2, seed=17),
-            p_zero=0.4,
-            lam=1.5,
-            y_split=0.5,
-        )
-
-    x = tm(16)
-
-    assert x.shape == (16, 2)
-    assert np.all(x[:, 0] >= 0.0)
-    assert np.all((0.0 <= x[:, 1]) & (x[:, 1] <= 1.0))
-    assert np.all(x[x[:, 0] == 0.0, 1] <= 0.5)
-    assert np.all(x[x[:, 0] > 0.0, 1] >= 0.5)
-
-
-def test_zero_inflated_y_split_preserves_replicated_two_dimensional_usage():
-    with pytest.warns(DeprecationWarning, match="2D zero-inflated"):
-        tm = ZeroInflatedExpUniform(
-            DigitalNetB2(2, seed=17, replications=2),
-            p_zero=0.4,
-            lam=1.5,
-            y_split=0.5,
-        )
-
-    x = tm(16)
-
-    assert x.shape == (2, 16, 2)
-    assert np.all(x[..., 0] >= 0.0)
-    assert np.all((0.0 <= x[..., 1]) & (x[..., 1] <= 1.0))
-    assert np.all(x[..., 1][x[..., 0] == 0.0] <= 0.5)
-    assert np.all(x[..., 1][x[..., 0] > 0.0] >= 0.5)
-
-
-def test_student_t_marginals_shape():
-    tm = SciPyWrapper(
-        sampler=DigitalNetB2(2, seed=5),
-        scipy_distribs=stats.t(df=5),
-    )
-    x = tm(8)
-    assert x.shape == (8, 2)
-
-
-def test_multivariate_student_t_joint_corr_and_cov():
-    if not hasattr(stats, "multivariate_t"):
-        pytest.skip("scipy.stats.multivariate_t not available in this SciPy version")
-
-    df = 5.0
-    rho = 0.8
-    loc = np.array([0.0, 0.0])
-    shape = np.array([[1.0, rho], [rho, 1.0]])
-
-    tm = StudentT(DigitalNetB2(2, seed=123), loc=loc, shape=shape, df=df)
-
-    n = 4096
-    x = tm(n)
-    emp_corr = np.corrcoef(x.T)[0, 1]
-
-    assert abs(emp_corr - rho) < 0.05
diff --git a/test/test_sr_check_links.py b/test/test_sr_check_links.py
new file mode 100644
index 000000000..49c148256
--- /dev/null
+++ b/test/test_sr_check_links.py
@@ -0,0 +1,198 @@
+import contextlib
+import io
+import shutil
+import ssl
+import sys
+import tempfile
+import unittest
+import urllib.error
+from pathlib import Path
+from unittest.mock import patch
+
+from scripts import check_links
+
+
+def _http_error(url, code):
+    return urllib.error.HTTPError(url, code, "test response", {}, None)
+
+
+class TestCheckLinks(unittest.TestCase):
+
+    def setUp(self):
+        self.tmp_path = Path(tempfile.mkdtemp())
+        self.addCleanup(shutil.rmtree, self.tmp_path, ignore_errors=True)
+
+    def _patch(self, target, name, value):
+        """monkeypatch.setattr equivalent: set now, auto-restore at test end."""
+        patcher = patch.object(target, name, value)
+        patcher.start()
+        self.addCleanup(patcher.stop)
+
+    def test_head_success_is_reachable(self):
+        with patch.object(check_links.urllib.request, "urlopen", return_value=object()) as urlopen:
+            self.assertIsNone(check_links._check_one("https://example.test", timeout=1))
+
+        self.assertEqual(urlopen.call_count, 1)
+        self.assertEqual(urlopen.call_args.args[0].get_method(), "HEAD")
+
+    def test_get_success_after_head_failure_is_reachable(self):
+        url = "https://example.test"
+        with patch.object(
+            check_links.urllib.request,
+            "urlopen",
+            side_effect=[_http_error(url, 405), object()],
+        ) as urlopen:
+            self.assertIsNone(check_links._check_one(url, timeout=1))
+
+        self.assertEqual(urlopen.call_count, 2)
+        self.assertEqual(urlopen.call_args_list[1].args[0].get_method(), "GET")
+
+    def test_not_found_and_gone_gets_are_broken(self):
+        for code in (404, 410):
+            with self.subTest(code=code):
+                url = f"https://example.test/{code}"
+                with patch.object(
+                    check_links.urllib.request,
+                    "urlopen",
+                    side_effect=[_http_error(url, code), _http_error(url, code)],
+                ):
+                    self.assertEqual(
+                        check_links._check_one(url, timeout=1),
+                        ("broken", f"{url} -- HTTP {code}"),
+                    )
+
+    def test_bot_block_and_rate_limit_are_warnings(self):
+        for code in (403, 429):
+            with self.subTest(code=code):
+                url = f"https://example.test/{code}"
+                with patch.object(
+                    check_links.urllib.request,
+                    "urlopen",
+                    side_effect=[_http_error(url, code), _http_error(url, code)],
+                ):
+                    severity, message = check_links._check_one(url, timeout=1)
+
+                self.assertEqual(severity, "warning")
+                self.assertIn(f"HTTP {code}", message)
+
+    def test_tls_and_timeout_failures_are_warnings(self):
+        failures = (
+            ssl.SSLCertVerificationError("certificate verify failed"),
+            TimeoutError("timed out"),
+        )
+        for failure in failures:
+            with self.subTest(failure=type(failure).__name__):
+                with patch.object(
+                    check_links.urllib.request,
+                    "urlopen",
+                    side_effect=[failure, failure],
+                ):
+                    severity, message = check_links._check_one(
+                        "https://example.test", timeout=1
+                    )
+
+                self.assertEqual(severity, "warning")
+                self.assertIn(str(failure), message)
+
+    def test_external_results_are_separated_and_duplicate_urls_checked_once(self):
+        (self.tmp_path / "page.html").write_text(
+            '<a href="https://example.test/missing">missing</a>'
+            '<a href="https://example.test/missing">duplicate</a>'
+            '<a href="https://example.test/blocked">blocked</a>',
+            encoding="utf-8",
+        )
+
+        def result_for(url, _timeout):
+            if url.endswith("/missing"):
+                return "broken", f"{url} -- HTTP 404"
+            return "warning", f"{url} -- HTTP 403"
+
+        with patch.object(check_links, "_check_one", side_effect=result_for) as check_one:
+            broken, warnings = check_links.check_external(self.tmp_path, workers=1)
+
+        self.assertEqual(check_one.call_count, 2)
+        self.assertEqual(
+            broken,
+            ["https://example.test/missing -- HTTP 404 (seen on page.html)"],
+        )
+        self.assertEqual(
+            warnings,
+            ["https://example.test/blocked -- HTTP 403 (seen on page.html)"],
+        )
+
+    def test_internal_links_strip_site_url_deployment_path(self):
+        target = self.tmp_path / "target"
+        target.mkdir()
+        (target / "index.html").write_text(
+            '<h2 id="section">Target</h2>', encoding="utf-8"
+        )
+        (self.tmp_path / "index.html").write_text(
+            '<a href="/QMCSoftware/target/">root-relative</a>'
+            '<a href="https://qmcsoftware.github.io/QMCSoftware/target/#section">absolute</a>',
+            encoding="utf-8",
+        )
+
+        self.assertEqual(
+            check_links.check_internal(
+                self.tmp_path, site_url="https://qmcsoftware.github.io/QMCSoftware/"
+            ),
+            [],
+        )
+
+    def test_external_check_skips_same_site_urls(self):
+        (self.tmp_path / "page.html").write_text(
+            '<a href="https://qmcsoftware.github.io/QMCSoftware/target/">same</a>'
+            '<a href="https://example.test/target/">external</a>',
+            encoding="utf-8",
+        )
+
+        with patch.object(check_links, "_check_one", return_value=None) as check_one:
+            broken, warnings = check_links.check_external(
+                self.tmp_path,
+                workers=1,
+                site_url="https://qmcsoftware.github.io/QMCSoftware/",
+            )
+
+        self.assertEqual(broken, [])
+        self.assertEqual(warnings, [])
+        self.assertEqual(check_one.call_count, 1)
+        self.assertEqual(check_one.call_args.args[0], "https://example.test/target/")
+
+    def test_external_warnings_do_not_make_main_fail(self):
+        self._patch(sys, "argv", ["check_links.py", str(self.tmp_path), "--external"])
+        self._patch(
+            check_links, "check_internal", lambda _site_dir, site_url=None: []
+        )
+        self._patch(
+            check_links,
+            "check_external",
+            lambda _site_dir, site_url=None: (
+                [],
+                ["https://example.test -- HTTP 403"],
+            ),
+        )
+
+        buf = io.StringIO()
+        with contextlib.redirect_stdout(buf):
+            self.assertEqual(check_links.main(), 0)
+        self.assertIn("0 broken link(s), 1 warning(s)", buf.getvalue())
+
+    def test_confirmed_external_breakage_makes_main_fail(self):
+        self._patch(sys, "argv", ["check_links.py", str(self.tmp_path), "--external"])
+        self._patch(
+            check_links, "check_internal", lambda _site_dir, site_url=None: []
+        )
+        self._patch(
+            check_links,
+            "check_external",
+            lambda _site_dir, site_url=None: (
+                ["https://example.test -- HTTP 404"],
+                [],
+            ),
+        )
+
+        self.assertEqual(check_links.main(), 1)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/test_sr_check_removed_urls.py b/test/test_sr_check_removed_urls.py
new file mode 100644
index 000000000..57a30223d
--- /dev/null
+++ b/test/test_sr_check_removed_urls.py
@@ -0,0 +1,165 @@
+import contextlib
+import io
+import shutil
+import sys
+import tempfile
+import unittest
+import urllib.error
+from pathlib import Path
+from unittest.mock import patch
+
+from scripts import check_removed_urls as cru
+
+SITE = "https://qmcsoftware.github.io/QMCSoftware/"
+
+
+def _sitemap(*paths):
+    locs = "".join(f"<loc>{SITE}{path}</loc>" for path in paths)
+    return f'<?xml version="1.0" encoding="UTF-8"?><urlset>{locs}</urlset>'
+
+
+def _config(redirect_maps=None):
+    plugins = ["material/search", {"mkdocs-jupyter": {"execute": False}}]
+    if redirect_maps is not None:
+        plugins.append({"redirects": {"redirect_maps": redirect_maps}})
+    return {"site_url": SITE, "plugins": plugins}
+
+
+class TestCheckRemovedUrls(unittest.TestCase):
+
+    def setUp(self):
+        self.tmp_path = Path(tempfile.mkdtemp())
+        self.addCleanup(shutil.rmtree, self.tmp_path, ignore_errors=True)
+        self._last_out = ""
+
+    def _patch(self, target, name, value):
+        """monkeypatch.setattr equivalent: set now, auto-restore at test end."""
+        patcher = patch.object(target, name, value)
+        patcher.start()
+        self.addCleanup(patcher.stop)
+
+    def _run(self, sitemap_paths, redirect_maps=None, extra_argv=(), base=None):
+        """Run main() offline against a temp sitemap and a temp docs/ tree."""
+        base = self.tmp_path if base is None else base
+        docs = base / "docs"
+        docs.mkdir(parents=True)
+        (docs / "README.md").write_text("home", encoding="utf-8")
+        (docs / "good_practices.md").write_text("page", encoding="utf-8")
+        sitemap = base / "sitemap.xml"
+        sitemap.write_text(_sitemap(*sitemap_paths), encoding="utf-8")
+
+        self._patch(cru, "read_config", lambda *a, **k: _config(redirect_maps))
+        self._patch(sys, "argv", [
+            "check_removed_urls.py", "--sitemap", str(sitemap), "--docs-dir", str(docs),
+            *extra_argv,
+        ])
+        buf = io.StringIO()
+        with contextlib.redirect_stdout(buf):
+            code = cru.main()
+        self._last_out = buf.getvalue()
+        return code
+
+    def test_url_path_and_source_round_trip(self):
+        for source, url_path in [("blogs/scipywrapper/index.md", "blogs/scipywrapper/"),
+                                 ("good_practices.md", "good_practices/"),
+                                 ("demos/quickstart.ipynb", "demos/quickstart/"),
+                                 ("index.md", ""), ("README.md", "")]:
+            self.assertEqual(cru.url_path_for_source(source), url_path)
+
+        for source in ("README.md", "good_practices.md", "demos/quickstart.ipynb",
+                       "api/index.md"):
+            path = self.tmp_path / source
+            path.parent.mkdir(parents=True, exist_ok=True)
+            path.write_text("page", encoding="utf-8")
+            self.assertTrue(
+                cru.source_exists(cru.url_path_for_source(source), self.tmp_path)
+            )
+        self.assertFalse(cru.source_exists("blogs/scipywrapper/", self.tmp_path))
+
+    def test_redirect_maps_reads_the_plugin_and_tolerates_its_absence(self):
+        entry = {"blogs/x/index.md": "https://qmcsoftware.org/blogs/x/"}
+        self.assertEqual(cru.redirect_maps(_config(entry)), entry)
+        self.assertEqual(cru.redirect_maps(_config()), {})
+        self.assertEqual(cru.redirect_maps({}), {})
+
+    def test_published_paths_separates_foreign_urls(self):
+        sitemap = _sitemap("", "good_practices/").replace(
+            "</urlset>", "<loc>https://example.test/other/</loc></urlset>")
+
+        self.assertEqual(
+            cru.published_paths(sitemap, SITE),
+            (["", "good_practices/"], ["https://example.test/other/"]),
+        )
+
+    def test_http_status_falls_back_to_get_when_head_is_unsupported(self):
+        url = "https://example.test"
+        error = urllib.error.HTTPError(url, 405, "test response", {}, None)
+        response = type("Response", (), {"status": 200, "__enter__": lambda s: s,
+                                         "__exit__": lambda s, *a: False})()
+        with patch.object(cru.urllib.request, "urlopen",
+                          side_effect=[error, response]) as urlopen:
+            self.assertEqual(cru.http_status(url, timeout=1), "200")
+
+        self.assertEqual(urlopen.call_count, 2)
+        self.assertEqual(urlopen.call_args_list[1].args[0].get_method(), "GET")
+
+    def test_removed_page_without_redirect_is_flagged(self):
+        code = self._run(["", "good_practices/", "blogs/scipywrapper/"])
+        out = self._last_out
+
+        self.assertEqual(code, 1)
+        self.assertIn("1 removed with no redirect", out)
+        self.assertIn(f"[ORPHAN]   {SITE}blogs/scipywrapper/", out)
+        self.assertIn("blogs/scipywrapper/index.md: <new URL or page>", out)
+
+    def test_removed_page_covered_by_a_redirect_passes(self):
+        code = self._run(
+            ["", "good_practices/", "blogs/scipywrapper/"],
+            redirect_maps={
+                "blogs/scipywrapper/index.md": "https://qmcsoftware.org/blogs/scipywrapper/"},
+        )
+        out = self._last_out
+
+        self.assertEqual(code, 0)
+        self.assertIn("0 removed with no redirect", out)
+        self.assertIn("[redirect]", out)
+        self.assertNotIn("[ORPHAN]", out)
+
+    def test_intact_site_passes(self):
+        self.assertEqual(self._run(["", "good_practices/"]), 0)
+        self.assertIn("2 still have a page source", self._last_out)
+
+    def test_verify_redirects_follows_the_target_status(self):
+        redirects = {"blogs/x/index.md": "https://qmcsoftware.org/blogs/x/"}
+        for status, expected_code in [("200", 0), ("404", 1)]:
+            with self.subTest(status=status):
+                self._patch(cru, "http_status", lambda *a, **k: status)
+                code = self._run(
+                    ["", "blogs/x/"],
+                    redirect_maps=redirects,
+                    extra_argv=("--verify-redirects",),
+                    base=self.tmp_path / status,
+                )
+                out = self._last_out
+
+                self.assertEqual(code, expected_code)
+                self.assertIn(status, out)
+                # The URL itself is covered, so a failure is the target, not an orphan.
+                self.assertNotIn("[ORPHAN]", out)
+
+    def test_unreachable_sitemap_fails_unless_offline_is_allowed(self):
+        self._patch(cru, "read_config", lambda *a, **k: _config())
+        argv = ["check_removed_urls.py", "--sitemap", str(self.tmp_path / "absent.xml")]
+
+        self._patch(sys, "argv", argv)
+        self.assertEqual(cru.main(), 1)
+
+        self._patch(sys, "argv", argv + ["--allow-offline"])
+        buf = io.StringIO()
+        with contextlib.redirect_stdout(buf):
+            self.assertEqual(cru.main(), 0)
+        self.assertIn("skipping the check", buf.getvalue())
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/test_sr_flatten_qmcpy_imports.py b/test/test_sr_flatten_qmcpy_imports.py
new file mode 100644
index 000000000..5a295121a
--- /dev/null
+++ b/test/test_sr_flatten_qmcpy_imports.py
@@ -0,0 +1,351 @@
+import json
+import shutil
+import tempfile
+import unittest
+from pathlib import Path
+
+from scripts.flatten_qmcpy_imports import (
+    _load_qmcpy_public_names,
+    flatten_imports,
+    main,
+)
+
+
+def _nested_import(module, imported):
+    return f"from {'qmcpy.' + module} import {imported}"
+
+
+class TestFlattenQmcpyImports(unittest.TestCase):
+
+    def setUp(self):
+        self.tmp_path = Path(tempfile.mkdtemp())
+        self.addCleanup(shutil.rmtree, self.tmp_path, ignore_errors=True)
+
+    def test_flatten_imports_basic(self):
+        source = (
+            _nested_import("integrand", "Keister")
+            + "\n"
+            + _nested_import("discrete_distribution.lattice", "Lattice as LD")
+            + "\nfrom qmcpy import DigitalNetB2\nimport qmcpy.util\n"
+        ).encode()
+
+        updated, count = flatten_imports(
+            source, frozenset({"DigitalNetB2", "Keister", "Lattice"})
+        )
+
+        self.assertEqual(count, 3)
+        self.assertEqual(
+            updated,
+            (
+                b"from qmcpy import DigitalNetB2, Keister, Lattice as LD\n"
+                b"import qmcpy.util\n"
+            ),
+        )
+
+    def test_flatten_preserves_private(self):
+        source = (
+            _nested_import("_internal._helpers", "PublicHelper")
+            + "\n"
+            + _nested_import(
+                "true_measure.uniform_triangle",
+                "UniformTriangle, _UniformTriangleAdapter",
+            )
+            + "\n"
+            + _nested_import(
+                "true_measure.copula",
+                "(\n    AbstractCopula,\n    _validate_dimension,\n)",
+            )
+            + "\n"
+            + _nested_import("integrand", "Keister")
+            + "\n"
+        ).encode()
+
+        updated, count = flatten_imports(source, frozenset({"Keister"}))
+
+        self.assertEqual(count, 1)
+        self.assertEqual(
+            updated,
+            source.replace(
+                _nested_import("integrand", "Keister").encode(),
+                b"from qmcpy import Keister",
+            ),
+        )
+
+    def test_private_module_splits_groups(self):
+        source = (
+            b"from qmcpy import Zeta\n"
+            b"from qmcpy._internal._helpers import PublicHelper\n"
+            b"from qmcpy import Alpha\n"
+        )
+
+        updated, count = flatten_imports(source)
+
+        self.assertEqual(count, 0)
+        self.assertEqual(updated, source)
+
+    def test_flatten_preserves_util_imports(self):
+        source = (
+            b"from qmcpy.util import ParameterError\n"
+            b"from qmcpy.util.transforms import tf_exp\n"
+        )
+
+        updated, count = flatten_imports(source, frozenset({"ParameterError", "tf_exp"}))
+
+        self.assertEqual((updated, count), (source, 0))
+
+    def test_flatten_keeps_nonpublic_names(self):
+        source = b"from qmcpy.stopping_criterion.pf_gp_ci import PFGPCIData\n"
+
+        updated, count = flatten_imports(source, frozenset({"PFGPCI"}))
+
+        self.assertEqual((updated, count), (source, 0))
+
+    def test_flatten_no_public_api_noop(self):
+        source = (_nested_import("integrand", "Keister") + "\n").encode()
+
+        updated, count = flatten_imports(source)
+
+        self.assertEqual((updated, count), (source, 0))
+
+    def test_flatten_preserve_str_literals(self):
+        source = b'text = """\nfrom qmcpy.integrand import Keister\n"""\n'
+
+        updated, count = flatten_imports(source, frozenset({"Keister"}))
+
+        self.assertEqual(count, 0)
+        self.assertEqual(updated, source)
+
+    def test_python_string_protection_applies_to_every_rewrite_stage(self):
+        string_body = (
+            b'text = """\n'
+            b"from qmcpy.integrand import Keister\n"
+            b"from qmcpy import Zeta,Beta\n"
+            b"from qmcpy import Alpha\n"
+            b"from qmcpy import *\n"
+            b"from qmcpy import *\n"
+            b'"""\n'
+        )
+        source = string_body + b"from qmcpy.integrand import Keister\n"
+
+        updated, count = flatten_imports(source, frozenset({"Keister"}))
+
+        self.assertEqual(count, 1)
+        self.assertEqual(updated, string_body + b"from qmcpy import Keister\n")
+
+    def test_python_tokenize_failure_is_fail_closed(self):
+        source = b'"""unterminated\nfrom qmcpy.integrand import Keister\n'
+
+        self.assertEqual(
+            flatten_imports(source, frozenset({"Keister"})), (source, 0)
+        )
+
+    def test_flatten_skip_star_expansion(self):
+        source = (
+            b"from qmcpy import *\n\n"
+            b"def f(Lattice):\n"
+            b"    return Lattice\n\n"
+            b"y = Keister(dimension=2)\n"
+            b"x = Lattice(dimension=2)\n"
+        )
+
+        updated, count = flatten_imports(source, frozenset({"Keister", "Lattice"}))
+
+        self.assertEqual(count, 0)
+        self.assertEqual(updated, source)
+
+    def test_notebook_star_dedup(self):
+        notebook = {
+            "cells": [
+                {
+                    "cell_type": "code",
+                    "source": [
+                        _nested_import("integrand", "*") + "\n",
+                        _nested_import("true_measure", "*"),
+                    ],
+                }
+            ]
+        }
+        source = json.dumps(notebook, indent=1).encode()
+
+        updated, count = flatten_imports(source, frozenset({"Keister"}))
+
+        self.assertEqual(count, 3)
+        self.assertEqual(
+            json.loads(updated)["cells"][0]["source"], ["from qmcpy import *"]
+        )
+
+    def test_named_imports_merge_sort(self):
+        source = (
+            b"from qmcpy import Zeta,Beta\n"
+            b"from qmcpy import Alpha\n"
+            b"\n"
+            b"from qmcpy import Gamma\n"
+        )
+
+        updated, count = flatten_imports(source)
+
+        self.assertEqual(count, 1)
+        self.assertEqual(
+            updated,
+            (
+                b"from qmcpy import Alpha, Beta, Zeta\n"
+                b"\n"
+                b"from qmcpy import Gamma\n"
+            ),
+        )
+        self.assertEqual(flatten_imports(updated), (updated, 0))
+
+    def test_merge_paren_and_single_line(self):
+        source = b"""from qmcpy import (
+    KernelDigShiftInvar,
+    KernelDigShiftInvarAdaptiveAlpha,
+    KernelDigShiftInvarCombined,
+    KernelShiftInvar,
+    KernelShiftInvarCombined,
+)
+from qmcpy import tf_exp_eps, tf_exp_eps_inv
+"""
+
+        updated, count = flatten_imports(source)
+
+        self.assertEqual(count, 1)
+        self.assertEqual(
+            updated,
+            b"""from qmcpy import (
+    KernelDigShiftInvar,
+    KernelDigShiftInvarAdaptiveAlpha,
+    KernelDigShiftInvarCombined,
+    KernelShiftInvar,
+    KernelShiftInvarCombined,
+    tf_exp_eps,
+    tf_exp_eps_inv,
+)
+""",
+        )
+        self.assertEqual(flatten_imports(updated), (updated, 0))
+
+    def test_merge_same_scope_only(self):
+        source = (
+            b"if enabled:\n"
+            b"    from qmcpy import Zeta\n"
+            b"    from qmcpy import Alpha as First\n"
+            b"else:\n"
+            b"    from qmcpy import Beta\n"
+            b"from qmcpy import _Private\n"
+            b"from qmcpy import Gamma  # keep this comment\n"
+        )
+
+        updated, count = flatten_imports(source)
+
+        self.assertEqual(count, 1)
+        self.assertEqual(
+            updated,
+            (
+                b"if enabled:\n"
+                b"    from qmcpy import Alpha as First, Zeta\n"
+                b"else:\n"
+                b"    from qmcpy import Beta\n"
+                b"from qmcpy import _Private\n"
+                b"from qmcpy import Gamma  # keep this comment\n"
+            ),
+        )
+
+    def test_notebook_named_merge(self):
+        notebook = {
+            "cells": [
+                {
+                    "cell_type": "code",
+                    "source": [
+                        "from qmcpy import Zeta\n",
+                        "from qmcpy import Alpha,Beta\n",
+                        "print(Alpha)\n",
+                    ],
+                }
+            ]
+        }
+        source = json.dumps(notebook, indent=1).encode()
+
+        updated, count = flatten_imports(source)
+
+        self.assertEqual(count, 1)
+        self.assertEqual(
+            json.loads(updated)["cells"][0]["source"],
+            [
+                "from qmcpy import Alpha, Beta, Zeta\n",
+                "print(Alpha)\n",
+            ],
+        )
+        self.assertEqual(flatten_imports(updated), (updated, 0))
+
+    def test_notebook_flattens_nested_imports_only_in_code_cells(self):
+        nested_import = _nested_import("integrand", "Keister") + "\n"
+        metadata_import = _nested_import("true_measure", "Gaussian") + "\n"
+        string_literal = f'text = "{nested_import.rstrip()}"\n'
+        multiline_string = ['text = """\n', nested_import, '"""\n']
+        notebook = {
+            "metadata": {"source": [metadata_import]},
+            "cells": [
+                {"cell_type": "markdown", "source": [nested_import]},
+                {"cell_type": "code", "source": [nested_import]},
+                {"cell_type": "code", "source": [string_literal]},
+                {"cell_type": "code", "source": multiline_string},
+            ]
+        }
+        source = json.dumps(notebook, indent=1).encode()
+
+        updated, count = flatten_imports(source, frozenset({"Keister"}))
+
+        cells = json.loads(updated)["cells"]
+        self.assertEqual(count, 1)
+        self.assertEqual(
+            json.loads(updated)["metadata"]["source"], [metadata_import]
+        )
+        self.assertEqual(cells[0]["source"], [nested_import])
+        self.assertEqual(cells[1]["source"], ["from qmcpy import Keister\n"])
+        self.assertEqual(cells[2]["source"], [string_literal])
+        self.assertEqual(cells[3]["source"], multiline_string)
+        self.assertEqual(
+            flatten_imports(updated, frozenset({"Keister"})), (updated, 0)
+        )
+
+    def test_markdown_import_examples_are_flattened(self):
+        path = self.tmp_path / "example.md"
+        path.write_bytes(
+            b'Example with unmatched prose delimiter: """\n\n'
+            b"```python\n"
+            b"from qmcpy.integrand import Keister\n"
+            b"```\n"
+        )
+
+        self.assertEqual(main([str(path)]), 0)
+        self.assertIn(b"from qmcpy import Keister", path.read_bytes())
+
+    def test_check_mode_no_write(self):
+        path = self.tmp_path / "example.py"
+        original = (_nested_import("true_measure", "Gaussian") + "\n").encode()
+        path.write_bytes(original)
+
+        self.assertEqual(main(["--check", str(path)]), 1)
+        self.assertEqual(path.read_bytes(), original)
+
+        self.assertEqual(main([str(path)]), 0)
+        self.assertEqual(path.read_bytes(), b"from qmcpy import Gaussian\n")
+        self.assertEqual(main(["--check", str(path)]), 0)
+
+    def test_public_names_optional_free_stable(self):
+        repository_root = Path(__file__).resolve().parent.parent
+        names = _load_qmcpy_public_names(repository_root)
+
+        self.assertIsNotNone(names)
+        self.assertIn("Gaussian", names)
+        self.assertIn("Keister", names)
+        # Optional dependencies are blocked in the probe context, so fallback
+        # exports are part of the deterministic name set.
+        self.assertIn("PFGPCI", names)
+        # Helpers that are deliberately not part of the top-level API.
+        self.assertNotIn("PFGPCIData", names)
+        self.assertNotIn("TriangularDistribution", names)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/test_sr_install_mpmc_pyg.py b/test/test_sr_install_mpmc_pyg.py
new file mode 100644
index 000000000..25e810ec1
--- /dev/null
+++ b/test/test_sr_install_mpmc_pyg.py
@@ -0,0 +1,92 @@
+"""Tests for the platform-specific MPMC dependency installer."""
+
+import subprocess
+import unittest
+from types import SimpleNamespace
+from unittest.mock import patch
+
+from qmcpy.util import install_mpmc_pyg
+
+
+def _torch(version="2.12.1+cpu", cuda=None, hip=None):
+    return SimpleNamespace(
+        __version__=version,
+        version=SimpleNamespace(cuda=cuda, hip=hip),
+    )
+
+
+class TestInstallMPMCPyG(unittest.TestCase):
+
+    def test_torch_versions_include_baseline_fallback(self):
+        """Wheel lookup tries an exact patch release, then its minor baseline."""
+        self.assertEqual(
+            install_mpmc_pyg.torch_versions("2.12.1+cpu"), ["2.12.1", "2.12.0"]
+        )
+        self.assertEqual(install_mpmc_pyg.torch_versions("2.12.0"), ["2.12.0"])
+
+        with self.assertRaisesRegex(RuntimeError, "Unable to parse torch version"):
+            install_mpmc_pyg.torch_versions("development")
+
+    def test_accelerator_tag(self):
+        """PyTorch build metadata maps to the expected PyG wheel tag."""
+        cases = [
+            (_torch(), "cpu"),
+            (_torch(cuda="12.6"), "cu126"),
+            (_torch(cuda="13.0.1"), "cu130"),
+        ]
+        for torch_module, expected in cases:
+            with self.subTest(expected=expected):
+                self.assertEqual(
+                    install_mpmc_pyg.accelerator_tag(torch_module), expected
+                )
+
+    def test_accelerator_tag_rejects_rocm(self):
+        """The installer directs unsupported ROCm users to upstream guidance."""
+        with self.assertRaisesRegex(RuntimeError, "does not currently support ROCm"):
+            install_mpmc_pyg.accelerator_tag(_torch(hip="6.3"))
+
+    def test_main_retries_with_torch_minor_baseline(self):
+        """A missing exact wheel page falls back to the minor baseline page."""
+        calls = []
+
+        def fake_run(*args):
+            calls.append(args)
+            if args[-1].endswith("torch-2.12.1+cpu.html"):
+                raise subprocess.CalledProcessError(1, args)
+
+        with patch.object(install_mpmc_pyg, "run", fake_run):
+            install_mpmc_pyg.main(_torch())
+
+        self.assertEqual(calls[0][-1], "torch-geometric>=2.6.1")
+        self.assertEqual(
+            calls[1][-1], "https://data.pyg.org/whl/torch-2.12.1+cpu.html"
+        )
+        self.assertEqual(
+            calls[2][-1], "https://data.pyg.org/whl/torch-2.12.0+cpu.html"
+        )
+        self.assertIn("--only-binary", calls[1])
+
+    def test_main_explains_that_torch_must_be_installed(self):
+        """Running the helper before installing the extra gives a useful error."""
+        def missing_torch(_name):
+            raise ModuleNotFoundError("No module named 'torch'", name="torch")
+
+        with patch.object(
+            install_mpmc_pyg.importlib, "import_module", missing_torch
+        ):
+            with self.assertRaisesRegex(RuntimeError, r"install 'qmcpy\[mpmc\]'"):
+                install_mpmc_pyg.main()
+
+    def test_main_reports_missing_wheel(self):
+        """Exhausting candidate wheel pages reports the build that failed."""
+        def fail_pyg_lib(*args):
+            if "pyg_lib>=0.6.0" in args:
+                raise subprocess.CalledProcessError(1, args)
+
+        with patch.object(install_mpmc_pyg, "run", fail_pyg_lib):
+            with self.assertRaisesRegex(RuntimeError, r"torch 2\.12\.1\+cpu \(cpu\)"):
+                install_mpmc_pyg.main(_torch())
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/test_sr_mpmc_optional_imports.py b/test/test_sr_mpmc_optional_imports.py
new file mode 100644
index 000000000..4628f2c93
--- /dev/null
+++ b/test/test_sr_mpmc_optional_imports.py
@@ -0,0 +1,114 @@
+import ast
+import builtins
+import unittest
+from pathlib import Path
+
+
+def _execute_optional_import(blocked_import):
+    repository_root = Path(__file__).resolve().parent.parent
+    init_path = repository_root / "qmcpy" / "__init__.py"
+    init_tree = ast.parse(init_path.read_text())
+    optional_import = next(
+        node
+        for node in init_tree.body
+        if isinstance(node, ast.Try)
+        and any(
+            isinstance(statement, ast.ImportFrom)
+            and statement.module == "discrete_distribution.mpmc"
+            for statement in node.body
+        )
+    )
+
+    import qmcpy
+
+    real_import = builtins.__import__
+
+    def guarded_import(name, globals=None, locals=None, fromlist=(), level=0):
+        missing_module = blocked_import(name, fromlist, level)
+        if missing_module is not None:
+            raise ModuleNotFoundError(
+                "blocked optional dependency",
+                name=missing_module,
+            )
+        return real_import(name, globals, locals, fromlist, level)
+
+    test_builtins = vars(builtins).copy()
+    test_builtins["__import__"] = guarded_import
+    namespace = {"__builtins__": test_builtins, "__package__": "qmcpy"}
+    module = ast.Module(body=[optional_import], type_ignores=[])
+    exec(compile(module, str(init_path), "exec"), namespace)
+    return namespace
+
+
+class TestMPMCOptionalImports(unittest.TestCase):
+
+    def test_mpmc_utils_remain_available_without_pyg(self):
+        try:
+            import torch  # noqa: F401
+        except ImportError:
+            self.skipTest("torch not available")
+
+        def block_pyg_models(name, fromlist, level):
+            if level == 1 and name == "discrete_distribution.mpmc.models":
+                return "torch_geometric"
+            return None
+
+        namespace = _execute_optional_import(block_pyg_models)
+
+        import qmcpy
+
+        self.assertIs(namespace["mpmc_utils"], qmcpy.mpmc_utils)
+        self.assertEqual(
+            namespace["mpmc_utils"].__name__,
+            "qmcpy.discrete_distribution.mpmc.utils",
+        )
+        self.assertNotIn("utils", namespace)
+
+        with self.assertRaisesRegex(
+            ModuleNotFoundError, "MPMC_net.*torch_geometric"
+        ) as cm:
+            namespace["MPMC_net"]()
+        self.assertEqual(cm.exception.name, "torch_geometric")
+
+    def test_mpmc_placeholders_report_missing_torch(self):
+        def block_torch_utils(name, fromlist, level):
+            if (
+                level == 1
+                and name == "discrete_distribution.mpmc"
+                and "utils" in fromlist
+            ):
+                return "torch"
+            return None
+
+        namespace = _execute_optional_import(block_torch_utils)
+
+        with self.assertRaisesRegex(ModuleNotFoundError, "mpmc_utils.*torch") as cm:
+            namespace["mpmc_utils"].L2star
+        self.assertEqual(cm.exception.name, "torch")
+
+        with self.assertRaisesRegex(ModuleNotFoundError, "MPMC_net.*torch") as cm:
+            namespace["MPMC_net"]()
+        self.assertEqual(cm.exception.name, "torch")
+
+    def test_mpmc_placeholder_missing_torch_scatter(self):
+        try:
+            import torch  # noqa: F401
+        except ImportError:
+            self.skipTest("torch not available")
+
+        def block_torch_scatter(name, fromlist, level):
+            if level == 1 and name == "discrete_distribution.mpmc.models":
+                return "torch_scatter"
+            return None
+
+        namespace = _execute_optional_import(block_torch_scatter)
+
+        with self.assertRaisesRegex(
+            ModuleNotFoundError, "MPMC_net.*torch_scatter"
+        ) as cm:
+            namespace["MPMC_net"]()
+        self.assertEqual(cm.exception.name, "torch_scatter")
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/test_sr_unwrap_markdown.py b/test/test_sr_unwrap_markdown.py
new file mode 100644
index 000000000..4115d9ad2
--- /dev/null
+++ b/test/test_sr_unwrap_markdown.py
@@ -0,0 +1,99 @@
+import unittest
+
+from scripts.unwrap_markdown import unwrap_markdown_text
+
+
+class TestUnwrapMarkdown(unittest.TestCase):
+
+    def test_unwraps_list_item_continuations(self):
+        cases = [
+            (
+                "- unordered first\n  unordered second\n",
+                "- unordered first unordered second\n",
+            ),
+            (
+                "- [ ] task first\n  task second\n",
+                "- [ ] task first task second\n",
+            ),
+            (
+                "10. ordered first\n    ordered second\n",
+                "10. ordered first ordered second\n",
+            ),
+        ]
+        for source, expected in cases:
+            with self.subTest(source=source):
+                updated = unwrap_markdown_text(source)
+
+                self.assertEqual(updated, expected)
+                self.assertEqual(unwrap_markdown_text(updated), updated)
+
+    def test_unwraps_adjacent_and_nested_list_items_separately(self):
+        source = (
+            "- parent first\n"
+            "  parent second\n"
+            "  - child first\n"
+            "    child second\n"
+            "- sibling first\n"
+            "  sibling second\n"
+        )
+
+        self.assertEqual(
+            unwrap_markdown_text(source),
+            (
+                "- parent first parent second\n"
+                "  - child first child second\n"
+                "- sibling first sibling second\n"
+            ),
+        )
+
+    def test_preserves_list_item_blocks_and_explicit_hard_breaks(self):
+        source = (
+            "- first paragraph\n"
+            "  continuation\n"
+            "\n"
+            "  second paragraph\n"
+            "  continuation\n"
+            "\n"
+            "- item before code\n"
+            "      indented code\n"
+            "\n"
+            "- explicit hard break  \n"
+            "  remains separate\n"
+        )
+
+        self.assertEqual(
+            unwrap_markdown_text(source),
+            (
+                "- first paragraph continuation\n"
+                "\n"
+                "  second paragraph continuation\n"
+                "\n"
+                "- item before code\n"
+                "      indented code\n"
+                "\n"
+                "- explicit hard break  \n"
+                "  remains separate\n"
+            ),
+        )
+
+    def test_unwraps_ordinary_paragraphs(self):
+        self.assertEqual(
+            unwrap_markdown_text("first line\nsecond line\n"),
+            "first line second line\n",
+        )
+
+    def test_preserves_horizontal_rules(self):
+        for rule in ["- - -", "* * *", "_ _ _"]:
+            with self.subTest(rule=rule):
+                source = f"{rule}\nfollowing paragraph\n"
+
+                self.assertEqual(unwrap_markdown_text(source), source)
+
+    def test_preserves_indented_code_that_looks_like_a_list(self):
+        source = "    - code first\n      code second\n"
+
+        self.assertEqual(unwrap_markdown_text(source), source)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/test_tm_copulas.py b/test/test_tm_copulas.py
new file mode 100644
index 000000000..ed29969fb
--- /dev/null
+++ b/test/test_tm_copulas.py
@@ -0,0 +1,1271 @@
+import unittest
+import warnings
+
+import numpy as np
+import scipy.stats as stats
+
+from qmcpy import (
+    AbstractCopula,
+    ClaytonCopula,
+    DigitalNetB2,
+    FrankCopula,
+    GaussianCopula,
+    GumbelCopula,
+    StudentTCopula,
+)
+
+from qmcpy.true_measure.copula import (
+    AbstractCopula as ModuleAbstractCopula,
+    _apply_marginal_ppfs,
+    _build_marginal_range,
+    _clip_unit_interval,
+    _marginal_cdfs_and_logpdf,
+    _validate_correlation_matrix,
+    _validate_dimension,
+    _validate_marginals,
+)
+
+from qmcpy.util import DimensionError, MethodImplementationError, ParameterError
+
+
+class PPFOnlyMarginal:
+    def ppf(self, u):
+        return np.asarray(u, dtype=float)
+
+
+class NonCallablePPFMarginal:
+    ppf = 1.0
+
+
+class UnitPDFMarginal:
+    def ppf(self, u):
+        return np.asarray(u, dtype=float)
+
+    def cdf(self, x):
+        return np.asarray(x, dtype=float)
+
+    def pdf(self, x):
+        return np.ones_like(np.asarray(x, dtype=float))
+
+
+class CDFOnlyMarginal(PPFOnlyMarginal):
+    def cdf(self, x):
+        return np.asarray(x, dtype=float)
+
+
+class BadIntervalMarginal(PPFOnlyMarginal):
+    def interval(self, confidence):
+        raise ValueError("interval unavailable")
+
+
+class BadRangeMarginal:
+    def ppf(self, u):
+        raise ValueError("ppf unavailable")
+
+
+def _equicorrelation(d, rho):
+    corr = np.full((d, d), rho, dtype=float)
+    np.fill_diagonal(corr, 1.0)
+    return corr
+
+
+def _make_copula(copula_cls, dimension=2, marginals=None, correlation=None, seed=7):
+    if marginals is None:
+        marginals = [stats.norm()] * dimension
+    if correlation is None:
+        correlation = np.eye(dimension)
+
+    kwargs = {}
+    if copula_cls is StudentTCopula:
+        kwargs["df"] = 4
+    if copula_cls is ClaytonCopula:
+        kwargs["theta"] = 2.0
+    if copula_cls is FrankCopula:
+        kwargs["theta"] = 5.0
+    if copula_cls is GumbelCopula:
+        kwargs["theta"] = 2.0
+
+    common = {
+        "sampler": DigitalNetB2(dimension, seed=seed),
+        "marginals": marginals,
+        **kwargs,
+    }
+    if copula_cls in [ClaytonCopula, FrankCopula, GumbelCopula]:
+        return copula_cls(**common)
+    return copula_cls(correlation=correlation, **common)
+
+
+class TestAbstractCopulaAndHelpers(unittest.TestCase):
+
+    def test_abstract_copula_is_importable_from_public_module_path(self):
+        self.assertIs(ModuleAbstractCopula, AbstractCopula)
+
+    def test_public_api_imports_and_normal_usage(self):
+        for copula_cls in [
+            GaussianCopula,
+            StudentTCopula,
+            ClaytonCopula,
+            FrankCopula,
+            GumbelCopula,
+        ]:
+            with self.subTest(copula_cls=copula_cls.__name__):
+                self.assertTrue(issubclass(copula_cls, AbstractCopula))
+
+                tm = _make_copula(copula_cls)
+                x = tm(8)
+                x_gen = tm.gen_samples(8)
+                v = tm.gen_copula_samples(8)
+
+                self.assertEqual(x.shape, (8, 2))
+                self.assertEqual(x_gen.shape, (8, 2))
+                self.assertEqual(v.shape, (8, 2))
+                self.assertTrue(np.all(np.isfinite(x)))
+                self.assertTrue(np.all(np.isfinite(x_gen)))
+                self.assertTrue(np.all((0 <= v) & (v <= 1)))
+
+    def test_abstract_copula_rejects_unimplemented_transform(self):
+        tm = AbstractCopula(
+            DigitalNetB2(2, seed=101),
+            marginals=[stats.uniform(), stats.uniform()],
+        )
+
+        with self.assertRaises(MethodImplementationError):
+            tm.copula_transform(np.full((3, 2), 0.5))
+
+    def test_abstract_copula_rejects_invalid_sampler(self):
+        with self.assertRaisesRegex(ParameterError, "sampler"):
+            AbstractCopula(object(), marginals=[stats.uniform()])
+
+    def test_validate_marginals_error_branches(self):
+        with self.assertRaisesRegex(ParameterError, "marginals"):
+            _validate_marginals(None)
+
+        with self.assertRaisesRegex(ParameterError, "at least one"):
+            _validate_marginals([])
+
+        with self.assertRaisesRegex(ParameterError, "ppf"):
+            _validate_marginals([NonCallablePPFMarginal()])
+
+    def test_validate_dimension_error_branches(self):
+        with self.assertRaisesRegex(DimensionError, "integer dimension"):
+            _validate_dimension(object(), [stats.uniform()])
+
+        with self.assertRaisesRegex(DimensionError, "marginals"):
+            _validate_dimension(3, [stats.uniform(), stats.uniform()])
+
+    def test_apply_marginal_ppfs_clips_endpoints_and_checks_dimension(self):
+        transformed = _apply_marginal_ppfs(
+            np.array([[0.0, 1.0], [1.0, 0.0]]),
+            [stats.norm(), stats.norm()],
+        )
+
+        self.assertEqual(transformed.shape, (2, 2))
+        self.assertTrue(np.all(np.isfinite(transformed)))
+
+        with self.assertRaisesRegex(DimensionError, "marginals"):
+            _apply_marginal_ppfs(np.full((2, 3), 0.5), [stats.uniform(), stats.uniform()])
+
+    def test_marginal_range_falls_back_when_interval_or_ppf_fails(self):
+        ranges = _build_marginal_range([BadIntervalMarginal(), BadRangeMarginal()])
+
+        self.assertEqual(ranges.shape, (2, 2))
+        self.assertTrue(np.all(np.isfinite(ranges[0])))
+        np.testing.assert_allclose(ranges[1], [-np.inf, np.inf])
+
+    def test_marginal_cdfs_and_logpdf_pdf_branch_and_errors(self):
+        x = np.array([[0.25, 0.75], [0.4, 0.6]])
+        u, log_density = _marginal_cdfs_and_logpdf(
+            x,
+            [UnitPDFMarginal(), UnitPDFMarginal()],
+        )
+
+        np.testing.assert_allclose(u, x)
+        np.testing.assert_allclose(log_density, np.zeros(2))
+
+        with self.assertRaisesRegex(ParameterError, "cdf"):
+            _marginal_cdfs_and_logpdf(x, [PPFOnlyMarginal(), UnitPDFMarginal()])
+
+        with self.assertRaisesRegex(ParameterError, "pdf"):
+            _marginal_cdfs_and_logpdf(x, [CDFOnlyMarginal(), UnitPDFMarginal()])
+
+    def test_validate_correlation_matrix_rejects_nonfinite_values(self):
+        with self.assertRaisesRegex(ValueError, "finite"):
+            _validate_correlation_matrix([[1.0, np.nan], [np.nan, 1.0]], 2)
+
+    def test_clip_unit_interval_uses_machine_epsilon(self):
+        clipped = _clip_unit_interval(np.array([0.0, 0.5, 1.0]))
+        eps = np.finfo(float).eps
+
+        np.testing.assert_allclose(clipped, [eps, 0.5, 1.0 - eps])
+
+    def test_copula_transform_outputs_dependent_uniforms_in_unit_cube(self):
+        for copula_cls in [GaussianCopula, StudentTCopula, ClaytonCopula, GumbelCopula, FrankCopula]:
+            with self.subTest(copula_cls=copula_cls.__name__):
+                tm = _make_copula(copula_cls, dimension=3)
+                u = np.array(
+                    [
+                        [0.1, 0.3, 0.7],
+                        [0.5, 0.5, 0.5],
+                        [0.9, 0.8, 0.2],
+                    ]
+                )
+
+                v = tm.copula_transform(u)
+
+                self.assertEqual(v.shape, u.shape)
+                self.assertTrue(np.all(np.isfinite(v)))
+                self.assertTrue(np.all((0.0 <= v) & (v <= 1.0)))
+
+    def test_copula_sample_shapes_are_preserved(self):
+        for copula_cls, dimension in [
+            (GaussianCopula, 3),
+            (StudentTCopula, 3),
+            (ClaytonCopula, 3),
+            (FrankCopula, 3),
+            (GumbelCopula, 3),
+        ]:
+            with self.subTest(copula_cls=copula_cls.__name__, dimension=dimension):
+                tm = _make_copula(copula_cls, dimension=dimension, seed=9)
+
+                one = tm(1)
+                many = tm(8)
+                batched_transform = tm._transform(np.full((2, 3, dimension), 0.5))
+
+                self.assertEqual(one.shape, (1, dimension))
+                self.assertEqual(many.shape, (8, dimension))
+                self.assertEqual(batched_transform.shape, (2, 3, dimension))
+                self.assertTrue(np.all(np.isfinite(one)))
+                self.assertTrue(np.all(np.isfinite(many)))
+                self.assertTrue(np.all(np.isfinite(batched_transform)))
+
+
+class TestEllipticalCopulas(unittest.TestCase):
+
+    def test_output_shape_with_nonnormal_marginals(self):
+        tm = GaussianCopula(
+            sampler=DigitalNetB2(2, seed=7),
+            marginals=[stats.beta(a=2, b=5), stats.gamma(a=3, scale=2)],
+            correlation=[[1.0, 0.4], [0.4, 1.0]],
+        )
+
+        x = tm(16)
+
+        self.assertEqual(x.shape, (16, 2))
+
+    def test_finite_output_for_normal_marginals(self):
+        tm = GaussianCopula(
+            sampler=DigitalNetB2(2, seed=11),
+            marginals=[stats.norm(), stats.norm(loc=1.0, scale=2.0)],
+            correlation=[[1.0, -0.3], [-0.3, 1.0]],
+        )
+
+        x = tm(128)
+
+        self.assertTrue(np.all(np.isfinite(x)))
+
+    def test_return_weights_shape_when_marginal_densities_available(self):
+        tm = GaussianCopula(
+            sampler=DigitalNetB2(2, seed=12),
+            marginals=[stats.norm(), stats.gamma(a=2.0)],
+            correlation=[[1.0, 0.25], [0.25, 1.0]],
+        )
+
+        x, weights = tm(32, return_weights=True)
+
+        self.assertEqual(x.shape, (32, 2))
+        self.assertEqual(weights.shape, (32,))
+        self.assertTrue(np.all(np.isfinite(weights)))
+        self.assertTrue(np.all(weights > 0.0))
+
+    def test_identity_correlation_matches_independent_marginal_transforms(self):
+        marginals = [stats.norm(loc=-1.0, scale=2.0), stats.gamma(a=2.0, scale=3.0)]
+        tm = GaussianCopula(
+            sampler=DigitalNetB2(2, seed=13),
+            marginals=marginals,
+            correlation=np.eye(2),
+        )
+        u = np.array([[0.2, 0.7], [0.4, 0.8], [0.9, 0.1]])
+
+        x = tm._transform(u)
+        expected = np.column_stack(
+            [marginals[j].ppf(u[:, j]) for j in range(len(marginals))]
+        )
+
+        np.testing.assert_allclose(x, expected, rtol=1e-12, atol=1e-12)
+
+    def test_positive_correlation_produces_positive_dependence(self):
+        rho = 0.75
+        tm = GaussianCopula(
+            sampler=DigitalNetB2(2, seed=17),
+            marginals=[stats.norm(), stats.norm()],
+            correlation=[[1.0, rho], [rho, 1.0]],
+        )
+
+        x = tm(4096)
+        empirical_corr = np.corrcoef(x.T)[0, 1]
+
+        self.assertGreater(empirical_corr, 0.5)
+        self.assertLess(abs(empirical_corr - rho), 0.2)
+
+    def test_elliptical_copulas_support_general_dimensions(self):
+        for copula_cls in [GaussianCopula, StudentTCopula]:
+            for dimension in [1, 3, 5]:
+                with self.subTest(copula_cls=copula_cls.__name__, dimension=dimension):
+                    correlation = _equicorrelation(dimension, 0.25)
+                    tm = _make_copula(
+                        copula_cls,
+                        dimension=dimension,
+                        marginals=[stats.norm()] * dimension,
+                        correlation=correlation,
+                        seed=19,
+                    )
+
+                    x = tm(16)
+                    one = tm(1)
+
+                    self.assertEqual(x.shape, (16, dimension))
+                    self.assertEqual(one.shape, (1, dimension))
+                    self.assertTrue(np.all(np.isfinite(x)))
+                    self.assertTrue(np.all(np.isfinite(one)))
+
+    def test_elliptical_copulas_handle_valid_near_singular_correlation(self):
+        for copula_cls in [GaussianCopula, StudentTCopula]:
+            with self.subTest(copula_cls=copula_cls.__name__):
+                dimension = 5
+                tm = _make_copula(
+                    copula_cls,
+                    dimension=dimension,
+                    marginals=[stats.norm()] * dimension,
+                    correlation=_equicorrelation(dimension, 0.999),
+                    seed=20,
+                )
+
+                x = tm(32)
+
+                self.assertEqual(x.shape, (32, dimension))
+                self.assertTrue(np.all(np.isfinite(x)))
+
+    def test_elliptical_copulas_reject_singular_correlation(self):
+        for copula_cls in [GaussianCopula, StudentTCopula]:
+            with self.subTest(copula_cls=copula_cls.__name__):
+                with self.assertRaisesRegex(ValueError, "positive definite"):
+                    _make_copula(
+                        copula_cls,
+                        dimension=3,
+                        marginals=[stats.norm(), stats.norm(), stats.norm()],
+                        correlation=np.ones((3, 3)),
+                        seed=22,
+                    )
+
+    def test_distribution_dimension_matches_number_of_marginals(self):
+        for copula_cls in [GaussianCopula, StudentTCopula, ClaytonCopula, FrankCopula, GumbelCopula]:
+            with self.subTest(copula_cls=copula_cls.__name__):
+                tm = _make_copula(
+                    copula_cls,
+                    dimension=5,
+                    marginals=[
+                        stats.norm(),
+                        stats.beta(a=2, b=5),
+                        stats.gamma(a=3),
+                        stats.expon(),
+                        stats.lognorm(s=0.5),
+                    ],
+                    correlation=np.eye(5),
+                )
+
+                x = tm(32)
+
+                self.assertEqual(x.shape, (32, 5))
+                self.assertTrue(np.all(np.isfinite(x)))
+
+    def test_invalid_dimension_mismatches_raise(self):
+        for copula_cls in [GaussianCopula, StudentTCopula]:
+            with self.subTest(copula_cls=copula_cls.__name__):
+                with self.assertRaisesRegex(DimensionError, "marginals"):
+                    _make_copula(
+                        copula_cls,
+                        dimension=2,
+                        marginals=[stats.norm(), stats.norm(), stats.norm()],
+                        correlation=np.eye(2),
+                    )
+
+                with self.assertRaisesRegex(ValueError, "shape"):
+                    _make_copula(
+                        copula_cls,
+                        dimension=2,
+                        marginals=[stats.norm(), stats.norm()],
+                        correlation=np.eye(3),
+                    )
+
+                with self.assertRaisesRegex(ValueError, "square"):
+                    _make_copula(
+                        copula_cls,
+                        dimension=2,
+                        marginals=[stats.norm(), stats.norm()],
+                        correlation=[[1.0, 0.2, 0.3], [0.2, 1.0, 0.4]],
+                    )
+
+    def test_archimedean_dimension_mismatch_raises_dimension_error(self):
+        for copula_cls in [ClaytonCopula, FrankCopula, GumbelCopula]:
+            with self.subTest(copula_cls=copula_cls.__name__):
+                with self.assertRaisesRegex(DimensionError, "marginals"):
+                    _make_copula(
+                        copula_cls,
+                        dimension=2,
+                        marginals=[stats.norm(), stats.norm(), stats.norm()],
+                    )
+
+    def test_invalid_correlation_matrices_raise_value_error(self):
+        correlations = [
+            [[1.0, 0.2], [0.3, 1.0]],
+            [[1.0, 0.2], [0.2, 0.9]],
+            [[1.0, 1.2], [1.2, 1.0]],
+        ]
+        for copula_cls in [GaussianCopula, StudentTCopula]:
+            for correlation in correlations:
+                with self.subTest(copula_cls=copula_cls.__name__, correlation=correlation):
+                    with self.assertRaises(ValueError):
+                        _make_copula(
+                            copula_cls,
+                            dimension=2,
+                            marginals=[stats.norm(), stats.norm()],
+                            correlation=correlation,
+                        )
+
+    def test_marginal_length_mismatch_raises_dimension_error(self):
+        with self.assertRaisesRegex(DimensionError, "marginals"):
+            GaussianCopula(
+                sampler=DigitalNetB2(2, seed=21),
+                marginals=[stats.norm()],
+                correlation=np.eye(2),
+            )
+
+    def test_marginal_without_ppf_raises_clear_error(self):
+        class NoPPF:
+            pass
+
+        with self.assertRaisesRegex(ParameterError, "ppf"):
+            GaussianCopula(
+                sampler=DigitalNetB2(1, seed=23),
+                marginals=[NoPPF()],
+                correlation=[[1.0]],
+            )
+
+    def test_common_scipy_frozen_marginals_work(self):
+        for copula_cls in [GaussianCopula, StudentTCopula, ClaytonCopula, FrankCopula, GumbelCopula]:
+            with self.subTest(copula_cls=copula_cls.__name__):
+                tm = _make_copula(
+                    copula_cls,
+                    dimension=5,
+                    marginals=[
+                        stats.norm(),
+                        stats.beta(a=2, b=5),
+                        stats.gamma(a=3),
+                        stats.expon(),
+                        stats.lognorm(s=0.5),
+                    ],
+                    correlation=np.eye(5),
+                    seed=47,
+                )
+
+                x = tm(128)
+
+                self.assertEqual(x.shape, (128, 5))
+                self.assertTrue(np.all(np.isfinite(x)))
+
+    def test_endpoint_uniforms_are_clipped_to_finite_outputs(self):
+        for copula_cls in [GaussianCopula, StudentTCopula, ClaytonCopula, FrankCopula, GumbelCopula]:
+            with self.subTest(copula_cls=copula_cls.__name__):
+                tm = _make_copula(
+                    copula_cls,
+                    dimension=5,
+                    marginals=[
+                        stats.norm(),
+                        stats.beta(a=2, b=5),
+                        stats.gamma(a=3),
+                        stats.expon(),
+                        stats.lognorm(s=0.5),
+                    ],
+                    correlation=np.eye(5),
+                    seed=53,
+                )
+                u = np.array(
+                    [
+                        [0.0, 1.0, 0.0, 1.0, 0.5],
+                        [1.0, 0.0, 1.0, 0.0, 0.5],
+                    ]
+                )
+
+                x = tm._transform(u)
+
+                self.assertEqual(x.shape, (2, 5))
+                self.assertTrue(np.all(np.isfinite(x)))
+
+    def test_student_t_copula_output_shape_and_finite_values(self):
+        tm = StudentTCopula(
+            sampler=DigitalNetB2(2, seed=29),
+            marginals=[stats.norm(), stats.gamma(a=3.0, scale=2.0)],
+            correlation=[[1.0, 0.5], [0.5, 1.0]],
+            df=4,
+        )
+
+        x = tm(128)
+
+        self.assertEqual(x.shape, (128, 2))
+        self.assertTrue(np.all(np.isfinite(x)))
+
+    def test_student_t_copula_positive_correlation_produces_positive_dependence(self):
+        tm = StudentTCopula(
+            sampler=DigitalNetB2(2, seed=31),
+            marginals=[stats.norm(), stats.norm()],
+            correlation=[[1.0, 0.7], [0.7, 1.0]],
+            df=5,
+        )
+
+        x = tm(4096)
+        empirical_corr = np.corrcoef(x.T)[0, 1]
+
+        self.assertGreater(empirical_corr, 0.45)
+
+    def test_student_t_copula_has_stronger_joint_tail_than_gaussian_copula(self):
+        rho = 0.7
+        df = 4
+        n = 2**12
+        marginals = [stats.norm(), stats.norm()]
+        correlation = [[1.0, rho], [rho, 1.0]]
+
+        gaussian = GaussianCopula(
+            sampler=DigitalNetB2(2, seed=101),
+            marginals=marginals,
+            correlation=correlation,
+        )
+        student_t = StudentTCopula(
+            sampler=DigitalNetB2(2, seed=101),
+            marginals=marginals,
+            correlation=correlation,
+            df=df,
+        )
+
+        x_gaussian = gaussian(n)
+        x_student_t = student_t(n)
+        threshold = stats.norm.ppf(0.99)
+
+        def joint_tail_rate(x):
+            tail_0 = x[:, 0] > threshold
+            return np.mean(x[tail_0, 1] > threshold)
+
+        gaussian_tail = joint_tail_rate(x_gaussian)
+        student_t_tail = joint_tail_rate(x_student_t)
+
+        self.assertGreater(student_t_tail, gaussian_tail + 0.08)
+
+    def test_student_t_copula_return_weights_shape_when_density_available(self):
+        tm = StudentTCopula(
+            sampler=DigitalNetB2(2, seed=37),
+            marginals=[stats.norm(), stats.gamma(a=2.0)],
+            correlation=[[1.0, 0.3], [0.3, 1.0]],
+            df=6,
+        )
+
+        x, weights = tm(32, return_weights=True)
+
+        self.assertEqual(x.shape, (32, 2))
+        self.assertEqual(weights.shape, (32,))
+        self.assertTrue(np.all(np.isfinite(weights)))
+        self.assertTrue(np.all(weights > 0.0))
+
+    def test_student_t_copula_boundary_df_values_are_finite(self):
+        for df in [1.0, 100.0]:
+            with self.subTest(df=df):
+                dimension = 3
+                tm = StudentTCopula(
+                    sampler=DigitalNetB2(dimension, seed=39),
+                    marginals=[stats.norm()] * dimension,
+                    correlation=_equicorrelation(dimension, 0.4),
+                    df=df,
+                )
+
+                x = tm(128)
+
+                self.assertEqual(x.shape, (128, dimension))
+                self.assertTrue(np.all(np.isfinite(x)))
+
+    def test_student_t_copula_large_df_is_close_to_gaussian_copula(self):
+        rho = 0.6
+        correlation = [[1.0, rho], [rho, 1.0]]
+        marginals = [stats.norm(), stats.norm()]
+        gaussian = GaussianCopula(
+            sampler=DigitalNetB2(2, seed=40),
+            marginals=marginals,
+            correlation=correlation,
+        )
+        student_t = StudentTCopula(
+            sampler=DigitalNetB2(2, seed=40),
+            marginals=marginals,
+            correlation=correlation,
+            df=100,
+        )
+
+        x_gaussian = gaussian(4096)
+        x_student_t = student_t(4096)
+        corr_gaussian = np.corrcoef(x_gaussian.T)[0, 1]
+        corr_student_t = np.corrcoef(x_student_t.T)[0, 1]
+
+        self.assertLess(abs(corr_student_t - corr_gaussian), 0.02)
+
+    def test_student_t_copula_invalid_df_raises_parameter_error(self):
+        for df in [0, -1, np.inf, "not-a-number"]:
+            with self.subTest(df=df):
+                with self.assertRaisesRegex(ParameterError, "df"):
+                    StudentTCopula(
+                        sampler=DigitalNetB2(2, seed=41),
+                        marginals=[stats.norm(), stats.norm()],
+                        correlation=np.eye(2),
+                        df=df,
+                    )
+
+    def test_student_t_copula_marginal_without_ppf_raises_clear_error(self):
+        class NoPPF:
+            pass
+
+        with self.assertRaisesRegex(ParameterError, "ppf"):
+            StudentTCopula(
+                sampler=DigitalNetB2(1, seed=43),
+                marginals=[NoPPF()],
+                correlation=[[1.0]],
+                df=4,
+            )
+
+
+class TestArchimedeanCopulas(unittest.TestCase):
+
+    def test_clayton_copula_output_shape_and_finite_values(self):
+        tm = ClaytonCopula(
+            sampler=DigitalNetB2(2, seed=57),
+            marginals=[stats.norm(), stats.gamma(a=3.0, scale=2.0)],
+            theta=2.0,
+        )
+
+        x = tm(128)
+
+        self.assertEqual(x.shape, (128, 2))
+        self.assertTrue(np.all(np.isfinite(x)))
+
+    def test_clayton_copula_return_weights_shape_when_density_available(self):
+        tm = ClaytonCopula(
+            sampler=DigitalNetB2(3, seed=59),
+            marginals=[stats.norm(), stats.gamma(a=2.0), stats.expon()],
+            theta=1.5,
+        )
+
+        x, weights = tm(32, return_weights=True)
+
+        self.assertEqual(x.shape, (32, 3))
+        self.assertEqual(weights.shape, (32,))
+        self.assertTrue(np.all(np.isfinite(weights)))
+        self.assertTrue(np.all(weights > 0.0))
+
+    def test_clayton_copula_invalid_theta_raises_parameter_error(self):
+        for theta in [0, -1, np.inf, "not-a-number"]:
+            with self.subTest(theta=theta):
+                with self.assertRaisesRegex(ParameterError, "theta"):
+                    ClaytonCopula(
+                        sampler=DigitalNetB2(2, seed=61),
+                        marginals=[stats.norm(), stats.norm()],
+                        theta=theta,
+                    )
+
+    def test_clayton_copula_supports_general_dimension(self):
+        for dimension in [2, 3, 5]:
+            with self.subTest(dimension=dimension):
+                tm = ClaytonCopula(
+                    sampler=DigitalNetB2(dimension, seed=63),
+                    marginals=[stats.norm()] * dimension,
+                    theta=2.0,
+                )
+
+                x = tm(128)
+
+                self.assertEqual(x.shape, (128, dimension))
+                self.assertTrue(np.all(np.isfinite(x)))
+
+    def test_clayton_copula_marginal_without_ppf_raises_clear_error(self):
+        class NoPPF:
+            pass
+
+        with self.assertRaisesRegex(ParameterError, "ppf"):
+            ClaytonCopula(
+                sampler=DigitalNetB2(2, seed=67),
+                marginals=[stats.norm(), NoPPF()],
+                theta=2.0,
+            )
+
+    def test_clayton_copula_common_scipy_frozen_marginals_work(self):
+        for marginals in [
+            [stats.norm(), stats.beta(a=2, b=5)],
+            [stats.gamma(a=3), stats.expon()],
+            [stats.lognorm(s=0.5), stats.norm()],
+        ]:
+            with self.subTest(marginals=[type(m.dist).__name__ for m in marginals]):
+                tm = ClaytonCopula(
+                    sampler=DigitalNetB2(2, seed=69),
+                    marginals=marginals,
+                    theta=2.0,
+                )
+
+                x = tm(128)
+
+                self.assertEqual(x.shape, (128, 2))
+                self.assertTrue(np.all(np.isfinite(x)))
+
+    def test_clayton_copula_endpoint_uniforms_are_clipped_to_finite_outputs(self):
+        tm = ClaytonCopula(
+            sampler=DigitalNetB2(2, seed=70),
+            marginals=[stats.norm(), stats.lognorm(s=0.5)],
+            theta=2.0,
+        )
+        u = np.array([[0.0, 1.0], [1.0, 0.0]])
+
+        x = tm._transform(u)
+
+        self.assertEqual(x.shape, (2, 2))
+        self.assertTrue(np.all(np.isfinite(x)))
+
+    def test_clayton_copula_tiny_theta_is_near_independent(self):
+        for dimension in [2, 3, 5]:
+            with self.subTest(dimension=dimension):
+                marginals = [stats.uniform()] * dimension
+                tm = ClaytonCopula(
+                    sampler=DigitalNetB2(dimension, seed=70),
+                    marginals=marginals,
+                    theta=1e-8,
+                )
+                u = np.array(
+                    [
+                        [0.2, 0.7, 0.4, 0.6, 0.8],
+                        [0.4, 0.8, 0.9, 0.3, 0.2],
+                        [0.9, 0.1, 0.3, 0.7, 0.5],
+                    ]
+                )[:, :dimension]
+
+                x = tm._transform(u)
+
+                self.assertEqual(x.shape, (3, dimension))
+                self.assertTrue(np.all(np.isfinite(x)))
+                np.testing.assert_allclose(x, u, atol=5e-6)
+
+    def test_clayton_copula_large_theta_is_finite(self):
+        for dimension in [2, 3, 5]:
+            for theta in [20.0, 50.0]:
+                with self.subTest(dimension=dimension, theta=theta):
+                    tm = ClaytonCopula(
+                        sampler=DigitalNetB2(dimension, seed=70),
+                        marginals=[stats.norm()] * dimension,
+                        theta=theta,
+                    )
+
+                    x = tm(128)
+
+                    self.assertEqual(x.shape, (128, dimension))
+                    self.assertTrue(np.all(np.isfinite(x)))
+
+    def test_clayton_copula_positive_dependence_behavior(self):
+        tm = ClaytonCopula(
+            sampler=DigitalNetB2(2, seed=71),
+            marginals=[stats.uniform(), stats.uniform()],
+            theta=2.0,
+        )
+
+        x = tm(4096)
+        empirical_corr = np.corrcoef(x.T)[0, 1]
+
+        self.assertGreater(empirical_corr, 0.45)
+
+    def test_clayton_copula_has_stronger_lower_tail_than_gaussian_copula(self):
+        theta = 2.0
+        n = 2**12
+        marginals = [stats.uniform(), stats.uniform()]
+        # Clayton Kendall tau is theta/(theta+2); convert to Gaussian rho.
+        rho = np.sin(np.pi * (theta / (theta + 2.0)) / 2.0)
+
+        clayton = ClaytonCopula(
+            sampler=DigitalNetB2(2, seed=73),
+            marginals=marginals,
+            theta=theta,
+        )
+        gaussian = GaussianCopula(
+            sampler=DigitalNetB2(2, seed=73),
+            marginals=marginals,
+            correlation=[[1.0, rho], [rho, 1.0]],
+        )
+
+        x_clayton = clayton(n)
+        x_gaussian = gaussian(n)
+        threshold = 0.05
+
+        def lower_tail_rate(x):
+            tail_0 = x[:, 0] < threshold
+            return np.mean(x[tail_0, 1] < threshold)
+
+        clayton_tail = lower_tail_rate(x_clayton)
+        gaussian_tail = lower_tail_rate(x_gaussian)
+
+        self.assertGreater(clayton_tail, gaussian_tail + 0.2)
+
+    def test_frank_copula_output_shape_for_two_dimensions(self):
+        tm = FrankCopula(
+            sampler=DigitalNetB2(2, seed=75),
+            marginals=[stats.norm(), stats.gamma(a=3.0, scale=2.0)],
+            theta=5.0,
+        )
+
+        x = tm(128)
+
+        self.assertEqual(x.shape, (128, 2))
+        self.assertTrue(np.all(np.isfinite(x)))
+
+    def test_frank_copula_positive_theta_supports_higher_dimensions(self):
+        for dimension in [3, 5]:
+            with self.subTest(dimension=dimension):
+                tm = FrankCopula(
+                    sampler=DigitalNetB2(dimension, seed=76),
+                    marginals=[stats.norm()] * dimension,
+                    theta=5.0,
+                )
+
+                x = tm(128)
+
+                self.assertEqual(x.shape, (128, dimension))
+                self.assertTrue(np.all(np.isfinite(x)))
+
+    def test_frank_copula_return_weights_shape_when_density_available(self):
+        tm = FrankCopula(
+            sampler=DigitalNetB2(3, seed=77),
+            marginals=[stats.norm(), stats.gamma(a=2.0), stats.expon()],
+            theta=4.0,
+        )
+
+        x, weights = tm(32, return_weights=True)
+
+        self.assertEqual(x.shape, (32, 3))
+        self.assertEqual(weights.shape, (32,))
+        self.assertTrue(np.all(np.isfinite(weights)))
+        self.assertTrue(np.all(weights > 0.0))
+
+    def test_frank_copula_invalid_theta_raises_parameter_error(self):
+        for theta in [0, np.inf, -np.inf, "not-a-number"]:
+            with self.subTest(theta=theta):
+                with self.assertRaisesRegex(ParameterError, "theta"):
+                    FrankCopula(
+                        sampler=DigitalNetB2(2, seed=78),
+                        marginals=[stats.norm(), stats.norm()],
+                        theta=theta,
+                    )
+
+    def test_frank_copula_negative_theta_rejected_above_two_dimensions(self):
+        with self.assertRaisesRegex(ParameterError, "d=2"):
+            FrankCopula(
+                sampler=DigitalNetB2(3, seed=79),
+                marginals=[stats.norm(), stats.norm(), stats.norm()],
+                theta=-2.0,
+            )
+
+    def test_frank_copula_dimension_mismatch_raises_dimension_error(self):
+        with self.assertRaisesRegex(DimensionError, "marginals"):
+            FrankCopula(
+                sampler=DigitalNetB2(2, seed=80),
+                marginals=[stats.norm(), stats.norm(), stats.norm()],
+                theta=5.0,
+            )
+
+    def test_frank_copula_marginal_without_ppf_raises_clear_error(self):
+        class NoPPF:
+            pass
+
+        with self.assertRaisesRegex(ParameterError, "ppf"):
+            FrankCopula(
+                sampler=DigitalNetB2(2, seed=82),
+                marginals=[stats.norm(), NoPPF()],
+                theta=5.0,
+            )
+
+    def test_frank_copula_positive_dependence_behavior(self):
+        tm = FrankCopula(
+            sampler=DigitalNetB2(2, seed=84),
+            marginals=[stats.uniform(), stats.uniform()],
+            theta=6.0,
+        )
+
+        x = tm(4096)
+        empirical_corr = np.corrcoef(x.T)[0, 1]
+
+        self.assertGreater(empirical_corr, 0.45)
+
+    def test_frank_copula_tiny_theta_is_close_to_independence(self):
+        for theta, dimension in [(1e-8, 3), (-1e-8, 2)]:
+            with self.subTest(theta=theta, dimension=dimension):
+                marginals = [stats.uniform()] * dimension
+                tm = FrankCopula(
+                    sampler=DigitalNetB2(dimension, seed=86),
+                    marginals=marginals,
+                    theta=theta,
+                )
+                u = np.array(
+                    [
+                        [0.2, 0.7, 0.4, 0.6, 0.8],
+                        [0.4, 0.8, 0.9, 0.3, 0.2],
+                        [0.9, 0.1, 0.3, 0.7, 0.5],
+                    ]
+                )[:, :dimension]
+
+                x = tm._transform(u)
+
+                self.assertEqual(x.shape, (3, dimension))
+                self.assertTrue(np.all(np.isfinite(x)))
+                np.testing.assert_allclose(x, u, atol=5e-6)
+
+    def test_frank_copula_large_theta_is_finite(self):
+        for theta, dimension in [(50.0, 5), (-50.0, 2)]:
+            with self.subTest(theta=theta, dimension=dimension):
+                tm = FrankCopula(
+                    sampler=DigitalNetB2(dimension, seed=87),
+                    marginals=[stats.norm()] * dimension,
+                    theta=theta,
+                )
+
+                x = tm(128)
+
+                self.assertEqual(x.shape, (128, dimension))
+                self.assertTrue(np.all(np.isfinite(x)))
+
+    def test_frank_copula_negative_theta_produces_negative_dependence_in_2d(self):
+        tm = FrankCopula(
+            sampler=DigitalNetB2(2, seed=88),
+            marginals=[stats.uniform(), stats.uniform()],
+            theta=-6.0,
+        )
+
+        x = tm(4096)
+        empirical_corr = np.corrcoef(x.T)[0, 1]
+
+        self.assertLess(empirical_corr, -0.35)
+
+    def test_gumbel_copula_output_shape_and_finite_values(self):
+        tm = GumbelCopula(
+            sampler=DigitalNetB2(2, seed=79),
+            marginals=[stats.norm(), stats.gamma(a=3.0, scale=2.0)],
+            theta=2.0,
+        )
+
+        x = tm(128)
+
+        self.assertEqual(x.shape, (128, 2))
+        self.assertTrue(np.all(np.isfinite(x)))
+
+    def test_gumbel_copula_return_weights_shape_when_density_available(self):
+        tm = GumbelCopula(
+            sampler=DigitalNetB2(3, seed=81),
+            marginals=[stats.norm(), stats.gamma(a=2.0), stats.expon()],
+            theta=1.5,
+        )
+
+        x, weights = tm(32, return_weights=True)
+
+        self.assertEqual(x.shape, (32, 3))
+        self.assertEqual(weights.shape, (32,))
+        self.assertTrue(np.all(np.isfinite(weights)))
+        self.assertTrue(np.all(weights > 0.0))
+
+    def test_gumbel_copula_invalid_theta_raises_parameter_error(self):
+        for theta in [0, 0.5, -1, np.inf, "not-a-number"]:
+            with self.subTest(theta=theta):
+                with self.assertRaisesRegex(ParameterError, "theta"):
+                    GumbelCopula(
+                        sampler=DigitalNetB2(2, seed=83),
+                        marginals=[stats.norm(), stats.norm()],
+                        theta=theta,
+                    )
+
+    def test_gumbel_copula_theta_one_is_independent_marginal_transform(self):
+        marginals = [stats.norm(loc=-1.0, scale=2.0), stats.gamma(a=2.0, scale=3.0)]
+        tm = GumbelCopula(
+            sampler=DigitalNetB2(2, seed=85),
+            marginals=marginals,
+            theta=1.0,
+        )
+        u = np.array([[0.2, 0.7], [0.4, 0.8], [0.9, 0.1]])
+
+        x = tm._transform(u)
+        expected = np.column_stack(
+            [marginals[j].ppf(u[:, j]) for j in range(len(marginals))]
+        )
+
+        np.testing.assert_allclose(x, expected, rtol=1e-12, atol=1e-12)
+
+    def test_gumbel_copula_theta_close_to_one_is_near_independent(self):
+        for dimension in [2, 3, 5]:
+            with self.subTest(dimension=dimension):
+                marginals = [stats.uniform()] * dimension
+                tm = GumbelCopula(
+                    sampler=DigitalNetB2(dimension, seed=85),
+                    marginals=marginals,
+                    theta=1.000001,
+                )
+                u = np.array(
+                    [
+                        [0.2, 0.7, 0.4, 0.6, 0.8],
+                        [0.4, 0.8, 0.9, 0.3, 0.2],
+                        [0.9, 0.1, 0.3, 0.7, 0.5],
+                    ]
+                )[:, :dimension]
+
+                x = tm._transform(u)
+
+                self.assertEqual(x.shape, (3, dimension))
+                self.assertTrue(np.all(np.isfinite(x)))
+                np.testing.assert_allclose(x, u, atol=5e-5)
+
+    def test_gumbel_copula_large_theta_is_finite(self):
+        for dimension in [2, 3, 5]:
+            for theta in [20.0, 50.0]:
+                with self.subTest(dimension=dimension, theta=theta):
+                    tm = GumbelCopula(
+                        sampler=DigitalNetB2(dimension, seed=86),
+                        marginals=[stats.norm()] * dimension,
+                        theta=theta,
+                    )
+
+                    x = tm(128)
+
+                    self.assertEqual(x.shape, (128, dimension))
+                    self.assertTrue(np.all(np.isfinite(x)))
+
+    def test_gumbel_copula_supports_general_dimension(self):
+        for dimension in [2, 3, 5]:
+            with self.subTest(dimension=dimension):
+                tm = GumbelCopula(
+                    sampler=DigitalNetB2(dimension, seed=87),
+                    marginals=[stats.norm()] * dimension,
+                    theta=2.0,
+                )
+
+                x = tm(128)
+
+                self.assertEqual(x.shape, (128, dimension))
+                self.assertTrue(np.all(np.isfinite(x)))
+
+    def test_gumbel_copula_marginal_without_ppf_raises_clear_error(self):
+        class NoPPF:
+            pass
+
+        with self.assertRaisesRegex(ParameterError, "ppf"):
+            GumbelCopula(
+                sampler=DigitalNetB2(2, seed=89),
+                marginals=[stats.norm(), NoPPF()],
+                theta=2.0,
+            )
+
+    def test_gumbel_copula_common_scipy_frozen_marginals_work(self):
+        for marginals in [
+            [stats.norm(), stats.beta(a=2, b=5)],
+            [stats.gamma(a=3), stats.expon()],
+            [stats.lognorm(s=0.5), stats.norm()],
+        ]:
+            with self.subTest(marginals=[type(m.dist).__name__ for m in marginals]):
+                tm = GumbelCopula(
+                    sampler=DigitalNetB2(2, seed=91),
+                    marginals=marginals,
+                    theta=2.0,
+                )
+
+                x = tm(128)
+
+                self.assertEqual(x.shape, (128, 2))
+                self.assertTrue(np.all(np.isfinite(x)))
+
+    def test_gumbel_copula_endpoint_uniforms_are_clipped_to_finite_outputs(self):
+        tm = GumbelCopula(
+            sampler=DigitalNetB2(2, seed=93),
+            marginals=[stats.norm(), stats.lognorm(s=0.5)],
+            theta=2.0,
+        )
+        u = np.array([[0.0, 1.0], [1.0, 0.0]])
+
+        x = tm._transform(u)
+
+        self.assertEqual(x.shape, (2, 2))
+        self.assertTrue(np.all(np.isfinite(x)))
+
+    def test_gumbel_copula_positive_dependence_behavior(self):
+        tm = GumbelCopula(
+            sampler=DigitalNetB2(2, seed=95),
+            marginals=[stats.uniform(), stats.uniform()],
+            theta=2.0,
+        )
+
+        x = tm(4096)
+        empirical_corr = np.corrcoef(x.T)[0, 1]
+
+        self.assertGreater(empirical_corr, 0.45)
+
+    def test_gumbel_copula_has_stronger_upper_tail_than_gaussian_copula(self):
+        theta = 2.0
+        n = 2**12
+        marginals = [stats.uniform(), stats.uniform()]
+        # Gumbel Kendall tau is 1 - 1/theta; convert to Gaussian rho.
+        rho = np.sin(np.pi * (1.0 - 1.0 / theta) / 2.0)
+
+        gumbel = GumbelCopula(
+            sampler=DigitalNetB2(2, seed=97),
+            marginals=marginals,
+            theta=theta,
+        )
+        gaussian = GaussianCopula(
+            sampler=DigitalNetB2(2, seed=97),
+            marginals=marginals,
+            correlation=[[1.0, rho], [rho, 1.0]],
+        )
+
+        x_gumbel = gumbel(n)
+        x_gaussian = gaussian(n)
+        threshold = 0.95
+
+        def upper_tail_rate(x):
+            tail_0 = x[:, 0] > threshold
+            return np.mean(x[tail_0, 1] > threshold)
+
+        gumbel_tail = upper_tail_rate(x_gumbel)
+        gaussian_tail = upper_tail_rate(x_gaussian)
+
+        self.assertGreater(gumbel_tail, gaussian_tail + 0.15)
+
+
+class TestCopulaWeightsFallbackAndSpawn(unittest.TestCase):
+
+    def test_copula_weight_fallback_warns_once_when_density_methods_are_missing(self):
+        for copula_cls in [GaussianCopula, StudentTCopula, ClaytonCopula, GumbelCopula, FrankCopula]:
+            with self.subTest(copula_cls=copula_cls.__name__):
+                tm = _make_copula(
+                    copula_cls,
+                    dimension=2,
+                    marginals=[PPFOnlyMarginal(), PPFOnlyMarginal()],
+                )
+                x = np.full((4, 2), 0.5)
+                expected_message = getattr(
+                    tm,
+                    "_missing_weight_warning_message",
+                    f"{copula_cls.__name__} marginals must implement 'cdf' and "
+                    "'pdf' or 'logpdf' to compute density weights. "
+                    "Weights will be treated as 1.",
+                )
+
+                self.assertNotIn("_unit_weight_with_warning", copula_cls.__dict__)
+                self.assertIs(
+                    tm._unit_weight_with_warning.__func__,
+                    AbstractCopula._unit_weight_with_warning,
+                )
+
+                with self.assertWarns(UserWarning) as wcm:
+                    weights = tm._weight(x)
+
+                with warnings.catch_warnings(record=True) as caught:
+                    warnings.simplefilter("always")
+                    second_weights = tm._weight(x)
+
+                np.testing.assert_allclose(weights, np.ones(4))
+                np.testing.assert_allclose(second_weights, np.ones(4))
+                self.assertEqual(str(wcm.warning), expected_message)
+                self.assertEqual(caught, [])
+
+    def test_student_t_weight_falls_back_when_multivariate_t_is_unavailable(self):
+        tm = StudentTCopula(
+            DigitalNetB2(2, seed=115),
+            marginals=[stats.norm(), stats.norm()],
+            correlation=np.eye(2),
+            df=4,
+        )
+        tm._mvt_scipy = None
+
+        with self.assertWarnsRegex(UserWarning, "Weights will be treated as 1"):
+            weights = tm._weight(np.full((3, 2), 0.25))
+
+        np.testing.assert_allclose(weights, np.ones(3))
+
+    def test_gaussian_weight_uses_pdf_branch_when_logpdf_is_unavailable(self):
+        tm = GaussianCopula(
+            DigitalNetB2(2, seed=117),
+            marginals=[UnitPDFMarginal(), UnitPDFMarginal()],
+            correlation=[[1.0, 0.4], [0.4, 1.0]],
+        )
+
+        weights = tm._weight(np.array([[0.25, 0.5], [0.75, 0.5]]))
+
+        self.assertEqual(weights.shape, (2,))
+        self.assertTrue(np.all(np.isfinite(weights)))
+        self.assertTrue(np.all(weights > 0.0))
+
+    def test_gumbel_theta_one_weight_is_independent_marginal_density(self):
+        tm = GumbelCopula(
+            DigitalNetB2(2, seed=119),
+            marginals=[stats.gamma(a=2.0), stats.expon()],
+            theta=1.0,
+        )
+        x = np.array([[1.0, 0.5], [2.0, 1.5]])
+        expected = stats.gamma(a=2.0).pdf(x[:, 0]) * stats.expon().pdf(x[:, 1])
+
+        weights = tm._weight(x)
+
+        np.testing.assert_allclose(weights, expected)
+
+    def test_gen_copula_samples_composed_transform_branch(self):
+        inner = GaussianCopula(
+            DigitalNetB2(2, seed=121),
+            marginals=[stats.uniform(), stats.uniform()],
+            correlation=[[1.0, 0.3], [0.3, 1.0]],
+        )
+        outer = ClaytonCopula(inner, marginals=[stats.uniform(), stats.uniform()], theta=1.5)
+
+        v = outer.gen_copula_samples(n_min=4, n_max=8)
+
+        self.assertEqual(v.shape, (4, 2))
+        self.assertTrue(np.all(np.isfinite(v)))
+        self.assertTrue(np.all((0.0 <= v) & (v <= 1.0)))
+
+    def test_copula_spawn_same_dimension_and_reject_different_dimension(self):
+        for copula_cls in [GaussianCopula, StudentTCopula, ClaytonCopula, GumbelCopula, FrankCopula]:
+            with self.subTest(copula_cls=copula_cls.__name__):
+                tm = _make_copula(copula_cls, dimension=2)
+
+                spawned = tm.spawn(s=1, dimensions=[2])
+                self.assertEqual(len(spawned), 1)
+                self.assertIsInstance(spawned[0], copula_cls)
+                self.assertEqual(spawned[0](4).shape, (4, 2))
+
+                with self.assertRaises(DimensionError):
+                    tm._spawn(DigitalNetB2(3, seed=123), 3)
+
+    def test_frank_one_dimensional_weight_covers_zero_order_eulerian_term(self):
+        tm = FrankCopula(
+            DigitalNetB2(1, seed=125),
+            marginals=[UnitPDFMarginal()],
+            theta=3.0,
+        )
+
+        weights = tm._weight(np.array([[0.25], [0.75]]))
+
+        self.assertEqual(weights.shape, (2,))
+        self.assertTrue(np.all(np.isfinite(weights)))
+        self.assertTrue(np.all(weights > 0.0))
+
+    def test_frank_rejects_large_negative_theta_when_exponential_overflows(self):
+        with np.errstate(over="ignore"):
+            with self.assertRaisesRegex(ParameterError, "too close to 0 or too large"):
+                FrankCopula(
+                    DigitalNetB2(2, seed=127),
+                    marginals=[stats.uniform(), stats.uniform()],
+                    theta=-1000.0,
+                )
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/test_tm_product_measure.py b/test/test_tm_product_measure.py
new file mode 100644
index 000000000..470356860
--- /dev/null
+++ b/test/test_tm_product_measure.py
@@ -0,0 +1,271 @@
+import unittest
+
+import numpy as np
+import scipy.stats as stats
+
+from qmcpy import (
+    AcceptanceRejection,
+    DigitalNetB2,
+    DummySampler,
+    Gaussian,
+    GaussianCopula,
+    ProductMeasure,
+    SciPyWrapper,
+    Uniform,
+    ZeroInflatedExpUniform,
+)
+from qmcpy.util import DimensionError, ParameterError
+
+
+class TestProductMeasure(unittest.TestCase):
+
+    def test_product_measure_zero_inflated_with_scipy_uniform_shape(self):
+        n = 32
+        marginals = [
+            ZeroInflatedExpUniform(DummySampler(1), p_zero=0.4, lam=1.5),
+            SciPyWrapper(DummySampler(1), stats.uniform(loc=2.0, scale=3.0)),
+        ]
+        tm = ProductMeasure(sampler=DigitalNetB2(2, seed=23), marginals=marginals)
+
+        x = tm(n)
+
+        self.assertEqual(x.shape, (n, 2))
+        self.assertTrue(np.any(x[:, 0] == 0.0))
+        self.assertTrue(np.all((2.0 <= x[:, 1]) & (x[:, 1] <= 5.0)))
+
+    def test_product_measure_replication_shape(self):
+        n = 16
+        r = 3
+        marginals = [
+            ZeroInflatedExpUniform(DummySampler(1), p_zero=0.4, lam=1.5),
+            Uniform(DummySampler(1), lower_bound=2.0, upper_bound=5.0),
+        ]
+        tm = ProductMeasure(
+            sampler=DigitalNetB2(2, seed=23, replications=r),
+            marginals=marginals,
+        )
+
+        x = tm(n)
+
+        self.assertEqual(x.shape, (r, n, 2))
+
+    def test_product_measure_marginals_with_different_dimensions(self):
+        n = 32
+        marginals = [
+            Gaussian(
+                DummySampler(2),
+                mean=[1.0, -1.0],
+                covariance=[[2.0, 0.25], [0.25, 1.0]],
+            ),
+            ZeroInflatedExpUniform(DummySampler(1), p_zero=0.4, lam=1.5),
+        ]
+        tm = ProductMeasure(sampler=DigitalNetB2(3, seed=31), marginals=marginals)
+
+        x = tm(n)
+
+        self.assertEqual(tm.d, 3)
+        self.assertTrue(np.array_equal(tm.marginal_dimensions, np.array([2, 1])))
+        self.assertEqual(x.shape, (n, 3))
+        self.assertTrue(np.all(np.isfinite(x[:, :2])))
+        self.assertTrue(np.all(x[:, 2] >= 0.0))
+
+    def test_product_measure_block_split_range_and_weight_product(self):
+        n = 16
+        marginals = [
+            Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0),
+            Uniform(
+                DummySampler(2),
+                lower_bound=[20.0, 30.0],
+                upper_bound=[24.0, 36.0],
+            ),
+        ]
+        tm = ProductMeasure(sampler=DigitalNetB2(3, seed=41), marginals=marginals)
+
+        u = tm.discrete_distrib.gen_samples(n)
+        x = tm._transform(u)
+        x_call, jac = tm(n, return_weights=True)
+        expected = np.concatenate(
+            [
+                marginals[0]._jacobian_transform_r(u[..., :1], return_weights=False),
+                marginals[1]._jacobian_transform_r(u[..., 1:], return_weights=False),
+            ],
+            axis=-1,
+        )
+        expected_range = np.array([[10.0, 12.0], [20.0, 24.0], [30.0, 36.0]])
+
+        self.assertEqual(x.shape, (n, 3))
+        self.assertTrue(np.allclose(tm.range, expected_range))
+        self.assertTrue(np.allclose(x, expected))
+        self.assertTrue(np.all((10.0 <= x[:, 0]) & (x[:, 0] <= 12.0)))
+        self.assertTrue(np.all((20.0 <= x[:, 1]) & (x[:, 1] <= 24.0)))
+        self.assertTrue(np.all((30.0 <= x[:, 2]) & (x[:, 2] <= 36.0)))
+        self.assertTrue(np.allclose(tm._weight(x), 1.0 / (2.0 * 4.0 * 6.0)))
+        self.assertEqual(x_call.shape, (n, 3))
+        self.assertTrue(np.allclose(jac, 2.0 * 4.0 * 6.0))
+
+    def test_product_measure_invalid_inputs(self):
+        with self.assertRaisesRegex(ParameterError, "nonempty list of marginals"):
+            ProductMeasure(sampler=DigitalNetB2(1, seed=7), marginals=[])
+
+        with self.assertRaisesRegex(ParameterError, "marginal"):
+            ProductMeasure(sampler=DigitalNetB2(1, seed=7), marginals=[object()])
+
+        with self.assertRaisesRegex(ParameterError, "AbstractDiscreteDistribution"):
+            ProductMeasure(sampler=object(), marginals=[Uniform(DummySampler(1))])
+
+        marginals = [Uniform(DummySampler(1))]
+        with self.assertRaisesRegex(DimensionError, "sum of marginal dimensions"):
+            ProductMeasure(sampler=DigitalNetB2(2, seed=7), marginals=marginals)
+
+    def test_product_measure_rejects_non_dimension_preserving_marginal(self):
+        marginal = AcceptanceRejection(
+            DigitalNetB2(2, seed=7),
+            lambda x: np.ones(len(x)),
+            1.0,
+            1.0,
+        )
+
+        with self.assertRaisesRegex(DimensionError, "dimension-preserving"):
+            ProductMeasure(DigitalNetB2(2, seed=11), [marginal])
+
+    def test_product_measure_spawn_preserves_marginal_blocks_and_replaces_outer_sampler(self):
+        marginals = [
+            Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0),
+            Uniform(
+                DummySampler(2),
+                lower_bound=[20.0, 30.0],
+                upper_bound=[24.0, 36.0],
+            ),
+        ]
+        tm = ProductMeasure(sampler=DigitalNetB2(3, seed=41), marginals=marginals)
+
+        spawn = tm.spawn(s=1)[0]
+
+        self.assertIsInstance(spawn, ProductMeasure)
+        self.assertEqual(spawn.d, 3)
+        self.assertEqual(spawn.marginals, tm.marginals)
+        self.assertIsNot(spawn.discrete_distrib, tm.discrete_distrib)
+        self.assertTrue(np.array_equal(spawn.marginal_dimensions, np.array([1, 2])))
+
+        with self.assertRaises(DimensionError):
+            tm.spawn(s=1, dimensions=4)
+
+    def test_product_measure_does_not_use_marginal_dummy_sampler_values(self):
+        marginals = [
+            Uniform(DummySampler(1), lower_bound=0.0, upper_bound=2.0),
+            Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0),
+        ]
+
+        with self.assertRaisesRegex(ParameterError, "construction placeholder"):
+            marginals[0].discrete_distrib(4)
+
+        tm = ProductMeasure(sampler=DigitalNetB2(2, seed=19), marginals=marginals)
+        x = tm(8)
+
+        self.assertEqual(x.shape, (8, 2))
+        self.assertTrue(np.all((0.0 <= x[:, 0]) & (x[:, 0] <= 2.0)))
+        self.assertTrue(np.all((10.0 <= x[:, 1]) & (x[:, 1] <= 12.0)))
+
+    def test_product_measure_same_outer_seed_matches_different_outer_seed_changes(self):
+        marginals = [
+            Uniform(DummySampler(1), lower_bound=0.0, upper_bound=2.0),
+            Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0),
+        ]
+
+        first = ProductMeasure(sampler=DigitalNetB2(2, seed=101), marginals=marginals)(16)
+        same_outer = ProductMeasure(sampler=DigitalNetB2(2, seed=101), marginals=marginals)(16)
+        different_outer = ProductMeasure(sampler=DigitalNetB2(2, seed=102), marginals=marginals)(16)
+
+        self.assertTrue(np.array_equal(first, same_outer))
+        self.assertFalse(np.array_equal(first, different_outer))
+
+    def test_product_measure_replication_means_close_to_uniform_targets(self):
+        n = 1024
+        r = 4
+        marginals = [
+            Uniform(DummySampler(1), lower_bound=0.0, upper_bound=2.0),
+            Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0),
+        ]
+        tm = ProductMeasure(
+            sampler=DigitalNetB2(2, seed=101, replications=r),
+            marginals=marginals,
+        )
+
+        x = tm(n)
+        replication_means = x.mean(axis=1)
+
+        self.assertEqual(x.shape, (r, n, 2))
+        self.assertTrue(np.allclose(replication_means[:, 0], 1.0, atol=0.03))
+        self.assertTrue(np.allclose(replication_means[:, 1], 11.0, atol=0.03))
+
+    def test_product_measure_with_scipywrapper_beta_marginal(self):
+        n = 64
+        marginals = [
+            Uniform(DummySampler(1), lower_bound=-1.0, upper_bound=1.0),
+            SciPyWrapper(DummySampler(1), stats.beta(a=2.0, b=5.0)),
+        ]
+        tm = ProductMeasure(sampler=DigitalNetB2(2, seed=71), marginals=marginals)
+
+        x = tm(n)
+
+        self.assertEqual(x.shape, (n, 2))
+        self.assertTrue(np.all((-1.0 <= x[:, 0]) & (x[:, 0] <= 1.0)))
+        self.assertTrue(np.all((0.0 <= x[:, 1]) & (x[:, 1] <= 1.0)))
+
+    def test_product_measure_matches_equivalent_scipywrapper(self):
+        n = 128
+        seed = 55
+        scipy_marginals = [stats.norm(loc=0.0, scale=1.0), stats.gamma(a=2.0, scale=1.0)]
+        product_marginals = [
+            SciPyWrapper(DummySampler(1), scipy_marginals[0]),
+            SciPyWrapper(DummySampler(1), scipy_marginals[1]),
+        ]
+
+        product_samples = ProductMeasure(
+            sampler=DigitalNetB2(2, seed=seed),
+            marginals=product_marginals,
+        )(n)
+        scipy_samples = SciPyWrapper(DigitalNetB2(2, seed=seed), scipy_marginals)(n)
+
+        self.assertTrue(np.array_equal(product_samples, scipy_samples))
+
+    def test_product_measure_with_gaussian_copula_marginal(self):
+        n = 64
+        copula = GaussianCopula(
+            DummySampler(2),
+            marginals=[stats.beta(a=2.0, b=5.0), stats.gamma(a=3.0, scale=1.0)],
+            correlation=[[1.0, 0.5], [0.5, 1.0]],
+        )
+        marginals = [copula, Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0)]
+        tm = ProductMeasure(sampler=DigitalNetB2(3, seed=81), marginals=marginals)
+
+        x = tm(n)
+
+        self.assertEqual(x.shape, (n, 3))
+        self.assertTrue(np.all((0.0 <= x[:, 0]) & (x[:, 0] <= 1.0)))
+        self.assertTrue(np.all(x[:, 1] >= 0.0))
+        self.assertTrue(np.all((10.0 <= x[:, 2]) & (x[:, 2] <= 12.0)))
+
+    def test_product_measure_recursive_transform_sampling_supported_but_weights_restricted(self):
+        recursive_marginal = Uniform(
+            Uniform(DummySampler(1), lower_bound=0.0, upper_bound=1.0),
+            lower_bound=2.0,
+            upper_bound=4.0,
+        )
+        direct_marginal = Uniform(DummySampler(1), lower_bound=10.0, upper_bound=12.0)
+        tm = ProductMeasure(
+            sampler=DigitalNetB2(2, seed=91),
+            marginals=[recursive_marginal, direct_marginal],
+        )
+
+        x = tm(16)
+
+        self.assertEqual(x.shape, (16, 2))
+        self.assertTrue(np.all((2.0 <= x[:, 0]) & (x[:, 0] <= 4.0)))
+        self.assertTrue(np.all((10.0 <= x[:, 1]) & (x[:, 1] <= 12.0)))
+        with self.assertRaisesRegex(ParameterError, "direct marginal"):
+            tm(16, return_weights=True)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/test_tm_scipy_wrapper_custom.py b/test/test_tm_scipy_wrapper_custom.py
new file mode 100644
index 000000000..105b984e5
--- /dev/null
+++ b/test/test_tm_scipy_wrapper_custom.py
@@ -0,0 +1,292 @@
+import unittest
+import warnings
+
+import numpy as np
+import scipy.stats as stats
+
+from qmcpy import DigitalNetB2, SciPyWrapper, StudentT, ZeroInflatedExpUniform
+
+from qmcpy.true_measure.triangular import TriangularDistribution
+from qmcpy.util import DimensionError, ParameterError
+
+
+MISSING_PDF_WARNING = "no 'pdf' or 'logpdf'"
+
+
+def _missing_pdf_warnings(caught):
+    return [
+        warning
+        for warning in caught
+        if issubclass(warning.category, UserWarning)
+        and MISSING_PDF_WARNING in str(warning.message)
+    ]
+
+
+class TestSciPyWrapperCustom(unittest.TestCase):
+
+    def test_mvn_dependence_correlation_and_moment(self):
+        """
+        Check that passing a SciPy multivariate normal through SciPyWrapper
+        preserves correlation and the mixed moment E[X1 X2].
+        """
+        sampler = DigitalNetB2(2, seed=5)
+        rho_target = 0.7
+        cov = [[1.0, rho_target], [rho_target, 1.0]]
+        mvn = stats.multivariate_normal(mean=[0.0, 0.0], cov=cov)
+        tm_mvn = SciPyWrapper(sampler, scipy_distribs=mvn)
+
+        n = 4096
+        x = tm_mvn(n)
+
+        rho_hat = np.corrcoef(x.T)[0, 1]
+        est_moment = np.mean(x[:, 0] * x[:, 1])
+
+        self.assertTrue(np.isfinite(rho_hat))
+        self.assertTrue(np.isfinite(est_moment))
+
+        self.assertLess(abs(rho_hat - rho_target), 0.05)
+        self.assertLess(abs(est_moment - rho_target), 0.05)
+
+    def test_triangular_custom_marginal_range_and_shape(self):
+        """
+        Make sure our custom triangular marginal behaves sensibly:
+        samples stay in the right interval and the empirical mean is close
+        to the analytic mean.
+        """
+        tri = TriangularDistribution(c=0.3, loc=-1.0, scale=2.0)
+        tm = SciPyWrapper(DigitalNetB2(1, seed=11), scipy_distribs=tri)
+
+        n = 4096
+        x = tm(n).ravel()
+
+        self.assertGreaterEqual(x.min(), -1.1)
+        self.assertLessEqual(x.max(), 1.1)
+
+        a = -1.0
+        b = 1.0
+        m = -1.0 + 0.3 * 2.0
+        true_mean = (a + b + m) / 3.0
+        emp_mean = x.mean()
+        self.assertLess(abs(emp_mean - true_mean), 0.05)
+
+    def test_zero_inflated_zero_rate(self):
+        """
+        Check that the zero-inflated exponential distribution preserves the
+        specified probability mass at X = 0.
+        """
+        p_zero = 0.4
+        sampler = DigitalNetB2(1, seed=17)
+        tm = ZeroInflatedExpUniform(sampler, p_zero=p_zero, lam=1.5)
+
+        n = 4096
+        samples = tm(n)
+        x = samples.ravel()
+        zero_rate = np.mean(x == 0.0)
+
+        self.assertEqual(samples.shape, (n, 1))
+        self.assertLess(abs(zero_rate - p_zero), 0.05)
+
+    def test_zero_inflated_replications_shape(self):
+        tm = ZeroInflatedExpUniform(
+            DigitalNetB2(1, seed=17, replications=2),
+            p_zero=0.4,
+            lam=1.5,
+        )
+
+        x = tm(8)
+
+        self.assertEqual(x.shape, (2, 8, 1))
+        self.assertTrue(np.all(x >= 0.0))
+
+    def test_zero_inflated_rejects_invalid_p_zero(self):
+        for p_zero in [0.0, 1.0, -0.1, 1.1]:
+            with self.subTest(p_zero=p_zero):
+                with self.assertRaisesRegex(ParameterError, "p_zero must be in"):
+                    ZeroInflatedExpUniform(
+                        DigitalNetB2(1, seed=17),
+                        p_zero=p_zero,
+                        lam=1.5,
+                    )
+
+    def test_zero_inflated_rejects_nonpositive_lam(self):
+        for lam in [0.0, -1.0]:
+            with self.subTest(lam=lam):
+                with self.assertRaisesRegex(ParameterError, "lam must be positive"):
+                    ZeroInflatedExpUniform(
+                        DigitalNetB2(1, seed=17),
+                        p_zero=0.4,
+                        lam=lam,
+                    )
+
+    def test_zero_inflated_requires_one_dimensional_sampler(self):
+        with self.assertRaisesRegex(
+            DimensionError, "requires a one-dimensional sampler"
+        ):
+            ZeroInflatedExpUniform(
+                DigitalNetB2(2, seed=17),
+                p_zero=0.4,
+                lam=1.5,
+            )
+
+    def test_zero_inflated_inverse_transform_exact_values(self):
+        tm = ZeroInflatedExpUniform(
+            DigitalNetB2(1, seed=17),
+            p_zero=0.4,
+            lam=2.0,
+        )
+        u = np.array([[0.0], [0.2], [0.4], [0.7], [0.9]])
+
+        x = tm._transform(u)
+
+        self.assertEqual(x.shape, (5, 1))
+        self.assertTrue(np.array_equal(x[:3], np.zeros((3, 1))))
+        self.assertTrue(np.all(x[3:] > 0.0))
+
+        u_positive = u[3:, 0]
+        u_rescaled = (u_positive - 0.4) / 0.6
+        expected = -np.log1p(-u_rescaled) / 2.0
+        self.assertTrue(np.allclose(x[3:, 0], expected))
+
+    def test_zero_inflated_inverse_transform_all_zero_branch(self):
+        tm = ZeroInflatedExpUniform(
+            DigitalNetB2(1, seed=17),
+            p_zero=0.4,
+            lam=2.0,
+        )
+        u = np.array([[0.0], [0.1], [0.4]])
+
+        x = tm._transform(u)
+
+        self.assertEqual(x.shape, (3, 1))
+        self.assertTrue(np.array_equal(x, np.zeros((3, 1))))
+
+    def test_zero_inflated_inverse_transform_clips_one(self):
+        tm = ZeroInflatedExpUniform(
+            DigitalNetB2(1, seed=17),
+            p_zero=0.4,
+            lam=2.0,
+        )
+        u = np.array([[1.0]])
+
+        x = tm._transform(u)
+
+        self.assertEqual(x.shape, (1, 1))
+        self.assertTrue(np.isfinite(x).all())
+        self.assertGreater(x[0, 0], 0.0)
+
+    def test_zero_inflated_construction_does_not_warn_about_missing_pdf(self):
+        with warnings.catch_warnings(record=True) as caught:
+            warnings.simplefilter("always")
+            tm = ZeroInflatedExpUniform(
+                DigitalNetB2(1, seed=17),
+                p_zero=0.4,
+                lam=1.5,
+            )
+
+        self.assertEqual(tm.d, 1)
+        self.assertEqual(_missing_pdf_warnings(caught), [])
+
+    def test_zero_inflated_sampling_does_not_warn_about_missing_pdf(self):
+        tm = ZeroInflatedExpUniform(DigitalNetB2(1, seed=17), p_zero=0.4, lam=1.5)
+
+        with warnings.catch_warnings(record=True) as caught:
+            warnings.simplefilter("always")
+            x = tm(8)
+
+        self.assertEqual(x.shape, (8, 1))
+        self.assertEqual(_missing_pdf_warnings(caught), [])
+
+    def test_zero_inflated_return_weights_warns_once_for_missing_pdf(self):
+        tm = ZeroInflatedExpUniform(DigitalNetB2(1, seed=17), p_zero=0.4, lam=1.5)
+
+        with self.assertWarnsRegex(UserWarning, MISSING_PDF_WARNING):
+            x, jac = tm(8, return_weights=True)
+
+        self.assertEqual(x.shape, (8, 1))
+        self.assertEqual(jac.shape, (8,))
+        self.assertTrue(np.allclose(jac, 1.0))
+
+        with warnings.catch_warnings(record=True) as caught:
+            warnings.simplefilter("always")
+            x_second, jac_second = tm(8, return_weights=True)
+
+        self.assertEqual(x_second.shape, (8, 1))
+        self.assertTrue(np.allclose(jac_second, 1.0))
+        self.assertEqual(_missing_pdf_warnings(caught), [])
+
+    def test_zero_inflated_y_split_warns_and_uses_one_dimensional_interface(self):
+        with self.assertWarnsRegex(DeprecationWarning, "y_split"):
+            tm = ZeroInflatedExpUniform(
+                DigitalNetB2(1, seed=17),
+                p_zero=0.4,
+                lam=1.5,
+                y_split=0.5,
+            )
+
+        x = tm(4)
+
+        self.assertEqual(x.shape, (4, 1))
+        self.assertTrue(np.all(x >= 0.0))
+
+    def test_zero_inflated_y_split_preserves_deprecated_two_dimensional_usage(self):
+        with self.assertWarnsRegex(DeprecationWarning, "2D zero-inflated"):
+            tm = ZeroInflatedExpUniform(
+                DigitalNetB2(2, seed=17),
+                p_zero=0.4,
+                lam=1.5,
+                y_split=0.5,
+            )
+
+        x = tm(16)
+
+        self.assertEqual(x.shape, (16, 2))
+        self.assertTrue(np.all(x[:, 0] >= 0.0))
+        self.assertTrue(np.all((0.0 <= x[:, 1]) & (x[:, 1] <= 1.0)))
+        self.assertTrue(np.all(x[x[:, 0] == 0.0, 1] <= 0.5))
+        self.assertTrue(np.all(x[x[:, 0] > 0.0, 1] >= 0.5))
+
+    def test_zero_inflated_y_split_preserves_replicated_two_dimensional_usage(self):
+        with self.assertWarnsRegex(DeprecationWarning, "2D zero-inflated"):
+            tm = ZeroInflatedExpUniform(
+                DigitalNetB2(2, seed=17, replications=2),
+                p_zero=0.4,
+                lam=1.5,
+                y_split=0.5,
+            )
+
+        x = tm(16)
+
+        self.assertEqual(x.shape, (2, 16, 2))
+        self.assertTrue(np.all(x[..., 0] >= 0.0))
+        self.assertTrue(np.all((0.0 <= x[..., 1]) & (x[..., 1] <= 1.0)))
+        self.assertTrue(np.all(x[..., 1][x[..., 0] == 0.0] <= 0.5))
+        self.assertTrue(np.all(x[..., 1][x[..., 0] > 0.0] >= 0.5))
+
+    def test_student_t_marginals_shape(self):
+        tm = SciPyWrapper(
+            sampler=DigitalNetB2(2, seed=5),
+            scipy_distribs=stats.t(df=5),
+        )
+        x = tm(8)
+        self.assertEqual(x.shape, (8, 2))
+
+    def test_multivariate_student_t_joint_corr_and_cov(self):
+        if not hasattr(stats, "multivariate_t"):
+            self.skipTest("scipy.stats.multivariate_t not available in this SciPy version")
+
+        df = 5.0
+        rho = 0.8
+        loc = np.array([0.0, 0.0])
+        shape = np.array([[1.0, rho], [rho, 1.0]])
+
+        tm = StudentT(DigitalNetB2(2, seed=123), loc=loc, shape=shape, df=df)
+
+        n = 4096
+        x = tm(n)
+        emp_corr = np.corrcoef(x.T)[0, 1]
+
+        self.assertLess(abs(emp_corr - rho), 0.05)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/test_true_measures.py b/test/test_tm_true_measures.py
similarity index 100%
rename from test/test_true_measures.py
rename to test/test_tm_true_measures.py
diff --git a/test/test_unwrap_markdown.py b/test/test_unwrap_markdown.py
deleted file mode 100644
index e227b1807..000000000
--- a/test/test_unwrap_markdown.py
+++ /dev/null
@@ -1,89 +0,0 @@
-import pytest
-
-from scripts.unwrap_markdown import unwrap_markdown_text
-
-
-@pytest.mark.parametrize(
-    ("source", "expected"),
-    [
-        (
-            "- unordered first\n  unordered second\n",
-            "- unordered first unordered second\n",
-        ),
-        (
-            "- [ ] task first\n  task second\n",
-            "- [ ] task first task second\n",
-        ),
-        (
-            "10. ordered first\n    ordered second\n",
-            "10. ordered first ordered second\n",
-        ),
-    ],
-)
-def test_unwraps_list_item_continuations(source, expected):
-    updated = unwrap_markdown_text(source)
-
-    assert updated == expected
-    assert unwrap_markdown_text(updated) == updated
-
-
-def test_unwraps_adjacent_and_nested_list_items_separately():
-    source = (
-        "- parent first\n"
-        "  parent second\n"
-        "  - child first\n"
-        "    child second\n"
-        "- sibling first\n"
-        "  sibling second\n"
-    )
-
-    assert unwrap_markdown_text(source) == (
-        "- parent first parent second\n"
-        "  - child first child second\n"
-        "- sibling first sibling second\n"
-    )
-
-
-def test_preserves_list_item_blocks_and_explicit_hard_breaks():
-    source = (
-        "- first paragraph\n"
-        "  continuation\n"
-        "\n"
-        "  second paragraph\n"
-        "  continuation\n"
-        "\n"
-        "- item before code\n"
-        "      indented code\n"
-        "\n"
-        "- explicit hard break  \n"
-        "  remains separate\n"
-    )
-
-    assert unwrap_markdown_text(source) == (
-        "- first paragraph continuation\n"
-        "\n"
-        "  second paragraph continuation\n"
-        "\n"
-        "- item before code\n"
-        "      indented code\n"
-        "\n"
-        "- explicit hard break  \n"
-        "  remains separate\n"
-    )
-
-
-def test_unwraps_ordinary_paragraphs():
-    assert unwrap_markdown_text("first line\nsecond line\n") == "first line second line\n"
-
-
-@pytest.mark.parametrize("rule", ["- - -", "* * *", "_ _ _"])
-def test_preserves_horizontal_rules(rule):
-    source = f"{rule}\nfollowing paragraph\n"
-
-    assert unwrap_markdown_text(source) == source
-
-
-def test_preserves_indented_code_that_looks_like_a_list():
-    source = "    - code first\n      code second\n"
-
-    assert unwrap_markdown_text(source) == source
diff --git a/test/test_ut_plot_and_stop.py b/test/test_ut_plot_and_stop.py
new file mode 100644
index 000000000..97ee233b8
--- /dev/null
+++ b/test/test_ut_plot_and_stop.py
@@ -0,0 +1,167 @@
+import builtins
+import sys
+import types
+import unittest
+from unittest.mock import patch
+
+import numpy as np
+
+import qmcpy
+from qmcpy import plot_proj
+from qmcpy.util import stop_notebook
+
+
+class FakeAxes:
+    def __init__(self):
+        self.removed = False
+        self.calls = []
+
+    def remove(self):
+        self.removed = True
+
+    def set_xlim(self, *a, **k):
+        self.calls.append(("set_xlim", a))
+
+    def set_ylim(self, *a, **k):
+        self.calls.append(("set_ylim", a))
+
+    def set_xticks(self, *a, **k):
+        self.calls.append(("set_xticks", a))
+
+    def set_yticks(self, *a, **k):
+        self.calls.append(("set_yticks", a))
+
+    def set_aspect(self, *a, **k):
+        self.calls.append(("set_aspect", a))
+
+    def grid(self, *a, **k):
+        self.calls.append(("grid", a))
+
+    def tick_params(self, *a, **k):
+        self.calls.append(("tick_params", a))
+
+    def set_xlabel(self, *a, **k):
+        self.calls.append(("set_xlabel", a))
+
+    def set_ylabel(self, *a, **k):
+        self.calls.append(("set_ylabel", a))
+
+    def scatter(self, *a, **k):
+        self.calls.append(("scatter", a))
+
+
+class FakeFig:
+    def __init__(self):
+        self.tl = False
+
+    def tight_layout(self, *a, **k):
+        self.tl = True
+
+
+def make_fake_matplotlib(nrows, ncols):
+    plt = types.ModuleType("matplotlib.pyplot")
+    plt.style = types.SimpleNamespace()
+    plt.style.use = lambda *a, **k: None
+    plt.rcParams = {
+        "font.family": "sans-serif",
+        "axes.prop_cycle": types.SimpleNamespace(
+            by_key=lambda: {"color": ["k", "b", "r"]}
+        ),
+    }
+
+    def subplots(nrows=1, ncols=1, figsize=None, squeeze=False):
+        fig = FakeFig()
+        ax = np.empty((nrows, ncols), dtype=object)
+        for i in range(nrows):
+            for j in range(ncols):
+                ax[i, j] = FakeAxes()
+        return fig, ax
+
+    plt.subplots = subplots
+    plt.suptitle = lambda *a, **k: None
+    return plt
+
+
+class DummySampler(qmcpy.AbstractDiscreteDistribution):
+    def __init__(self, d=2):
+        super().__init__(dimension=d, replications=1, seed=1, d_limit=10, n_limit=100)
+
+    def _gen_samples(self, n_min, n_max, return_binary=False, warn=True):
+        n = n_max - n_min
+        return np.tile(np.arange(n)[:, None] / max(1, n - 1), (1, 1, self.d)).reshape(
+            self.replications, n, self.d
+        )
+
+    def __repr__(self):
+        return "DummySampler"
+
+
+class TestPlotProjAndStopNotebook(unittest.TestCase):
+
+    def test_plot_proj_with_fake_matplotlib_and_sampler(self):
+        # Inject fake matplotlib.pyplot
+        fake_plt = make_fake_matplotlib(1, 1)
+        # Create a proper matplotlib package module with colors submodule
+        fake_matplotlib = types.ModuleType("matplotlib")
+        fake_matplotlib.pyplot = fake_plt
+        fake_matplotlib.colors = types.SimpleNamespace()
+
+        with patch.dict(
+            sys.modules,
+            {"matplotlib.pyplot": fake_plt, "matplotlib": fake_matplotlib},
+        ):
+            sampler = DummySampler(d=3)
+            fig, ax = plot_proj(
+                sampler,
+                n=4,
+                d_horizontal=1,
+                d_vertical=2,
+                math_ind=True,
+                marker_size=1,
+                figfac=1,
+            )
+
+        self.assertIsInstance(fig, FakeFig)
+        self.assertIsInstance(ax, np.ndarray)
+        # At least one axes should have scatter calls or be removed
+        found = False
+        for a in ax.flatten():
+            if getattr(a, "removed", False) or any(c[0] == "scatter" for c in a.calls):
+                found = True
+                break
+        self.assertTrue(found)
+
+    def test_plot_proj_with_callable_sampler(self):
+        # sampler not instance of AbstractDiscreteDistribution -> uses t_i labels
+        fake_plt = make_fake_matplotlib(1, 1)
+        fake_matplotlib = types.ModuleType("matplotlib")
+        fake_matplotlib.pyplot = fake_plt
+        fake_matplotlib.colors = types.SimpleNamespace()
+
+        with patch.dict(
+            sys.modules,
+            {"matplotlib.pyplot": fake_plt, "matplotlib": fake_matplotlib},
+        ):
+            def sampler_callable(n):
+                return np.zeros((n, 1))
+
+            fig, ax = plot_proj(
+                sampler_callable, n=3, d_horizontal=0, d_vertical=0, math_ind=False
+            )
+
+        self.assertIsInstance(fig, FakeFig)
+
+    def test_stop_notebook_yes_and_no(self):
+        # When input is 'yes' nothing should happen
+        with patch.object(builtins, "input", lambda prompt="": "yes"):
+            # Should not raise
+            stop_notebook("prompt")
+
+        # When input is not 'yes' should exit
+        with patch.object(builtins, "input", lambda prompt="": "no"):
+            with self.assertRaises(SystemExit):
+                stop_notebook("prompt")
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/test_util.py b/test/test_ut_util.py
similarity index 100%
rename from test/test_util.py
rename to test/test_ut_util.py

From b7647cf92ae9c2148d99ab39f10d07dab874c2ed Mon Sep 17 00:00:00 2001
From: sou-cheng-choi <sou.cheng.choi@stanfordalumni.org>
Date: Mon, 31 Aug 2026 11:08:01 +0800
Subject: [PATCH 2/3] Respond to comments from gitHub-code-quality

---
 test/test_ig_financial_option_quick.py | 5 ++---
 test/test_ut_plot_and_stop.py          | 5 ++---
 2 files changed, 4 insertions(+), 6 deletions(-)

diff --git a/test/test_ig_financial_option_quick.py b/test/test_ig_financial_option_quick.py
index 921701712..b05b8708b 100644
--- a/test/test_ig_financial_option_quick.py
+++ b/test/test_ig_financial_option_quick.py
@@ -2,11 +2,10 @@
 
 import numpy as np
 
-from qmcpy import FinancialOption
-import qmcpy
+from qmcpy import AbstractDiscreteDistribution, FinancialOption
 
 
-class SmallSampler(qmcpy.AbstractDiscreteDistribution):
+class SmallSampler(AbstractDiscreteDistribution):
     def __init__(self, d=3):
         super().__init__(
             dimension=d, replications=1, seed=123, d_limit=100, n_limit=1024
diff --git a/test/test_ut_plot_and_stop.py b/test/test_ut_plot_and_stop.py
index 97ee233b8..4ce81545b 100644
--- a/test/test_ut_plot_and_stop.py
+++ b/test/test_ut_plot_and_stop.py
@@ -6,8 +6,7 @@
 
 import numpy as np
 
-import qmcpy
-from qmcpy import plot_proj
+from qmcpy import AbstractDiscreteDistribution, plot_proj
 from qmcpy.util import stop_notebook
 
 
@@ -82,7 +81,7 @@ def subplots(nrows=1, ncols=1, figsize=None, squeeze=False):
     return plt
 
 
-class DummySampler(qmcpy.AbstractDiscreteDistribution):
+class DummySampler(AbstractDiscreteDistribution):
     def __init__(self, d=2):
         super().__init__(dimension=d, replications=1, seed=1, d_limit=10, n_limit=100)
 

From f46f7f334b053e539e23273fcf830187d0cf2867 Mon Sep 17 00:00:00 2001
From: sou-cheng-choi <sou.cheng.choi@stanfordalumni.org>
Date: Mon, 31 Aug 2026 17:20:23 +0800
Subject: [PATCH 3/3] Fix Torch fallback test arguments

---
 test/test_ft_fast_transform_fallbacks.py | 38 +++++++++++++++++-------
 1 file changed, 27 insertions(+), 11 deletions(-)

diff --git a/test/test_ft_fast_transform_fallbacks.py b/test/test_ft_fast_transform_fallbacks.py
index d60266030..7bf9c3493 100644
--- a/test/test_ft_fast_transform_fallbacks.py
+++ b/test/test_ft_fast_transform_fallbacks.py
@@ -2,6 +2,11 @@
 
 import numpy as np
 
+try:
+    import torch
+except ImportError:
+    torch = None
+
 from qmcpy import (
     fftbr,
     fftbr_torch,
@@ -35,17 +40,28 @@ def test_non_torch_transforms_basic(self):
         omega2 = omega_fwht(3)
         self.assertEqual(omega2.shape[0], 2**3)
 
-    def test_torch_fallbacks_raise(self):
-        with self.assertRaises(Exception):
-            fftbr_torch()
-        with self.assertRaises(Exception):
-            ifftbr_torch()
-        with self.assertRaises(Exception):
-            fwht_torch()
-        with self.assertRaises(Exception):
-            omega_fftbr_torch()
-        with self.assertRaises(Exception):
-            omega_fwht_torch()
+    def test_torch_transforms_or_fallbacks(self):
+        if torch is None:
+            calls = (
+                (fftbr_torch, np.zeros(8, dtype=complex)),
+                (ifftbr_torch, np.zeros(8, dtype=complex)),
+                (fwht_torch, np.zeros(8)),
+                (omega_fftbr_torch, 3),
+                (omega_fwht_torch, 3),
+            )
+            for transform, argument in calls:
+                with self.subTest(transform=transform.__name__):
+                    with self.assertRaisesRegex(ModuleNotFoundError, "requires torch"):
+                        transform(argument)
+            return
+
+        complex_x = torch.zeros(8, dtype=torch.complex64)
+        real_x = torch.zeros(8)
+        self.assertEqual(fftbr_torch(complex_x).shape, complex_x.shape)
+        self.assertEqual(ifftbr_torch(complex_x).shape, complex_x.shape)
+        self.assertEqual(fwht_torch(real_x).shape, real_x.shape)
+        self.assertEqual(omega_fftbr_torch(3).shape[0], 2**3)
+        self.assertEqual(omega_fwht_torch(3).shape[0], 2**3)
 
 
 if __name__ == "__main__":