jeremymanning · Oct 29, 2026 · Oct 29, 2026
diff --git a/CLAUDE.md b/CLAUDE.md
@@ -12,7 +12,7 @@ Hierarchical Topographic Factor Analysis (HTFA) - A lightweight Python implement
 make install   # uv venv (.venv, Python 3.12) + pip install -e ".[dev,docs]"
 make test      # .venv/bin/pytest -q
 make lint      # ruff check, ruff format --check, mypy
-make docs      # sphinx-build -> docs/_build/html (needs docs/conf.py, added in a later task)
+make docs      # sphinx-build -> docs/_build/html
 
 # Run a specific test
 .venv/bin/pytest tests/test_package.py::test_version -xvs
@@ -23,12 +23,12 @@ Run `make lint` and `make test` before committing.
 
 ## Architecture
 
-Planned layout (modules land task by task during the rewrite; see docs/superpowers/plans/2026-10-23-htfa-rewrite.md):
+Package layout:
 
 ```
 htfa/
 ├── __init__.py       # public API, __version__
-├── _kernels.py       # numerical kernels (RBF factors, weights, objective)
+├── _kernels.py       # numerical kernels, factor mixtures, weights, objective
 ├── _engines.py       # compute engines (NumPy / torch)
 ├── _lazy.py          # lazy optional-dependency imports
 ├── tfa.py            # single-subject TFA
@@ -37,6 +37,8 @@ htfa/
 ├── fit.py            # main API entry point
 ├── results.py        # results container, plotting, export
 ├── connectivity.py   # factor-level connectivity
+├── groups.py         # nested templates and permutation comparisons
+├── stats.py          # reconstruction statistics and nominal parameter counts
 └── synthetic.py      # synthetic data generation
 tests/                # pytest suite (real implementations, no mocks)
 docs/                 # Sphinx documentation

diff --git a/README.md b/README.md
@@ -22,7 +22,7 @@ htfa is a standalone Python implementation of these algorithms. Its numerical va
 
 - No MPI or local compilation is needed when installing the prebuilt dependency wheels for a supported platform.
 - Analytic Jacobian for the center/width least-squares problem (`htfa/_kernels.py`).
-- A NumPy engine and a torch engine (CUDA, Apple Metal, CPU), chosen with the estimators' `backend` argument; `auto` uses a GPU only for large problems when torch is already installed.
+- A NumPy engine and a torch engine (CUDA, Apple Metal, CPU), chosen with the estimators' `backend` argument; `auto` uses a GPU only for large problems when torch is already installed. HTFA distributes independent subjects across visible CUDA GPUs; `gpu_devices` selects a subset.
 - Two solver modes: `mode='fast'` (the default) and `mode='exact'` (the reference for numerical equivalence and BrainIAK checks). The modes optimize the same objective but can reach different local minima; see [Performance](https://htfa.readthedocs.io/en/latest/performance.html) and [Validation](https://htfa.readthedocs.io/en/latest/validation.html).
 - Factor shapes: Gaussian (`kernel='rbf'`, the default), `'laplace'` and `'mexican_hat'`, and mixtures of several components per factor (`n_components`).
 - Nested subject groups with group-level templates and permutation comparisons.
@@ -32,6 +32,8 @@ htfa is a standalone Python implementation of these algorithms. Its numerical va
 
 ## Quick tour
 
+The [quick-start guide](https://htfa.readthedocs.io/en/latest/quickstart.html) walks through fitting, inspecting, and plotting a model, then loading your own images.
+
 Generate synthetic data for three subjects with known factor parameters, fit HTFA, and compare the recovered centers with the truth:
 
 ```python
@@ -66,6 +68,11 @@ For real data, pass a folder with one `.nii` / `.nii.gz` file per subject:
 res = htfa.fit("path/to/folder", K=50, random_state=0)
 ```
 
+Nested folders supply group names and hierarchical templates from their paths.
+Use `groups=False` to ignore folder-derived grouping, or pass explicit `groups`
+to define the experimental hierarchy yourself. Each file remains one subject;
+repeated runs from one person are not independent subjects.
+
 The inputs `htfa.fit` accepts, its parameters and the results object are described in the [documentation](https://htfa.readthedocs.io/en/latest/) and in [docs/API.md](https://github.com/ContextLab/htfa/blob/main/docs/API.md).
 
 ## Installation

diff --git a/docs/API.md b/docs/API.md
@@ -2,7 +2,7 @@
 
 The Sphinx documentation describes the API in these pages:
 
-- [Inputs](inputs.rst): how folders, BIDS lists, masks, arrays and images are read
+- [Inputs](inputs.rst): how recursive folders supply hierarchy names, and how BIDS lists, masks, arrays and images are read
 - [API reference](api.rst): `htfa.fit`, the estimators, results, connectivity and synthetic data, with model parameters
 - [Optional dependencies](optional_dependencies.rst): extras and installation on demand
 - [Performance](performance.rst): backends and benchmarks

diff --git a/docs/_static/custom.css b/docs/_static/custom.css
@@ -38,3 +38,29 @@ a:hover {
 pre {
     border-radius: 4px;
 }
+
+/* Tutorial cards link directly to the executed notebooks. */
+.tutorial-gallery {
+    display: grid;
+    grid-template-columns: repeat(auto-fit, minmax(min(100%, 240px), 1fr));
+    gap: 1rem;
+    margin: 1.5rem 0;
+}
+.tutorial-card {
+    border: 1px solid var(--color-background-border);
+    border-radius: 6px;
+    padding: 0.8rem;
+}
+.tutorial-card img {
+    width: 100%;
+    height: 170px;
+    object-fit: contain;
+    background: white;
+}
+.tutorial-card > p:nth-child(2) {
+    font-weight: 600;
+}
+.tutorial-card > p:last-child {
+    font-size: 0.9rem;
+    margin-bottom: 0;
+}
diff --git a/docs/_static/tutorials/01_quickstart.png b/docs/_static/tutorials/01_quickstart.png
diff --git a/docs/_static/tutorials/02_tfa_single_subject.png b/docs/_static/tutorials/02_tfa_single_subject.png
diff --git a/docs/_static/tutorials/03_htfa_multi_subject.png b/docs/_static/tutorials/03_htfa_multi_subject.png
diff --git a/docs/_static/tutorials/04_inputs_nifti_folders_nilearn_arrays.png b/docs/_static/tutorials/04_inputs_nifti_folders_nilearn_arrays.png
diff --git a/docs/_static/tutorials/05_2d_images.png b/docs/_static/tutorials/05_2d_images.png
diff --git a/docs/_static/tutorials/06_connectivity_and_dynamics.png b/docs/_static/tutorials/06_connectivity_and_dynamics.png
diff --git a/docs/_static/tutorials/09_word_decoding_mitchell2008.png b/docs/_static/tutorials/09_word_decoding_mitchell2008.png
diff --git a/docs/_static/tutorials/10_dynamic_isfc_animation.png b/docs/_static/tutorials/10_dynamic_isfc_animation.png
diff --git a/docs/_static/tutorials/11_plotting_and_appearance.png b/docs/_static/tutorials/11_plotting_and_appearance.png
diff --git a/docs/_static/tutorials/12_kernel_shapes.png b/docs/_static/tutorials/12_kernel_shapes.png
diff --git a/docs/_static/tutorials/13_group_hierarchy.png b/docs/_static/tutorials/13_group_hierarchy.png
diff --git a/docs/_static/tutorials/14_fit_statistics.png b/docs/_static/tutorials/14_fit_statistics.png
diff --git a/docs/_static/tutorials/15_real_fmri_kernel_comparison.png b/docs/_static/tutorials/15_real_fmri_kernel_comparison.png
diff --git a/docs/api.rst b/docs/api.rst
@@ -117,13 +117,16 @@ Common model parameters
      - Components per factor; more than 1 fits each factor as a mixture
    * - ``backend``
      - ``'auto'``
-     - ``'numpy'``, ``'cpu'``, ``'torch'``, ``'cuda'``, ``'mps'`` or ``'auto'``; see :doc:`performance`
+     - ``'numpy'``, ``'cpu'``, ``'torch'``, ``'cuda'``, ``'cuda:N'``, ``'mps'`` or ``'auto'``; see :doc:`performance`
    * - ``mode``
      - ``'fast'``
      - Step model of the center/width solver: ``'fast'`` keeps the Hessian where it is indefinite, ``'exact'`` (the reference, used by the equivalence tests and BrainIAK cross-checks) takes Gauss-Newton steps there; see :doc:`performance`
    * - ``n_jobs`` (HTFA)
      - ``None``
-     - Concurrent subject fits; NumPy fits use threads in one process
+     - Concurrent subject fits; NumPy fits use threads in one process. CUDA uses at most one worker per selected GPU; ``None`` uses all selected GPUs.
+   * - ``gpu_devices`` (HTFA)
+     - ``None``
+     - Visible CUDA indices, e.g. ``[0, 1]``; defaults to all visible GPUs when CUDA is selected. See :doc:`performance`.
    * - ``verbose``
      - ``False``
      - Print per-iteration progress

diff --git a/docs/changelog.rst b/docs/changelog.rst
@@ -25,7 +25,8 @@ Fitting
   Levenberg–Marquardt updates. ``'trf'`` and ``'dogbox'`` solvers use SciPy's
   ``least_squares`` with analytic, matrix-free Jacobians.
 - NumPy and torch backends support CPU, CUDA and Apple Metal. ``backend='auto'``
-  chooses a backend using the available devices and problem size.
+  chooses a backend using the available devices and problem size. HTFA can
+  distribute subjects across multiple CUDA GPUs.
 - ``mode='fast'`` is the default. ``mode='exact'`` is the numerical reference
   for equivalence and BrainIAK checks. Both modes optimize the same objective;
   their different step models can reach different local minima on real data.

diff --git a/docs/index.rst b/docs/index.rst
@@ -38,6 +38,7 @@ the factor count for your analysis:
    :caption: User guide
 
    installation
+   quickstart
    inputs
    optional_dependencies
    performance

diff --git a/docs/inputs.rst b/docs/inputs.rst
@@ -16,7 +16,7 @@ matrix and the ``(V, D)`` coordinates of its ``V`` locations.
    * - Path to a ``.nii`` / ``.nii.gz`` file
      - One subject
    * - Path to a folder
-     - One subject per NIfTI or image file (see below)
+     - One subject per NIfTI or image file, recursively; subfolders define groups
    * - nibabel image (for example from nilearn)
      - One subject
    * - ``(X, Y, Z, T)``, ``(X, Y, Z)`` or ``(H, W)`` numpy array, no ``coords``
@@ -26,7 +26,8 @@ matrix and the ``(V, D)`` coordinates of its ``V`` locations.
    * - ``.png``, ``.jpg``, ``.tif`` or other image file
      - One subject with one sample (TIFF pages become samples)
    * - List of any of the above
-     - One subject per element; ``coords`` may be a matching list
+     - Folders expand recursively; each file, array, or image is one subject.
+       ``coords`` may be a matching list.
 
 2D images need the ``images`` extra and nilearn masks need the ``nilearn`` extra;
 see :doc:`optional_dependencies`.
@@ -45,7 +46,10 @@ Folders of NIfTI files
 - Files are taken in natural filename order (``sub-2`` before ``sub-10``). The
   extension check ignores case. Dot files and other files (``.json`` sidecars,
   ``.tsv``) are skipped.
-- The folder is not searched recursively.
+- Subfolders are searched recursively. Hidden directories and directory symlinks
+  are skipped; empty subfolders contribute no groups. Each subject's
+  ``group_path`` records the folder names relative to the supplied root.
+  The root itself is the global template, not an extra group.
 - Voxel coordinates are world coordinates computed from each file's affine, so
   they are in mm for standard NIfTI files.
 - Subjects can have different numbers of voxels and of time points.
@@ -54,13 +58,77 @@ Folders of NIfTI files
 - ``stack=True`` treats every 3D volume in the folder as one sample of a single
   subject's time series. Stacked volumes must share shape and affine;
   explicit coordinate arrays must also have identical row order. Inputs with
-  and without world-coordinate affines cannot be stacked together.
+  and without world-coordinate affines cannot be stacked together. Stacking
+  combines every selected file, including nested folders, and clears group
+  metadata; it does not fit separate conditions.
+
+Folder-derived hierarchies
+--------------------------
+
+A nested folder tree supplies group templates automatically. This illustrative
+layout groups independent participants by acquisition site, then by condition:
+
+.. code-block:: text
+
+   participants/
+       site-a/
+           control/
+               subject-01.nii.gz
+           treatment/
+               subject-02.nii
+       site-b/
+           control/
+               subject-03.nii.gz
+
+.. code-block:: python
+
+   subjects = htfa.load("participants")
+   print(subjects[0].group_path)  # ('site-a', 'control')
+   res = htfa.fit("participants", K=5, random_state=0)
+   print(list(res.group_centers))
+   # ['site-a', 'site-b', ('site-a', 'control'),
+   #  ('site-a', 'treatment'), ('site-b', 'control')]
+
+Nodes are identified by the full folder path, so the two ``control`` folders
+are distinct groups. Files directly inside a group folder attach to that group;
+files directly inside the supplied root attach to the global template. This
+also supports trees with unequal depths without adding artificial groups.
+A single loaded subject still fits TFA, while retaining its ``group_path``.
+
+The inferred hierarchy uses ``template_update="fixed"`` by default. Set
+``template_update="empirical"`` to estimate hierarchy variances from the data;
+``"brainiak"`` does not support group templates and raises an error when groups
+are inferred. Pass explicit ``groups=`` to override the folder-derived groups,
+or ``groups=False`` to fit the same files with only a global template. Explicit
+label sequences retain the estimator's requirement to set a compatible
+``template_update``. To construct an arbitrary tree directly, pass a
+:class:`htfa.groups.GroupTree` made with ``GroupTree.from_paths`` as ``groups``.
+
+A list of folders includes each supplied folder's name as its first group
+level. A list of individual files has no inferred grouping. If different roots
+have the same name, explicit groups disambiguate them. Files are not recognized
+as independent participants from their names: repeated runs or conditions for
+the same people need an appropriate analysis design, not separate independent
+subject entries. ``subject_ids=`` can detect repeated identifiers in a grouped
+fit.
+
+Group fitting and plotting support unequal depths and repeated leaf names.
+Group permutation tests require equally deep paths and distinct labels for
+groups under different parents, because their level-wise shuffling treats
+subgroups as exchangeability units. To run such a test, select equally deep group paths
+and pass explicit distinct labels (for example ``site-a-control`` and
+``site-b-control``). The test rejects unsupported trees rather than assigning
+an invalid null distribution. Choosing labels does not establish independence
+or exchangeability of the underlying observations.
 
 BIDS datasets
 -------------
 
 The Brain Imaging Data Structure (BIDS) organizes imaging files in a nested
-folder tree. Pass a list of paths rather than the dataset root. For example,
+folder tree. Its structural folders (such as ``func`` and subject directories)
+are not automatically experimental groups. Pass a selected list of paths
+rather than the dataset root to avoid treating those folders as a hierarchy.
+For example,
 this pattern selects resting-state runs (``task-rest``); change the pattern
 for your task and preprocessed files:
 

diff --git a/docs/optional_dependencies.rst b/docs/optional_dependencies.rst
@@ -23,7 +23,7 @@ These extras are declared in ``pyproject.toml``:
    * - ``tutorials``
      - hypertools, sentence-transformers (only to re-execute the word-decoding tutorial)
    * - ``gpu``
-     - torch (``backend='torch'``, ``'cuda'`` or ``'mps'``)
+     - torch (``backend='torch'``, ``'cuda'``, ``'cuda:N'`` or ``'mps'``)
    * - ``dev``
      - pytest, pytest-cov, ruff, mypy, build, twine, plus plotting and image dependencies
    * - ``docs``
@@ -104,13 +104,13 @@ process. The function overrides the environment variable.
 torch
 -----
 
-``backend='torch'``, ``'cuda'`` or ``'mps'`` imports torch, installing
+``backend='torch'``, ``'cuda'``, ``'cuda:N'`` or ``'mps'`` imports torch, installing
 ``torch>=2.2`` when it is missing. On hosts other than macOS where ``nvidia-smi``
 is found and reports a CUDA driver, the installer picks the PyTorch wheel index
 for the newest supported CUDA version the driver can run (13.2, 13.0 or 12.6),
 or the CPU-only index when the driver is older than 12.6. Without an NVIDIA
 driver, and on macOS, it uses the default PyPI wheel.
-``backend='auto'`` never installs torch: it uses torch only if it is already
+Without an explicit ``gpu_devices`` selection, ``backend='auto'`` never installs torch: it uses torch only if it is already
 installed (see :doc:`performance`).
 
 Rendering without a display

diff --git a/docs/performance.rst b/docs/performance.rst
@@ -16,6 +16,41 @@ workload ``V * T * K`` reaches ``htfa._engines.GPU_MIN_WORK``. Otherwise it
 uses NumPy. This heuristic is not a guarantee that a GPU is faster; compare
 explicit backends on the workload of interest.
 
+Multiple GPUs
+-------------
+
+HTFA distributes independent subject fits across CUDA GPUs, with one serial
+queue per device. With ``backend='cuda'``, all visible devices are selected
+by default. ``n_jobs`` caps concurrent devices; ``None`` uses all selected
+devices and ``1`` uses the first. Select a subset explicitly:
+
+.. code-block:: python
+
+   import htfa
+
+   result = htfa.fit(
+       images, K=50, backend="cuda", gpu_devices=[0, 1], n_jobs=2
+   )
+   print(result.model.subject_engines_)
+
+Here ``images`` is a list of registered multi-subject images or input arrays.
+``gpu_devices`` also works with ``backend='auto'`` or ``'torch'`` and forces
+CUDA and raises an error when CUDA or a requested index is unavailable;
+there is no silent CPU fallback for an explicit device selection. It cannot
+be combined with ``backend='numpy'``, ``'cpu'``, ``'mps'``, or ``'cuda:N'``.
+Without that argument, ``'auto'`` uses the largest subsampled subject workload
+to decide whether to distribute HTFA subjects across CUDA devices.
+``backend='cuda:1'`` pins TFA or HTFA to one visible GPU. Indices follow
+PyTorch's ``CUDA_VISIBLE_DEVICES`` numbering; see `CUDA semantics
+<https://docs.pytorch.org/docs/stable/notes/cuda.html>`_.
+
+Subject order and seeds are preserved when work completes out of order.
+Different GPU models can produce floating-point differences; parallelism is
+not a guarantee of identical local minima or a speedup. The global template
+update remains on the CPU. Each subject's working arrays must fit on its
+assigned device: this does not shard a single-subject TFA fit across GPUs.
+MPS uses one device and remains serial across subjects.
+
 Solver modes
 ------------