disco-space a0eaa21cc0a6, brain.toml

README.md

@@ -8,13 +8,13 @@ sdk_version: 6.28.0
 app_file: app.py
 pinned: false
 startup_duration_timeout: 1h
-license: mit
+license: agpl-3.0
 short_description: A real brain replayed from packs, to an 84-region connectome
 ---
 
 # The brain Space on ZeroGPU
 
-The same application as `rfick/disco-zero` (dmrai-lab/disco-space), serving `space/brain.toml` (`DISCO_CONFIG`):
+The same application as `rfick/disco-replay` (dmrai-lab/disco-space), serving `space/brain.toml` (`DISCO_CONFIG`):
 a real brain as a compositional replay phantom. The input is the multi-tissue CSD of one scan, MASiVar `sub-cIs1`
 (OpenNeuro ds003416, Cai et al. MRM 2021, CC0): its white-matter FOD field, its WM / GM / CSF fractions and an
 84-region parcellation, made once offline and published as the dataset `SubstrateCommons/masivar-brain`. Each run
@@ -27,7 +27,7 @@ https://github.com/dmrai-lab/disco-space (`requirements-zero.txt`, `app.py`, `sp
 **The scanner.** The menu holds the ideal scanner and two catalogued machines, the Siemens Prisma 3 T and Terra 7 T,
 each played at every voxel as it delivers there with the head centre at isocentre (voxels binned into encoding
 classes at 1 % of b). The Hyperfine Swoop 64 mT is not on this menu, and the page says why under it: its own gradient
-is exact in dmipy-sim's closed form, but it bins this head into 9,752 classes at 2-3 minutes each
+is exact in dmipy-replay's closed form, but it bins this head into 9,752 classes at 2-3 minutes each
 (`docs/scanner.md`). The DiSCo Spaces play it exactly.
 
 **Where the work runs.** The packs' pose responses (the only physics that depends on the tissue and the field) are
@@ -49,7 +49,7 @@ reservation is killed. Only the device part of a run holds the GPU, with the pac
 cached (part of the reservation shown, which falls once a run has cached them); the page's files and figures after
 it do not.
 
-Not re-measured at the current pins (dmipy-sim 1b79d2e; the tables below name their own). Measured on an L40S with the BATMAN development fixture (96 x 96 x 60 at 2.5 mm, 90,205 brain voxels, 36,605 in the
+Not re-measured on the dmipy 3.0 line (the tables below name their own pins). Measured on an L40S with the BATMAN development fixture (96 x 96 x 60 at 2.5 mm, 90,205 brain voxels, 36,605 in the
 WM stop mask; dmipy-sim 7f6f1fa, dmipy-fit 0c7dde8 with the single-tissue reconstruction, dmipy-tract 7da22c3;
 `tools/measure_brain.py`, 2026-09-30), steady state, seconds:
 

app.py

@@ -4,14 +4,15 @@ the source's share a run needs that the page has not cached, a brain's pack resp
 as the request arrives); the cache lookup before it and the page's drawing and files after it stay in this process
 (disco-space#7). The configuration is ``DISCO_CONFIG`` (``space/config.toml``, the DiSCo Space, unless the
 Space sets another: ``brain.toml`` for the brain Space). The compute backend is torch (``DISCO_BACKEND=torch``), the
-layout's tiles are rebuilt on the device inside every call (``resident = false``: nothing survives between calls),
+layout's packed words are streamed to the device a chunk at a time inside every call (``resident = false``: nothing
+survives between calls),
 deterministic algorithms are on and TF32 off inside the call. The source's data come from the Hub at the config's
 revision, downloaded once when the container starts (the parent process loads it before the page serves; the forked
 GPU worker inherits it). ``spaces`` must be imported before torch."""
 import os
 
 os.environ.setdefault("DISCO_BACKEND", "torch")
-os.environ.setdefault("DISCO_RESIDENT", "0")             # the pool drops the device between calls: tiles rebuilt per call
+os.environ.setdefault("DISCO_RESIDENT", "0")             # the pool drops the device between calls: words streamed per call
 os.environ.setdefault("JAX_PLATFORMS", "cpu")            # jax is a library dependency here; it never sees the GPU
 
 from space.limits import cap_threads, cap_affinity      # noqa: E402  (no numpy behind it)
@@ -47,7 +48,7 @@ def probe():
     for the deployment check, `/probe`)."""
     import glob, json, time
     import numpy as np
-    st = A._load(); src = st["source"]; where = getattr(getattr(src, "moments", None), "path", "")
+    st = A._load(); src = st["source"]; where = getattr(getattr(src, "layout", None), "path", "")
     out = dict(torch=torch.__version__, cuda=torch.cuda.is_available(), device=torch.cuda.get_device_name(0) if torch.cuda.is_available() else None,
                gpu_memory_gb=round(torch.cuda.get_device_properties(0).total_memory / 1e9, 1) if torch.cuda.is_available() else None,
                layout=where, data=sorted(os.path.basename(p) for p in glob.glob(os.path.join(where, "*")))[:12], backend=os.environ.get("DISCO_BACKEND"),
@@ -59,24 +60,21 @@ def probe():
         out["host_ram_total_gb"] = round(int(mem["MemTotal"].split()[0]) / 1e6, 1); out["host_ram_available_gb"] = round(int(mem["MemAvailable"].split()[0]) / 1e6, 1)
     except Exception as e:
         out["meminfo"] = repr(e)[:100]
-    files = sorted(glob.glob(os.path.join(where, "m_*.npy")))
+    files = sorted(glob.glob(os.path.join(where, "p_*.npy")))
     if files:
         t0 = time.perf_counter(); m = np.load(files[0], mmap_mode="r"); chunk = np.array(m[:100000]); dt = time.perf_counter() - t0
         out["read_mb"] = round(chunk.nbytes / 1e6, 1); out["read_mb_per_s"] = round(chunk.nbytes / 1e6 / dt, 1)
-        t0 = time.perf_counter(); torch.as_tensor(chunk, device="cuda"); torch.cuda.synchronize(); out["upload_mb_per_s"] = round(chunk.nbytes / 1e6 / (time.perf_counter() - t0), 1)
-    # the moment image inside this worker: the first compiled call, a second one (the kernel alone), the eager kernel
+        t0 = time.perf_counter(); torch.as_tensor(chunk.view(np.int32), device="cuda"); torch.cuda.synchronize(); out["upload_mb_per_s"] = round(chunk.nbytes / 1e6 / (time.perf_counter() - t0), 1)
+    # the time-segment image inside this worker, streamed and (when it fits) resident: first and second call
     try:
-        from dmipy_sim.replay import shape_moments as SM
-        kernel = dict(SM._TORCH_KERNEL)
-        b = np.full(184, 1e9); u = np.tile([0.6, 0.0, 0.8], (184, 1)); name = "d10.2-D16.7"
-        try:
-            t0 = time.perf_counter(); src.moments.image(name, b, u, backend="torch", resident=True); torch.cuda.synchronize(); out["image_first_compiled_s"] = round(time.perf_counter() - t0, 2)
-            t0 = time.perf_counter(); src.moments.image(name, b, u, backend="torch", resident=True); torch.cuda.synchronize(); out["image_second_compiled_s"] = round(time.perf_counter() - t0, 2)
-            SM._TORCH_KERNEL["fn"] = SM._tile_sums_torch
-            t0 = time.perf_counter(); src.moments.image(name, b, u, backend="torch", resident=True); torch.cuda.synchronize(); out["image_eager_s"] = round(time.perf_counter() - t0, 2)
-            t0 = time.perf_counter(); src.moments.image(name, b, u, backend="torch", resident=False); torch.cuda.synchronize(); out["image_eager_nonresident_s"] = round(time.perf_counter() - t0, 2)
-        finally:
-            SM._TORCH_KERNEL.update(kernel)                  # the probe leaves the kernel as it found it
+        name = "d10.2-D16.7"
+        seq = src.sequence(name, np.full(184, 1000.0), np.tile([0.6, 0.0, 0.8], (184, 1)))
+        for resident in (False, True):
+            for k in ("first", "second"):
+                t0 = time.perf_counter(); src.layout.image(seq, partition=name, backend="torch", resident=resident); torch.cuda.synchronize()
+                out[f"image_{'resident' if resident else 'streamed'}_{k}_s"] = round(time.perf_counter() - t0, 2)
+        out["gpu_peak_gb"] = round(torch.cuda.max_memory_allocated() / 1e9, 2)
+        src.layout.release()
     except Exception as e:
         out["image_probe_error"] = repr(e)[:300]
     return json.dumps(out)

data/SOURCE.md

@@ -2,7 +2,7 @@
 
 From the Diffusion-Simulated Connectivity (DiSCo) dataset, Rafael-Patino J., Girard G., Truffet R., Pizzolato M.,
 Caruyer E., Thiran J.-P., "The Diffusion-Simulated Connectivity (DiSCo) dataset", Data in Brief 38 (2021) 107429,
-https://doi.org/10.1016/j.dib.2021.107429; data at https://doi.org/10.5281/zenodo.4634234. Licence CC BY 4.0.
+https://doi.org/10.1016/j.dib.2021.107429; data on Mendeley Data, https://doi.org/10.17632/fgf86jdfg6.1 (files checked by sha256 against the published release). Licence CC BY 4.0.
 
 | file | what |
 |---|---|
@@ -11,6 +11,6 @@ https://doi.org/10.1016/j.dib.2021.107429; data at https://doi.org/10.5281/zenod
 | `DiSCo_gradients.bvals`, `DiSCo_gradients_dipy.bvecs` | the dataset's 364-measurement gradient table (b in s/mm²), the "DiSCo 364" acquisition |
 | `DiSCo_Connectivity_Matrix_Strands_Count.txt` | ground truth: strands per region pair |
 | `DiSCo_Connectivity_Matrix_Cross-Sectional_Area.txt` | ground truth: total cross-sectional area per region pair |
-| `DiSCo_Strands_Trajectories.tck`, `DiSCo_Strands_Diameters.txt` | the 12,196 strands' centerlines (in units of the 25 µm voxel, i.e. the image grid) and inner diameters (mm): the ground-truth view |
+| `DiSCo_Strands_Diameters.txt` (the centrelines, `DiSCo_Strands_Trajectories.tck`, are read from `SubstrateCommons/disco-replay` `substrate/` at the pinned revision) | the 12,196 strands' centerlines (in units of the 25 µm voxel, i.e. the image grid) and inner diameters (mm): the ground-truth view |
 
 Unmodified copies of the dataset's DiSCo1 files, redistributed here under CC BY 4.0 so the Space is self-contained.

data/classes/qte-now-80.npy

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:6310e918172eb0086c9c31158307d8e4cfeca8c25b1a9fbab0a25e6838c71706
+size 24128

data/protocols/c2.bvals

@@ -0,0 +1 @@
+0.0 0.0 0.0 0.0 0.0 0.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0 6000.0

data/protocols/c2.bvecs

@@ -0,0 +1,3 @@
+0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 -0.32949755 0.33414940 -0.60771156 -0.82127342 0.01160086 -0.12962089 0.26498506 0.56912603 -0.40095140 0.11311349 -0.25030891 0.77356999 -0.52375456 0.94083501 -0.63920498 -0.67014535 -0.45412061 0.61922918 0.79604789 -0.78743880 -0.18826079 -0.94768325 -0.58407055 -0.02913436 -0.62545224 -0.33733820 -0.21403085 0.00421902 0.58345560 -0.35220267 0.45295742 0.91782422 0.94162595 0.28168108 0.79855761 0.99971654 0.48125673 0.36140462 0.24089400 0.26950471 -0.82478877 0.13564362 0.29332779 0.56938601 -0.08884351 -0.33494071 -0.85970235 0.73020584 0.92961722 -0.03034955 -0.60692305 0.80574615 0.08571867 -0.94986286 -0.03689412 -0.79313489 0.57162376 0.36836640 -0.49539615 0.74900912
+0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.62521613 0.94159591 -0.79319112 -0.55746095 0.64167929 -0.22505973 0.77269471 -0.07315939 0.80871935 0.37214629 -0.52159459 -0.63362406 0.30920269 -0.33860613 -0.01697085 0.68218268 -0.25517415 0.28358495 -0.51642751 0.22226796 0.96109710 0.09345700 -0.46706432 -0.91920621 0.52925384 -0.89663480 0.37609875 0.07555217 0.60185062 -0.73383382 -0.38312209 -0.21553952 0.10663176 -0.64667073 0.02326119 -0.01496704 0.81874023 0.15874085 -0.16398178 -0.95125392 -0.13803286 0.94259977 -0.84100598 -0.78380323 0.84364697 0.05931043 0.41440579 -0.31387397 0.35178020 -0.74332437 -0.70211555 0.37244222 -0.45885689 -0.22994878 -0.99750236 -0.42974786 -0.61279834 0.51158361 0.86507917 0.62230435
+0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.70748580 0.04172908 0.03917277 0.12143827 0.76688533 0.96568450 0.57682390 0.81898917 0.43036146 0.92125592 0.81564976 0.01048864 0.79377255 0.01324307 0.76884913 0.29245853 0.85361621 0.73221227 0.31561113 0.57492355 0.20211443 0.30522492 0.66386182 0.39269726 0.57332353 0.28679084 0.90152123 0.99713292 0.54530302 0.58089690 0.80501369 0.33337877 0.31932779 0.70885311 0.60146876 0.01851587 0.31313957 0.91879706 0.95659817 0.14994393 0.54833422 0.30513354 0.45460724 0.24789529 0.52949658 0.94037067 0.29862973 0.60686289 0.10982952 0.66824231 0.37239518 0.46049966 0.88436569 0.21185870 0.06023173 0.43157134 0.54564134 0.77626568 0.07887130 0.22742611

data/protocols/clinical.bvals

@@ -0,0 +1 @@
+0.0 0.0 0.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0

data/protocols/clinical.bvecs

@@ -0,0 +1,3 @@
+0.00000000 0.00000000 0.00000000 -0.21204524 -0.95684375 0.79698041 -0.41182571 -0.81011573 -0.11514776 -0.92695104 0.43532954 0.92492393 -0.42978830 -0.64150046 -0.39834740 -0.23078110 -0.72213785 -0.55205845 0.72558771 0.38087223 -0.63961904 0.96273433 0.35854872 0.01680837 0.04499735 0.64600654 0.69874109 -0.28489889 0.47941432 0.02769143 0.18831204 0.17906768 0.77676570
+0.00000000 0.00000000 0.00000000 -0.75054853 0.24143549 0.60309286 -0.90060083 0.12191004 0.46881409 -0.24577886 -0.12018932 -0.25417356 0.08442250 -0.30479278 0.91680538 -0.37028342 0.63418578 0.48863467 -0.66199141 0.35783255 -0.64027306 0.19605979 -0.79351726 -0.96157594 0.97637349 -0.44154424 0.48464294 0.81674577 0.81785461 0.02205862 -0.52460313 0.73817677 0.05109827
+0.00000000 0.00000000 0.00000000 0.62587036 0.16173726 0.03318481 0.13898823 0.57345483 0.87575929 0.28346871 0.89221226 0.28268628 0.89897434 0.70397338 0.02805805 0.89979457 0.27627039 0.67562388 0.18785595 0.85257974 0.42536795 0.18628785 0.49170436 0.27402405 0.21135290 0.62266703 0.52619589 0.50175589 0.31823852 0.99937310 0.83025908 0.65040743 0.62771332

data/protocols/ogse.bvals

@@ -0,0 +1 @@
+0.0 0.0 0.0 0.0 0.0 0.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0

data/protocols/ogse.bvecs

@@ -0,0 +1,3 @@
+0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 -0.32949755 0.33414940 -0.60771156 -0.82127342 0.01160086 -0.12962089 0.26498506 0.56912603 -0.40095140 0.11311349 -0.25030891 0.77356999 -0.52375456 0.94083501 -0.63920498 -0.67014535 -0.45412061 0.61922918 0.79604789 -0.78743880 -0.18826079 -0.94768325 -0.58407055 -0.02913436 -0.62545224 -0.33733820 -0.21403085 0.00421902 0.58345560 -0.35220267 0.45295742 0.91782422 0.94162595 0.28168108 0.79855761 0.99971654 0.48125673 0.36140462 0.24089400 0.26950471 -0.82478877 0.13564362 0.29332779 0.56938601 -0.08884351 -0.33494071 -0.85970235 0.73020584 0.92961722 -0.03034955 -0.60692305 0.80574615 0.08571867 -0.94986286 -0.03689412 -0.79313489 0.57162376 0.36836640 -0.49539615 0.74900912
+0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.62521613 0.94159591 -0.79319112 -0.55746095 0.64167929 -0.22505973 0.77269471 -0.07315939 0.80871935 0.37214629 -0.52159459 -0.63362406 0.30920269 -0.33860613 -0.01697085 0.68218268 -0.25517415 0.28358495 -0.51642751 0.22226796 0.96109710 0.09345700 -0.46706432 -0.91920621 0.52925384 -0.89663480 0.37609875 0.07555217 0.60185062 -0.73383382 -0.38312209 -0.21553952 0.10663176 -0.64667073 0.02326119 -0.01496704 0.81874023 0.15874085 -0.16398178 -0.95125392 -0.13803286 0.94259977 -0.84100598 -0.78380323 0.84364697 0.05931043 0.41440579 -0.31387397 0.35178020 -0.74332437 -0.70211555 0.37244222 -0.45885689 -0.22994878 -0.99750236 -0.42974786 -0.61279834 0.51158361 0.86507917 0.62230435
+0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.70748580 0.04172908 0.03917277 0.12143827 0.76688533 0.96568450 0.57682390 0.81898917 0.43036146 0.92125592 0.81564976 0.01048864 0.79377255 0.01324307 0.76884913 0.29245853 0.85361621 0.73221227 0.31561113 0.57492355 0.20211443 0.30522492 0.66386182 0.39269726 0.57332353 0.28679084 0.90152123 0.99713292 0.54530302 0.58089690 0.80501369 0.33337877 0.31932779 0.70885311 0.60146876 0.01851587 0.31313957 0.91879706 0.95659817 0.14994393 0.54833422 0.30513354 0.45460724 0.24789529 0.52949658 0.94037067 0.29862973 0.60686289 0.10982952 0.66824231 0.37239518 0.46049966 0.88436569 0.21185870 0.06023173 0.43157134 0.54564134 0.77626568 0.07887130 0.22742611

data/protocols/research.bvals

@@ -0,0 +1 @@
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 1000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 2000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0 3000.0

data/protocols/research.bvecs

@@ -0,0 +1,3 @@
+0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.59576860 0.93816635 -0.98905967 0.62451169 0.60683766 0.46290349 0.52148610 0.67795329 0.57111780 -0.17229806 -0.60572887 0.30809534 0.09534500 -0.77316638 0.09213804 -0.93519076 -0.55432683 0.81317227 0.87438895 -0.32853522 -0.56029699 0.97598393 -0.23544716 0.73124810 -0.00744073 -0.91070654 -0.16119654 0.09450063 0.39933516 0.57567786 -0.50070270 -0.83540161 0.98983934 -0.60281077 0.26703580 0.14275761 0.76160413 0.22068804 0.46004974 -0.16201257 -0.76341115 -0.91346498 0.19530566 -0.83054568 -0.07599847 -0.32379331 0.20157094 -0.60199623 -0.96376113 -0.43532233 -0.37974241 -0.67767550 0.33952724 0.44205976 -0.70961734 0.80187651 0.74817156 0.36770201 -0.86826836 -0.41544466 0.88698753 -0.58844532 0.67213557 -0.48615765 0.05142893 -0.36298452 0.90965478 -0.05848180 0.42943351 -0.78648992 -0.65587245 -0.39711485 -0.74887479 0.36635901 -0.05166445 -0.35496186 0.04103047 -0.12970923 0.06803766 -0.79174685 -0.18015876 0.82118333 -0.24712053 0.33019542 0.62372888 0.18643620 -0.40296738 0.19068873 0.92264726 -0.12540350 -0.33939275 0.16937678 0.89443500 0.34132607 -0.11320959 -0.81549102 0.64271660 -0.56923596 -0.76110906 -0.18178130 -0.28013937 -0.91360323 0.62079660 -0.46983928 0.38745233 -0.38416406 0.86822315 0.38120650 0.97531028 0.14051826 0.09884027 0.91267749 0.12291292 0.09417199 0.03506780 0.75051710 -0.63253219 -0.93656587 -0.74230411 -0.17556086 -0.90079207 0.56642162 -0.07126760 -0.38884372 -0.81568460 -0.07870329 0.77216408 0.34610728 -0.14906497 0.49412932 -0.41089422 -0.70913609 0.73202853 0.04978052 0.25106660 0.14193245 -0.75510751 0.66482405 -0.54595570 -0.64863191 0.86777392 -0.61767795 0.38730006 0.53672683 -0.96571430 -0.08597049 -0.19297778 -0.01677332 -0.42439546 0.58300144 0.95287282 0.78830792 -0.15891744 0.44046530 0.43636594 -0.12828867 0.62440566 -0.62827991 0.05265537 0.30402380 0.72392095 0.89347943 -0.22828902 -0.35081051 -0.86751866 0.30659065 -0.53668306 0.50998312 0.30891840 0.56165136 0.19374169 -0.58996776 -0.53446692 -0.79846068 -0.35390888 -0.97788739 0.75071745 -0.99547768 0.84048957 -0.31808885 0.94484484 -0.84507777 -0.27683998 0.88646476 -0.53324067 -0.41203317 -0.89513761 0.95935617 0.49252736 0.02044253 -0.75474992 -0.27425308 -0.60748804 -0.14677795 0.74445785 0.65184168 0.07207788 0.13286927 -0.77145863 0.12854416 0.84895556 0.38505005 -0.46946387 -0.32491845 0.56555779 0.24922406 0.10881076 0.60604096 0.99616741 -0.67208291 0.73654120 0.00020812 -0.46032509 -0.37727868 -0.61581083 0.35902896 0.36400089 -0.77937346 0.11745982 -0.65124185 0.34467616 -0.59758113 0.51582304 0.71983068 0.52860193 -0.89761540 -0.08842218 0.24092043 -0.82485704 0.33566474 0.53588812 0.86327729 0.39890986 -0.29144269 0.27586801 -0.20370148 -0.90224938 0.24743918 0.68573160 0.92545984 0.79372111 0.00141729 -0.47546969 -0.52426327 -0.43561588 0.43586461 -0.96073423 -0.73996389 -0.35118151 -0.19461791 0.59799415 0.16707138 0.27516140 -0.49623143 -0.64003755 0.84463988 -0.02563607 -0.16777918 0.75162129 -0.24697176 0.06620066 -0.71414354 -0.10757112 -0.15534669 -0.87945643 0.53522975 -0.96915272 -0.01847400 -0.98015195 0.72625438
+0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 -0.34328993 -0.02894166 -0.07537980 -0.77822261 0.16522689 -0.09717810 0.72628440 -0.09066748 -0.57192510 -0.42570880 -0.06645117 0.95007460 -0.54106777 0.62788656 0.60148275 -0.34745037 0.23970275 -0.22071348 0.45331689 -0.92347344 -0.82719563 0.16614636 -0.86119434 0.32668414 -0.83868609 0.03894678 0.12418008 0.15327271 -0.76705287 0.53968238 0.69338938 0.29984505 -0.11459020 -0.33206066 0.83397117 0.95050176 -0.46440436 0.39494465 0.39803948 0.97133778 -0.39143607 -0.22386741 -0.90417263 -0.48421087 0.90476887 0.33809415 -0.97354662 -0.58847011 0.18757384 0.53221513 -0.23412573 -0.68869468 -0.32800806 0.86362704 0.58171516 -0.56702420 0.53449434 -0.54890950 0.43490556 0.90722617 0.26762841 0.77273332 0.71376276 -0.79077834 0.99864256 -0.69682954 -0.30236529 0.38688077 -0.87328527 -0.61725755 0.42880755 0.03881415 0.13169272 0.15159132 -0.96429629 0.86683029 0.78388462 -0.13470650 -0.30775216 -0.13527179 0.59493900 0.07050457 0.76497525 0.64939030 -0.70278082 -0.72691934 -0.49299344 -0.08314931 -0.38487194 -0.66612704 -0.76522024 0.28827414 -0.44369383 0.07146179 0.94669326 0.42041533 0.31113019 -0.73033240 0.24207340 -0.97875668 -0.58813087 0.14427350 -0.61873144 -0.38421131 -0.77768481 0.84994249 0.49572154 0.86362573 0.05213835 0.97147518 0.72495971 0.33554201 0.88218359 -0.02003787 -0.29239537 0.45134487 0.66200677 0.30802605 -0.65774550 -0.08149437 -0.37795979 -0.77430176 -0.75072201 0.51297984 -0.02603652 0.18323856 0.58426463 0.74015658 0.82866461 -0.86769512 0.70423284 -0.38430215 -0.16766046 -0.96939287 -0.45494868 -0.86386005 -0.53601743 -0.41311978 -0.82943151 -0.18515687 0.20623515 0.06647726 0.91913452 0.55014454 -0.07461521 0.99616771 0.39808667 -0.54237055 -0.11408583 0.70047299 -0.23070901 -0.54778626 0.65259125 -0.56524431 0.30300002 -0.89603466 0.76911011 0.49328279 0.50839647 -0.19727400 -0.68756584 -0.07298230 -0.34776838 -0.89275958 -0.23690190 -0.92239684 0.30665494 -0.28982189 0.53553365 -0.00057833 -0.68829445 0.78718234 -0.58433993 0.57476084 0.15904028 -0.20890286 0.09035692 0.08881355 -0.33727115 0.94161979 -0.07278939 0.42737516 0.50462476 0.45055703 0.43525339 0.66745109 -0.11119547 0.18188953 0.85647022 0.87327104 0.31801063 0.93861836 -0.74126614 0.71700485 -0.47129467 -0.30718999 0.06162447 0.73130663 -0.30652513 -0.95586512 0.37808679 0.71043815 -0.01998251 -0.59424799 -0.80730057 -0.17827675 -0.82406971 0.12498971 -0.06188978 -0.74004371 0.64270367 0.53716192 -0.36470688 0.25424118 0.17272700 -0.86935053 -0.74307542 0.02510924 -0.42406819 0.56495917 0.07704918 -0.49651760 -0.12763604 0.30301213 -0.54601361 -0.33753347 -0.91631661 0.96247567 -0.54996550 -0.93995524 0.50041170 0.11962207 -0.38316851 -0.21390578 0.87976651 0.04508451 0.16816510 -0.60994601 0.55493582 -0.32361836 -0.56216559 0.96810208 -0.70511581 0.84300859 -0.89431241 0.31580809 0.24785562 -0.54227393 -0.86733432 -0.78135457 -0.69126622 0.31346996 0.53654334 0.78059548 -0.17313862 -0.23968231 -0.19664312 -0.97660760 -0.05244464 0.85364998 -0.99780619 0.65146757 0.30182409 -0.44708130 0.47460028 0.75394102 -0.24229357 -0.65640243 -0.02491093 -0.68741560
+0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.72609352 0.34497286 0.12680240 0.06599022 0.77746262 0.88106571 0.44784285 0.72949212 0.58883457 0.88830479 0.79289140 0.04935098 0.83555670 0.08928721 0.79355471 0.06853090 0.79703472 0.53855029 0.17305420 0.19814494 0.04259901 0.14089267 0.45045416 0.59879353 0.54456431 0.41121374 0.97907864 0.98365497 0.50215667 0.61427839 0.51817755 0.46064857 0.08418529 0.72550321 0.48289127 0.27598311 0.45198201 0.89180464 0.79367424 0.17393922 0.51379103 0.33980158 0.37990465 0.27519758 0.41906746 0.88365731 0.10759245 0.53972536 0.18968538 0.72611399 0.89497533 0.25778975 0.88155089 0.24234584 0.39755591 0.18835502 0.39313498 0.75066210 0.23867804 0.06592729 0.37633517 0.23793974 0.19691747 0.37191446 0.00825285 0.61860394 0.28478662 0.92027341 0.23012974 0.02065724 0.62125309 0.91694779 0.64949486 0.91804202 0.25973726 0.35015328 0.61954935 0.98235924 0.94903081 0.59568320 0.78332012 0.56629237 0.59476408 0.68502790 0.34214237 0.66093102 0.77108674 0.97812269 0.02440558 0.73522019 0.54703798 0.94244869 0.05587325 0.93722448 0.30158822 0.39777551 0.70008098 0.37760430 0.60175864 0.09482047 0.75869890 0.38015035 0.48143845 0.79475324 0.49506264 0.36057695 0.02118273 0.32989726 0.21459601 0.19102511 0.68166269 0.23330527 0.45458171 0.99535427 0.95565433 0.48271305 0.40205729 0.16722536 0.12788851 0.98108973 0.21382246 0.28217607 0.65676281 0.76527920 0.57791066 0.97991297 0.24979487 0.57652233 0.53953183 0.05423453 0.57898363 0.59113270 0.66032129 0.24041500 0.85439292 0.48332287 0.37748902 0.62236727 0.11821902 0.73823681 0.45214543 0.78361644 0.07204445 0.63973845 0.24865328 0.01609243 0.89682026 0.83997193 0.89826105 0.41163929 0.19699424 0.28018000 0.74085755 0.69748778 0.84721648 0.42504576 0.13633491 0.60161155 0.85951175 0.93201529 0.05649487 0.44313440 0.90936312 0.28268766 0.43736583 0.23491750 0.78608781 0.80988918 0.78598545 0.82737381 0.69908148 0.17967194 0.61064880 0.17919408 0.92165866 0.00979013 0.65441496 0.03370965 0.42405833 0.11032524 0.31932731 0.32123829 0.81774903 0.10572887 0.72540256 0.62027229 0.43169923 0.21575899 0.15451722 0.48680571 0.57377853 0.20923867 0.28545189 0.68143985 0.47293111 0.69335181 0.99549345 0.66898159 0.55757862 0.26419361 0.36922194 0.58907902 0.88272554 0.73572911 0.16855317 0.95189536 0.55594004 0.78555199 0.06180726 0.02529581 0.21080574 0.84347912 0.80937612 0.89051795 0.76872778 0.33959955 0.56155344 0.62605634 0.89798071 0.50666083 0.93555426 0.62957690 0.84713377 0.62452176 0.64996087 0.28347443 0.39057060 0.12489169 0.13095346 0.06175207 0.68000886 0.49034976 0.83309832 0.93236553 0.38717920 0.97799442 0.39707249 0.75282117 0.47097601 0.19696510 0.23232877 0.25055210 0.52606108 0.12035172 0.10219646 0.84278544 0.12472902 0.39798545 0.35270770 0.59296620 0.40565257 0.93478540 0.79775150 0.38003298 0.74858197 0.47868139 0.98013993 0.13449069 0.65750665 0.45857023 0.00052158 0.25610353 0.94727539 0.88090052 0.03620729 0.38092263 0.04513119 0.75418470 0.19667640 0.00379049

requirements.txt

@@ -1,13 +1,13 @@
-# The ZeroGPU (Gradio SDK) deployment: PyTorch on the shared pool; jax stays as a CPU library dependency of dmipy-sim
-# and never touches the GPU. The dmipy packages by commit (the torch backends: dmipy-sim#511, dmipy-fit#38,
-# dmipy-tract#5).
+# The ZeroGPU (Gradio SDK) deployment: PyTorch on the shared pool; jax stays as a CPU library dependency of the dmipy
+# packages (the brain page's pose expansion runs on it) and never touches the GPU. The dmipy packages are the released
+# 3.0 line from PyPI; dmipy-replay's extras: `pose` (the pose expansion and the CPU thread cap), `jax`, `torch`.
 torch>=2.4
 tomli; python_version < "3.11"          # the ZeroGPU image is Python 3.10 (its README python_version is ignored)
 spaces>=0.30
 gradio==6.28.0
-jax>=0.4.30
+jax>=0.6.2
 jaxopt==0.8.5
-numpy>=1.26,<3
+numpy>=2.0,<3
 scipy>=1.11
 nibabel==5.4.2
 safetensors>=0.4
@@ -16,8 +16,8 @@ requests>=2.31
 matplotlib>=3.8
 plotly>=5.20
 pillow>=10
-dipy>=1.9
-# dmipy-sim: main at the merge of dmipy-sim#598 (the field declared by the spec's pools)
-dmipy-sim @ git+https://github.com/dmrai-lab/dmipy-sim@1b79d2e06a958578bc523c84bacf54bae86fbb95
-dmipy-fit @ git+https://github.com/dmrai-lab/dmipy-fit@c3dcc546595cbde8e8e7c7bb2ff4e87f099d5dab
-dmipy-tract @ git+https://github.com/dmrai-lab/dmipy-tract@7da22c3a9473c50a08cc7c6246ec16b68a64d877
+dmipy-replay[pose,jax,torch]==0.2.0
+dmipy-analytic==3.0.0
+dmipy-sim==3.0.0
+dmipy-fit==3.0.0
+dmipy-tract==0.1.0

space/app.py

@@ -31,11 +31,14 @@ from dataclasses import replace
 
 from . import pipeline as P
 from . import sources
+from . import tensor as T
 from . import viewers as V
 
 MAX_SHELLS = 4
 CUSTOM = "custom shells"
 UPLOADED = "uploaded scheme"
+SNR_NOTE = ("µFA is noise-limited: the b-tensor preset starts at SNR 100, because at SNR 30 with the tissue at 3 T the "
+            "spherical encoding's contrast in the densest strand voxels is below the noise and µFA reads 0 there.")
 PER_SHELL = 7                                            # on, timing class, b, directions, delta, Delta, TE
 NO_KNOB = P.NO_KNOB
 IDEAL = P.IDEAL
@@ -44,9 +47,12 @@ SAMPLE = 10_000                                          # streamlines kept in t
 OUTPUTS = ("result", "headline", "dwi_view", "tract_view", "mats", "timings", "tck", "volumes", "z_slider", "m_slider",
            "tract_view_b", "mats_b", "b_row", "floor_view", "accuracy", "spread_view",
            "explore_view", "metric_view", "ingredient_pool", "ingredient_contact", "ingredient_field", "layer_table", "layer_choice",
-           "ez_slider", "em_slider", "truth_view", "fractions_view", "lobar_view", "roundtrip", "response_view")   # the run button's outputs, in order
+           "ez_slider", "em_slider", "truth_view", "fractions_view", "lobar_view", "roundtrip", "response_view",
+           "tensor_row", "tensor_view", "tensor_pick", "powder_view", "tensor_view_b", "connectome_row")   # the run button's outputs, in order
 EXPLORE_MODES = ("signal", "minus the previous layer", "B − A", "(B − A) ÷ floor")
-METRICS = ("DWI", "MD (µm²/ms)", "FA")
+METRICS = ("DWI", "MD (µm²/ms)", "FA", "µFA")                 # µFA on a b-tensor run (the linear vs spherical powder signal)
+ROWS = ("every row", "LTE", "STE")                             # a b-tensor run's DWI rows: all, or one encoding's (the measurement slider counts within it)
+PICK = 8                                                       # pixels per voxel of the clickable µFA slice
 _state = {"source": None, "error": None}
 _lock = threading.Lock()
 _RUNS = tempfile.mkdtemp(prefix="disco-runs-")          # one directory per run tag, replaced on every run
@@ -104,9 +110,11 @@ def gradient_text(cfg, protocol, shapes, scanner):
     """The per-shell gradient table as Markdown, and whether every shell is playable on the menu's ``scanner``."""
     rows = P.playable(protocol, shapes, P.machine(cfg, scanner))
     lines = ["| shell | b (s/mm²) | timing | needs | limit |", "|---|---|---|---|---|"]
+    from . import classes as C
     for name, b, delta, Delta, G, G_max, ok in rows:
         limit = "any" if G_max is None else f"{G_max * 1e3:.0f} mT/m {'✓' if ok else '✗ cannot play'}"
-        lines.append(f"| {name} | {b:g} | δ {delta * 1e3:g} / Δ {Delta * 1e3:g} ms | {G * 1e3:.0f} mT/m | {limit} |")
+        timing = C.describe(shapes[name]) if shapes.get(name, {}).get("kind") == "ogse" else f"δ {delta * 1e3:g} / Δ {Delta * 1e3:g} ms"
+        lines.append(f"| {name} | {b:g} | {timing} | {G * 1e3:.0f} mT/m | {limit} |")
     return "\n".join(lines), all(r[-1] for r in rows)
 
 
@@ -196,8 +204,20 @@ def _fill(source, runs, ladder_on, prepared):
     return out, computed, time.perf_counter() - t0
 
 
+def _tensor_state(tm):
+    """A b-tensor run's maps for the page: FA, MD and µFA in float32, the powder averages per (encoding, b), the b-values."""
+    return dict(fa=np.asarray(tm["fa"], np.float32), md=np.asarray(tm["md"], np.float32), ufa=np.asarray(tm["ufa"], np.float32),
+                powder={k: np.asarray(v, np.float32) for k, v in tm["powder"].items()}, bvals=list(tm["bvals"]))
+
+
 def _result_state(res, sample, load_seconds):
-    """What the results tab needs, kept per session: float32 volumes, the peaks, the streamline sample, the numbers."""
+    """What the results tab needs, kept per session: float32 volumes, the peaks, the streamline sample, the numbers; on
+    a b-tensor run the maps instead of the peaks, the connectome and the streamlines."""
+    if res.protocol.tensor:
+        return dict(dwi=res.dwi.astype(np.float32), floor=res.floor.astype(np.float32), floor_median=P.floor_stats(res)["median"], meas=res.meas,
+                    peaks=None, peak_amp=None, matrix=None, score=None, seconds=res.seconds, load_seconds=load_seconds, tractogram=None,
+                    n_streamlines=0, shape=res.dwi.shape[:3], name=res.protocol.name, extras=dict(fractions=None),
+                    tensor=_tensor_state(res.extras["tensor"]))
     pk, amp = P.peaks(res.sh)
     fr = res.extras.get("fractions")
     return dict(dwi=res.dwi.astype(np.float32), floor=res.floor.astype(np.float32), floor_median=P.floor_stats(res)["median"], meas=res.meas,
@@ -222,11 +242,15 @@ def _explorer_state(results, ladder, ingredients, knob, source):
     layers = list(ladder) + [(name_a, ra.clean)] + ([(name_b, results["B"].clean)] if "B" in results else [])
     metrics = {}
     for label, vol in layers:
-        try:
-            md, fa = P.dti(vol, ra.meas, mask)
-        except ValueError:                                  # too few low-b rows for a tensor: no metric maps
-            md = fa = None
-        metrics[label] = (None if md is None else md.astype(np.float16), None if fa is None else fa.astype(np.float16))
+        ufa = None
+        if ra.protocol.tensor:                               # the b-tensor maps of the noise-free layer: DTI on the linear rows, the powder fit
+            tm = T.maps(vol, ra.meas, mask); md, fa, ufa = tm["md"], tm["fa"], tm["ufa"]
+        else:
+            try:
+                md, fa = P.dti(vol, ra.meas, mask)
+            except ValueError:                              # too few low-b rows for a tensor: no metric maps
+                md = fa = None
+        metrics[label] = tuple(None if x is None else np.asarray(x).astype(np.float16) for x in (md, fa, ufa))
     return dict(layers=[(label, vol.astype(np.float16)) for label, vol in layers], metrics=metrics, floor=ra.floor.astype(np.float32),
                 meas=ra.meas, mask=mask, ingredients=source.ingredient_layers(ingredients), differences=P.layer_differences(layers, ra.meas, mask),
                 snr=ra.snr, physics=ra.physics)
@@ -255,6 +279,10 @@ def _write_files(tag, res, source):
     sample as .tck, the DWI and FOD volumes, the source's own files; ``(sample, {"tck": [...], "volumes": [...]})``."""
     out = os.path.join(_RUNS, tag); shutil.rmtree(out, ignore_errors=True); os.makedirs(out)
     stem = f"{source.describe(source.cfg)['files']}_{tag}_{res.protocol.name.replace(' ', '_')}"
+    if res.protocol.tensor:                              # a b-tensor run: the DWI with its table and encodings, and the maps
+        import re
+        stem = re.sub(r"[^A-Za-z0-9_.+-]", "_", stem)
+        return None, dict(tck=[], volumes=P.write_tensor_volumes(res, out, prefix=stem, affine=source.affine))
     tck = os.path.join(out, f"{stem}.tck"); res.tractogram.to_tck(tck)
     sample = P.sample_tractogram(res.tractogram, SAMPLE)
     sample_path = os.path.join(out, f"{stem}_sample{SAMPLE // 1000}k.tck"); sample.to_tck(sample_path)
@@ -262,24 +290,36 @@ def _write_files(tag, res, source):
     return sample, dict(tck=[tck, sample_path], volumes=[vols["dwi"], vols["bvals"], vols["bvecs"], vols["fod"]] + source.extra_files(res, out, stem))
 
 
-def explore(ex, layer, mode, metric, z, m):
+def row_of(meas, rows, m):
+    """The measurement index the slider's ``m`` names: the ``m``-th row of the encoding ``rows`` (``LTE`` / ``STE``) on
+    a b-tensor run, else row ``m``; clipped to the rows there are."""
+    enc = getattr(meas, "encoding", None)
+    if enc is None or rows not in ROWS[1:]:
+        return min(int(m), len(meas.bvals) - 1)
+    idx = np.flatnonzero(np.asarray(enc) == rows)
+    return int(idx[min(int(m), len(idx) - 1)])
+
+
+def explore(ex, layer, mode, metric, z, m, rows=ROWS[0]):
     """The explorer's map: ``layer`` of the state's layers, in ``mode`` (the signal, its difference to the previous
-    layer, B minus A, or that over the replay floor), as the DWI at measurement ``m`` or a tensor metric, slice ``z``."""
+    layer, B minus A, or that over the replay floor), as the DWI at measurement ``m`` (on a b-tensor run the ``m``-th
+    row of the encoding ``rows``) or a metric (MD, FA, and on a b-tensor run µFA), slice ``z``."""
     if not ex:
         return None
     labels = [l for l, _ in ex["layers"]]; vols = dict(ex["layers"])
     if layer not in vols:
         layer = labels[-1]
-    z = int(z); m = min(int(m), len(ex["meas"].bvals) - 1)
+    z = int(z); m = row_of(ex["meas"], rows, m)
     k = METRICS.index(metric) if metric in METRICS else 0
 
     def field(label):
         if k == 0:
             return np.asarray(vols[label], np.float32)[..., m]
-        md, fa = ex["metrics"][label]
-        v = md if k == 1 else fa
+        v = ex["metrics"][label][k - 1]
         return None if v is None else np.asarray(v, np.float32)
-    what = f"{metric} of {layer}" if k else f"{layer}: b = {ex['meas'].bvals[m]:g} s/mm², measurement {m}"
+    enc = getattr(ex["meas"], "encoding", None)
+    tag = "" if enc is None else f" ({enc[m]})"
+    what = f"{metric} of {layer}" if k else f"{layer}: b = {ex['meas'].bvals[m]:g} s/mm²{tag}, measurement {m}"
     x = field(layer)
     if x is None:
         return None
@@ -355,7 +395,7 @@ def compute(prepared, preset, n_b0, snr_on, snr, density, max_angle, step_mm, ke
         yield (f"**Replay DWI Explorer · the ingredient maps** … ({time.perf_counter() - t0:.0f} s so far)", 0.85)
         t = time.perf_counter(); ingredients = source.ingredients(meas_a, physics_a, prepared["runs"]["A"]); seconds["explorer · ingredients"] = time.perf_counter() - t
         spread = None
-        if int(n_keys) > 1:                                 # A's tracking repeated over further keys: the tractogram's own spread
+        if int(n_keys) > 1 and not results["A"].protocol.tensor:                                 # A's tracking repeated over further keys: the tractogram's own spread
             keys = [int(key) + 1 + i for i in range(int(n_keys) - 1)]
             mats = [results["A"].matrix]; scores = [results["A"].score]
             for i, (k, M, sc, secs) in enumerate(P.repeat_tracking(results["A"], source, tracking, keys)):
@@ -397,9 +437,11 @@ def present(payload, state):
     post["page · states"] = time.perf_counter() - t; t = time.perf_counter()
     labels = [l for l, _ in ex["layers"]]; name_a = labels[-2] if "B" in results else labels[-1]
     ra = rs["A"]; rb = rs.get("B")
-    headline = source.score_text("A", results["A"])
+    tensor = results["A"].protocol.tensor
+    text = (lambda tag, r: tensor_text(tag, r)) if tensor else source.score_text
+    headline = text("A", results["A"])
     if rb:
-        headline += "<br>" + source.score_text("B", results["B"]) + "<br>" + source.compare_text(source.compare(results["A"], results["B"]), knob)
+        headline += "<br>" + text("B", results["B"]) + ("" if tensor else "<br>" + source.compare_text(source.compare(results["A"], results["B"]), knob))
     if spread:
         headline += (f"<br>**A over {spread['n']} tracker keys: {spread['key']} {spread['pearson_mean']:.3f} ± {spread['pearson_std']:.3f}**, "
                      f"median pair count CV {spread['cv_median']:.2f}, {spread['pairs_always']} pairs in every run, {spread['pairs_any']} in any.")
@@ -413,26 +455,49 @@ def present(payload, state):
     rs["explorer"] = ex                                  # the page state: A, B and the explorer's layers
     box = V.grid_box(ra["shape"], source.affine) if d["views"]["truth"] else ra["shape"]
     pa = source.cached(results["A"].meas, results["A"].physics)
-    out = (rs, headline, V.dwi_slice(ra["dwi"], ra["meas"], z0, m0, peaks=ra["peaks"], peak_amp=ra["peak_amp"], overlay=True, label="A: "),
-           V.tractogram3d(ra["tractogram"], regions, box, total=ra["n_streamlines"]),
-           source.matrices(results["A"].matrix, results["A"].score, results["A"].reference),
+    connectome = not tensor
+    t_a, t_pick, t_b = tensor_views(rs, z0)
+    out = (rs, headline, V.dwi_slice(ra["dwi"], ra["meas"], z0, m0, peaks=ra["peaks"], peak_amp=ra["peak_amp"], overlay=connectome, label="A: "),
+           V.tractogram3d(ra["tractogram"], regions, box, total=ra["n_streamlines"]) if connectome else None,
+           source.matrices(results["A"].matrix, results["A"].score, results["A"].reference) if connectome else None,
            timings, [t_ for f in files.values() for t_ in f["tck"]], [v for f in files.values() for v in f["volumes"]],
            _slider_update(z0, ra["dwi"].shape[2] - 1), _slider_update(m0, results["A"].protocol.n_meas - 1),
-           V.tractogram3d(rb["tractogram"], regions, box, total=rb["n_streamlines"]) if rb else None,
-           source.matrices(results["B"].matrix, results["B"].score, results["B"].reference) if rb else None,
-           gr.update(visible=rb is not None),
+           V.tractogram3d(rb["tractogram"], regions, box, total=rb["n_streamlines"]) if rb and connectome else None,
+           source.matrices(results["B"].matrix, results["B"].score, results["B"].reference) if rb and connectome else None,
+           gr.update(visible=rb is not None and connectome),
            V.floor_slice(ra["floor"], z0, ra["floor_median"], label="A: ") if d["views"]["floor"] else None, source.accuracy(results["A"]),
            V.spread_matrices(spread) if spread else None,
            explore(ex, name_a, "minus the previous layer" if len(labels) > 1 else "signal", METRICS[0], z0, m0),
-           explore(ex, name_a, "signal", METRICS[2], z0, m0), *ingredient_views(ex, z0), layer_table(ex),
+           explore(ex, name_a, "signal", METRICS[3] if tensor else METRICS[2], z0, m0), *ingredient_views(ex, z0), layer_table(ex),
            gr.update(choices=labels, value=name_a), _slider_update(z0, ra["dwi"].shape[2] - 1), _slider_update(m0, results["A"].protocol.n_meas - 1),
            source.truth_view(ra["dwi"], ra["meas"], z0, m0), source.fractions_view(ra["extras"], z0),
-           source.lobar(results["A"].matrix, results["A"].score, results["A"].reference), source.roundtrip_rows(results["A"]),
-           source.response_view(pa, results["A"].extras))
+           source.lobar(results["A"].matrix, results["A"].score, results["A"].reference) if connectome else None,
+           source.roundtrip_rows(results["A"]) if connectome else None,
+           source.response_view(pa, results["A"].extras) if connectome else None,
+           gr.update(visible=tensor), t_a, t_pick, None, t_b, gr.update(visible=connectome))
     timings.append(["page · figures", f"{time.perf_counter() - t:.2f}"])
     return out
 
 
+def tensor_text(tag, res):
+    """The headline of a b-tensor run: the maps' medians over the voxels the ground truth's strands cross, and over the
+    densest quarter of them (strand length in the upper quartile) by their crossing count, the run's settings."""
+    tm = res.extras["tensor"]; physics = res.physics
+    n_dirs, length = _strands()
+    sig = np.isfinite(res.dwi[..., 0]); strand = sig & (length > 0)
+    dense = strand & (length >= np.quantile(length[strand], 0.75)) if strand.any() else strand
+    med = lambda x, m: float(np.nanmedian(np.asarray(x)[m])) if np.any(m) else float("nan")
+    snr = "" if res.snr is None else f", SNR {res.snr:g}"
+    by = ", ".join(f"{k}{'+' if k == 3 else ''} directions {med(tm['ufa'], dense & ((n_dirs >= 3) if k == 3 else (n_dirs == k))):.2f} / "
+                   f"{med(tm['fa'], dense & ((n_dirs >= 3) if k == 3 else (n_dirs == k))):.2f}" for k in (1, 2, 3))
+    return (f"**{tag}: in the voxels the strands cross, µFA median {med(tm['ufa'], strand):.3f}, FA {med(tm['fa'], strand):.3f}, MD "
+            f"{med(tm['md'], strand):.3f} µm²/ms**; the densest quarter by its crossing count, µFA / FA: {by} "
+            f"({res.protocol.n_meas} measurements: linear and spherical b-tensors at b {', '.join(f'{b:g}' for b in tm['bvals'])} s/mm², "
+            f"{physics.label() if physics else 'bare diffusion'}{snr}, {res.seconds['total']:.1f} s; replay floor median "
+            f"{P.floor_stats(res)['median']:.4f}). µFA is the anisotropy the spherical encoding removes and the orientation "
+            "dispersion does not: high where FA is low in crossings. Click the µFA slice for a voxel's powder curves.")
+
+
 def run_pipeline(compute_fn, *args, progress=None):
     """The run button, a generator: the source's share of the runs the page has cached (:func:`prepare_runs`, a
     lookup), then the stage texts of ``compute_fn`` (:func:`compute`, or it wrapped for a GPU pool) into the headline
@@ -473,13 +538,14 @@ def estimated_seconds(preset, n_b0, snr_on, snr, density, max_angle, step_mm, ke
     (:meth:`~space.pipeline.Source.estimated_seconds`) with the source's share the call computes because the page has
     it not cached (:func:`uncached_responses`). The pool refuses a request above the visitor's daily quota
     (:data:`GPU_TIERS`) and kills a run that outlives its reservation, so this is the measured cost with its margin,
-    not a generous one; 480 when the inputs make no protocol."""
+    not a generous one; 480 when the inputs make no protocol (the refusal is a ``ValueError``; any other error is a
+    defect and raises)."""
     try:
         cfg = P.config()
         S = sources.source_class(cfg)
         n = len(S.panel(cfg).fields)
         protocol = _protocol_from_inputs(S, cfg, preset, n_b0, *rest[n:], scheme=scheme_file, full=S.mode == "full")
-    except Exception:
+    except ValueError:
         return 480
     responses = uncached_responses(preset, n_b0, snr_on, snr, density, max_angle, step_mm, key, scheme_file, knob, scanner, n_keys, ladder_on, *rest)
     try:
@@ -523,19 +589,63 @@ def _slider_update(value, maximum):
     return gr.update(value=int(value), maximum=int(maximum))
 
 
-def redraw_explorer(rs, layer, mode, metric, z, m):
+def redraw_explorer(rs, layer, mode, metric, z, m, rows=ROWS[0]):
     ex = rs and rs.get("explorer")
-    return (explore(ex, layer, mode, metric, z, m), explore(ex, layer, mode, METRICS[2] if metric == METRICS[0] else metric, z, m), *ingredient_views(ex, z))
+    second = (METRICS[3] if ex and ex["meas"].encoding is not None else METRICS[2]) if metric == METRICS[0] else metric
+    return (explore(ex, layer, mode, metric, z, m, rows), explore(ex, layer, mode, second, z, m, rows), *ingredient_views(ex, z))
 
 
-def redraw_slice(rs, z, m, overlay, which):
+def redraw_slice(rs, z, m, overlay, which, rows=ROWS[0]):
     if not rs or which not in rs:
-        return None, None, None, None
+        return None, None, None, None, None, None, None
     r = rs[which]; source = _load()["source"]
     views = source.describe(source.cfg)["views"]
-    return (V.dwi_slice(r["dwi"], r["meas"], int(z), int(m), peaks=r["peaks"], peak_amp=r["peak_amp"], overlay=bool(overlay), label=f"{which}: "),
+    m = row_of(r["meas"], rows, m)
+    return (V.dwi_slice(r["dwi"], r["meas"], int(z), m, peaks=r["peaks"], peak_amp=r["peak_amp"], overlay=bool(overlay) and r["peaks"] is not None, label=f"{which}: "),
             V.floor_slice(r["floor"], int(z), r["floor_median"], label=f"{which}: ") if views["floor"] else None,
-            source.truth_view(r["dwi"], r["meas"], int(z), int(m)), source.fractions_view(r["extras"], int(z)))
+            source.truth_view(r["dwi"], r["meas"], int(z), m), source.fractions_view(r["extras"], int(z)), *tensor_views(rs, int(z)))
+
+
+def tensor_views(rs, z):
+    """A b-tensor run's maps at slice ``z`` for A (with the ground truth's crossings outlined) and B, and the clickable
+    µFA slice of A; Nones on any other run."""
+    if not rs or "tensor" not in rs.get("A", {}) or rs["A"]["tensor"] is None:
+        return None, None, None
+    crossing = _crossings()
+    a = V.tensor_maps(rs["A"]["tensor"], crossing, z, label="A: ")
+    b = V.tensor_maps(rs["B"]["tensor"], crossing, z, label="B: ") if rs.get("B") and rs["B"].get("tensor") else None
+    return a, V.map_picker(rs["A"]["tensor"]["ufa"][:, :, z], PICK, vmin=0.0, vmax=1.0), b
+
+
+def pick_voxel(rs, z, evt):
+    """The powder curve (linear and spherical, A and dashed B) of the voxel clicked on the µFA slice: ``evt.index`` is
+    the click's pixel ``(x, y)``, a voxel ``PICK`` pixels wide, the slice drawn with y up."""
+    if not rs or not rs.get("A", {}).get("tensor"):
+        return None
+    tm = rs["A"]["tensor"]; nx, ny = tm["ufa"].shape[:2]
+    x, y = (int(v) for v in evt.index)
+    i, j = min(x // PICK, nx - 1), min(ny - 1 - y // PICK, ny - 1)
+    other = rs["B"]["tensor"] if rs.get("B") and rs["B"].get("tensor") else None
+    return V.powder_curve(tm, (i, j, int(z)), other=other)
+
+
+_CROSSINGS = {}
+
+
+def _crossings():
+    """The voxels the ground truth's strands cross in two or more directions (the source's strand configuration),
+    computed once per process."""
+    return _strands()[0] >= 2
+
+
+def _strands():
+    """``(n_dirs, length)`` per voxel of the ground truth's strands (the source's strand configuration), computed once
+    per process."""
+    if "c" not in _CROSSINGS:
+        from .sources.disco import strand_configuration
+        n_dirs, _, length = strand_configuration(shape=_load()["source"].mask.shape)
+        _CROSSINGS["c"] = (np.asarray(n_dirs), np.asarray(length))
+    return _CROSSINGS["c"]
 
 
 def ground_truth_views():
@@ -563,13 +673,26 @@ def _widget(c):
     raise ValueError(f"unknown control kind {c.kind!r}")
 
 
+def unavailable_page(title, text):
+    """The Blocks of a source whose inputs cannot be read: its text (:meth:`space.pipeline.Source.unavailable`) and nothing
+    to run."""
+    import gradio as gr
+    with gr.Blocks(title=title) as demo:
+        gr.Markdown(text)
+    return demo
+
+
 def build(runner=None, cfg=None):
-    """The Blocks of the configured source (``cfg``, else :func:`space.pipeline.config`). ``runner`` wraps
+    """The Blocks of the configured source (``cfg``, else :func:`space.pipeline.config`), or the page that says why the
+    source cannot be served (:func:`unavailable_page`). ``runner`` wraps
     :func:`compute` for the run button (the ZeroGPU entry passes ``spaces.GPU(...)``): the device part of a run; the
     page looks up its cache before it and draws from its payload after it, in this process."""
     import gradio as gr
     cfg = cfg or P.config()
     S = sources.source_class(cfg)
+    refusal = S.unavailable(cfg)
+    if refusal is not None:
+        return unavailable_page(*refusal)
     d = S.describe(cfg)
     full = S.mode == "full"
     shapes = S.shapes_of(cfg); presets = S.presets(cfg) + [CUSTOM] + ([UPLOADED] if full else [])
@@ -624,7 +747,9 @@ def build(runner=None, cfg=None):
                         physics_inputs = [widgets[name] for name in panel.fields]
                         with gr.Row():
                             snr_on = gr.Checkbox(value=True, label="add Rician noise")
-                            snr = gr.Slider(5, 100, value=30, step=1, label="SNR at M0 (a full water voxel before relaxation; each voxel's b = 0 SNR follows its tissue)")
+                            snr = gr.Slider(5, 100, value=P.preset_snr(cfg, presets[0]), step=1, label="SNR at M0 (a full water voxel before relaxation; each voxel's b = 0 SNR follows its tissue)")
+                        if any("tensor" in p for p in (cfg.get("presets") or {}).values()):
+                            gr.Markdown(SNR_NOTE)
                         density = gr.Slider(tc["density"][0], tc["density"][1], value=tc["density"][2], step=tc["density"][3], label=tc["density"][4])
                         max_angle = gr.Slider(tc["max_angle"][0], tc["max_angle"][1], value=tc["max_angle"][2], step=tc["max_angle"][3], label=tc["max_angle"][4])
                         step_mm = gr.Slider(tc["step"][0], tc["step"][1], value=tc["step"][2], step=tc["step"][3], label=tc["step"][4])
@@ -636,7 +761,7 @@ def build(runner=None, cfg=None):
                         scanner_terms = gr.Markdown(P.scanner_text(cfg, IDEAL, S.refused_machines(cfg)))
                         gradients = gr.Markdown()
                         n_keys = gr.Slider(1, 8, value=1, step=1, label="repeat A's tracking over N keys (the tractogram's own spread)")
-                        ladder_on = gr.Checkbox(value=not full, visible=not full, label=d["ladder"])
+                        ladder_on = gr.Checkbox(value=False, visible=not full, label=d["ladder"])
                         go = gr.Button(d["run_label"], variant="primary")
                         gpu_text = gr.Markdown(visible=runner is not None)
                         headline = gr.Markdown()
@@ -657,6 +782,7 @@ def build(runner=None, cfg=None):
                 with gr.Row():
                     ez_slider = gr.Slider(0, 39, value=20, step=1, label="axial slice z")
                     em_slider = gr.Slider(0, 363, value=0, step=1, label="measurement (DWI only)")
+                    erows = gr.Radio(list(ROWS), value=ROWS[0], label="rows (a b-tensor run: the measurement counts within the encoding)")
                 with gr.Row():
                     explore_view = gr.Image(label="the chosen layer and mode", type="pil")
                     metric_view = gr.Image(label="the same for FA (or the chosen metric)", type="pil")
@@ -671,17 +797,25 @@ def build(runner=None, cfg=None):
                             m_slider = gr.Slider(0, 363, value=0, step=1, label="measurement")
                             overlay = gr.Checkbox(value=True, label="FOD principal directions")
                             which = gr.Radio(["A", "B"], value="A", label="run")
+                            rows_radio = gr.Radio(list(ROWS), value=ROWS[0], label="rows (b-tensor runs)")
                     with gr.Column(scale=1, visible=d["views"]["truth"]):
                         truth_view = gr.Image(label="the same slice with the input FOD's principal directions", type="pil")
                     with gr.Column(scale=1):
                         tract_view = gr.Plot(label="tractogram A")
                 fractions_view = gr.Image(label="recovered vs input fractions", type="pil", visible=d["views"]["fractions"])
-                mats = gr.Image(label="connectome A vs its truth", type="pil")
-                lobar_view = gr.Image(label="connectome A vs its truth over the lobar groups", type="pil", visible=d["views"]["lobar"])
-                with gr.Row(visible=False) as b_row:
-                    tract_view_b = gr.Plot(label="tractogram B")
-                    mats_b = gr.Image(label="connectome B vs its truth", type="pil")
-                spread_view = gr.Image(label="A over N tracker keys: mean and spread per pair", type="pil")
+                with gr.Column(visible=False) as tensor_row:            # a b-tensor run: the maps and a voxel's powder curves
+                    tensor_view = gr.Image(label="b-tensor maps of A: FA, µFA, MD (the slice above)", type="pil")
+                    with gr.Row():
+                        tensor_pick = gr.Image(label="µFA of A: click a voxel", type="pil", interactive=False)
+                        powder_view = gr.Image(label="the clicked voxel's powder signal: linear and spherical, A and B (dashed)", type="pil")
+                    tensor_view_b = gr.Image(label="b-tensor maps of B", type="pil")
+                with gr.Column() as connectome_row:
+                    mats = gr.Image(label="connectome A vs its truth", type="pil")
+                    lobar_view = gr.Image(label="connectome A vs its truth over the lobar groups", type="pil", visible=d["views"]["lobar"])
+                    with gr.Row(visible=False) as b_row:
+                        tract_view_b = gr.Plot(label="tractogram B")
+                        mats_b = gr.Image(label="connectome B vs its truth", type="pil")
+                    spread_view = gr.Image(label="A over N tracker keys: mean and spread per pair", type="pil")
                 roundtrip = gr.Dataframe(headers=["what", "value"], label="the round trip: the reconstruction and the connectome against the input",
                                          interactive=False, wrap=True, visible=d["views"]["roundtrip"])
                 with gr.Accordion("accuracy: what this replay is an approximation of", open=False):
@@ -699,6 +833,7 @@ def build(runner=None, cfg=None):
         def run_with_progress(*args, progress=gr.Progress()):
             yield from run_pipeline(compute_fn, *args, progress=progress)
         tissue_numbers = [widgets[name] for name in panel.catalogue] + [catalogue_note]
+        preset.change(lambda name: P.preset_snr(cfg, name), inputs=preset, outputs=snr, show_progress="hidden")
         field_preset.change(lambda name: [field_presets[name]] + S.catalogue_numbers(cfg, field_presets[name]), inputs=field_preset,
                             outputs=[widgets["field_T"]] + tissue_numbers, show_progress="hidden")
         reset.click(lambda f: S.catalogue_numbers(cfg, f), inputs=widgets["field_T"], outputs=tissue_numbers, show_progress="hidden")
@@ -725,9 +860,11 @@ def build(runner=None, cfg=None):
                        floor_view=floor_view, accuracy=accuracy, spread_view=spread_view, explore_view=explore_view, metric_view=metric_view,
                        ingredient_pool=ingredient_pool, ingredient_contact=ingredient_contact, ingredient_field=ingredient_field,
                        layer_table=layer_tbl, layer_choice=layer_choice, ez_slider=ez_slider, em_slider=em_slider, truth_view=truth_view,
-                       fractions_view=fractions_view, lobar_view=lobar_view, roundtrip=roundtrip, response_view=response_view)
-        for ctl in (layer_choice, mode, metric, ez_slider, em_slider):
-            ctl.change(redraw_explorer, inputs=[result, layer_choice, mode, metric, ez_slider, em_slider],
+                       fractions_view=fractions_view, lobar_view=lobar_view, roundtrip=roundtrip, response_view=response_view,
+                       tensor_row=tensor_row, tensor_view=tensor_view, tensor_pick=tensor_pick, powder_view=powder_view,
+                       tensor_view_b=tensor_view_b, connectome_row=connectome_row)
+        for ctl in (layer_choice, mode, metric, ez_slider, em_slider, erows):
+            ctl.change(redraw_explorer, inputs=[result, layer_choice, mode, metric, ez_slider, em_slider, erows],
                        outputs=[explore_view, metric_view, ingredient_pool, ingredient_contact, ingredient_field], show_progress="hidden")
         run_inputs = [preset, n_b0, snr_on, snr, density, max_angle, step_mm, key, scheme_file, knob, scanner, n_keys, ladder_on, *physics_inputs, *shell_inputs]
         if runner is not None:                                   # the pool: what the run reserves, live as the inputs change
@@ -736,8 +873,13 @@ def build(runner=None, cfg=None):
                              outputs=[outputs[name] for name in OUTPUTS], concurrency_limit=1, api_name="run_pipeline")   # the endpoint tools/live.py drives
         if runner is not None:                                   # the run's responses are now cached on the page: the reservation falls
             run_event.then(gpu_seconds_text, inputs=run_inputs, outputs=gpu_text, show_progress="hidden")
-        for ctl in (z_slider, m_slider, overlay, which):
-            ctl.change(redraw_slice, inputs=[result, z_slider, m_slider, overlay, which], outputs=[dwi_view, floor_view, truth_view, fractions_view], show_progress="hidden")
+        for ctl in (z_slider, m_slider, overlay, which, rows_radio):
+            ctl.change(redraw_slice, inputs=[result, z_slider, m_slider, overlay, which, rows_radio],
+                       outputs=[dwi_view, floor_view, truth_view, fractions_view, tensor_view, tensor_pick, tensor_view_b], show_progress="hidden")
+
+        def on_pick(rs, z, evt: gr.SelectData):
+            return pick_voxel(rs, z, evt)
+        tensor_pick.select(on_pick, inputs=[result, z_slider], outputs=powder_view, show_progress="hidden")
         gt_tab.select(ground_truth_views, outputs=[strands_view, gt_matrix])
         demo.load(lambda: (_load().get("error") and f"**the source did not load:** {_load()['error']}") or "", outputs=headline)
     return demo

space/brain.toml

@@ -89,7 +89,7 @@ m0_gm_knob = [0.7, 1.0]
 
 # The scanner menu after "ideal": catalogue machines whose every catalogued term the brain applies, the head centre at
 # isocentre (docs/scanner.md). Each voxel is binned into an encoding class by what the machine delivers there, to
-# `tolerance` of b (dmipy_sim.phantom.bore.encoding_classes), and every class's played acquisition is expanded by the
+# `tolerance` of b (dmipy_replay.phantom.bore.encoding_classes), and every class's played acquisition is expanded by the
 # packs. `refused` are the machines kept off the menu, each shown under it with the reason (Brain.refused_machines):
 # the Swoop, whose own gradient makes the head thousands of encoding classes (docs/scanner.md).
 [scanners]
@@ -103,7 +103,7 @@ refused = [
 tolerance = 0.01
 
 # The pack menu per tissue (the first is the default): a label and a pack URI (hf:// or a local path, as
-# dmipy_sim.phantom.PackSubstrate resolves it). A Hub pack declaring several windows loads the windows a run reaches
+# dmipy_replay.phantom.PackSubstrate resolves it). A Hub pack declaring several windows loads the windows a run reaches
 # (Brain.pack); a local path or a one-window pack loads whole. SubstrateCommons family records feed this later (disco-space#8, decision 5).
 [packs]
 wm = [
@@ -159,7 +159,7 @@ margin = 1.3
 # of 100 / 204 ms (ReplayPack.prefix per shard: 1,303 / 2,656 saves), three runs, the third adding the 85 ms prefix
 # (1,107 saves) and playing the d17-D30 class instead of the scan's. With the field
 # tier off the slope is 0.0008-0.0016 s per save; the field states' slopes (-0.0015 to +0.0062) are within the
-# shared box's noise, and their save-dependent work is the same expansion, so the three take the largest no_field
+# the host's noise, and their save-dependent work is the same expansion, so the three take the largest no_field
 # slope. fixed and per_meas were fit on the 1,002-save pack, so its saves are priced twice: 1.6 s of margin.
 response_cold = { fixed = 22.0, per_meas = 0.07, per_save = 0.0016 }
 response_warm = { fixed = 2.5, per_meas = 0.04, per_save = 0.0016 }
@@ -173,3 +173,7 @@ response_per_class = 5.6
 1 = 1.5                       # 36,605 seeds: 0.7 s
 2 = 3.5                       # 292,840 seeds: 1.9-2.1 s
 4 = 35.0                      # 2,342,720 seeds: 13.6-23.2 s
+
+# The noise a run starts with: the SNR at M0 every preset starts at.
+[noise]
+snr = 30

space/classes.py

@@ -0,0 +1,119 @@
+"""The timing classes as waveforms: each stored class of the time-segment layout is ONE waveform shape (a square-pulse
+PGSE, a stimulated echo, a cosine OGSE with its edges ramped), and a shell on it plays that shape at the amplitude
+its b-value needs, along each of its directions. The shape is built once from the configuration's ``[shapes]`` entry
+by dmipy-replay's sequence builders (:func:`template`), on the layout's grid (TE 53.5 ms, 1000 samples), and checked
+against the layout's stored partition when the layout loads; :func:`play` scales it to the rows. A class has a fixed
+shape: an OGSE built at b = 1000 s/mm^2 with 200 T/m/s edges plays the same edge durations at every other b (the
+layout stores that shape, nothing else)."""
+from __future__ import annotations
+
+import functools
+
+import numpy as np
+
+KINDS = ("pgse", "pgste", "ogse", "tensor")
+N_T = 1000                                           # the layout's grid: 1000 samples over the echo time
+
+
+def kind(sh):
+    """The class's pulse kind: ``pgse`` (the default), ``pgste`` or ``ogse``."""
+    k = sh.get("kind", "pgse")
+    if k not in KINDS:
+        raise ValueError(f"a timing class's kind is one of {KINDS}, got {k!r}")
+    return k
+
+
+def _key(sh):
+    return tuple(sorted((k, v) for k, v in sh.items() if k != "label"))
+
+
+@functools.lru_cache(maxsize=None)
+def _template(key):
+    from dmipy_replay.sequences import ogse, pgse, pgste
+    sh = dict(key)
+    k = kind(sh)
+    if k == "pgse":
+        return pgse([[0.0, 0.0, 1.0]], float(sh["delta"]), float(sh["Delta"]), gradient_strengths=0.1, TE=float(sh["TE"]),
+                    n_t=N_T, slew_rate=np.inf)
+    if k == "pgste":                                  # TM = Delta - delta between the 90s; TE = 2 delta + TM
+        return pgste([[0.0, 0.0, 1.0]], float(sh["delta"]), float(sh["Delta"]) - float(sh["delta"]), gradient_strengths=0.1,
+                     n_t=N_T, slew_rate=np.inf, ste_flip_angles=(90.0, 90.0, 90.0))
+    return ogse([[0.0, 0.0, 1.0]], float(sh["frequency"]), float(sh["duration"]), shape="cosine",
+                bvalues=[float(sh["b_built"]) * 1e6], TE=float(sh["TE"]), n_t=N_T, slew_rate=float(sh["slew_rate"]))
+
+
+def template(sh):
+    """The class's waveform along z, one row (a ScannerSequence); a ``tensor`` class is its three shapes, one row each
+    along its own axis (:func:`space.tensor.template`): ``pgse`` square lobes ``delta`` wide ``Delta`` apart
+    at ``TE``; ``pgste`` square lobes storing for ``TM = Delta - delta``; ``ogse`` a cosine of ``frequency`` (Hz) over
+    ``duration`` (s) on each side of the 180, built at ``b_built`` (s/mm^2) with ``slew_rate`` (T/m/s) edges."""
+    if kind(sh) == "tensor":
+        from . import tensor as T
+        return T.template(sh)
+    return _template(_key(sh))
+
+
+def tensor_rows(sh, encodings, bvals, dirs):
+    """The rows of a ``tensor`` class: each row's encoding (``b0``, ``LTE``, ``STE``) at its b (s/mm^2) about its axis,
+    the physical gradient ``(n, 1000, 3)`` (T/m) -- the class's waveform times the row's level matrix
+    (:func:`space.tensor.levels`); a b = 0 row plays nothing."""
+    from . import tensor as T
+    W = T.waveform(sh)
+    G = np.zeros((len(encodings), N_T, 3))
+    for i, (e, b, d) in enumerate(zip(encodings, bvals, dirs)):
+        if e != "b0" and float(b) > 0:
+            G[i] = W @ T.levels(sh, str(e), float(b) * 1e6, d)
+    return G
+
+
+def tensor_peak(sh, encoding, b, n_dirs):
+    """The largest per-axis gradient (T/m) of a tensor shell's rows (``n_dirs`` axes of :func:`space.tensor.fibonacci`)."""
+    from . import tensor as T
+    return float(np.abs(tensor_rows(sh, [encoding] * int(n_dirs), [b] * int(n_dirs), T.fibonacci(int(n_dirs)))).max())
+
+
+def unit_b(sh):
+    """``(profile (n_t,), b per (T/m)^2)``: the class's physical gradient along its axis normalised to a peak of 1, and
+    the b-value (s/m^2) of that shape at a peak amplitude of 1 T/m."""
+    t = template(sh)
+    G = np.asarray(t.G, np.float64)[0, :, 2]
+    peak = float(np.abs(G).max())
+    return G / peak, float(t.b()[0]) / peak ** 2
+
+
+def amplitude(sh, b):
+    """The peak gradient (T/m) the class plays to reach ``b`` (s/mm^2)."""
+    return float(np.sqrt(max(float(b), 0.0) * 1e6 / unit_b(sh)[1]))
+
+
+def play(sh, bvals_SI, dirs):
+    """The ScannerSequence of rows on the class: its shape at each row's amplitude along each row's unit direction
+    (``bvals_SI`` in s/m^2, ``dirs`` ``(n, 3)``); a b = 0 row plays nothing."""
+    from dmipy_replay.acquisition.scanner_sequence import ScannerSequence
+    t = template(sh)
+    prof, unit = unit_b(sh)
+    b = np.asarray(bvals_SI, np.float64).ravel()
+    u = np.asarray(dirs, np.float64).reshape(-1, 3)
+    n = np.linalg.norm(u, axis=1)
+    u = np.where(n[:, None] > 0, u / np.where(n > 0, n, 1.0)[:, None], 0.0)
+    amp = np.sqrt(np.maximum(b, 0.0) / unit)
+    G = prof[None, :, None] * (amp[:, None] * u)[:, None, :]
+    return ScannerSequence(G=G.astype(np.float32), dt=t.dt, rf=t.rf, readout=t.readout)
+
+
+def describe(sh):
+    """The class's timing in a few words, for the gradient table."""
+    k = kind(sh)
+    if k == "tensor":
+        return "b-tensor class (three shapes: linear, planar, spherical at any rotation)"
+    if k == "ogse":
+        return f"cosine {float(sh['frequency']):g} Hz, {float(sh['duration']) * 1e3:g} ms per side"
+    if k == "pgste":
+        return f"STE δ {float(sh['delta']) * 1e3:g} / Δ {float(sh['Delta']) * 1e3:g} ms"
+    return f"δ {float(sh['delta']) * 1e3:g} / Δ {float(sh['Delta']) * 1e3:g} ms"
+
+
+def n_terms(sh):
+    """The class's stored terms: a PGSE, stimulated-echo or OGSE class two (one per encoding window), the b-tensor class
+    six (its three gradient shapes in both windows)."""
+    return 6 if kind(sh) == "tensor" else 2

space/config.toml

@@ -3,14 +3,20 @@
 source = "disco"                                 # space/sources/disco.py: the DiSCo phantom (demo: its layout, full: its columnar pack)
 
 [data]
-repo = "SubstrateCommons/disco-replay"          # the public replay dataset
-moments = "disco/moments"                        # the shape-moment layout written by dmipy-sim's write_shape_moments
-revision = "2db4e7839bc21c52c49d2aa7d399dde13be4cca7"   # the dataset revision the layout is read at (with the background moments, disco-space#29)
-source_manifest_sha256 = "499e411b7f7412ffa72a8885ee22a7bed9198e992a517b81c273ceeaeafd629e"   # the columnar manifest the layout was contracted from (its manifest.json names it)
+repo = "SubstrateCommons/disco-time-segments"   # the public time-segment layout of the DiSCo v2 walk
+layout = "disco"                                 # the layout's prefix in the dataset (dmipy_replay.replay.time_segments, format time_segments/2)
+revision = "a38679f061f1487b6688ebb37b7ee949da000ea3"   # the dataset revision the layout is read at (the OGSE class cos-f50-T20 and the b-tensor class qte-now-80)
+manifest_sha256 = "0fe044fd6963de74c6174276e35dd2ddc219d9d4ef2a476df0150ada48e23ef9"   # the layout manifest the page was accepted on (checked at load)
+pack = "hf://SubstrateCommons/disco-replay/disco"         # the replay pack the layout was contracted from: the acceptance test's reference
+pack_revision = "b06eb4179a99ef885af73de1ddfbc2f3c230de37"
+strands = "substrate/DiSCo_Strands_Trajectories.tck"   # the ground truth's strand centrelines (CC BY 4.0) in the replay dataset, read at pack_revision
 
-# The timing classes ("shapes") the layout holds: square pulses (slew_rate = inf) at TE = 53.5 ms, PGSE unless
-# kind = "pgste" (a stimulated echo storing for TM = Delta - delta). delta / Delta / TE in seconds, the same numbers the
-# layout's manifest records (checked at load); label is what the app shows.
+# The timing classes ("shapes"), each ONE waveform the layout stores (space/classes.py builds it with dmipy-replay's
+# sequence builders on the layout's grid, TE = 53.5 ms, 1000 samples; checked against the layout's partition at load): square
+# pulses (slew_rate = inf), PGSE unless kind = "pgste" (a stimulated echo storing for TM = Delta - delta) or "ogse"
+# (a cosine of `frequency` Hz over `duration` s on each side of the 180, its edges ramped at `slew_rate` T/m/s as built
+# at `b_built` s/mm^2: the class stores that one shape at every b). delta / Delta / TE in seconds (an OGSE's delta /
+# Delta are the scheme's nominal block and block separation); label is what the app shows.
 [shapes."d10.2-D16.7"]
 label = "DiSCo δ 10.2 / Δ 16.7 ms (shells b = 0, 1000, 1925)"
 delta = 0.0102
@@ -49,6 +55,31 @@ delta = 0.0076
 Delta = 0.0459
 TE = 0.0535
 kind = "pgste"
+# The paper's cosine OGSE (experiments/disco/disco_acq.py "ogse"): 50 Hz, one 20 ms period each side of the 180.
+[shapes."cos-f50-T20"]
+label = "cosine OGSE 50 Hz, one 20 ms period each side of the 180 (edges 200 T/m/s, as built at b = 1000)"
+kind = "ogse"
+frequency = 50.0
+duration = 0.020
+slew_rate = 200.0
+b_built = 1000.0
+delta = 0.019975
+Delta = 0.020029
+TE = 0.0535
+
+# The b-tensor class: three gradient shapes on the spin echo, each in both encoding windows (six stored terms), so one
+# class plays every linear, planar and spherical encoding of them at any rotation and b (space/tensor.py). The shapes
+# are the three axes of NOW's maximum-b spherical encoding (Sjolund et al. 2015) at TE 53.5 ms designed at 80 mT/m and
+# 200 T/m/s per axis (a 4 ms window for the 180), the physical gradient in `waveform` (data/); rotated, a linear or
+# spherical row at b 2000 needs about 124 mT/m on one axis: Connectome-class hardware. delta / Delta are the scheme's
+# nominal encoding windows (the CSD and the scheme never read them: a b-tensor run fits maps, not FODs).
+[shapes."qte-now-80"]
+label = "b-tensor class qte-now-80: NOW spherical encoding at TE 53.5 ms, 80 mT/m and 200 T/m/s per axis (Connectome-class when rotated)"
+kind = "tensor"
+waveform = "classes/qte-now-80.npy"
+delta = 0.0229
+Delta = 0.0253
+TE = 0.0535
 
 # DiSCo's own protocol: which timing class each shell of the dataset's gradient table plays (b in s/mm^2, the
 # table's rows are assigned to the nearest b here).
@@ -61,20 +92,35 @@ shells = [
   { b = 13183, shape = "d17.7-D35.8" },
 ]
 
-# Presets the app offers besides "DiSCo 364" (which comes from the dataset's table). b in s/mm^2.
+# Presets the app offers besides "DiSCo 364" (which comes from the dataset's table): the paper's flagship
+# acquisitions (experiments/disco/disco_acq.py), each its gradient table in data/protocols/ (the same files the paper
+# replays: b = 0 rows first, then every shell's directions spread by electrostatic repulsion) on one stored class.
 [presets."clinical b1000 x 30"]
-n_b0 = 1
-shells = [{ shape = "d12-D24", b = 1000, n_dirs = 30 }]
+table = "protocols/clinical"
+shape = "d12-D24"
 [presets."research 3-shell x 90"]
-n_b0 = 3
-shells = [
-  { shape = "d17-D30", b = 1000, n_dirs = 90 },
-  { shape = "d17-D30", b = 2000, n_dirs = 90 },
-  { shape = "d17-D30", b = 3000, n_dirs = 90 },
-]
+table = "protocols/research"
+shape = "d17-D30"
 [presets."Connectome 2.0 b6000 x 60"]
-n_b0 = 2
-shells = [{ shape = "d8-D20", b = 6000, n_dirs = 60 }]
+table = "protocols/c2"
+shape = "d8-D20"
+[presets."cosine OGSE 50 Hz b1000 x 60"]
+table = "protocols/ogse"
+shape = "cos-f50-T20"
+# Tensor-valued encoding: linear and spherical b-tensors interleaved at four b-values on the b-tensor class, the
+# b-values and direction counts of Szczepankiewicz et al. 2019's clinical protocol (6, 6, 10, 16 per encoding); the
+# same axes for both encodings at each b; one b = 0 row. Its run fits FA (DTI on the linear rows), MD and the
+# microscopic FA (the powder averages), not a connectome.
+[presets."b-tensor: linear + spherical (Szczepankiewicz 2019), b 100-2000 s/mm²"]
+shape = "qte-now-80"
+n_b0 = 1
+snr = 100                                 # its own starting SNR: µFA is noise-limited (the other presets start at [noise] snr)
+tensor = { bvalues = [100, 700, 1400, 2000], directions = [6, 6, 10, 16], encodings = ["LTE", "STE"] }
+
+# The noise a run starts with: the SNR at M0 (a full water voxel before relaxation) every preset starts at unless it
+# states its own.
+[noise]
+snr = 30
 
 # The tissue and the scanner: the field presets the page offers (T), the default field,
 # and whether a run starts with the tiers on. Tissue defaults come from dmipy-sim's white-matter catalogue at the
@@ -85,10 +131,10 @@ default_field = 3.0
 default_on = true
 pools = ["intra", "extra"]          # the pack's seeded pools (its spec also names a myelin pool nobody was seeded in)
 
-# The scanner menu after "ideal": catalogue machines of dmipy-sim (ScannerLimits.of(key)), each applied with every term
+# The scanner menu after "ideal": catalogue machines of dmipy-replay (ScannerLimits.of(key)), each applied with every term
 # its entry carries (docs/scanner.md); the phantom sits at a distance from isocentre along a scanner axis (the panel's
 # placement, at most 7.9 cm: the whole 1 mm cube inside the Swoop's 8 cm anchor). The Swoop's own gradient needs the
-# layout's background moments (shape_moments.stamp_background); a layout without them refuses it by name.
+# layout's background vectors (segment_rows(background=True)); a layout without them refuses it by name.
 [scanners]
 machines = [
   { label = "Hyperfine Swoop 64 mT (permanent magnet, field across the patient)", key = "hyperfine_swoop_64mT" },
@@ -99,8 +145,8 @@ default_offset_cm = 7.9
 max_offset_cm = 7.9
 delivery_per_meas = 0.0036   # s per measurement per replay: bore.delivered_moments over DiSCo's 31,802 live voxels, 1.3 s at 364 (a CPU host)
 
-# The mode: demo = the prebaked shape-moment layout (the hosted Spaces; seconds; the stored pulse
-# timings), full = the columnar replay pack itself (minutes per volume; any timing; the scheme upload), read from
+# The mode: demo = the time-segment layout (the hosted Spaces; seconds; the stored waveforms), full = the columnar
+# replay pack itself (minutes per volume; any timing; the scheme upload), read from
 # DISCO_COLUMNS (a directory with the columns, or the Hub at 45 MB/s). DISCO_MODE in the environment overrides.
 [mode]
 mode = "demo"
@@ -108,11 +154,12 @@ columns = "hf://SubstrateCommons/disco-replay/disco"
 workers = 16
 
 # The compute backend: jax (the reference; the Docker Space) or torch (Hugging Face's shared GPU pool, disco-space#2);
-# DISCO_BACKEND in the environment overrides. resident = false rebuilds the layout's device tensors on every call
-# (a pool that drops the device between calls).
+# DISCO_BACKEND in the environment overrides. resident = false streams the layout's packed words to the device a
+# chunk at a time on every image (nothing kept: a pool that drops the device between calls, a GPU smaller than one
+# class's columns, a CPU host); true keeps the columns an image reads on the device across calls.
 [compute]
 backend = "jax"
-resident = true
+resident = false
 
 [tracking]
 density = 4          # seeds per region voxel = density^3 (4 -> 659,840 seeds, the DiSCo score's setting)
@@ -124,7 +171,8 @@ max_steps = 500
 [gate]
 pearson_count_min = 0.90
 missed_pairs_max = 0
-reference_max_abs_diff = 1e-3     # the replay of DiSCo 364 against disco/reference/disco_replay_bare.nii.gz
+reference_max_abs_diff = 1e-3     # the layout's signal against the replay pack's own replay of the same rows at the gate's voxels, over the bare b = 0 signal
+reference_voxels = 32             # voxels spread over the pack's index the reference is replayed at
 
 # The reconstruction: DiSCo is single-fibre white matter, so its CSD is the single-tissue Tournier 2007 one on the
 # response white_matter_response_tournier07 estimates from the volume. DISCO_RECONSTRUCTION in the environment overrides.

space/limits.py

@@ -52,7 +52,7 @@ def cap_threads(quota=None):
 def cap_affinity(quota=None):
     """Restricts this process to the first ``quota`` CPUs it may run on (``os.sched_setaffinity``; ``quota`` overrides
     the cgroup read), so that every pool that sizes itself on the schedulable CPUs (XLA's, and torch's) matches the
-    quota; returns the CPUs it now runs on, or None when there is no quota or the process already runs on no more
+    quota; returns the CPUs it runs on, or None when there is no quota or the process already runs on no more
     than it. Children (the warm-up's process, the forked GPU worker) inherit it."""
     n = cpu_quota() if quota is None else quota
     allowed = sorted(os.sched_getaffinity(0))

space/pipeline.py

@@ -97,13 +97,15 @@ def backend(cfg):
 class Shell:
     """One shell: the timing class it plays, its b-value (s/mm^2) and its direction count. A shell made by
     :meth:`free` carries its own pulse timing (``delta``, ``Delta``, ``TE`` in seconds) instead of a stored class,
-    and ``shape`` names it."""
+    and ``shape`` names it. On a b-tensor class (:mod:`space.tensor`) ``encoding`` is the b-tensor's shape, ``LTE``
+    (linear) or ``STE`` (spherical), and a direction is the encoding's axis."""
     shape: str
     b: float
     n_dirs: int
     delta: Optional[float] = None
     Delta: Optional[float] = None
     TE: Optional[float] = None
+    encoding: Optional[str] = None
 
     @classmethod
     def free(cls, b, n_dirs, delta, Delta, TE):
@@ -145,16 +147,26 @@ class Protocol:
     def n_meas(self):
         return self.n_b0 + sum(s.n_dirs for s in self.shells)
 
+    @property
+    def tensor(self):
+        """Whether the shells are b-tensor encodings (every shell names one; a protocol does not mix them)."""
+        enc = [s.encoding is not None for s in self.shells]
+        if any(enc) and not all(enc):
+            raise ValueError("a protocol's shells are all b-tensor encodings or none")
+        return all(enc)
+
 
 @dataclass(frozen=True)
 class Measurements:
-    """Per row: b (s/mm^2), unit direction, timing-class name, delta / Delta / TE (s)."""
+    """Per row: b (s/mm^2), unit direction, timing-class name, delta / Delta / TE (s), and on a b-tensor protocol the
+    encoding (``b0``, ``LTE``, ``STE``; None otherwise)."""
     bvals: np.ndarray
     dirs: np.ndarray
     shape: np.ndarray
     delta: np.ndarray
     Delta: np.ndarray
     TE: np.ndarray
+    encoding: Optional[np.ndarray] = None
 
     @property
     def b0(self):
@@ -185,7 +197,10 @@ def measurements(protocol, shapes):
         b = np.asarray(protocol.bvals, np.float64)
     delta = np.array([timing[n]["delta"] for n in names]); Delta = np.array([timing[n]["Delta"] for n in names])
     TE = np.array([timing[n]["TE"] for n in names])
-    return Measurements(b, dirs, np.asarray(names), delta, Delta, TE)
+    enc = None
+    if protocol.tensor:
+        enc = np.asarray(["b0"] * protocol.n_b0 + [s.encoding for s in protocol.shells for _ in range(s.n_dirs)])
+    return Measurements(b, dirs, np.asarray(names), delta, Delta, TE, enc)
 
 
 def protocol_from_rows(bvals, dirs, group, make_shell, *, name):
@@ -230,15 +245,41 @@ def protocol_from_scheme(path, *, name="uploaded scheme"):
     return protocol_from_rows(b, np.asarray(enc.gradient_directions, np.float64), group.ravel(), make_shell, name=name)[0]
 
 
+def preset_snr(cfg, name):
+    """The SNR a preset starts at: its own ``snr`` when the preset states one (the b-tensor preset: µFA is
+    noise-limited), else ``[noise] snr``."""
+    p = cfg.get("presets", {}).get(name) or {}
+    return float(p.get("snr", cfg["noise"]["snr"]))
+
+
 def preset_protocol(cfg, name):
-    """A ``[presets]`` entry of the config as a Protocol."""
+    """A ``[presets]`` entry of the config as a Protocol: a stored gradient table (``table``: ``data/<table>.bvals`` and
+    ``.bvecs``, every row's own b and direction, each shell its rounded b, all on the class ``shape``), or ``shells``
+    of Fibonacci directions after ``n_b0`` b = 0 rows."""
     p = cfg["presets"][name]
+    if "tensor" in p:                                  # b-tensor encodings interleaved per b, the same axes for each
+        from . import tensor as T
+        t = p["tensor"]
+        shells, dirs = [], [np.tile([0.0, 0.0, 1.0], (int(p.get("n_b0", 1)), 1))]
+        for b, n in zip(t["bvalues"], t["directions"]):
+            for e in t["encodings"]:
+                shells.append(Shell(p["shape"], float(b), int(n), encoding=e)); dirs.append(T.fibonacci(int(n)))
+        return Protocol(tuple(shells), n_b0=int(p.get("n_b0", 1)), directions=np.concatenate(dirs), name=name)
+    if "table" in p:
+        bv = np.loadtxt(os.path.join(DATA_DIR, p["table"] + ".bvals")).ravel()
+        dirs = np.loadtxt(os.path.join(DATA_DIR, p["table"] + ".bvecs"))
+        dirs = dirs.T if dirs.shape[0] == 3 else dirs
+        shell = np.round(bv, -1)
+        return protocol_from_rows(bv, dirs, shell, lambda g, rows: Shell(p["shape"], float(g), int(len(rows))), name=name)[0]
     return Protocol(tuple(Shell(s["shape"], float(s["b"]), int(s["n_dirs"])) for s in p["shells"]), n_b0=int(p["n_b0"]), name=name)
 
 
 def retime(protocol, shape):
     """``protocol`` with every shell on the timing class ``shape`` (b-values, directions and counts kept): the A/B
-    knob that swaps a PGSE class for its stimulated-echo twin, or one δ / Δ for another."""
+    knob that swaps a PGSE class for its stimulated-echo twin, or one δ / Δ for another. A b-tensor protocol's shells
+    are encodings of its class, not timings: it is refused by name."""
+    if protocol.tensor:
+        raise ValueError("a b-tensor acquisition's shells are encodings of its one class: a pulse-timing knob does not apply")
     return Protocol(tuple(Shell(shape, s.b, s.n_dirs) for s in protocol.shells), n_b0=protocol.n_b0,
                     directions=protocol.directions, bvals=protocol.bvals, name=f"{protocol.name} on {shape}")
 
@@ -258,6 +299,26 @@ def write_volumes(res, out_dir, *, prefix, affine):
     return paths
 
 
+def write_tensor_volumes(res, out_dir, *, prefix, affine):
+    """A b-tensor run's files: the DWI as NIfTI with its bvals (s/mm^2), bvecs (each row's encoding axis, 3 x N) and
+    ``.btensor`` (each row's encoding, ``b0`` / ``LTE`` / ``STE``, one per line), and the FA, MD (µm²/ms) and µFA maps as
+    NIfTI (float32 on ``affine``, NaN as 0); returns the paths."""
+    import nibabel as nib
+    os.makedirs(out_dir, exist_ok=True)
+    aff = np.asarray(affine, np.float64); tm = res.extras["tensor"]; out = []
+    p = os.path.join(out_dir, f"{prefix}_dwi.nii.gz"); nib.save(nib.Nifti1Image(np.nan_to_num(res.dwi).astype(np.float32), aff), p); out.append(p)
+    p = os.path.join(out_dir, f"{prefix}.bvals"); np.savetxt(p, res.meas.bvals[None, :], fmt="%.1f"); out.append(p)
+    p = os.path.join(out_dir, f"{prefix}.bvecs"); np.savetxt(p, res.meas.dirs.T, fmt="%.8f"); out.append(p)
+    p = os.path.join(out_dir, f"{prefix}.btensor")
+    with open(p, "w") as f:
+        f.write("\n".join(str(e) for e in res.meas.encoding) + "\n")
+    out.append(p)
+    for key in ("fa", "md", "ufa"):
+        p = os.path.join(out_dir, f"{prefix}_{key}.nii.gz")
+        nib.save(nib.Nifti1Image(np.nan_to_num(np.asarray(tm[key], np.float64)).astype(np.float32), aff), p); out.append(p)
+    return out
+
+
 # ---- the physics: tissue and scanner --------------------------------------------------------------------------------
 CATALOGUE_FIELDS = (1.5, 3.0, 7.0)                       # where dmipy-sim's white-matter catalogue has cited relaxation
 CATALOGUE_POOLS = ("intra", "extra", "myelin")           # the pools the catalogue has T2 and T1 for
@@ -295,7 +356,7 @@ class Physics:
     """The tissue and the scanner a replay is evaluated at: the field (T) and its direction in the substrate frame,
     T2 and T1 per seeded pool (s), the walls' surface relaxivity ``rho2`` (m/s), the sheath's ``chi_iso`` and
     ``chi_aniso`` (SI, the field source); ``relaxation`` / ``contact`` / ``field`` switch the three tiers, so a tier
-    is a knob the page can turn off one at a time. :meth:`tissue` is the :class:`dmipy_sim.spec.tissue.Tissue` for
+    is a knob the page can turn off one at a time. :meth:`tissue` is the :class:`dmipy_replay.spec.tissue.Tissue` for
     the replay (None when every tier is off: bare diffusion). ``scanner`` is the machine (a catalogue key; None: the
     ideal scanner), whose field the physics is at, and ``offset_m`` where the phantom's centre sits from isocentre
     (metres, the scanner's frame): together they say what gradient every voxel is delivered (docs/scanner.md)."""
@@ -353,7 +414,7 @@ class Physics:
     def tissue(self, unseeded=()):
         """The Tissue for the replay; ``unseeded`` names the substrate's pools no walker was seeded in, which the
         library's mapping still wants a value for: they take the first seeded pool's (nothing reads it)."""
-        from dmipy_sim.spec.tissue import Tissue
+        from dmipy_replay.spec.tissue import Tissue
         if self.bare:
             return None
         first = self.pools[0]
@@ -462,7 +523,8 @@ class Source:
     The page's side: :meth:`describe` (the title and every text that names the source), :meth:`presets` and
     :meth:`protocol` (the acquisitions offered), :meth:`panel` (the tissue-and-scanner inputs),
     :meth:`catalogue_numbers`, :meth:`physics_from`, :meth:`knobs` and :meth:`apply_knob` (B, one knob away from A),
-    :meth:`refused_machines` (the catalogued machines kept off the scanner menu, and why),
+    :meth:`refused_machines` (the catalogued machines kept off the scanner menu, and why), :meth:`unavailable` (the
+    page's text when the inputs the configuration names cannot be read),
     :meth:`tracking_controls` and :meth:`estimated_seconds` (the GPU seconds a run reserves on a shared pool).
 
     The run's side: :attr:`regions` and :attr:`affine`, :meth:`validate`, :meth:`prepare` (the share of a replay
@@ -488,6 +550,12 @@ class Source:
         page shows under the menu. None here."""
         return {}
 
+    @classmethod
+    def unavailable(cls, cfg):
+        """``(title, text)`` when the inputs the configuration names cannot be read, the text naming each and why, else
+        None: the Space then serves that text as its page. Read at startup, before the page is built. None here."""
+        return None
+
     @classmethod
     def describe(cls, cfg):
         """The texts that name the source: ``title``, ``heading`` (the page's first paragraph), ``acquisition``
@@ -785,7 +853,7 @@ def csd(dwi, meas, mask, backend="jax"):
     t0 = time.perf_counter()
     data = np.nan_to_num(dwi, nan=0.0)
     sch = scheme(meas)
-    S0, response, _ = white_matter_response_tournier07(sch, data[mask])
+    S0, response, _ = white_matter_response_tournier07(sch, data[mask], backend=backend)
     mc = MultiCompartmentSphericalHarmonicsModel(models=[response], sh_order=SH_ORDER)
     fitted = mc.fit(sch, data, mask=mask, solver=f"csd_tournier07_{backend}", verbose=False)
     return np.asarray(fitted.fitted_parameters["sh_coeff"], np.float64), time.perf_counter() - t0, (S0, response)
@@ -983,7 +1051,7 @@ def b0_snr(res):
 def gradient_needed(b, delta, Delta):
     """The square-pulse PGSE amplitude (T/m) that gives ``b`` (s/mm^2) at ``delta`` / ``Delta`` (s):
     b = γ² G² δ² (Δ − δ/3)."""
-    from dmipy_sim.acquisition.scanners import GAMMA
+    from dmipy_analytic.constants import GAMMA
     return float(np.sqrt(b * 1e6 / (GAMMA ** 2 * delta ** 2 * (Delta - delta / 3.0))))
 
 
@@ -992,7 +1060,7 @@ IDEAL = "ideal: no scanner terms, any gradient amplitude"
 
 def machines(cfg):
     """The configuration's scanner menu after :data:`IDEAL`: ``{label: catalogue key}`` (``[scanners] machines``), each a
-    machine of dmipy-sim's catalogue whose every catalogued term the page applies (docs/scanner.md)."""
+    machine of dmipy-replay's catalogue whose every catalogued term the page applies (docs/scanner.md)."""
     return {m["label"]: m["key"] for m in (cfg.get("scanners") or {}).get("machines", [])}
 
 
@@ -1007,21 +1075,26 @@ def machine(cfg, label):
 
 
 def limits(key):
-    """The :class:`dmipy_sim.acquisition.scanners.ScannerLimits` of a catalogue key."""
-    from dmipy_sim.acquisition.scanners import ScannerLimits
+    """The :class:`dmipy_replay.acquisition.scanners.ScannerLimits` of a catalogue key."""
+    from dmipy_replay.acquisition.scanners import ScannerLimits
     return ScannerLimits.of(key)
 
 
 def playable(protocol, shapes, scanner=None):
     """Per shell, the gradient it needs and whether ``scanner`` (a catalogue key, or None for no limit) can play it:
     rows of ``(shell, b, delta, Delta, G_needed, G_max or None, ok)`` in SI. A stimulated-echo class needs the same
-    amplitude as its PGSE twin (the same δ and diffusion time)."""
+    amplitude as its PGSE twin (the same δ and diffusion time); an OGSE class the peak of its own waveform
+    (:func:`space.classes.amplitude`; its delta / Delta are the scheme's nominal block and separation)."""
+    from . import classes as C
     G_max = float(limits(scanner).G_max) if scanner else None
     rows = []
     for s in protocol.shells:
         t = dict(delta=s.delta, Delta=s.Delta) if s.free_timing else shapes[s.shape]
         delta, Delta = float(t["delta"]), float(t["Delta"])
-        G = gradient_needed(s.b, delta, Delta)
+        if t.get("kind") == "tensor":
+            G = C.tensor_peak(t, s.encoding, s.b, s.n_dirs)                 # per axis, over the shell's rotations
+        else:
+            G = C.amplitude(t, s.b) if t.get("kind") == "ogse" else gradient_needed(s.b, delta, Delta)
         rows.append((s.shape, s.b, delta, Delta, G, G_max, G_max is None or G <= G_max))
     return rows
 
@@ -1041,7 +1114,7 @@ def scanner_terms(key):
     """The rows of a machine's term table: ``(term, applied, where the number comes from)``. ``applied`` is whether
     the catalogue carries the term (the page applies every one it carries); the source is the catalogue leaves'
     confidence as the page names it (:data:`CONFIDENCE`) with their citation keys, or why the term is absent."""
-    from dmipy_sim.acquisition import scanner_constants as scc
+    from dmipy_replay.acquisition import scanner_constants as scc
     L = limits(key)
     _key, entry, _kind = scc.resolve(key)
     rows = []
@@ -1120,6 +1193,9 @@ class Result:
     extras: dict = field(default_factory=dict)
 
 
+TENSOR_STAGES = ("replay", "noise", "maps")              # a b-tensor protocol's stages
+
+
 def run_stages(source, protocol, *, snr=None, tracking=Tracking(), noise_seed=0, physics=None, prepared=None):
     """The pipeline as a generator: before each of the source's ``STAGES`` it yields ``(stage, k, n)`` (so a page
     can show where the run is, from inside any worker), and last the :class:`Result`. ``physics`` is the tissue and
@@ -1131,10 +1207,20 @@ def run_stages(source, protocol, *, snr=None, tracking=Tracking(), noise_seed=0,
     def at(stage):
         print(f"[pipeline] {stage} at +{time.perf_counter() - t_run:.1f} s", flush=True)      # the server log shows where a run is
         return stage, stages.index(stage), len(stages)
+    if protocol.tensor:                                  # a b-tensor protocol: microstructure maps, not a connectome
+        stages = TENSOR_STAGES
     meas = source.validate(protocol, physics)
     yield at("replay"); dwi, floor, s0_factor, t_replay = source.replay(meas, physics, prepared)
     yield at("noise"); t0 = time.perf_counter(); noisy = add_noise(dwi, snr, s0_factor, seed=noise_seed, backend=source.backend); t_noise = time.perf_counter() - t0
     signal = np.isfinite(dwi[..., 0])
+    if protocol.tensor:
+        from . import tensor as T
+        yield at("maps"); t0 = time.perf_counter(); tmaps = T.maps(noisy, meas, signal)
+        seconds = dict(replay=t_replay, noise=t_noise, maps=time.perf_counter() - t0)
+        seconds["total"] = sum(seconds.values())
+        yield Result(protocol, meas, noisy, floor, s0_factor, None, None, None, None, None, snr=snr, physics=physics, clean=dwi,
+                     seconds=seconds, extras=dict(tensor=tmaps))
+        return
     yield at("csd"); sh, t_csd, extras = source.reconstruct(noisy, s0_factor, meas, signal)
     yield at("track"); t0 = time.perf_counter(); fld, seeds = source.tracking_inputs(sh, tracking.density); tg, _ = track(fld, seeds, tracking, source.backend); t_track = time.perf_counter() - t0
     seconds = dict(replay=t_replay, noise=t_noise, **t_csd, track=t_track)

space/reference.py

@@ -0,0 +1,42 @@
+"""The acceptance gate's reference: the replay pack the layout was contracted from, replayed walker by walker at a set
+of voxels (only those voxels' rows are read, by range, from the Hub or a local copy). This is the full replay, not a
+rendering: every stored band, every field mode, the contact, the pools; the layout's image of the same rows must equal
+it to ``[gate] reference_max_abs_diff`` (S/S0)."""
+from __future__ import annotations
+
+import numpy as np
+
+
+def open_pack(uri, workers=8):
+    """The columnar replay pack at ``uri`` (``hf://owner/name/prefix`` or a directory)."""
+    from dmipy_replay.replay.columnar import ColumnarPack
+    return ColumnarPack(uri, workers=workers)
+
+
+def gate_voxels(pack, n):
+    """``n`` voxels spread evenly over the pack's index (its (voxel, pool) rows in order), as ijk triples: the spot
+    check's rule."""
+    keys = sorted({tuple(r["ijk"]) for r in pack.index["rows"]})
+    pick = np.unique(np.linspace(0, len(keys) - 1, min(int(n), len(keys))).round().astype(int))
+    return [keys[i] for i in pick]
+
+
+def replay_at(pack, voxels, sequence, tissue=None, scanner=None):
+    """``(n_voxels, n_meas)``: each voxel's weighted ensemble magnitude under ``sequence`` (normalised by its walkers'
+    weight, as the layout's image is), from the pack's own rows at ``voxels`` (every band; the field modes and the
+    contact read when the tissue needs them)."""
+    from dmipy_replay.phantom.grid import Grid
+    grid = Grid.from_meta(pack.index["grid"])
+    need_field = tissue is not None and (getattr(tissue, "chi_iso", None) or getattr(tissue, "chi_aniso", None))
+    view = pack.view(voxels=voxels, contact=tissue is not None,
+                     modes=(pack.path_groups[-1] if need_field and pack.path_groups else 0))
+    w, ew, E = view.walker_signals(sequence, tissue=tissue, scanner=scanner, tol=None)
+    ijk, _ = grid.bin(view.r0)
+    v = np.ravel_multi_index(ijk.T, tuple(grid.shape))
+    flat = np.array([np.ravel_multi_index(k, tuple(grid.shape)) for k in voxels])
+    pos = {int(f): i for i, f in enumerate(flat)}
+    idx = np.array([pos[int(x)] for x in v])
+    num = np.zeros((len(voxels), E.shape[1]), complex)
+    np.add.at(num, idx, np.asarray(ew)[:, None] * np.asarray(E))
+    den = np.bincount(idx, np.asarray(w, np.float64), minlength=len(voxels))
+    return np.abs(num / den[:, None])

space/sources/brain.py

@@ -4,14 +4,14 @@ replay packs, replayed on the device, reconstructed, tracked and scored against
 The input is an **asset** (the directory layout of ``SubstrateCommons/masivar-brain``: ``manifest.json``,
 ``fod_wm.npy``, ``fractions.npy``, ``mask.npy``, ``labels.npy``, ``regions.json``, ``stop_mask.npy``,
 ``mean_b0.npy``): the multi-tissue CSD of a real scan, made once offline, with its parcellation. The replay is the
-kernel route of dmipy-sim's phantom (RPH.md 6): per voxel and measurement
+kernel route of dmipy-replay's phantom (RPH.md 6): per voxel and measurement
 
     S = | f_wm m0_wm FOD . C_wm  +  f_gm m0_gm E_gm  +  f_csf m0_csf E_csf |
 
 where ``C_wm`` is the WM pack's pose response contracted to the FOD's harmonics
-(:meth:`dmipy_sim.replay.ReplayPack.pose_responses` at ``keep=(8, 0)`` through :func:`dmipy_sim.replay.so3.axis_density_coeffs`),
+(:meth:`dmipy_replay.ReplayPack.pose_responses` at ``keep=(8, 0)`` through :func:`dmipy_replay.replay.so3.axis_density_coeffs`),
 ``E_gm`` the GM pack's response to an isotropic orientation distribution, ``E_csf`` free water in closed form
-(:class:`dmipy_sim.phantom.FreeWater`), each at the run's tissue (T2 / T1 per pool at the echo, the WM walls' surface
+(:class:`dmipy_replay.phantom.FreeWater`), each at the run's tissue (T2 / T1 per pool at the echo, the WM walls' surface
 relaxivity, the myelin sheath's susceptibility at the scanner's field and its direction), and ``m0`` the proton
 density per tissue. The pose responses (:meth:`Brain.prepare`: seconds on the CPU) are cached in the page's process per
 pack, protocol, tissue and field; a run looks them up on the page and computes the missing ones inside the GPU call,
@@ -46,15 +46,22 @@ N_COEF = 45
 REACH_MM = 250.0                                          # the longest streamline the tracker must be able to reach
 PROFILE_ANGLES = np.linspace(0.0, 90.0, 19)               # degrees between the gradient and a fibre, for the response plot
 C00 = 1.0 / (2.0 * np.sqrt(np.pi))                        # the l = 0 coefficient of a unit-integral density
+POSE_BACKEND = "jax"                                      # where the pose expansion's walker work runs: JAX on the CPU, which the [budget] is measured on
 #: why a magnet with its own gradient (a field law: the Swoop) is not on the brain's scanner menu, as the page says it
 FIELD_LAW_REFUSAL = ("Its own gradient (the permanent magnet's field law, 1.4 mT/m at 8 cm) is on through every pulse and dead "
-                     "time, so every voxel plays a two-direction waveform of its own. dmipy-sim's closed form expands that "
+                     "time, so every voxel plays a two-direction waveform of its own. dmipy-replay's closed form expands that "
                      "waveform exactly; the cost keeps it off this menu. Binned to the menu's 1 % of b, this head needs 9,752 "
                      "encoding classes for the scan's protocol (665 at 3 %), and one class's pose responses take 2-3 minutes "
                      "(WM 104-166 s, GM 4-21 s at 485 measurements on a GH200), so one protocol at one tissue and field is about "
                      "two weeks of one device. The DiSCo Space plays the Swoop exactly.")
 
 
+#: what the page says when the brain cannot be composed from the inputs its configuration names
+UNAVAILABLE = ("**The brain page is not available: its packs are being re-walked for dmipy 3.0.** It composes a brain from one "
+               "white-matter pack and one grey-matter pack and refuses to run until every input its configuration names is one "
+               "this release serves. What is absent or refused now:")
+
+
 # ---- the asset ---------------------------------------------------------------------------------------------------------
 
 def asset_location(cfg):
@@ -94,7 +101,8 @@ def manifest(where):
 def declared_segments(uri):
     """The segment table ``{n, n_t, T}`` (RPK.md 4.3) a pack declares, read without its arrays: a Hub pack's
     record in its dataset's ``manifest.json`` (``hf://owner/name/path``), a local pack's header."""
-    from dmipy_sim.replay.publish import MANIFEST, header_of, is_hub_uri, parse_uri
+    from dmipy_replay.replay.hub import MANIFEST, is_hub_uri, parse_uri
+    from dmipy_sim.replay.publish import header_of
     if is_hub_uri(uri):
         from huggingface_hub import hf_hub_download
         repo, path = parse_uri(uri)
@@ -111,9 +119,61 @@ def declared_segments(uri):
     return dict(n=int(seg["n"]), n_t=int(seg["n_t"]), T=float(seg["T"]))
 
 
+def _is_404(e):
+    """Whether ``e`` is an HTTP client error whose response says 404 -- the Hub's answer that a repository, a revision or
+    a file is not there (``urllib`` errors carry ``code``, the Hub client's carry ``response``)."""
+    return getattr(e, "code", None) == 404 or getattr(getattr(e, "response", None), "status_code", None) == 404
+
+
+def asset_refusal(where):
+    """Why the asset at ``where`` cannot be read, or None when its manifest reads: a Hub revision (or dataset) the Hub does
+    not have, or a directory without ``manifest.json``. An error that says nothing about the asset itself (the Hub
+    unreachable) is raised, not named."""
+    hub = where.startswith("hf://")
+    try:
+        manifest(where)
+    except Exception as e:
+        if hub and _is_404(e):
+            repo, _, rev = where[len("hf://"):].partition("@")
+            return f"the Hub has no manifest.json of {repo} at revision {rev or 'main'}"
+        if not hub and isinstance(e, FileNotFoundError):
+            return f"{where} has no manifest.json"
+        raise
+    return None
+
+
+def pack_refusal(uri):
+    """Why this page cannot serve the pack at ``uri``, or None: the pack is absent from the Hub or from the disk, the reader
+    refuses its format (the reader's own message names what it does not read), or its window record is not the segment
+    table this page reads (:func:`declared_segments`). A Hub pack is read for its header and its dataset's record only;
+    an error that says nothing about the pack (the Hub unreachable) is raised, not named."""
+    from dmipy_replay import ReplayPack
+    from dmipy_replay.replay.hub import is_hub_uri
+    if is_hub_uri(uri):
+        try:
+            ReplayPack.open(uri)
+        except ValueError as e:
+            return f"the reader refuses it: {e}"
+        except Exception as e:
+            if _is_404(e):
+                return "it is not in its dataset on the Hub"
+            raise
+    elif not os.path.exists(uri):
+        return "there is no such file"
+    try:
+        declared_segments(uri)
+    except ValueError as e:
+        return f"this page cannot read its window record: {e}"
+    except Exception as e:
+        if _is_404(e):
+            return "its dataset has no manifest.json to read the window record from"
+        raise
+    return None
+
+
 def image_rotation(man):
-    """``R`` (image -> scanner) of the asset's grid (:meth:`dmipy_sim.phantom.Grid.from_oblique_affine`)."""
-    from dmipy_sim.phantom import Grid
+    """``R`` (image -> scanner) of the asset's grid (:meth:`dmipy_replay.phantom.Grid.from_oblique_affine`)."""
+    from dmipy_replay.phantom import Grid
     _, R = Grid.from_oblique_affine(np.asarray(man["grid"]["affine"], np.float64), tuple(man["grid"]["shape"]))
     return np.asarray(R, np.float64)
 
@@ -169,11 +229,11 @@ class Asset:
         mask = np.load(os.path.join(d, "mask.npy")).astype(bool)
         if fod.shape != shape + (N_COEF,) or mask.shape != shape:
             raise ValueError(f"the asset's FOD {fod.shape} and mask {mask.shape} are not on its grid {shape} with {N_COEF} coefficients")
-        from dmipy_sim.phantom.orientation import SH_BASES
+        from dmipy_replay.phantom.orientation import SH_BASES
         basis = str(man["fod"]["basis"]).split()[0]
         if SH_BASES.get(basis) != ("tournier07", False) or int(man["fod"]["lmax"]) != LMAX or man["fod"]["frame"] != "image":
-            raise ValueError(f"the asset's FOD is {man['fod']}: this source reads dmipy-sim's required basis (a basis whose "
-                             f"name, the first word, dmipy_sim.phantom.orientation.SH_BASES maps to orthonormal tournier07: "
+            raise ValueError(f"the asset's FOD is {man['fod']}: this source reads dmipy-replay's required basis (a basis whose "
+                             f"name, the first word, dmipy_replay.phantom.orientation.SH_BASES maps to orthonormal tournier07: "
                              f"tournier07, mrtrix3, dmipy-fit) at lmax {LMAX} in the image frame")
         if list(man["fractions"]) != list(TISSUES):
             raise ValueError(f"the asset's fractions are {man['fractions']}, this source reads {list(TISSUES)}")
@@ -231,7 +291,7 @@ class BrainPhysics:
     ``chi_iso`` / ``chi_aniso`` the myelin sheath's susceptibility (SI, the WM pack's field source), ``D_csf`` free
     water's diffusivity (m^2/s), ``m0`` the proton density per tissue relative to CSF, ``wm_pack`` / ``gm_pack`` the
     packs by their configuration label; ``relaxation``, ``contact`` and ``field`` switch the three tiers; ``scanner``
-    the machine (a key of dmipy-sim's catalogue): its slew rate shapes every pulse's ramps, its field and direction are
+    the machine (a key of dmipy-replay's catalogue): its slew rate shapes every pulse's ramps, its field and direction are
     the physics', and every voxel is played the gradient it delivers there with the head centre at isocentre
     (:meth:`Brain.prepare`); None is the ideal scanner, ideal pulses with vertical ramps. All off is bare diffusion at
     the same proton densities and packs."""
@@ -284,7 +344,7 @@ class BrainPhysics:
         (vertical ramps) on the ideal scanner."""
         if self.scanner is None:
             return np.inf
-        from dmipy_sim.acquisition.scanners import ScannerLimits
+        from dmipy_replay.acquisition.scanners import ScannerLimits
         return float(ScannerLimits.of(self.scanner).slew_max)
 
     def label(self):
@@ -296,13 +356,13 @@ class BrainPhysics:
         return f"{self.field_T:g} T along ({u[0]:.2f}, {u[1]:.2f}, {u[2]:.2f}), tiers {tiers}{grad}"
 
     def tissues(self, wm_pack, gm_pack):
-        """``(wm, gm, csf)`` :class:`dmipy_sim.spec.Tissue` for the packs ``wm_pack`` and ``gm_pack``: the WM pack's
+        """``(wm, gm, csf)`` :class:`dmipy_replay.spec.tissue.Tissue` for the packs ``wm_pack`` and ``gm_pack``: the WM pack's
         pools with their T2 / T1, the walls' rho2 and, on a WM pack whose spec declares its susceptibility field
         present (a magnetic myelin pool), the sheath's susceptibility; the GM pack's pools at grey matter's T2 / T1;
         free water's D, T2 and T1. A tier off leaves its values out (``None`` for a pack with nothing on: its bare
         diffusion). A WM pack whose field is absent has no sheath, so its tissue carries no susceptibility and its
         field tier is zero at any field."""
-        from dmipy_sim.spec import Tissue
+        from dmipy_replay.spec.tissue import Tissue
         wm_pools, gm_pools = Brain.pools(wm_pack), Brain.pools(gm_pack)
         relax, contact = self.relaxation, self.contact
         field = self.field and wm_pack.susceptibility_field == "present"
@@ -335,7 +395,7 @@ def specimen_pose(R, b0_direction):
     image -> scanner rotation and ``W`` the rotation taking the field's scanner-frame direction onto z (the bore's
     field axis), so the packs see the field at ``R.T b0_direction`` in the image frame. A sequence composed at this
     pose carries the gradients ``pose @ g`` (:meth:`Brain.sequence`)."""
-    from dmipy_sim.replay.so3 import rotation_of
+    from dmipy_replay.replay.so3 import rotation_of
     return rotation_of(np.asarray(b0_direction, np.float64)).T @ np.asarray(R, np.float64)
 
 
@@ -503,6 +563,27 @@ class Brain(P.Source):
             out[m["label"]] = FIELD_LAW_REFUSAL
         return out
 
+    @classmethod
+    def unavailable(cls, cfg):
+        """``(title, text)`` when the brain cannot be composed from the inputs the configuration names -- the asset's
+        manifest, or a pack of either menu, is absent, refused by the reader or not readable by this page
+        (:func:`asset_refusal`, :func:`pack_refusal`) -- the text naming each input and why; None when every input reads.
+        The Space serves this text as its page instead of failing to start."""
+        found = []
+        where = asset_location(cfg)
+        why = asset_refusal(where)
+        if why:
+            found.append(f"- the asset `{where}`: {why}")
+        for tissue in ("wm", "gm"):
+            for p in cfg["packs"][tissue]:
+                why = pack_refusal(p["uri"])
+                if why:
+                    found.append(f"- the {tissue.upper()} pack \"{p['label']}\" `{p['uri']}`: {why}")
+        if not found:
+            return None
+        title = cfg["describe"]["title"]
+        return title, "\n\n".join([f"# {title}", UNAVAILABLE, "\n".join(found)])
+
     @classmethod
     def pack_menu(cls, cfg, tissue):
         """The configuration's packs for ``tissue`` (``wm`` / ``gm``): ``{label: uri}`` in the menu's order."""
@@ -520,7 +601,7 @@ class Brain(P.Source):
         """``{class: k}``: per timing class of ``meas``, the number of the pack's windows its acquisition reaches.
         ``segments`` is a declared table ``{n, n_t, T}`` (:func:`declared_segments`) or a pack's meta (its
         ``walk_params.segments``). The acquisition lasts ``T_acq = (G.shape[1] - 1) dt`` of the class's sequence,
-        as :meth:`dmipy_sim.replay.ReplayPack._compile` reads it, and window ``i`` starts at ``i T``; a window whose
+        as :meth:`dmipy_replay.ReplayPack._span` reads it, and window ``i`` starts at ``i T``; a window whose
         first save sits at or beyond ``T_acq`` is not reached, so the save at a window's end belongs to that window
         alone. A pack of one window (``n = 1``) is reached as one window whatever the class; ``k``
         above ``n`` is an acquisition beyond the walk (:meth:`windows_needed` refuses it)."""
@@ -554,12 +635,12 @@ class Brain(P.Source):
 
     def pack(self, tissue, label, k):
         """The pack ``label`` of ``tissue``'s menu holding at least its first ``k`` windows: a Hub URI whose
-        declared table has ``n > 1`` loads ``windows=range(k)`` (:class:`dmipy_sim.phantom.PackSubstrate`; the pack
+        declared table has ``n > 1`` loads ``windows=range(k)`` (:class:`dmipy_replay.phantom.PackSubstrate`; the pack
         records them as ``windows_present``); a local path or a one-window pack loads whole. A later request for
         more windows than it holds reloads with the larger range (the library fetches only the new ranges) and
         replaces the cached pack."""
-        from dmipy_sim.phantom import PackSubstrate
-        from dmipy_sim.replay.publish import is_hub_uri
+        from dmipy_replay.phantom import PackSubstrate
+        from dmipy_replay.replay.hub import is_hub_uri
         key = (tissue, label)
         seg = self.segments(tissue, label)
         if key in self.packs and self.windows_held(self.packs[key], seg) >= int(k):
@@ -640,7 +721,7 @@ class Brain(P.Source):
                    "the head's own susceptibility and RF effects are not modelled",
                    "the head centre at the scanner's isocentre (a machine of the scanner menu plays every voxel the gradient it delivers there; docs/scanner.md)",
                    "no per-voxel split-half floor: the packs' certified floors are the accuracy"],
-            accuracy=("A brain replay is exact composition (the tests check it against dmipy-sim's own `Phantom.compose(...).replay`) of "
+            accuracy=("A brain replay is exact composition (the tests check it against dmipy-replay's own `Phantom.compose(...).replay`) of "
                       "the packs' responses, so its accuracy is the packs' certified Monte-Carlo floors; there is no per-voxel split-half "
                       "floor, because every voxel reads the same two walks."),
             ingredients=("proton density: Σ f·M0 (relative to CSF)", "relaxation: the b = 0 signal's T2 / T1 weight at the echo",
@@ -840,16 +921,16 @@ class Brain(P.Source):
         the image-frame directions turned by ``pose``, trapezoid pulses at the class's delta / Delta / TE with ramps at
         ``slew_rate`` (vertical at the default, infinite), a stimulated echo storing for TM = Delta - delta when the
         class is one."""
-        import dmipy_sim as d
+        from dmipy_replay.sequences import pgse, pgste
         names = np.unique(meas.shape[rows])
         if len(names) != 1:
             raise ValueError(f"one timing class per sequence, got {list(names)}")
         t = self.shapes[str(names[0])]
         dirs = (meas.dirs[rows] @ np.asarray(pose).T).tolist(); b = meas.bvals[rows] * 1e6
         if t.get("kind", "pgse") == "pgste":
-            return d.pgste(dirs, float(t["delta"]), float(t["Delta"]) - float(t["delta"]), bvalues=b, TE=float(t["TE"]), n_t=1000,
-                           slew_rate=float(slew_rate), ste_flip_angles=(90.0, 90.0, 90.0))
-        return d.pgse(dirs, float(t["delta"]), float(t["Delta"]), bvalues=b, TE=float(t["TE"]), n_t=1000, slew_rate=float(slew_rate))
+            return pgste(dirs, float(t["delta"]), float(t["Delta"]) - float(t["delta"]), bvalues=b, TE=float(t["TE"]), n_t=1000,
+                         slew_rate=float(slew_rate), ste_flip_angles=(90.0, 90.0, 90.0))
+        return pgse(dirs, float(t["delta"]), float(t["Delta"]), bvalues=b, TE=float(t["TE"]), n_t=1000, slew_rate=float(slew_rate))
 
     def _b0_sequence(self, name, pose, slew_rate=np.inf):
         """The class ``name``'s b = 0 measurement as a sequence at ``pose`` (:meth:`sequence`): what the b = 0 signal
@@ -895,10 +976,10 @@ class Brain(P.Source):
     @functools.cached_property
     def bore_grid(self):
         """``(grid, index)``: the brain's voxels as they sit in the bore, the head centre (the mask's centroid) at
-        isocentre: :meth:`dmipy_sim.phantom.Grid.from_oblique_affine` of the cropped asset's affine, whose rotation
+        isocentre: :meth:`dmipy_replay.phantom.Grid.from_oblique_affine` of the cropped asset's affine, whose rotation
         is :attr:`Asset.R` (the frame the page composes in), and ``index`` the brain voxels' own indices, each at
         ``affine @ [ijk, 1]`` in the scanner for a left- or right-handed affine."""
-        from dmipy_sim.phantom import Grid
+        from dmipy_replay.phantom import Grid
         A = np.asarray(self.affine, np.float64)
         centre = (A[:3, :3] @ self.vox.mean(0) + A[:3, 3]) * 1e-3
         grid, R = Grid.from_oblique_affine(A, self.mask.shape, isocenter_m=centre)
@@ -908,9 +989,9 @@ class Brain(P.Source):
 
     def encoding(self, meas, physics, rows):
         """``(class_of_voxel, played)``: the timing class ``rows`` of ``meas`` as the run's machine plays it at every brain
-        voxel, binned to ``[scanners] tolerance`` of b (:func:`dmipy_sim.phantom.bore.encoding_classes` on
+        voxel, binned to ``[scanners] tolerance`` of b (:func:`dmipy_replay.phantom.bore.encoding_classes` on
         :attr:`bore_grid`), the played acquisitions in the image frame; computed once per machine, slew and rows."""
-        from dmipy_sim.phantom.bore import encoding_classes
+        from dmipy_replay.phantom.bore import encoding_classes
         import hashlib
         h = hashlib.sha256()
         for x in (meas.bvals[rows], meas.dirs[rows], meas.shape[rows].astype(str)):
@@ -943,7 +1024,7 @@ class Brain(P.Source):
         or ``cold`` (a band to compile first); ``saves`` the WM pack's saves the entry's classes span
         (:meth:`saves_spanned` of :meth:`windows_needed`). ``[budget]`` prices them (``response_<state>``); on a machine
         an entry counts its rows once per encoding class times ``response_per_class`` (a class's rows each play their own
-        amplitude, so the shells no longer share the expansion's bodies)."""
+        amplitude, so the shells do not share the expansion's bodies)."""
         out = []; keys = set(); bands = set(self._compiled)
         for meas, ph in entries:
             key = self.response_key(meas, ph)
@@ -961,19 +1042,19 @@ class Brain(P.Source):
 
     def prepare(self, meas, physics=None):
         """The pose responses of a replay, on the CPU: per timing class of ``meas``, the WM and GM packs' pose
-        responses at the run's tissue, field and specimen pose (dmipy-sim's closed form) contracted to the FOD's
+        responses at the run's tissue, field and specimen pose (dmipy-replay's closed form) contracted to the FOD's
         harmonics, free water's closed form; the true response curves for the explorer and the b = 0 weight of each
         tier alone (both at the first class's shells, within the windows ``meas`` reaches). Loads the windows of the
         chosen packs the classes reach (:meth:`windows_needed`, :meth:`pack`); computed every call and kept in this
         process's cache under :meth:`response_key` (:meth:`cached` reads it)."""
-        from dmipy_sim.phantom import FreeWater
-        from dmipy_sim.replay import so3
+        from dmipy_replay.phantom import FreeWater
+        from dmipy_replay.replay import so3
         t0 = time.perf_counter(); stages = {}
         wm_pack = self.pack("wm", physics.wm_pack, self.windows_needed("wm", physics.wm_pack, meas))
         stages["load_wm"] = time.perf_counter() - t0
         gm_pack = self.pack("gm", physics.gm_pack, self.windows_needed("gm", physics.gm_pack, meas))
         stages["load_gm"] = time.perf_counter() - t0 - stages["load_wm"]
-        for name in np.unique(meas.shape):                                  # the walks and the WM pack's tiers, now loaded
+        for name in np.unique(meas.shape):                                  # the walks and the WM pack's tiers, loaded
             self.validate(P.Protocol((P.Shell(str(name), 1.0, 1),), name="check"), physics)
         t_wm, t_gm, t_csf = physics.tissues(wm_pack, gm_pack)
         pose = specimen_pose(self.asset.R, physics.b0_direction)
@@ -983,7 +1064,7 @@ class Brain(P.Source):
         # per timing class, the acquisitions the voxels are played: the commanded one, or one per encoding class of the
         # machine (its delivered gradient, in the image frame, turned to the specimen's pose)
         plan, seqs = [], []
-        from dmipy_sim.acquisition.waveforms import rotate_waveform
+        from dmipy_replay.acquisition.waveforms import rotate_waveform
         for rows in classes:
             if physics.scanner is None:
                 cls, played = np.zeros(len(self.vox), int), [self.sequence(meas, rows, pose, physics.slew_rate)]
@@ -994,9 +1075,9 @@ class Brain(P.Source):
         prof = self._profile_sequence(meas, classes[0], pose, physics.slew_rate)   # the explorer's curves ride in the same pass
         batch = seqs + ([prof[0]] if prof is not None else [])
         t = time.perf_counter()
-        r_wm = wm_pack.pose_responses(batch, tissue=t_wm, scanner=physics.scanner_field, pose=pose, keep=(LMAX, 0))
+        r_wm = wm_pack.pose_responses(batch, tissue=t_wm, scanner=physics.scanner_field, pose=pose, keep=(LMAX, 0), backend=POSE_BACKEND)
         stages["pose_wm"] = time.perf_counter() - t; t = time.perf_counter()
-        r_gm = gm_pack.pose_responses(batch, tissue=t_gm, scanner=physics.scanner_field, pose=pose, keep=(0, 0))
+        r_gm = gm_pack.pose_responses(batch, tissue=t_gm, scanner=physics.scanner_field, pose=pose, keep=(0, 0), backend=POSE_BACKEND)
         stages["pose_gm"] = time.perf_counter() - t
         free = FreeWater(m0=1.0, tissue=t_csf)
         K = len(seqs)
@@ -1044,7 +1125,7 @@ class Brain(P.Source):
         """The b = 0 signal per tissue with each tier alone, over bare: ``relaxation`` (T2 at the echo, T1 over a
         stimulated echo's storage) per tissue, ``contact`` (the WM walls' survival) for WM; a tier that is off is
         absent. The b = 0 signal has no orientation, so one pose of the pack answers it."""
-        from dmipy_sim.phantom import FreeWater
+        from dmipy_replay.phantom import FreeWater
         seq = self._b0_sequence(str(meas.shape[rows[0]]), pose, physics.slew_rate)
         out = {}
         alone = dict(relaxation=False, contact=False, field=False)
@@ -1148,7 +1229,7 @@ class Brain(P.Source):
     def tracking_inputs(self, sh, density):
         """``(field, seeds)``: the FOD turned from the image frame into the scanner frame on the asset's affine, the
         domain the white matter and the regions (leaving it ends a streamline), ``density^3`` seeds per WM voxel."""
-        from dmipy_sim.replay.so3 import rotate_sh
+        from dmipy_replay.replay.so3 import rotate_sh
         from dmipy_tract import FODField, seeds_from_mask
         a = self.asset
         world = np.zeros(np.shape(sh), np.float32)
@@ -1234,7 +1315,7 @@ class Brain(P.Source):
                 ["regions", f"{self.regions.n} ({len(self.regions.groups)} lobar groups)"],
                 ["the scan's timing class", self.shapes[SCAN]["label"]],
                 *[list(n) for n in self.asset.notes],
-                ["replay", "exact composition of the packs' pose responses (dmipy-sim's phantom route), float32 on the device"],
+                ["replay", "exact composition of the packs' pose responses (dmipy-replay's phantom route), float32 on the device"],
                 ["per-voxel split-half floor", "none: every voxel reads the same two walks; the packs' certified floors below are the accuracy"]]
         for (tissue, label), pk in self.packs.items():
             seg = self.segments(tissue, label); held = self.windows_held(pk, seg)
@@ -1345,7 +1426,7 @@ def _warm_worker(cfg, where, entries, queue):
     """The warm-up's process (:meth:`Brain.warm_in_background`): its own :class:`Brain` on the same configuration
     and asset, every ``(meas, physics)`` of ``entries`` computed on the CPU and put on ``queue`` as
     ``(response_key, Kernels)``, then ``None``. One line per entry on stdout (the Space's run log) with its index,
-    the seconds it took, :attr:`Kernels.stages`, the pose expansion's phases (dmipy-sim's run record, which the
+    the seconds it took, :attr:`Kernels.stages`, the pose expansion's phases (dmipy-replay's phase record, which the
     expansion joins), the process's user and system CPU seconds and page faults over the entry (a system-bound entry
     is paying the kernel, for pages or I/O, not computing) and the seconds since the start."""
     os.environ.setdefault("JAX_PLATFORMS", "cpu")
@@ -1357,11 +1438,11 @@ def _warm_worker(cfg, where, entries, queue):
     try:
         src = Brain(cfg, asset=where)
         print(f"warm-up: the source loaded, RSS {_rss_gb():.1f} GB ({time.perf_counter() - t0:.0f} s since the start)", flush=True)
-        from dmipy_sim.run import Run
+        from dmipy_replay.replay.phases import Phases
         import resource
         for i, (meas, ph) in enumerate(entries, 1):
             t = time.perf_counter(); r0 = resource.getrusage(resource.RUSAGE_SELF)
-            with Run("brain-warm", params=dict(entry=i)) as run:              # the pose expansion's phases join this record
+            with Phases("brain-warm", params=dict(entry=i)) as run:           # the pose expansion's phases join this record
                 kernels = src.prepare(meas, ph)
             r1 = resource.getrusage(resource.RUSAGE_SELF)
             phases = ", ".join(f"{n} {v:.1f}" for n, v in run.phase_seconds().items())
@@ -1441,9 +1522,9 @@ def _floor(pack):
 @functools.lru_cache(maxsize=1)
 def axis_map():
     """``(n_so3(8, 0), 45)``: a compact even-order density's harmonics to its SO(3) coefficients
-    (:func:`dmipy_sim.replay.so3.axis_density_coeffs`, which is linear), so ``resp.retained(8, 0) @ axis_map()`` is the
+    (:func:`dmipy_replay.replay.so3.axis_density_coeffs`, which is linear), so ``resp.retained(8, 0) @ axis_map()`` is the
     response contracted to the FOD's harmonics."""
-    from dmipy_sim.replay import so3
+    from dmipy_replay.replay import so3
     return np.stack([so3.axis_density_coeffs(e, LMAX, 0) for e in np.eye(N_COEF)], 1)
 
 

space/sources/disco.py

@@ -1,20 +1,22 @@
 """The DiSCo source: the DiSCo phantom's 40^3 grid replayed from its stored walk, tracked from the sixteen regions
 and scored against the dataset's strand-count and cross-sectional-area matrices (CC BY 4.0, ``data/SOURCE.md``).
 
-Two sources replay the same walk. :class:`Layout` (demo mode) is dmipy-sim's shape-moment layout: the DiSCo replay
-pack contracted once per stored pulse timing ("class"), so a protocol on those classes at any b-values, directions,
+Two sources replay the same walk. :class:`Layout` (demo mode) is dmipy-replay's time-segment layout: the DiSCo walk
+contracted once per stored waveform ("class"), so a protocol on those classes at any b-values, directions,
 tissue and field is an elementwise kernel on device-resident rows, in seconds. :class:`Columns` (full mode) is the
 columnar replay pack itself, read once per run for any pulse timing, in minutes. Both are a :class:`Disco`: the
 page's texts, the tissue panel of the white-matter catalogue, the regions, the ground truth and the GPU budget.
 """
 from __future__ import annotations
 
+import functools
 import os
 import time
 
 import numpy as np
 
 from .. import pipeline as P
+from .. import tensor as T
 
 N_REGIONS = 16
 VOXEL_M = 25e-6                                          # DiSCo's voxel; the strands' coordinate unit
@@ -119,7 +121,8 @@ class Disco(P.Source):
             truth_tab="2 · ground truth",
             truth_labels=("the strands", "the ground-truth matrices"),
             stages={"replay": "replaying the grid from the stored walk", "noise": "adding Rician noise", "csd": "fitting CSD (order 8)",
-                    "track": "tracking from the sixteen regions", "score": "scoring the connectome"},
+                    "track": "tracking from the sixteen regions", "score": "scoring the connectome",
+                    "maps": "fitting FA, MD and µFA (the b-tensor run's maps)"},
             fixed=["the DiSCo phantom's geometry and its one stored walk", "the sixteen regions and the ground-truth matrices",
                    "the pulse timing classes of the layout (demo mode)"],
             accuracy=("A replay is a measured approximation of the stored walk, not a rendering: the walk's two halves are replayed "
@@ -222,7 +225,8 @@ class Disco(P.Source):
             out[f"SNR → {v}"] = ("snr", float(v))
         out["noise → off"] = ("snr", None)
         for name, sh in cfg["shapes"].items():
-            out[f"every shell's pulse timing → {name} ({sh['label']})"] = ("shape", name)
+            if sh.get("kind") != "tensor":               # the b-tensor class is its own acquisition (encodings, not shells)
+                out[f"every shell's pulse timing → {name} ({sh['label']})"] = ("shape", name)
         return out
 
     @classmethod
@@ -262,16 +266,21 @@ class Disco(P.Source):
 
     @classmethod
     def estimated_seconds(cls, cfg, protocol, *, density, knob, n_keys, ladder, responses=(), scanner=None):
-        """Measured on the pool at DiSCo 364 with every tier (#7, dmipy-sim#522/#523): the worker's start and the
-        payload's handoff 8 s, noise to scoring 8 s, the first replay 0.045 s per measurement (the tiles uploaded
-        inside the call), each further replay on the resident tiles 0.032 s per measurement (the ladder is three, B
-        one plus its 8 s of stages, a new field direction included), 10 s per extra tracker key; on a machine every
-        replay adds its delivery, ``[scanners] delivery_per_meas`` s per measurement (the delivered vectors of the
-        layout's live voxels on the host and their upload); times 1.3, within 30 and 480 s. DiSCo prepares nothing,
-        so ``responses`` is empty."""
+        """Measured on the pool (RTX PRO 6000 MIG 2g.48gb, the v2 layout streamed, every tier at 3 T, 2026-10-30): a
+        replay costs per measurement per stored term of its class (a PGSE, stimulated-echo or OGSE class 2, the
+        b-tensor class 6: :func:`space.classes.n_terms`), 0.05 s the first in a call (DiSCo 364: 32.4 s; the b-tensor
+        preset: 23.4 s) and 0.045 s each further one (the ladder is three, B one); the worker's start and the
+        payload's handoff 8 s, the stages after the replay 10 s (noise to scoring, or the b-tensor maps), B's 10 s
+        more, 10 s per extra tracker key; on a machine every replay adds its delivery, ``[scanners]
+        delivery_per_meas`` s per measurement; times 1.3, within 30 and 480 s. DiSCo prepares nothing, so
+        ``responses`` is empty."""
+        from .. import classes as C
         n = protocol.n_meas
-        replays = 1 + (3 if ladder else 0) + (1 if knob != P.NO_KNOB else 0)
-        secs = 16 + 0.045 * n + (3 * 0.032 * n if ladder else 0) + ((8 + 0.032 * n) if knob != P.NO_KNOB else 0) + 10 * (int(n_keys) - 1)
+        terms = sum(C.n_terms(cfg["shapes"][s.shape]) * s.n_dirs if s.shape in cfg["shapes"] else 2 * s.n_dirs for s in protocol.shells)
+        terms += C.n_terms(cfg["shapes"][protocol.shells[0].shape]) * protocol.n_b0 if protocol.shells[0].shape in cfg["shapes"] else 2 * protocol.n_b0
+        b = knob != P.NO_KNOB
+        replays = 1 + (3 if ladder else 0) + (1 if b else 0)
+        secs = 8 + 10 + 0.05 * terms + 0.045 * terms * (replays - 1) + (10 if b else 0) + 10 * (int(n_keys) - 1)
         if scanner is not None:
             secs += replays * float(cfg["scanners"]["delivery_per_meas"]) * n
         return int(min(480, max(30, 1.3 * secs)))
@@ -299,7 +308,7 @@ class Disco(P.Source):
         """The ground truth tab: a sample of the strands in 3-D with the region markers, and the two matrices."""
         from dmipy_sim.io.strands import read_tck, read_diameters
         from .. import viewers as V
-        strands = read_tck(os.path.join(P.DATA_DIR, "DiSCo_Strands_Trajectories.tck"), coordinate_unit_m=VOXEL_M)
+        strands = read_tck(strands_path(self.cfg), coordinate_unit_m=VOXEL_M)
         diameters = read_diameters(os.path.join(P.DATA_DIR, "DiSCo_Strands_Diameters.txt"), diameter_unit_m=1e-3)
         return (V.strands3d([s / VOXEL_M for s in strands], diameters, V.region_markers(self.regions), self.mask.shape),
                 V.ground_truth_matrix(self.gt_count, self.gt_area, P.connected_pairs(self.gt_count)))
@@ -343,6 +352,52 @@ class Disco(P.Source):
         return [pool, contact, fld]
 
 
+@functools.lru_cache(maxsize=2)
+def strands_path(cfg=None):
+    """The DiSCo strand centrelines (``.tck``, voxel units; CC BY 4.0): ``[data] strands`` in the replay dataset
+    (``[data] pack``'s repository) at ``[data] pack_revision``, fetched into the Hub cache once."""
+    from huggingface_hub import hf_hub_download
+    d = (cfg or P.config())["data"]
+    repo = "/".join(d["pack"][len("hf://"):].split("/")[:2])
+    return hf_hub_download(repo, d["strands"], repo_type="dataset", revision=d["pack_revision"])
+
+
+def strand_configuration(tck=None, shape=(40, 40, 40), *, angle_deg=30.0, share=0.10):
+    """The ground truth's fibre configuration per voxel from the released strand trajectories (voxel units):
+    ``(n_dirs, dispersion, length)`` -- the distinct strand directions crossing the voxel (the strands' length-weighted
+    mean tangents in it, clustered at ``angle_deg``, a direction counted when it carries at least ``share`` of the voxel's
+    strand length), ``1 - lambda_1 / trace`` of the length-weighted dyadic tensor of the tangents, and the strand length
+    in the voxel (voxel units)."""
+    from dmipy_sim.io.strands import read_tck
+    strands = read_tck(tck or strands_path(), coordinate_unit_m=1.0)
+    n_vox = int(np.prod(shape))
+    Tm = np.zeros((n_vox, 3, 3)); per = {}
+    for c in strands:
+        c = np.asarray(c, float)
+        if len(c) < 2:
+            continue
+        seg = np.diff(c, axis=0); mid = 0.5 * (c[1:] + c[:-1]); L = np.linalg.norm(seg, axis=1); ok = L > 0
+        t = seg[ok] / L[ok, None]; ijk = np.floor(mid[ok]).astype(int); inside = np.all((ijk >= 0) & (ijk < np.asarray(shape)), axis=1)
+        v = np.ravel_multi_index(ijk[inside].T, shape); w = L[ok][inside]; t = t[inside]
+        np.add.at(Tm, v, (t[:, :, None] * t[:, None, :]) * w[:, None, None])
+        for vox in np.unique(v):
+            k = v == vox
+            per.setdefault(int(vox), []).append((np.linalg.eigh(np.einsum("n,ni,nj->ij", w[k], t[k], t[k]))[1][:, -1], float(w[k].sum())))
+    ev = np.linalg.eigvalsh(Tm); tr = ev.sum(-1)
+    disp = np.where(tr > 0, 1.0 - ev[:, -1] / np.maximum(tr, 1e-30), np.nan)
+    ndirs = np.zeros(n_vox, int); cos = np.cos(np.radians(angle_deg))
+    for vox, items in per.items():
+        items = sorted(items, key=lambda x: -x[1]); total = sum(w for _, w in items); clusters = []
+        for d, w in items:
+            for cl in clusters:
+                if abs(d @ cl[0]) >= cos:
+                    cl[1] += w; break
+            else:
+                clusters.append([d, w])
+        ndirs[vox] = sum(1 for _, w in clusters if w >= share * total)
+    return ndirs.reshape(shape), disp.reshape(shape), tr.reshape(shape)
+
+
 PAIRS = np.triu_indices(N_REGIONS, 1)                # the 120 region pairs
 
 
@@ -363,48 +418,65 @@ def compare(a, b):
 
 
 class Layout(Disco):
-    """Demo mode: dmipy-sim's shape-moment layout (``dmipy_sim.replay.shape_moments``) at the config's dataset
-    revision, or a local copy; every stored class's timing is the config's, and a tissue needs the layout's tiers."""
+    """Demo mode: dmipy-replay's time-segment layout (``dmipy_replay.replay.time_segments``) at the config's dataset revision,
+    or a local copy. Every stored class is one waveform (:mod:`space.classes`), checked at load to be the layout's
+    partition of that name; a shell plays it at its b-values and directions."""
     mode = "demo"
 
     @classmethod
     def load(cls, cfg, *, local=None):
-        """The layout at ``local``, ``DISCO_MOMENTS``, or the config's dataset revision on the Hub."""
-        return cls(cfg, local=local or os.environ.get("DISCO_MOMENTS"))
+        """The layout at ``local``, ``DISCO_LAYOUT``, or the config's dataset revision on the Hub."""
+        return cls(cfg, local=local or os.environ.get("DISCO_LAYOUT"))
 
     def __init__(self, cfg, *, local=None):
-        from dmipy_sim.replay.shape_moments import ShapeMoments
+        import hashlib
+        from dmipy_replay.replay.time_segments import MANIFEST, TimeSegmentLayout
+        from .. import classes as C
         super().__init__(cfg)
         d = cfg["data"]
-        uri = local or f"hf://{d['repo']}/{d['moments']}"
-        self.moments = ShapeMoments.open(uri, revision=d.get("revision") or None) if uri.startswith("hf://") else ShapeMoments(uri)
+        uri = local or f"hf://{d['repo']}/{d['layout']}"
+        names = list(self.shapes)
+        self.layout = TimeSegmentLayout.open(uri, revision=d.get("revision") or "main", partitions=names)
         self.resident = P.resident(cfg)
-        want = d.get("source_manifest_sha256")
-        got = self.moments.manifest["source"]["manifest_sha256"]
+        want = d.get("manifest_sha256")
+        with open(os.path.join(self.layout.path, MANIFEST), "rb") as f:
+            got = hashlib.sha256(f.read()).hexdigest()
         if want and got != want:
-            raise ValueError(f"the layout at {uri} was contracted from columnar manifest {got[:12]}, the config pins {want[:12]}")
-        missing = [s for s in self.shapes if s not in self.moments.shapes]
+            raise ValueError(f"the layout at {uri} has manifest {got[:12]}, the config pins {want[:12]}")
+        self.manifest_sha256 = got
+        missing = [s for s in names if s not in self.layout.partitions]
         if missing:
-            raise ValueError(f"the layout at {uri} lacks the timing classes {missing}; it holds {self.moments.shapes}")
-        for name, t in self.shapes.items():
-            enc = self.moments.manifest["shapes"][name].get("encoding") or {}
-            for k in ("delta", "Delta", "TE"):
-                v = enc.get(k)
-                v = v[0] if isinstance(v, list) and v else v
-                if v is not None and abs(float(v) - float(t[k])) > 1e-9:
-                    raise ValueError(f"class {name}: the layout's {k} is {v} s, the config's {t[k]} s")
-        self.tiers = bool(self.moments.manifest.get("tiers"))
-        self.M0 = None                                   # the bare b = 0 map (every voxel's walker weight): set by warm()
-        if self.tiers:                                   # the pools with walkers are the spec's first n_pools ids
-            t = self.moments.manifest["tiers"]
+            raise ValueError(f"the layout at {uri} lacks the timing classes {missing}; it holds {list(self.layout.partitions)}")
+        for name, sh in self.shapes.items():                 # each class's waveform is the stored partition's
+            got_name, _ = self.layout.fit(C.template(sh), partition=name)
+            if not self.layout.partitions[name].is_constant and C.kind(sh) not in ("ogse", "tensor"):
+                raise ValueError(f"class {name}: the layout stores shaped terms, the config a square-pulse {C.kind(sh)}")
+        m = self.layout.manifest
+        self.tiers = bool(m.get("tiers"))
+        self.M0 = self.layout.weights.sum(1).reshape(self.layout.grid.shape)   # every voxel's walker weight: the bare b = 0 signal
+        if self.tiers:                                   # the pools with walkers: the pool ids of the rows that weigh something
+            t = m["tiers"]
             spec_pools = sorted(t["substrate"]["pools"], key=lambda q: q["id"]) if t.get("substrate") else []
-            seeded = tuple(q["name"] for q in spec_pools[:int(t["pools"])])
+            ids = self.seeded_pool_ids()
+            seeded = tuple(q["name"] for q in spec_pools if q["id"] in ids)
             if set(seeded) != set(self.pools):
                 raise ValueError(f"the layout's seeded pools are {seeded}, the config's [physics] pools {self.pools}")
-            self.unseeded = tuple(q["name"] for q in spec_pools[int(t["pools"]):])
-        self.check_grid(self.moments.grid.shape, "the layout")
+            self.unseeded = tuple(q["name"] for q in spec_pools if q["id"] not in ids)
+        self.check_grid(self.layout.grid.shape, "the layout")
+
+    def seeded_pool_ids(self, chunk=1 << 16):
+        """The pool ids of the layout's rows with weight (a padding row weighs nothing), read once in chunks of tiles."""
+        w = self.layout._column("w"); pool = self.layout._column("pool"); ids = set()
+        for s in range(0, self.layout.n_tiles, chunk):
+            q = np.asarray(pool[s:s + chunk]); ids |= set(np.unique(q[np.asarray(w[s:s + chunk]) > 0]).tolist())
+        return ids
+
+    def background(self, name):
+        """Whether the class's sequence group carries the background vectors a magnet's own gradient needs."""
+        return bool(self.layout.manifest["groups"][self.layout.manifest["partitions"][name]["group"]].get("background"))
 
     def validate(self, protocol, physics):
+        from .. import classes as C
         meas = P.measurements(protocol, self.shapes)
         if physics and not physics.bare and not self.tiers:
             raise ValueError("this layout holds bare diffusion only: switch the tissue panel off")
@@ -412,21 +484,35 @@ class Layout(Disco):
             raise ValueError(f"the tissue names the pools {physics.pools}; this layout's are {self.pools}")
         if physics and physics.scanner is not None:
             L = P.limits(physics.scanner)
-            if L.has_field_law and not self.moments.background:
-                raise ValueError(f"the {physics.scanner}'s own gradient needs the layout's background moments, and this layout has "
-                                 "none (dmipy_sim.replay.shape_moments.stamp_background adds them in one pass)")
-            if L.has_transmit_profile and any("rf" not in self.moments.manifest["shapes"][n] for n in np.unique(meas.shape)):
-                raise ValueError(f"the {physics.scanner}'s transmit scale needs the classes' RF schedules, which this layout does not "
-                                 "record (stamp_background records them)")
+            shaped = [n for n in np.unique(meas.shape) if not self.layout.partitions[n].is_constant]
+            if shaped:
+                raise ValueError(f"a machine's delivered gradient is one vector per constant segment, and the classes {shaped} are "
+                                 "shaped waveforms (an OGSE): choose the ideal scanner for them")
+            if L.has_field_law and not all(self.background(n) for n in np.unique(meas.shape)):
+                raise ValueError(f"the {physics.scanner}'s own gradient needs the layout's background vectors, and this layout "
+                                 "lacks them for a class of this run")
+        for name in np.unique(meas.shape):
+            C.kind(self.shapes[name])
         return meas
 
-    def delivered(self, name, bvals, dirs, physics):
-        """What the run's machine delivers at every voxel with rows for the class ``name`` at ``bvals`` (s/m^2) along
-        ``dirs``: the layout's :meth:`~dmipy_sim.replay.shape_moments.ShapeMoments.delivery` with the grid's centre at
-        the physics' offset; None on the ideal scanner."""
+    def sequence(self, name, bvals, dirs, encodings=None):
+        """The rows on class ``name`` (b in s/mm^2) as the class's waveform at their amplitudes and directions; on a
+        b-tensor class each row's ``encodings`` entry (``b0``, ``LTE``, ``STE``) about its direction."""
+        from .. import classes as C
+        sh = self.shapes[name]
+        if C.kind(sh) == "tensor":
+            if encodings is None:
+                raise ValueError(f"class {name} is a b-tensor class: its rows name their encodings")
+            return T._sequence(C.tensor_rows(sh, encodings, bvals, dirs))
+        return C.play(sh, np.asarray(bvals, np.float64) * 1e6, dirs)
+
+    def delivered(self, name, seq, physics):
+        """What the run's machine delivers at every voxel with rows for the class ``name`` playing ``seq``: the layout's
+        :meth:`~dmipy_replay.replay.time_segments.TimeSegmentLayout.delivery` with the grid's centre at the physics'
+        offset; None on the ideal scanner."""
         if physics is None or physics.scanner is None:
             return None
-        return self.moments.delivery(name, bvals, dirs, P.limits(physics.scanner), self.moments.grid.centred_at(physics.offset_m))
+        return self.layout.delivery(seq, P.limits(physics.scanner), self.layout.grid.centred_at(physics.offset_m), partition=name)
 
     def replay(self, meas, physics=None, prepared=None):
         t0 = time.perf_counter()
@@ -436,36 +522,34 @@ class Layout(Disco):
         b0 = physics.b0_direction if physics else P.B0_ALONG_Z
         for name in np.unique(meas.shape):
             rows = meas.shape == name
-            S[..., rows], f = self.moments.image(name, meas.bvals[rows] * 1e6, meas.dirs[rows], backend=self.backend,
-                                                 resident=True, tissue=tissue, scanner=scanner, b0_direction=b0,
-                                                 delivered=self.delivered(name, meas.bvals[rows] * 1e6, meas.dirs[rows], physics))
-            floor = np.fmax(floor, f)                    # resident within a run: release() drops the device copies after it
-        if self.M0 is None:                              # the bare b = 0 map: the walker weight of every voxel
-            self.M0 = self.moments.image(meas.shape[0], np.zeros(1), np.array([P.B0_ALONG_Z]), backend=self.backend, resident=True)[0][..., 0]
+            seq = self.sequence(name, meas.bvals[rows], meas.dirs[rows], None if meas.encoding is None else meas.encoding[rows])
+            img, f = self.layout.image(seq, partition=str(name), backend=self.backend, resident=self.resident, tissue=tissue,
+                                       scanner=scanner, b0_direction=b0, delivered=self.delivered(str(name), seq, physics))
+            S[..., rows] = img * self.M0[..., None]       # the voxel's signal: its walkers' weight times the normalised one
+            floor = np.fmax(floor, f)
         dwi, factor = P.s0_normalised(S, meas.b0, self.M0)
         return dwi, floor, factor, time.perf_counter() - t0
 
     def release(self):
         if not self.resident:
-            self.moments.release()
+            self.layout.release()
 
     def ingredients(self, meas, physics=None, prepared=None):
         """What each tier multiplies into the sum, per voxel, for the run's first class (the layout's
-        :meth:`~dmipy_sim.replay.shape_moments.ShapeMoments.tier_maps` on the device columns the replay left
-        resident): ``intra_fraction`` (the walker weight in the intra-axonal pool over the voxel's),
-        ``wall_contact_um`` (the walkers' boundary local time l under the class's gate, a length; the layout stores
-        it signed as the exponent's term, -l, so the contact tier's weight is exp(rho2 c / D) = exp(-rho2 l / D)),
-        ``contact_survival`` (that factor at the run's rho2, None without the contact tier), ``field_rad`` (the spread
-        over the voxel's walkers of the dephasing phase the sheath's field gives them by the echo at the run's field
-        and direction, the exact per-walker phase the kernel applies, in radians; None without the field tier), and
-        ``D_walk`` (m^2/s)."""
+        :meth:`~dmipy_replay.replay.time_segments.TimeSegmentLayout.tier_maps`): ``intra_fraction`` (the walker weight in
+        the intra-axonal pool over the voxel's), ``wall_contact_um`` (the walkers' boundary local time l under the
+        class's gate, a length; the layout stores it signed as the exponent's term, -l, so the contact tier's weight is
+        exp(rho2 c / D) = exp(-rho2 l / D)), ``contact_survival`` (that factor at the run's rho2, None without the
+        contact tier), ``field_rad`` (the spread over the voxel's walkers of the dephasing phase the sheath's field
+        gives them by the echo at the run's field and direction, the exact per-walker phase the kernel applies, in
+        radians; None without the field tier), and ``D_walk`` (m^2/s)."""
         if not self.tiers:
             return None
         tissue, scanner = P.tissue_and_scanner(physics, self.unseeded)
         b0 = physics.b0_direction if physics else P.B0_ALONG_Z
-        maps = self.moments.tier_maps(str(meas.shape[0]), tissue, scanner, backend=self.backend, resident=True, b0_direction=b0)
+        maps = self.layout.tier_maps(str(meas.shape[0]), tissue, scanner, backend=self.backend, resident=self.resident, b0_direction=b0)
         contact = maps["contact"]
-        return dict(D_walk=float(self.moments.manifest["tiers"]["D_walk"]), intra_fraction=maps["pool"]["intra"],
+        return dict(D_walk=float(self.layout.manifest["tiers"]["D_walk"]), intra_fraction=maps["pool"]["intra"],
                     wall_contact_um=None if contact is None else -contact * 1e6,
                     contact_survival=maps["contact_weight"] if (physics and physics.contact) else None,
                     field_rad=maps["phase_std"] if (physics and physics.field) else None)
@@ -477,33 +561,41 @@ class Layout(Disco):
 
     def warm(self):
         """DiSCo's own protocol replayed once at the default physics (every class it plays compiled at the row
-        counts the default run uses), then one direction on every other class (its tier group's terms contracted
-        and on the device). When the tiles are not kept resident (a pool that drops the device between calls), the
-        padded host arrays are loaded into this process instead, so a forked worker inherits them and pays the
-        transfer alone."""
-        if not self.resident:
-            self.moments.preload(list(self.shapes))
+        counts the default run uses), then one direction on every other class (on a b-tensor class one row of each
+        encoding). When the columns are not kept resident
+        (a pool that drops the device between calls), the padded host arrays are loaded into this process instead, so
+        a forked worker inherits them and pays the transfer alone."""
+        if not self.resident and self.backend == "torch":
+            self.layout.preload(list(self.shapes))
             return
         physics = self.default_physics()
         disco = disco_protocol(self.cfg)[0]
         self.replay(P.measurements(disco, self.shapes), physics)
-        rest = [name for name in self.shapes if name not in {s.shape for s in disco.shells}]
+        played = {s.shape for s in disco.shells}
+        rest = [n for n, sh in self.shapes.items() if n not in played and sh.get("kind") != "tensor"]
         if rest:
             self.replay(P.measurements(P.Protocol(tuple(P.Shell(name, 1000.0, 1) for name in rest), n_b0=1, name="warm"), self.shapes), physics)
+        for name in [n for n, sh in self.shapes.items() if sh.get("kind") == "tensor"]:     # a b-tensor class: one row per encoding
+            shells = tuple(P.Shell(name, 1000.0, 1, encoding=e) for e in ("LTE", "STE"))
+            self.replay(P.measurements(P.Protocol(shells, n_b0=1, name="warm"), self.shapes), physics)
 
     def accuracy(self, res=None):
-        m = self.moments.manifest; src = m["source"]
-        rows = [["temporal bands kept (K)", str(m["K"])],
-                ["band error at the built amplitude (worst class)", f"{m['band_error']:.2e}"],
-                ["band tolerance (× the pack's floor)", f"{m['tol']:g}"],
-                ["source pack", str(src.get("pack"))],
-                ["source pack's certified median floor", f"{src['floor']:.4g}"],
-                ["source manifest sha256", src["manifest_sha256"][:16]],
-                ["layout written by", f"{m.get('code', {}).get('commit', '?')[:12]} on {m.get('created', '?')}"],
-                ["pulses", "square (slew rate ∞), one TE per class"],
+        from .. import classes as C
+        m = self.layout.manifest; src = m["source"]; cert = m["certificate"]; env = src["envelope"]
+        rows = [["source", f"the {src['kind']}'s coefficients ({len(src.get('blocks') or [])} blocks); classes appended from the "
+                           "quantised pack carry its quantisation in their own term" if any(p.get("source") for p in m["partitions"].values())
+                           else f"the {src['kind']}'s coefficients ({len(src.get('blocks') or [])} blocks)"],
+                ["floor target (the words' budget)", f"{float(src['floor_target']):.4g}"],
+                ["layout's worst signal error over its classes and tiers (err_max)", f"{float(cert['err_max']):.2e}"],
+                ["envelope: gradient per axis / B0 / ρ₂/D", f"{float(env['G_max']) * 1e3:.0f} mT/m ({env.get('scanner') or env.get('bound')}) / "
+                                                            f"{float(env['B0_max_T']):g} T / {float(env['rho2_over_D_max']):g} 1/m"],
+                ["layout manifest sha256", self.manifest_sha256[:16]],
+                ["layout written by", f"{(m.get('code') or {}).get('commit') or '?'} on {m.get('created', '?')}"],
+                ["pulses", "one waveform per class: square lobes (slew rate ∞) for the PGSE and stimulated-echo classes, "
+                           "the cosine OGSE with its 200 T/m/s edges as built at b = 1000"],
                 ["tiers stored", "pool, contact, field" if self.tiers else "none (bare diffusion only)"],
-                ["background moments (a magnet's own gradient)", "stored, certified to {:g} mT/m".format(self.moments.background["amplitude"] * 1e3)
-                 if self.moments.background else "not stored: a machine with its own gradient is refused"]]
+                ["background vectors (a magnet's own gradient)", "stored, certified to {:g} mT/m".format(float(src["background_amplitude"]) * 1e3)
+                 if any(g.get("background") for g in m["groups"].values()) else "not stored: a machine with its own gradient is refused"]]
         if res is not None and res.physics is not None and res.physics.scanner is not None:
             rows.append(["scanner", f"{res.physics.scanner}, the phantom's centre at "
                                     f"({', '.join(f'{x * 100:.1f}' for x in res.physics.offset_m)}) cm from isocentre: every catalogued term "
@@ -511,22 +603,22 @@ class Layout(Disco):
         rows += super().accuracy(res)
         if res is not None:
             for name in np.unique(res.meas.shape):
-                sh = m["shapes"][name]; enc = sh.get("encoding") or {}
-                first = {k: (v[0] if isinstance(v, list) and v else v) for k, v in enc.items()}
-                rows.append([f"class {name}", f"δ {first.get('delta')} / Δ {first.get('Delta')} / TE {first.get('TE')} s, "
-                                              f"built at {float(sh.get('amplitude_built') or 0) * 1e3:.0f} mT/m, pathway amplitude {sh.get('pathway') or 1}"])
+                p = m["partitions"][str(name)]
+                rows.append([f"class {name}", f"{C.describe(self.shapes[str(name)])}, TE {float(self.shapes[str(name)].get('TE', 0.0535)) * 1e3:g} ms; "
+                                              f"{4 * int(p['n_words'])} B per walker, rounding variance {float(p['rounding']['variance']):.2e} rad²"])
         return rows
 
 
 class Columns(Disco):
-    """Full mode: the columnar replay pack (:class:`dmipy_sim.replay.columnar.ColumnarPack`) at
-    :func:`space.pipeline.columns_uri`, read once per run: one ScannerSequence per run (one echo time, one pulse kind, square
-    pulses), the tissue and field of the run, the field along the pack's z."""
+    """Full mode: the columnar replay pack (:class:`dmipy_replay.replay.columnar.ColumnarPack`) at
+    :func:`space.pipeline.columns_uri`, read once per run: one ScannerSequence per run (one echo time, one RF schedule;
+    the stored classes' waveforms, square pulses for a free timing), the tissue and field of the run, the field along
+    the pack's z."""
     mode = "full"
     RATE = {True: 45e6, False: 500e6}                # bytes/s read from the Hub, from a local disk
 
     def __init__(self, cfg, *, uri=None):
-        from dmipy_sim.replay.columnar import ColumnarPack
+        from dmipy_replay.replay.columnar import ColumnarPack
         super().__init__(cfg)
         self.uri = uri or P.columns_uri(cfg)
         self.pack = ColumnarPack(self.uri, workers=int(cfg["mode"]["workers"]))
@@ -534,15 +626,17 @@ class Columns(Disco):
         self.remote = self.uri.startswith("hf://")
 
     def kinds(self, meas):
-        """The pulse kinds (``pgse`` / ``pgste``) of the classes the rows play; a free timing is PGSE."""
-        return sorted({self.shapes[n].get("kind", "pgse") if n in self.shapes else "pgse" for n in np.unique(meas.shape)})
+        """The pulse kinds (``pgse`` / ``pgste`` / ``ogse``) of the classes the rows play; a free timing is PGSE."""
+        from .. import classes as C
+        return sorted({C.kind(self.shapes[n]) if n in self.shapes else "pgse" for n in np.unique(meas.shape)})
 
     def validate(self, protocol, physics):
         meas = P.measurements(protocol, self.shapes)
         if len(np.unique(np.round(meas.TE, 9))) != 1:
             raise ValueError("full mode replays one echo time per run")
-        if len(self.kinds(meas)) != 1:
-            raise ValueError("full mode replays one pulse kind (PGSE or stimulated echo) per run")
+        kinds = self.kinds(meas)
+        if len(kinds) != 1 and not set(kinds) <= {"pgse", "ogse"}:
+            raise ValueError("full mode replays one RF schedule per run (spin echoes, or a stimulated echo)")
         if physics and not physics.along_z:
             raise ValueError("full mode replays the field along z only")
         if physics and physics.scanner is not None:
@@ -551,15 +645,27 @@ class Columns(Disco):
         return meas
 
     def sequence(self, meas):
-        """The one ScannerSequence of the measurements: square pulses, every row's own delta / Delta at the rows'
-        common TE; a stimulated echo stores for TM = Delta - delta between 90° pulses."""
-        import dmipy_sim as d
+        """The one ScannerSequence of the measurements, rows in order: a stored class's rows play its waveform
+        (:mod:`space.classes`, the layout's own), a free shell square PGSE lobes at its delta / Delta, all at the rows'
+        common TE on the layout's grid."""
+        from dmipy_replay.acquisition.scanner_sequence import ScannerSequence
+        from dmipy_replay.sequences import pgse
+        from .. import classes as C
         TE = float(np.unique(np.round(meas.TE, 9))[0])
-        kind, = self.kinds(meas)
-        if kind == "pgste":
-            return d.pgste(meas.dirs.tolist(), meas.delta, meas.Delta - meas.delta, bvalues=meas.bvals * 1e6, n_t=1000,
-                           slew_rate=np.inf, ste_flip_angles=(90.0, 90.0, 90.0))
-        return d.pgse(meas.dirs.tolist(), meas.delta, meas.Delta, bvalues=meas.bvals * 1e6, TE=TE, n_t=1000, slew_rate=np.inf)
+        G = np.zeros((len(meas.bvals), C.N_T, 3), np.float32); tmpl = None
+        for name in np.unique(meas.shape):
+            rows = meas.shape == name
+            if name in self.shapes and C.kind(self.shapes[name]) == "tensor":
+                seq = T._sequence(C.tensor_rows(self.shapes[name], meas.encoding[rows], meas.bvals[rows], meas.dirs[rows]))
+            elif name in self.shapes:
+                seq = C.play(self.shapes[name], meas.bvals[rows] * 1e6, meas.dirs[rows])
+            else:
+                seq = pgse(meas.dirs[rows].tolist(), meas.delta[rows], meas.Delta[rows], bvalues=meas.bvals[rows] * 1e6, TE=TE,
+                           n_t=C.N_T, slew_rate=np.inf)
+            if tmpl is not None and (seq.dt != tmpl.dt or seq.rf.to_dicts() != tmpl.rf.to_dicts()):
+                raise ValueError("full mode replays one RF schedule and grid per run")
+            tmpl = seq; G[rows] = np.asarray(seq.G, np.float32)
+        return ScannerSequence(G=G, dt=tmpl.dt, rf=tmpl.rf, readout=tmpl.readout)
 
     def plan(self, meas, physics=None):
         """What the replay will read (the pack's plan: bands, modes, tiers, rows, bytes) and ``estimated_seconds``
@@ -570,11 +676,11 @@ class Columns(Disco):
         return plan
 
     def replay(self, meas, physics=None, prepared=None, *, progress=None):
-        from dmipy_sim.replay.study import Study, Protocol as SProtocol, Acquisition
+        from dmipy_replay.replay.study import Acquisition, Study, StudyProtocol
         t0 = time.perf_counter()
         tissue, scanner = P.tissue_and_scanner(physics, self.unseeded)
-        study = Study(SProtocol([Acquisition(self.sequence(meas), name="run")]), tissues=[tissue], scanners=[scanner])
-        S, floor, plan = self.pack.image_study(study, tol=P.BAND_TOL, chunk_rows=1_000_000, progress=progress)
+        study = Study(StudyProtocol([Acquisition(self.sequence(meas), name="run")]), tissues=[tissue], scanners=[scanner])
+        S, floor, plan = self.pack.image_study(study, tol=P.BAND_TOL, chunk_rows=1_000_000, backend="jax", progress=progress)
         self.last_plan = plan
         dwi, factor = P.s0_normalised(S[0], meas.b0)      # the pack's bare weights are not read here: the factor is 1
         return dwi, floor[0], factor, time.perf_counter() - t0
@@ -583,7 +689,7 @@ class Columns(Disco):
         m = self.pack.meta
         rows = [["source pack", str(m.get("id"))], ["source pack's certified median floor", f"{self.pack.floor:.4g}"],
                 ["bands kept", str(self.pack.K)], ["band tolerance (× the pack's floor)", f"{P.BAND_TOL:g}"],
-                ["pulses", "square (slew rate ∞), one TE and one pulse kind per run"]]
+                ["pulses", "the stored classes' waveforms, square (slew rate ∞) for a free timing; one TE and one RF schedule per run"]]
         plan = getattr(self, "last_plan", None)
         if plan:
             rows.append(["last replay read", f"{plan.get('rows', 0):,} rows, {plan.get('bytes', 0) / 1e9:.1f} GB, bands {plan.get('K')}, modes {plan.get('modes')}"])

space/tensor.py

@@ -0,0 +1,245 @@
+"""The b-tensor encoding class: three gradient shapes psi_1..3 on the spin echo's two encoding windows, stored as one
+class of the layout (six terms: each shape before and after the 180), and the rows that play any b-tensor of them.
+
+A row plays ``G(t) = sum_k psi_k(t) a_k``, the 3x3 level matrix ``A`` (row k the vector shape k is played along), so
+its b-tensor is ``B = A^T M A`` with ``M`` the shapes' Gram matrix (``M_kl = int k_k k_l dt``, the moments to the
+readout). Every encoding is a choice of ``A``: spherical ``A = sqrt(b/3) M^-1/2 R^T`` (``B = b/3 I`` at any rotation
+``R``), planar ``A = sqrt(b/2) M^-1/2 [u; v; 0]`` (``B = b/2 (u u^T + v v^T)``), linear ``A = sqrt(b) c d^T`` with ``c``
+the shape mix of the most b per peak gradient (``c^T M c = 1``; ``B = b d d^T``). The layout reads every row from
+the same stored vectors, exactly as a pulse pair's rows are read: the encoding is exact in the substrate (the mean of
+the walkers' phase factors), not a Gaussian approximation.
+
+The waveform is ``[tensor] waveform`` in ``data/``: the physical gradient ``(1000, 3)`` (T/m) on the layout's
+spin-echo grid (TE 53.5 ms), axis k the shape psi_k; the page checks it against the stored partition at load.
+
+The maps: MD the first cumulant of the spherical powder average, and the microscopic fractional anisotropy from the
+linear against the spherical powder average (Lasic et al. 2014, Front Phys 2:11; Szczepankiewicz et al. 2016,
+NeuroImage 125:653): the gamma model fitted to the two powder signals jointly with one intercept (the b = 0 rows'
+noise shared by every ``ln S/S0``), ``V_LTE - V_STE = (4/5) (D_iso D_Delta)^2``, and muFA is the FA of the axisymmetric
+microscopic tensor with that ``D_iso`` and ``D_Delta``."""
+from __future__ import annotations
+
+import functools
+import os
+
+import numpy as np
+
+from . import pipeline as P
+
+ENCODINGS = ("LTE", "PTE", "STE")
+B_DELTA = {"LTE": 1.0, "PTE": -0.5, "STE": 0.0}
+N_T = 1000
+TE = 0.0535
+
+
+@functools.lru_cache(maxsize=None)
+def _waveform(path):
+    G = np.load(path)
+    if G.shape != (N_T, 3) or not np.isfinite(G).all():
+        raise ValueError(f"the b-tensor waveform {path} is not a finite (1000, 3) gradient")
+    return np.asarray(G, np.float64)
+
+
+def waveform(tc):
+    """``(1000, 3)``: the class's physical gradient (T/m), axis k the shape psi_k at its designed amplitude."""
+    return _waveform(os.path.join(P.DATA_DIR, tc["waveform"]))
+
+
+@functools.lru_cache(maxsize=None)
+def _grid():
+    from dmipy_replay.sequences import pgse
+    return pgse([[0.0, 0.0, 1.0]], 0.0102, 0.0167, gradient_strengths=0.1, TE=TE, n_t=N_T, slew_rate=np.inf)
+
+
+def _sequence(G):
+    """``G (n, 1000, 3)`` on the spin-echo grid, RF schedule and readout of the layout's group g0."""
+    from dmipy_replay.acquisition.scanner_sequence import ScannerSequence
+    g = _grid()
+    return ScannerSequence(G=np.asarray(G, np.float32), dt=g.dt, rf=g.rf, readout=g.readout)
+
+
+def template(tc):
+    """The class as the layout stores it: three rows, row k shape psi_k along axis k."""
+    P_ = waveform(tc)
+    G = np.zeros((3, N_T, 3))
+    for k in range(3):
+        G[k, :, k] = P_[:, k]
+    return _sequence(G)
+
+
+def gram(tc):
+    """``M (3, 3)`` (s/m^2 per (T/m)^2 of level): the b-tensor of ``sum_k psi_k a_k`` is ``A^T M A``."""
+    from dmipy_replay.acquisition.waveforms import btensor_from_gradient
+    g = _grid()
+    G = waveform(tc)[None] * np.asarray(g.effective_gate, np.float64)[None, :, None]
+    return btensor_from_gradient(G, g.dt)[0]
+
+
+def _inv_sqrt(M):
+    w, V = np.linalg.eigh(M)
+    return V @ np.diag(w ** -0.5) @ V.T
+
+
+def lte_mix(tc, n=4000):
+    """``c (3,)`` with ``c^T M c = 1``: the shape mix whose single waveform ``sum_k c_k psi_k`` gives the most b per
+    squared peak gradient (searched over ``n`` Fibonacci directions of the sphere of mixes)."""
+    M = gram(tc); W = waveform(tc)
+    c = fibonacci(n, hemisphere=True)
+    peak = np.abs(W @ c.T).max(axis=0)
+    gain = np.einsum("nk,kl,nl->n", c, M, c) / peak ** 2
+    best = c[int(np.argmax(gain))]
+    return best / np.sqrt(best @ M @ best)
+
+
+def fibonacci(n, *, hemisphere=True):
+    """``n`` unit vectors spread on the (upper half-) sphere, deterministic."""
+    i = np.arange(n) + 0.5
+    z = 1.0 - i / n if hemisphere else 1.0 - 2.0 * i / n
+    phi = np.pi * (3.0 - np.sqrt(5.0)) * i
+    r = np.sqrt(np.clip(1.0 - z ** 2, 0.0, None))
+    return np.stack([r * np.cos(phi), r * np.sin(phi), z], axis=1)
+
+
+def frame(d):
+    """``R (3, 3)`` orthonormal with first column ``d`` (the encoding's axis), the rest completing it."""
+    d = np.asarray(d, np.float64); d = d / np.linalg.norm(d)
+    ref = np.array([1.0, 0.0, 0.0]) if abs(d[0]) < 0.9 else np.array([0.0, 1.0, 0.0])
+    u = ref - (ref @ d) * d; u /= np.linalg.norm(u)
+    return np.stack([d, u, np.cross(d, u)], axis=1)
+
+
+def levels(tc, encoding, b_SI, d):
+    """``A (3, 3)``: the level matrix of one row of ``encoding`` at ``b_SI`` (s/m^2) about the axis ``d`` (the
+    direction of a linear encoding, the normal of a planar one, the first axis of a spherical one's rotation)."""
+    M = gram(tc); Mi = _inv_sqrt(M); R = frame(d)
+    if encoding == "LTE":
+        return np.sqrt(b_SI) * np.outer(lte_mix(tc), R[:, 0])
+    if encoding == "PTE":
+        return np.sqrt(b_SI / 2.0) * Mi @ np.stack([R[:, 1], R[:, 2], np.zeros(3)])
+    if encoding == "STE":
+        return np.sqrt(b_SI / 3.0) * Mi @ R.T
+    raise ValueError(f"an encoding is one of {ENCODINGS}, got {encoding!r}")
+
+
+def rows(tc, bvals, n_dirs):
+    """The acquisition: one b = 0 row, then for each b (s/mm^2) ``n_dirs`` rows of each encoding (their axes spread on
+    the hemisphere). Returns ``(sequence, encoding (n,), b (n,) s/mm^2)``."""
+    W = waveform(tc)
+    dirs = fibonacci(int(n_dirs))
+    G = [np.zeros((N_T, 3))]; enc = ["b0"]; bs = [0.0]
+    for b in bvals:
+        for e in ENCODINGS:
+            for d in dirs:
+                G.append(W @ levels(tc, e, float(b) * 1e6, d)); enc.append(e); bs.append(float(b))
+    return _sequence(np.stack(G)), np.array(enc), np.array(bs)
+
+
+def check(tc, seq, enc, bs):
+    """The rows' b-tensors as built: ``{encoding: (max |b - b_target| / b, max |b_delta - target|)}`` from
+    dmipy-replay's own b-tensor of the sequence, and the peak gradient per axis (T/m) of every row."""
+    from dmipy_replay.acquisition.waveforms import btensor_invariants
+    B = seq.btensor()
+    b, b_delta, _ = btensor_invariants(B)
+    out = {}
+    for e in ENCODINGS:
+        m = enc == e
+        if m.any():
+            out[e] = (float(np.max(np.abs(b[m] / 1e6 - bs[m]) / bs[m])), float(np.max(np.abs(b_delta[m] - B_DELTA[e]))))
+    return out, np.abs(np.asarray(seq.G, np.float64)).max(axis=(1, 2))
+
+
+def powder(S, enc, bs, bvals):
+    """``{(encoding, b): S/S0 averaged over the encoding's rows at b}`` per voxel, ``S`` ``(..., n)`` normalised to its
+    b = 0 row."""
+    return {(e, float(b)): np.nanmean(S[..., (enc == e) & (bs == float(b))], axis=-1) for e in ENCODINGS for b in bvals}
+
+
+def gamma_joint(b, yL, yS, wL, wS, w0, iters=80):
+    """``(MD, V_L, V_S, c)`` per voxel of the gamma model fitted jointly to the linear and the spherical powder signal
+    (``ln S/S0`` ``(..., n_b)`` at ``b`` ms/um^2): one mean diffusivity (the first cumulant, the same for every b-tensor
+    of one trace), one intercept ``c`` (the log of the voxel's b = 0 signal over its measured one) and a variance per
+    encoding, ``y_e(b) = c - ln(1 + r_e b MD) / r_e``, ``V_e = r_e MD^2``; weighted least squares, each point weighted by
+    the rows it averages (``wL``, ``wS`` ``(n_b,)``) and the b = 0 point ``y = 0`` by the b = 0 rows ``w0``. The shared MD
+    and intercept leave the curvature of each curve to its variance alone: the noise of the b = 0 rows (shared by every
+    ``ln S/S0``) is not read as curvature. A vectorised Gauss-Newton over every voxel, MD in [1e-4, 5] um^2/ms, r in
+    [1e-6, 50]."""
+    sh = yL.shape[:-1]; n_b = len(b)
+    Y = np.concatenate([yL.reshape(-1, n_b), yS.reshape(-1, n_b), np.zeros((int(np.prod(sh)), 1))], axis=1)
+    W = np.sqrt(np.concatenate([np.asarray(wL, float), np.asarray(wS, float), [float(w0)]]))
+    n = Y.shape[0]
+    p = np.stack([np.clip(-np.nan_to_num(0.5 * (Y[:, 0] + Y[:, n_b]) / b[0]), 1e-3, None), np.full(n, 0.1), np.full(n, 0.1),
+                  np.zeros(n)], axis=1)
+    lo = np.array([1e-4, 1e-6, 1e-6, -2.0]); hi = np.array([5.0, 50.0, 50.0, 2.0])
+
+    def f(p):
+        g = lambda r: -np.log1p(r[:, None] * b[None] * p[:, 0:1]) / r[:, None]
+        return np.concatenate([p[:, 3:4] + g(p[:, 1]), p[:, 3:4] + g(p[:, 2]), p[:, 3:4]], axis=1)
+    for _ in range(int(iters)):
+        F = f(p); res = (Y - F) * W
+        J = np.empty((n, Y.shape[1], 4))
+        for k in range(4):
+            e = 1e-6 * np.maximum(np.abs(p[:, k]), 1e-3); q = p.copy(); q[:, k] += e
+            J[:, :, k] = (f(q) - F) / e[:, None] * W
+        A = np.einsum("nik,nil->nkl", J, J) + 1e-12 * np.eye(4)
+        g = np.einsum("nik,ni->nk", J, res)
+        with np.errstate(invalid="ignore"):
+            step = np.linalg.solve(A, g[..., None])[..., 0]
+        p = np.clip(p + np.nan_to_num(step), lo, hi)
+    md, rL, rS, c = (p[:, k].reshape(sh) for k in range(4))
+    return md, rL * md ** 2, rS * md ** 2, c
+
+
+def md_ufa(pw, bvals, counts=None, n_b0=1):
+    """``(MD (um^2/ms), muFA)`` per voxel from the powder averages: the gamma model fitted jointly to the linear and the
+    spherical powder signal over the b-values with one MD and one intercept (:func:`gamma_joint`; ``counts``
+    ``{(encoding, b): rows averaged}`` and ``n_b0`` the b = 0 rows weight the points, equal when None), MD the shared first cumulant,
+    ``V_A = V_LTE - V_STE`` the anisotropic variance (the orientation dispersion of the microscopic tensors, ``(4/5)
+    (D_iso D_Delta)^2``), and muFA the FA of the axisymmetric microscopic tensor with ``D_iso = MD`` and that
+    ``D_Delta``; 0 where ``V_A`` is not positive (a spherical signal whose variance exceeds the linear one's: on
+    restricted tissue the two waveforms' spectra differ, and the spherical one sees more of the time dependence), NaN
+    where a signal is not positive."""
+    keys = sorted(float(x) for x in bvals)                                 # the powder averages' own keys, s/mm^2
+    b = np.asarray(keys) * 1e-3                                            # ms/um^2
+    cnt = (lambda e: [1.0] * len(keys)) if counts is None else (lambda e: [float(counts[(e, k)]) for k in keys])
+    with np.errstate(divide="ignore", invalid="ignore"):
+        yL = np.stack([np.log(pw[("LTE", k)]) for k in keys], axis=-1); yS = np.stack([np.log(pw[("STE", k)]) for k in keys], axis=-1)
+        bad = ~np.isfinite(np.concatenate([yL, yS], axis=-1)).all(-1)
+        md, vL, vS, _ = gamma_joint(b, np.nan_to_num(yL), np.nan_to_num(yS), cnt("LTE"), cnt("STE"), 1.0 if counts is None else n_b0)
+    x = np.clip(1.25 * (vL - vS), 0.0, None)                              # (D_iso D_Delta)^2, (um^2/ms)^2
+    with np.errstate(divide="ignore", invalid="ignore"):
+        ufa = np.sqrt(1.5) * np.sqrt(6.0 * x / (3.0 * md ** 2 + 6.0 * x))
+    return np.where(bad, np.nan, md), np.where(bad, np.nan, ufa)
+
+
+def dti_fa(S, enc, bs, axes, b_max=1400.0):
+    """FA per voxel of the diffusion tensor fitted log-linearly to the linear rows at b <= ``b_max`` (s/mm^2) and the
+    b = 0 rows (``S`` S/S0 ``(..., n)``, ``axes`` each row's encoding axis); NaN where a signal is not positive."""
+    m = (np.asarray(enc) == "LTE") & (np.asarray(bs) <= b_max)
+    if m.sum() < 6:
+        raise ValueError(f"a tensor fit needs six or more linear rows at b <= {b_max:g}, got {int(m.sum())}")
+    g = np.asarray(axes, np.float64)[m]; b = np.asarray(bs, np.float64)[m] * 1e-3
+    X = -b[:, None] * np.stack([g[:, 0] ** 2, g[:, 1] ** 2, g[:, 2] ** 2, 2 * g[:, 0] * g[:, 1], 2 * g[:, 0] * g[:, 2],
+                                2 * g[:, 1] * g[:, 2]], 1)
+    with np.errstate(divide="ignore", invalid="ignore"):
+        y = np.log(S[..., m])
+    d = np.einsum("pk,...k->...p", np.linalg.pinv(X), y)
+    D = np.stack([np.stack([d[..., 0], d[..., 3], d[..., 4]], -1), np.stack([d[..., 3], d[..., 1], d[..., 5]], -1),
+                  np.stack([d[..., 4], d[..., 5], d[..., 2]], -1)], -2)
+    ok = np.isfinite(D).all(axis=(-1, -2))
+    ev = np.full(D.shape[:-1], np.nan); ev[ok] = np.linalg.eigvalsh(D[ok])
+    md = ev.mean(-1)
+    with np.errstate(divide="ignore", invalid="ignore"):
+        return np.sqrt(1.5 * ((ev - md[..., None]) ** 2).sum(-1) / (ev ** 2).sum(-1))
+
+
+def maps(dwi, meas, signal):
+    """What a b-tensor run reconstructs, per voxel where ``signal``: ``fa`` (:func:`dti_fa` on the linear rows), ``md``
+    and ``ufa`` (:func:`md_ufa` on the powder averages), and ``powder`` ``{(encoding, b): S/S0}`` with ``bvals``, the
+    b-values the encodings share."""
+    enc, bs = np.asarray(meas.encoding), np.asarray(meas.bvals, np.float64)
+    bvals = sorted(set(float(b) for b, e in zip(bs, enc) if e == "LTE") & set(float(b) for b, e in zip(bs, enc) if e == "STE"))
+    S = np.where(signal[..., None], dwi, np.nan)
+    pw = {(e, b): np.nanmean(S[..., (enc == e) & (bs == b)], axis=-1) for e in ("LTE", "STE") for b in bvals}
+    counts = {(e, b): int(((enc == e) & (bs == b)).sum()) for e in ("LTE", "STE") for b in bvals}
+    md, ufa = md_ufa(pw, bvals, counts, n_b0=int((enc == "b0").sum()))
+    return dict(fa=dti_fa(S, enc, bs, meas.dirs), md=md, ufa=ufa, powder=pw, bvals=bvals)

space/viewers.py

@@ -318,3 +318,52 @@ def response_curves(true, estimated, note=""):
     axes[1].set_title("GM and CSF" + (f" ({note})" if note else ""), fontsize=9); axes[1].legend(fontsize=7)
     fig.tight_layout()
     return _image(fig)
+
+
+def tensor_maps(maps, crossing, z, label="A: "):
+    """A b-tensor run's maps at axial slice ``z``: FA (DTI on the linear rows) beside the microscopic FA, both 0 to 1
+    with the voxels the ground truth's strands cross in two or more directions (``crossing``, boolean on the grid)
+    outlined, and MD (um^2/ms)."""
+    plt = _mpl()
+    fig, axes = plt.subplots(1, 3, figsize=(11.0, 3.6))
+    cr = np.asarray(crossing, float)[:, :, z].T
+    for ax, key, title, kw in ((axes[0], "fa", "FA (DTI, linear rows)", dict(vmin=0, vmax=1, cmap="viridis")),
+                               (axes[1], "ufa", "µFA (linear vs spherical powder)", dict(vmin=0, vmax=1, cmap="viridis")),
+                               (axes[2], "md", "MD (µm²/ms)", dict(vmin=0, vmax=1.0, cmap="magma"))):
+        im = ax.imshow(np.asarray(maps[key], np.float64)[:, :, z].T, origin="lower", **kw)
+        if key != "md":
+            ax.contour(cr, levels=[0.5], colors="w", linewidths=0.6)
+        ax.set_xticks([]); ax.set_yticks([]); ax.set_title(label + title, fontsize=9)
+        fig.colorbar(im, ax=ax, fraction=0.046, pad=0.02)
+    fig.suptitle("white outline: voxels the ground truth's strands cross in two or more directions", fontsize=8)
+    return _image(fig)
+
+
+def powder_curve(maps, ijk, other=None, labels=("A", "B")):
+    """The powder-averaged signal of voxel ``ijk`` against b, linear vs spherical (``other``: B's maps, dashed)."""
+    plt = _mpl()
+    fig, ax = plt.subplots(figsize=(4.6, 3.4))
+    for m, ls, tag in ((maps, "-", labels[0]), (other, "--", labels[1])):
+        if m is None:
+            continue
+        b = [0.0] + list(m["bvals"])
+        for e, c in (("LTE", "C0"), ("STE", "C3")):
+            y = [1.0] + [float(m["powder"][(e, x)][tuple(ijk)]) for x in m["bvals"]]
+            ax.plot(b, y, ls + "o", color=c, ms=3, label=f"{tag} {e}")
+    ax.set_yscale("log"); ax.set_xlabel("b (s/mm²)"); ax.set_ylabel("powder S/S0")
+    ax.set_title(f"voxel {tuple(int(x) for x in ijk)}: µFA {float(maps['ufa'][tuple(ijk)]):.2f}, FA {float(maps['fa'][tuple(ijk)]):.2f}", fontsize=9)
+    ax.legend(fontsize=7)
+    fig.tight_layout()
+    return _image(fig)
+
+
+def map_picker(sl, scale, *, vmin=0.0, vmax=1.0, cmap="viridis"):
+    """A 2-D map ``sl`` (x, y) as an image ``scale`` pixels per voxel with y up and no margins, so a click's pixel names
+    its voxel exactly (``x // scale``, ``ny - 1 - y // scale``); NaN grey."""
+    from PIL import Image
+    import matplotlib
+    a = np.asarray(sl, np.float64).T[::-1]                       # rows top to bottom = y from high to low
+    rgba = matplotlib.colormaps[cmap]((np.clip(np.nan_to_num(a, nan=vmin), vmin, vmax) - vmin) / (vmax - vmin))
+    rgba[~np.isfinite(a)] = (0.5, 0.5, 0.5, 1.0)
+    img = (np.repeat(np.repeat(rgba[..., :3], int(scale), axis=0), int(scale), axis=1) * 255).astype(np.uint8)
+    return Image.fromarray(img)

tests/test_acceptance.py

@@ -1,9 +1,20 @@
-"""The gate: on a GPU with the layout reachable, DiSCo's own protocol replayed from the layout equals the published
-bare reference volume, and the noiseless pipeline scores the connectome within the reference's band with no ground-
-truth pair missed. Prints the per-stage timings (first call and steady) and writes them to ``gate.json``.
+"""The gate: the page's replay from the time-segment layout equals the replay pack's own replay of the same rows, and
+the noiseless pipeline scores the connectome within the reference's band with no ground-truth pair missed.
 
-    JAX_PLATFORMS=cuda XLA_FLAGS=--xla_gpu_deterministic_ops=true pytest tests/test_acceptance.py -s
-    DISCO_MOMENTS=/local/dir pytest ...        # a local copy of the layout instead of the Hub
+The reference (:mod:`space.reference`) is the replay pack the layout was contracted from (``[data] pack``), replayed
+walker by walker -- every band, every field mode, the contact -- at ``[gate] reference_voxels`` voxels spread over its
+index; the page's signal at those voxels, over the voxel's bare b = 0 signal (the quantity both certificates bound;
+S/S0 divides it by the relaxed b = 0 signal, down to 0.3 at 7 T, and is recorded beside it), must be within
+``[gate] reference_max_abs_diff`` of the reference's, for DiSCo's
+own protocol bare and in the catalogue's white matter at 3 T and 7 T, and for every preset (the paper's flagship
+acquisitions, the cosine OGSE among them) bare and at 3 T. It runs where the layout is (``DISCO_LAYOUT`` a local copy,
+else the Hub) and reads the pack's rows at the gate's voxels (``DISCO_PACK`` a local copy, else the Hub); on a CPU
+host every image streams the layout chunk by chunk (minutes per class), so it runs when asked:
+
+    DISCO_GATE=1 DISCO_LAYOUT=/local/layout DISCO_PACK=/local/coefficients pytest tests/test_acceptance.py -s
+    DISCO_GATE=1 JAX_PLATFORMS=cuda XLA_FLAGS=--xla_gpu_deterministic_ops=true pytest tests/test_acceptance.py -s
+
+The per-case worst difference and timings go to ``gate.json`` (``GATE_OUT``).
 """
 import json
 import os
@@ -13,95 +24,100 @@ import numpy as np
 import pytest
 
 from space import pipeline as P
+from space import reference as R
 from space import sources
 from space.sources import disco as D
 
 CFG = P.config()
 OUT = os.environ.get("GATE_OUT", os.path.join(P.ROOT, "gate.json"))
+pytestmark = pytest.mark.skipif(not os.environ.get("DISCO_GATE"), reason="the gate runs when asked (DISCO_GATE=1)")
+
+
+def _record(key, value):
+    table = json.load(open(OUT)) if os.path.exists(OUT) else {}
+    table[key] = value
+    with open(OUT, "w") as f:
+        json.dump(table, f, indent=1, default=float)
 
 
 @pytest.fixture(scope="module")
 def layout():
-    import jax
-    if jax.devices()[0].platform != "gpu":
-        pytest.skip("the gate runs on a GPU")
-    pytest.importorskip("huggingface_hub")
     t0 = time.perf_counter()
     lay = sources.source(CFG)
-    lay.warm()
     lay.load_seconds = time.perf_counter() - t0
     return lay
 
 
-@pytest.mark.parametrize("run", ["bare", "3T", "7T"])
-def test_discos_protocol_from_the_layout_is_the_published_reference_volume(layout, run):
-    """S0-normalised, on every voxel the reference has signal: the replay through the moments is the replay the
-    reference was made with, to the float32 arithmetic of the layout; bare, and in the catalogue's white matter at
-    3 T and 7 T with every tier on and the field along z (the references' setting)."""
-    import nibabel as nib
-    from huggingface_hub import hf_hub_download
-    p, order = D.disco_protocol(CFG)                                   # order: the table's rows in the protocol's order
+@pytest.fixture(scope="module")
+def pack():
+    pk = R.open_pack(os.environ.get("DISCO_PACK") or CFG["data"]["pack"])
+    return pk, R.gate_voxels(pk, CFG["gate"]["reference_voxels"])
+
+
+CASES = [(D.DISCO_364, run) for run in ("bare", "3T", "7T")] + [(name, run) for name in CFG["presets"] for run in ("bare", "3T")]
+
+
+@pytest.mark.parametrize("preset,run", CASES)
+def test_the_layout_is_the_replay_packs_replay_at_the_gates_voxels(layout, pack, preset, run):
+    """At the gate's voxels: the page's replay through the layout (its classes' waveforms, its tiers) is the pack's replay
+    of the same rows, to ``reference_max_abs_diff`` in units of the voxel's bare b = 0 signal."""
+    pk, voxels = pack
+    protocol = D.Disco.protocol(CFG, preset)
     physics = None if run == "bare" else P.Physics.at(float(run[:-1]))
-    m = layout.validate(p, physics)
-    dwi, floor, factor, secs = layout.replay(m, physics)
-    ref_path = hf_hub_download(CFG["data"]["repo"], f"disco/reference/disco_replay_{run}.nii.gz", repo_type="dataset", revision=CFG["data"]["revision"])
-    ref = np.asarray(nib.load(ref_path).dataobj, np.float64)[..., order]
-    both = np.isfinite(dwi) & (ref != 0)
-    d = np.abs(dwi - ref)[both]
-    print(f"\nreference check {run}: max |dS| {d.max():.2e}, 99.9 % {np.quantile(d, 0.999):.2e}, median {np.median(d):.2e} "
-          f"over {both.any(-1).sum()} voxels; replay {secs:.2f} s")
-    assert (~np.isfinite(dwi) & (ref != 0)).any(-1).sum() == 0, "voxels the reference has signal in and the layout has no rows"
-    assert d.max() < CFG["gate"]["reference_max_abs_diff"]
+    meas = layout.validate(protocol, physics)
+    t0 = time.perf_counter()
+    dwi, floor, factor, secs = layout.replay(meas, physics)
+    tissue, scanner = P.tissue_and_scanner(physics, layout.unseeded)
+    ref = np.zeros((len(voxels), len(meas.bvals)))
+    t1 = time.perf_counter()
+    for name in np.unique(meas.shape):
+        rows = meas.shape == name
+        ref[:, rows] = R.replay_at(pk, voxels, layout.sequence(name, meas.bvals[rows], meas.dirs[rows], None if meas.encoding is None else meas.encoding[rows]), tissue, scanner)
+    t_ref = time.perf_counter() - t1
+    ref_S0 = ref[:, meas.b0].mean(1, keepdims=True)
+    got = np.stack([dwi[tuple(v)] for v in voxels])
+    d = np.abs(got - ref / ref_S0)                                   # S/S0
+    d_S = np.abs(got * ref_S0 - ref)                                 # the signal over the voxel's bare b = 0 signal: what the certificates bound
+    worst, worst_S = float(np.nanmax(d)), float(np.nanmax(d_S))
+    print(f"\n{preset} {run}: max |dS| {worst_S:.2e} (S over the bare b = 0 signal), max |d(S/S0)| {worst:.2e} over "
+          f"{len(voxels)} voxels x {len(meas.bvals)} measurements (S0 down to {float(ref_S0.min()):.3f}); "
+          f"layout replay {secs:.1f} s, reference {t_ref:.1f} s")
+    _record(f"{preset} | {run}", dict(worst_S=worst_S, worst_S_over_S0=worst, median_S=float(np.nanmedian(d_S)), S0_min=float(ref_S0.min()),
+                                      voxels=len(voxels), n_meas=int(len(meas.bvals)), replay_seconds=secs, reference_seconds=t_ref,
+                                      bound=CFG["gate"]["reference_max_abs_diff"], manifest_sha256=layout.manifest_sha256))
+    assert np.isfinite(got).all(), "a gate voxel the pack has rows in has no signal from the layout"
+    assert worst_S < CFG["gate"]["reference_max_abs_diff"]
 
 
 def test_the_noiseless_pipeline_scores_within_the_references_band(layout):
-    """Pearson vs strand count at or above the gate's floor (the replay reference's band is 0.912 +- 0.003; the
-    tracker on the reference volume measured 0.927) and no ground-truth pair missed; the timings recorded."""
+    """Pearson vs strand count at or above the gate's floor and no ground-truth pair missed on DiSCo's own protocol;
+    the timings recorded (first call with every compile, then steady)."""
     p, _ = D.disco_protocol(CFG)
     tr = P.Tracking(**{k: CFG["tracking"][k] for k in ("density", "step_mm", "max_angle", "max_steps")})
-    first = P.run(layout, p, snr=None, tracking=tr)                    # includes every compile
+    first = P.run(layout, p, snr=None, tracking=tr)
     steady = P.run(layout, p, snr=None, tracking=tr)
     s = steady.score
-    table = dict(device=str(__import__("jax").devices()[0]), layout_load_seconds=layout.load_seconds,
-                 first_call=first.seconds, steady=steady.seconds, score=s, streamlines=len(steady.tractogram),
-                 seeds=int(len(steady.seeds)), n_meas=p.n_meas)
-    print("\ngate: " + json.dumps(table, indent=1))
-    with open(OUT, "w") as f:
-        json.dump(table, f, indent=1)
+    _record("pipeline", dict(layout_load_seconds=layout.load_seconds, first_call=first.seconds, steady=steady.seconds, score=s,
+                             streamlines=len(steady.tractogram), seeds=int(len(steady.seeds)), n_meas=p.n_meas))
     assert s["pearson_count"] >= CFG["gate"]["pearson_count_min"], s
     assert s["missed_pairs"] <= CFG["gate"]["missed_pairs_max"], s
-    assert first.score["pearson_count"] == s["pearson_count"]   # the same key: the same tractogram (deterministic XLA ops)
+    assert first.score["pearson_count"] == s["pearson_count"]
 
 
 @pytest.mark.parametrize("machine", [m["key"] for m in CFG["scanners"]["machines"]])
-def test_a_machine_on_discos_protocol_moves_the_image_and_is_timed(layout, machine):
-    """Each machine of the menu at 7.9 cm from isocentre along R-L on DiSCo's own protocol (none of them can play it:
-    the gradient limit is the page's, the replay's cost is measured here at the gate's protocol), against the ideal
-    scanner at the machine's field and direction: the delivered image is finite where the ideal one is and moves it;
-    the first and steady replay seconds of both go to ``gate_scanners.json`` beside ``gate.json``. A machine with its
-    own gradient needs the layout's background moments and is skipped by name on a layout without them."""
+def test_a_machine_on_discos_protocol_moves_the_image(layout, machine):
+    """Each machine of the menu at 7.9 cm from isocentre along R-L on DiSCo's own protocol against the ideal scanner at
+    the machine's field and direction: the delivered image is finite where the ideal one is and moves it."""
     from dataclasses import replace
     L = P.limits(machine)
-    if L.has_field_law and not layout.moments.background:
-        pytest.skip(f"{machine}'s own gradient needs the layout's background moments (shape_moments.stamp_background)")
     p, _ = D.disco_protocol(CFG)
     ideal = P.Physics.at(float(L.field_T), b0_direction=tuple(L.b0_axis or (0.0, 0.0, 1.0)))
     far = replace(ideal, scanner=machine, offset_m=(0.079, 0.0, 0.0))
     m = layout.validate(p, far)
-    secs = {}
-    for tag, ph in (("ideal", ideal), ("machine", far)):
-        for k in ("first", "steady"):
-            dwi, _, factor, s = layout.replay(m, ph)
-            secs[f"{tag}_{k}"] = s
-        if tag == "ideal":
-            ref = dwi
+    ref, *_ = layout.replay(m, ideal)
+    dwi, _, _, secs = layout.replay(m, far)
     both = np.isfinite(ref)
     assert (np.isfinite(dwi) == both).all()
     moved = float(np.max(np.abs(dwi - ref)[both]))
-    print(f"\n{machine} at 7.9 cm on DiSCo 364: max |dS| vs ideal {moved:.3e}; seconds {json.dumps(secs)}")
+    _record(f"machine | {machine}", dict(max_abs_dS=moved, seconds=secs))
     assert moved > 1e-3
-    path = os.path.join(os.path.dirname(OUT), "gate_scanners.json")
-    table = json.load(open(path)) if os.path.exists(path) else {}
-    table[machine] = dict(seconds=secs, max_abs_dS=moved, n_meas=p.n_meas, offset_m=list(far.offset_m), device=str(__import__("jax").devices()[0]))
-    with open(path, "w") as f:
-        json.dump(table, f, indent=1)

tests/test_app.py

@@ -8,6 +8,7 @@ import pytest
 
 from space import app as A
 from space import pipeline as P
+from space import viewers as V
 from space.sources import disco as D
 
 CFG = P.config()
@@ -21,8 +22,8 @@ def test_the_three_kinds_of_acquisition():
     assert disco.name == "DiSCo 364" and disco.n_meas == 364 and disco.n_b0 == 4        # the table's own, not the slider's
     name = next(iter(CFG["presets"]))
     preset = A._protocol_from_inputs(S, CFG, name, 5, *ROWS)
-    assert preset.name == name and preset.n_b0 == CFG["presets"][name]["n_b0"]
-    assert [s.shape for s in preset.shells] == [s["shape"] for s in CFG["presets"][name]["shells"]]
+    assert preset.name == name and preset.n_b0 == P.preset_protocol(CFG, name).n_b0 == 3        # the table's own
+    assert {s.shape for s in preset.shells} == {CFG["presets"][name]["shape"]}
     custom = A._protocol_from_inputs(S, CFG, A.CUSTOM, 2, *ROWS)
     assert custom.n_b0 == 2 and [(s.shape, s.b, s.n_dirs) for s in custom.shells] == [("d12-D24", 1000.0, 30), ("d8-D20", 3000.0, 45)]
     assert not any(s.free_timing for s in custom.shells)
@@ -68,7 +69,8 @@ def test_a_knob_changes_one_thing_of_a():
     v = _values()
     prot = A._protocol_from_inputs(S, CFG, A.CUSTOM, 2, *ROWS)
     ks = S.knobs(CFG)
-    assert ks[A.NO_KNOB] is None and len(ks) == 1 + 5 + 2 + 3 + 1 + 3 + len(CFG["shapes"])
+    timing = [n for n, sh in CFG["shapes"].items() if sh.get("kind") != "tensor"]
+    assert ks[A.NO_KNOB] is None and len(ks) == 1 + 5 + 2 + 3 + 1 + 3 + len(timing)
     for name, change in ks.items():
         if change is None:
             continue
@@ -215,11 +217,14 @@ def test_the_estimated_seconds_grow_with_the_run_and_stay_in_the_pools_window():
     small = A.estimated_seconds("clinical b1000 x 30", 1, True, 30, 2, 30.0, 0.5, 0, None, A.NO_KNOB, A.IDEAL, 1, False, *v, *ROWS)
     assert 30 <= small < one < two <= 480 and one < keys <= 480 and one < ladder <= 480
     both = A.estimated_seconds("DiSCo 364", 1, True, 30, 2, 30.0, 0.5, 0, None, "SNR → 10", A.IDEAL, 1, True, *v, *ROWS)
-    assert one < ladder < both <= 120          # DiSCo alone, with the ladder, and A + B + ladder all fit a logged-out visitor's quota
+    assert one <= 120 < ladder < both <= 480   # DiSCo alone fits a logged-out visitor's quota; the ladder and B need an account
+    bt = next(n for n, p in CFG["presets"].items() if "tensor" in p)
+    tensor = A.estimated_seconds(bt, 1, True, 100, 2, 30.0, 0.5, 0, None, A.NO_KNOB, A.IDEAL, 1, False, *v, *ROWS)
+    assert small < tensor <= 120               # 77 rows on six terms: 23 s of replay measured on the pool
     assert A.estimated_seconds("nonsense", 1, True, 30, 2, 30.0, 0.5, 0, None, A.NO_KNOB, A.IDEAL, 1, True, *v, *ROWS) == 480
     text = A.gpu_seconds_text("DiSCo 364", 1, True, 30, 2, 30.0, 0.5, 0, None, A.NO_KNOB, A.IDEAL, 1, True, *v, *ROWS)
-    assert f"reserves {ladder} s" in text and "logged out" in text.split("started by")[1].split(".")[0]
-    assert len(A.OUTPUTS) == 30
+    assert f"reserves {ladder} s" in text and "logged out" not in text.split("started by")[1].split(".")[0]
+    assert len(A.OUTPUTS) == 36                 # the 30 of a connectome run and the b-tensor run's six
 
 
 def test_the_explorer_draws_layers_differences_and_metrics_from_its_state():
@@ -230,15 +235,16 @@ def test_the_explorer_draws_layers_differences_and_metrics_from_its_state():
     rng = np.random.default_rng(0)
     bare = rng.uniform(0.3, 1.0, (4, 4, 2, 13)).astype(np.float16); a = (bare * 0.98).astype(np.float16); b = (bare * 0.97).astype(np.float16)
     md = rng.random((4, 4, 2)).astype(np.float16); fa = rng.random((4, 4, 2)).astype(np.float16)
-    ex = dict(layers=[("bare diffusion", bare), ("A = bare + field", a), ("B = A with field → 7 T", b)], metrics={k: (md, fa) for k in ("bare diffusion", "A = bare + field", "B = A with field → 7 T")},
+    ex = dict(layers=[("bare diffusion", bare), ("A = bare + field", a), ("B = A with field → 7 T", b)], metrics={k: (md, fa, None) for k in ("bare diffusion", "A = bare + field", "B = A with field → 7 T")},
               floor=np.full((4, 4, 2), 0.01, np.float32), meas=m, mask=np.ones((4, 4, 2), bool),
               ingredients=S.ingredient_layers(dict(intra_fraction=rng.random((4, 4, 2)), wall_contact_um=rng.random((4, 4, 2)), contact_survival=None, field_rad=None, D_walk=6e-10)),
               differences=P.layer_differences([("bare diffusion", bare.astype(np.float32)), ("A", a.astype(np.float32)), ("B", b.astype(np.float32))], m, np.ones((4, 4, 2), bool)),
               snr=None, physics=None)
     for mode in A.EXPLORE_MODES:
-        for metric in A.METRICS:
+        for metric in A.METRICS[:3]:
             img = A.explore(ex, "A = bare + field", mode, metric, 1, 5)
             assert img is not None and img.size[0] > 100, (mode, metric)
+    assert A.explore(ex, "A = bare + field", "signal", A.METRICS[3], 1, 5) is None        # no µFA off a b-tensor run
     assert A.explore(ex, "bare diffusion", "minus the previous layer", A.METRICS[0], 0, 0) is not None
     pool, contact, fld = A.ingredient_views(ex, 1)
     assert pool is not None and contact is not None and fld is None
@@ -283,3 +289,82 @@ def test_the_payload_carries_float32_volumes():
     assert light.dwi.dtype == np.float32 and light.clean.dtype == np.float32 and res.dwi.dtype == np.float64
     np.testing.assert_allclose(light.dwi, dwi, rtol=1e-6); assert light.snr == 30.0 and light.meas is m
     assert A._light(P.Result(p, m, dwi, floor, None, None, None, None, None, {})).clean is None
+
+
+def _tensor_rs():
+    """A results state of a b-tensor run on a 6 x 5 x 2 grid: A and B, their maps, the measurement rows."""
+    p = P.preset_protocol(CFG, next(n for n, v in CFG["presets"].items() if "tensor" in v))
+    m = P.measurements(p, CFG["shapes"])
+    rng = np.random.default_rng(1)
+    bvals = sorted(set(float(b) for b in m.bvals if b > 0))
+
+    def maps(scale):
+        return dict(fa=rng.random((6, 5, 2)) * scale, md=rng.random((6, 5, 2)), ufa=rng.random((6, 5, 2)) * scale, bvals=bvals,
+                    powder={(e, b): np.exp(-b * 1e-3 * (0.5 if e == "LTE" else 0.7)) * np.ones((6, 5, 2)) for e in ("LTE", "STE") for b in bvals})
+    dwi = rng.uniform(0.2, 1.0, (6, 5, 2, len(m.bvals))).astype(np.float32)
+    run = lambda s: dict(dwi=dwi, floor=np.full((6, 5, 2), 0.01, np.float32), floor_median=0.01, meas=m, peaks=None, peak_amp=None,
+                         extras=dict(fractions=None), tensor=A._tensor_state(maps(s)))
+    return dict(A=run(1.0), B=run(0.5)), m
+
+
+def test_a_b_tensor_runs_rows_maps_and_voxel_click():
+    """On a b-tensor run the measurement slider counts within the chosen encoding; the µFA slice is a picker whose
+    click pixel names its voxel (y up); the clicked voxel's powder curves draw for A and B; the maps draw for both."""
+    rs, m = _tensor_rs()
+    enc = np.asarray(m.encoding)
+    assert A.row_of(m, "every row", 3) == 3
+    assert enc[A.row_of(m, "STE", 0)] == "STE" and enc[A.row_of(m, "LTE", 10 ** 6)] == "LTE"
+    assert A.row_of(m, "STE", 0) == int(np.flatnonzero(enc == "STE")[0])
+    pick = V.map_picker(rs["A"]["tensor"]["ufa"][:, :, 1], A.PICK)
+    assert pick.size == (6 * A.PICK, 5 * A.PICK)
+    seen = []
+    import space.viewers as VV
+    orig = VV.powder_curve
+    try:
+        VV.powder_curve = lambda tm, ijk, other=None: seen.append((tuple(ijk), other is not None)) or "img"
+        evt = types.SimpleNamespace(index=(2 * A.PICK + 3, 0))                   # x in voxel 2, the top row = y 4
+        assert A.pick_voxel(rs, 1, evt) == "img" and seen[-1] == ((2, 4, 1), True)
+    finally:
+        VV.powder_curve = orig
+    assert V.powder_curve(rs["A"]["tensor"], (1, 1, 0), other=rs["B"]["tensor"]).size[0] > 100
+    n_dirs = np.ones((6, 5, 2), int); n_dirs[2:4, 1:3] = 2
+    A._CROSSINGS["c"] = (n_dirs, np.ones((6, 5, 2)))
+    try:
+        a, picker, b = A.tensor_views(rs, 0)
+        assert a is not None and picker is not None and b is not None
+    finally:
+        A._CROSSINGS.clear()
+    assert A.tensor_views(dict(A=dict(tensor=None)), 0) == (None, None, None)
+
+
+def test_the_explorer_on_a_b_tensor_run_draws_the_encodings_rows_and_b_minus_a_of_mufa():
+    rs, m = _tensor_rs()
+    a = rs["A"]["tensor"]; b = rs["B"]["tensor"]
+    ex = dict(layers=[("A = bare diffusion", rs["A"]["dwi"].astype(np.float16)), ("B = A with SNR → 10", rs["B"]["dwi"].astype(np.float16))],
+              metrics={"A = bare diffusion": (a["md"], a["fa"], a["ufa"]), "B = A with SNR → 10": (b["md"], b["fa"], b["ufa"])},
+              floor=np.full((6, 5, 2), 0.01, np.float32), meas=m, mask=np.ones((6, 5, 2), bool), ingredients=[None, None, None],
+              differences=[], snr=None, physics=None)
+    for metric in A.METRICS:
+        for mode in A.EXPLORE_MODES:
+            assert A.explore(ex, "A = bare diffusion", mode, metric, 1, 2, "STE") is not None, (metric, mode)
+    assert A.explore(ex, "A = bare diffusion", "B − A", A.METRICS[3], 0, 0) is not None
+    with pytest.raises(ValueError, match="encodings of its one class"):
+        P.retime(P.preset_protocol(CFG, next(n for n, v in CFG["presets"].items() if "tensor" in v)), "d12-D24")
+
+
+def test_every_stage_of_a_run_has_its_text():
+    """The progress text names every stage either kind of run passes (the connectome's and the b-tensor's)."""
+    texts = S.describe(CFG)["stages"]
+    assert set(P.TENSOR_STAGES) <= set(texts) and set(D.Disco.STAGES) <= set(texts)
+
+
+def test_a_b_tensor_headline_reads_the_strand_voxels_by_crossing_count(monkeypatch):
+    rs, m = _tensor_rs()
+    n_dirs = np.zeros((6, 5, 2), int); n_dirs[:2] = 1; n_dirs[2:4] = 2; n_dirs[4:] = 3
+    length = np.ones((6, 5, 2)); length[0] = 0
+    monkeypatch.setitem(A._CROSSINGS, "c", (n_dirs, length))
+    tm = rs["A"]["tensor"]
+    res = types.SimpleNamespace(extras=dict(tensor=tm), physics=None, dwi=rs["A"]["dwi"], snr=30.0, seconds=dict(total=1.0),
+                                protocol=types.SimpleNamespace(n_meas=len(m.bvals)), floor=rs["A"]["floor"])
+    text = A.tensor_text("A", res)
+    assert "voxels the strands cross" in text and "2 directions" in text and "3+ directions" in text and "SNR 30" in text

tests/test_brain.py

@@ -1,5 +1,5 @@
 """The brain source: the asset's layout (the BATMAN fixture of ``tools/build_brain_fixture.py``, local only), the
-kernel-route replay against dmipy-sim's own ``Phantom.compose(...).replay`` (the oracle), the connectome of a known
+kernel-route replay against dmipy-replay's own ``Phantom.compose(...).replay`` (the oracle), the connectome of a known
 tractogram, the reservation model, the panel and its knobs, the refusals, and the page's chain on a small crop of the
 fixture. Skipped where the fixture or the packs are not on this machine; ``DISCO_BRAIN_FIXTURE`` and
 ``DISCO_BRAIN_PACKS`` point at them."""
@@ -12,12 +12,12 @@ import pytest
 from space import pipeline as P
 from space.sources import brain as B
 
-FIXTURE = os.environ.get("DISCO_BRAIN_FIXTURE", os.path.expanduser("~/dmrai-ws/data/batman/brain_fixture"))
-PACKS = os.environ.get("DISCO_BRAIN_PACKS", os.path.expanduser("~/dmrai-ws/packs"))
+FIXTURE = os.environ.get("DISCO_BRAIN_FIXTURE", "")          # the directory tools/build_brain_fixture.py wrote
+PACKS = os.environ.get("DISCO_BRAIN_PACKS", "")              # the directory holding the WM and GM replay packs
 WM_PACK = os.path.join(PACKS, "single_bundle_1s_c3_seg125ms.rpk")
 GM_PACK = os.path.join(PACKS, "gm_spheres_100ms.rpk")
-has_fixture = pytest.mark.skipif(not os.path.exists(os.path.join(FIXTURE, "manifest.json")), reason="the BATMAN brain fixture is not on this machine")
-has_packs = pytest.mark.skipif(not (os.path.exists(WM_PACK) and os.path.exists(GM_PACK)), reason="the WM and GM packs are not on this machine")
+has_fixture = pytest.mark.skipif(not os.path.exists(os.path.join(FIXTURE, "manifest.json")), reason="DISCO_BRAIN_FIXTURE does not name a BATMAN brain fixture")
+has_packs = pytest.mark.skipif(not (os.path.exists(WM_PACK) and os.path.exists(GM_PACK)), reason="DISCO_BRAIN_PACKS does not name a directory with the WM and GM packs")
 
 
 def brain_cfg(asset=FIXTURE):
@@ -63,9 +63,9 @@ def test_the_fixture_loads_in_the_assets_layout(brain):
 
 def compose_reference(brain, meas, ph, *, scanner, n=200):
     """``(got, ref)`` on ``n`` voxels drawn with seed 0, in M0 units: the page's contraction (:meth:`Brain.prepare`,
-    :meth:`Brain.replay`) and dmipy-sim's ``Phantom.compose(...).replay(seq, scanner=scanner, pose=...)`` of the
+    :meth:`Brain.replay`) and dmipy-replay's ``Phantom.compose(...).replay(seq, scanner=scanner, pose=...)`` of the
     same packs as loaded, tissue, fractions and FODs, NaN where a voxel carries no signal."""
-    from dmipy_sim.phantom import FreeWater, Grid, Inert, ODF, PackSubstrate, Phantom
+    from dmipy_replay.phantom import FreeWater, Grid, Inert, ODF, PackSubstrate, Phantom
     dwi, floor, factor, _ = brain.replay(meas, ph, brain.prepare(meas, ph))
     assert floor.max() == 0 and np.nanmax(dwi) <= 1 + 1e-5
     pick = np.sort(np.random.default_rng(0).choice(len(brain.vox), n, replace=False))
@@ -90,7 +90,7 @@ def compose_reference(brain, meas, ph, *, scanner, n=200):
 
 @has_fixture
 @has_packs
-def test_the_replay_is_dmipy_sims_phantom_composition(brain):
+def test_the_replay_is_dmipy_replays_phantom_composition(brain):
     """At the scan's protocol, 3 T with every tier on (relaxation, contact, the sheath's field along the bore), noise
     off: the Space's float32 contraction over 200 voxels equals ``Phantom.compose(...).replay(seq, pose=...)`` of the
     same packs, tissue, fractions and FODs to 1e-5 in M0 units, the field given as ``scanner=ph.scanner_field``."""
@@ -106,13 +106,13 @@ def test_the_replay_is_dmipy_sims_phantom_composition(brain):
 
 @pytest.mark.skipif(not os.environ.get("DISCO_HUB_TESTS"), reason="reads the Hub: set DISCO_HUB_TESTS=1")
 @has_fixture
-def test_a_windowed_multi_seed_pack_replays_as_dmipy_sims_phantom_composition():
+def test_a_windowed_multi_seed_pack_replays_as_dmipy_replays_phantom_composition():
     """The page's default packs (brain.toml's first menu entries: the WM walk of 1 s in 125 ms windows from several
     seeded walks, the GM walk of 250 ms in 125 ms windows), read from the Hub by window, every tier on at 3 T: a
     class within window 0 loads ``windows=range(1)`` and the page's contraction over 200 voxels equals
     ``Phantom.compose(...).replay`` of the same packs to float precision (1e-6 in M0 units; measured 7e-8). The
     sheath's field is the composition's ``scanner=ph.scanner_field``; the same composition without it
-    (``scanner=None``, the ideal machine's key ``ph.scanner``) is refused by dmipy-sim, since the tissue's
+    (``scanner=None``, the ideal machine's key ``ph.scanner``) is refused by dmipy-replay, since the tissue's
     susceptibility needs a field to act in (disco-space#36)."""
     cfg = P.config(os.path.join(P.HERE, "brain.toml"))
     cfg["asset"] = {"local": FIXTURE}
@@ -428,16 +428,19 @@ def test_the_asset_as_the_masivar_build_writes_it_loads(tmp_path):
 
 @pytest.fixture(scope="module")
 def two_windows(tmp_path_factory):
-    """A cylinder walk of 200 ms stored in two 100 ms windows (300 walkers, every tier, built with the bank) and
-    the same pack truncated to its first window (no re-encode), both on disk: ``(two, one)`` paths."""
+    """A cylinder walk of 200 ms stored in two 100 ms windows (300 walkers, built with the bank for the Connectome 2.0
+    class) and the same pack truncated to its first window (no re-encode), both on disk: ``(two, one)`` paths."""
     import dmipy_sim as d
-    from dmipy_sim.replay import read_rpk
+    from dmipy_replay import envelope_of, read_rpk
     from dmipy_sim.replay.bank import build_replay_pack
+    from dmipy_sim.replay.reencode import truncate_pack
     out = tmp_path_factory.mktemp("windows")
-    walk = d.simulate_trajectories(300, 2e-9, d.Cylinder(radius=3e-6, orientation=(0.0, 0.0, 1.0)), 0.2, 2e-4, seed=0, require_gpu=False)
+    envelope = envelope_of("connectome_2", B0_list=(3.0,), rho2_max=1e-4)
+    walk = d.walk(d.spec.spec_of(d.Cylinder(radius=3e-6, orientation=(0.0, 0.0, 1.0))), 300, 0.2, diffusivity=2e-9, seed=0,
+                  K=48, window_T=0.1, require_gpu=False, envelope=envelope)
     two = str(out / "two.rpk"); one = str(out / "one.rpk")
-    build_replay_pack(walk, id="test/two-windows", license="x", citation="x", K=48, out_path=two, segment_T=0.1)
-    read_rpk(two).truncate(1, out_path=one)
+    build_replay_pack(walk, id="test/two-windows", license="x", citation="x", K=48, window_T=0.1, out_path=two)
+    truncate_pack(read_rpk(two), 1, out_path=one)
     return two, one
 
 
@@ -447,6 +450,10 @@ def windowed(two_windows):
     if not os.path.exists(os.path.join(FIXTURE, "manifest.json")):
         pytest.skip("the BATMAN brain fixture is not on this machine")
     two, one = two_windows
+    try:
+        B.declared_segments(two)
+    except ValueError as e:                  # a pack of this release records its windows {n, T}; the source reads the table {n, n_t, T}
+        pytest.skip(f"the brain source does not read this release's window record yet: {e}"[:300])
     cfg = brain_cfg()
     menu = [{"label": "two", "uri": two}, {"label": "one", "uri": one}]
     cfg["packs"] = {"wm": list(menu), "gm": list(menu)}
@@ -495,8 +502,8 @@ def test_a_hub_pack_loads_the_windows_reached_and_reloads_for_more(windowed, two
     loads nothing, a request for two reloads with ``range(2)`` and replaces the cached pack; a local path loads whole.
     (The library's window-ranged read is stood in for by the fixture stamped with the windows asked for;
     ``test_the_hub_fixture_loads_by_window`` reads the Hub.)"""
-    import dmipy_sim.phantom as phantom
-    from dmipy_sim.replay import read_rpk
+    import dmipy_replay.phantom as phantom
+    from dmipy_replay import read_rpk
     two, _ = two_windows
     uri = "hf://SubstrateCommons/test/packs/two.rpk"
     calls = []
@@ -532,8 +539,7 @@ def test_a_hub_pack_loads_by_window(windowed, tmp_path, monkeypatch):
     byte range into an empty Hub cache (the library keeps a window read there and serves a later, smaller request
     from it): ``windows=range(1)`` for a class within window 0, then ``range(2)`` when a class reaching window 1 asks,
     which replaces the cached pack."""
-    import huggingface_hub.constants
-    monkeypatch.setattr(huggingface_hub.constants, "HF_HUB_CACHE", str(tmp_path))
+    monkeypatch.setenv("DMIPY_REPLAY_CACHE", str(tmp_path))
     cfg = P.config(os.path.join(P.HERE, "brain.toml"))
     label, uri = next(iter(B.Brain.pack_menu(cfg, "gm").items()))
     src = B.Brain.__new__(B.Brain)
@@ -643,7 +649,7 @@ def test_the_warm_set_has_one_entry_per_response_and_the_high_fields_last():
 @has_fixture
 @has_packs
 def test_a_machine_plays_every_pulse_at_its_slew_and_field():
-    """``BrainPhysics.scanner`` names a machine of dmipy-sim's catalogue: the pulses are trapezoids with ramps at its
+    """``BrainPhysics.scanner`` names a machine of dmipy-replay's catalogue: the pulses are trapezoids with ramps at its
     slew (the same b), its field and direction replace the panel's, the response key tells the machines apart, the label
     names it, a magnet with its own gradient (the Swoop) is refused by name, and a name outside the catalogue too."""
     cfg = brain_cfg()
@@ -723,7 +729,7 @@ def test_a_left_handed_brains_bore_maps_are_the_catalogues_field_at_its_affines_
     its transmit scale and its own gradient, the Prisma's gradient nonlinearity -- at every sampled brain voxel equal
     the catalogue's law at ``affine @ [ijk, 1]`` from the head centre, on both sides of the midline (a mirrored head
     would read the law at the wrong side's places); its left-handed and right-handed copies give every voxel the same values."""
-    from dmipy_sim.phantom.bore import b0_offset_map, b1_scale_map, background_gradient_map, gradient_tensor_map
+    from dmipy_replay.phantom.bore import b0_offset_map, b1_scale_map, background_gradient_map, gradient_tensor_map
     asset = left_handed_asset(tmp_path)
     left = B.Brain(brain_cfg(asset), asset=asset)
     right = B.Brain(brain_cfg(), asset=FIXTURE)
@@ -760,9 +766,9 @@ def test_a_machine_plays_each_voxel_its_encoding_class_and_a_representative_exac
     within the configured tolerance of its class's), and a class's first member -- the voxel whose exact values the
     class was composed from -- gets exactly what its own delivered acquisition, expanded by the packs and composed with
     its FOD and fractions, gives: the replay to 1e-5 in M0 units; the machine moves the image off the ideal scanner's."""
-    from dmipy_sim.phantom import FreeWater
-    from dmipy_sim.phantom.bore import encoding_classes
-    from dmipy_sim.acquisition.waveforms import rotate_waveform
+    from dmipy_replay.acquisition.waveforms import rotate_waveform
+    from dmipy_replay.phantom import FreeWater
+    from dmipy_replay.phantom.bore import encoding_classes
     cfg = brain_cfg()
     src = B.Brain(cfg, asset=FIXTURE)
     values = default_values(cfg)
@@ -793,8 +799,8 @@ def test_a_machine_plays_each_voxel_its_encoding_class_and_a_representative_exac
     for v in reps:
         _, own = encoding_classes(P.limits(ph.scanner), grid, src.sequence(meas, rows, np.eye(3), ph.slew_rate), idx[v:v + 1], tolerance=None)
         seq = rotate_waveform(own[0], pose)
-        Cw = wm_pack.pose_responses([seq], tissue=t_wm, scanner=ph.scanner_field, pose=pose, keep=(B.LMAX, 0))[0].retained(B.LMAX, 0) @ B.axis_map()
-        eg = gm_pack.pose_responses([seq], tissue=t_gm, scanner=ph.scanner_field, pose=pose, keep=(0, 0))[0].retained(0, 0)[:, 0] * B.axis_map()[0, 0] * B.C00
+        Cw = wm_pack.pose_responses([seq], tissue=t_wm, scanner=ph.scanner_field, pose=pose, keep=(B.LMAX, 0), backend=B.POSE_BACKEND)[0].retained(B.LMAX, 0) @ B.axis_map()
+        eg = gm_pack.pose_responses([seq], tissue=t_gm, scanner=ph.scanner_field, pose=pose, keep=(0, 0), backend=B.POSE_BACKEND)[0].retained(0, 0)[:, 0] * B.axis_map()[0, 0] * B.C00
         ec = FreeWater(m0=1.0, tissue=t_csf).response(seq)
         f = src.f[v].astype(np.float64); m0 = ph.m0
         ref = np.abs(f[0] * m0["wm"] * (src.fod_unit[v].astype(np.float64) @ Cw.T) + f[1] * m0["gm"] * eg + f[2] * m0["csf"] * ec)
@@ -817,4 +823,113 @@ def test_the_warm_set_holds_the_scan_on_every_machine_that_can_play_it():
     got = {ph.scanner: ph for m, ph in src.warm_entries() if ph.scanner is not None}
     assert set(got) == can and can
     for key, ph in got.items():
-        assert ph.field_T == pytest.approx(P.limits(key).field_T) and ph.slew_rate == pytest.approx(P.limits(key).slew_max)
+        assert ph.field_T == pytest.approx(P.limits(key).field_T) and ph.slew_rate == pytest.approx(P.limits(key).slew_max)
+
+# ---- the page when the Hub does not serve what the configuration names -------------------------------------------------
+
+BRAIN_TOML = P.config(os.path.join(P.HERE, "brain.toml"))
+
+
+def _not_found():
+    import urllib.error
+    return urllib.error.HTTPError("https://huggingface.co/datasets/o/n/resolve/main/p.rpk", 404, "Not Found", {}, None)
+
+
+def test_a_pack_is_refused_by_name_when_it_is_absent_refused_or_unreadable(monkeypatch, tmp_path):
+    """The three ways a pack of the menu cannot be served, each named: absent from its dataset (HTTP 404), refused by
+    the reader (its message, which names what it does not read), and a window record this page does not read; a pack
+    that opens and has its table is served; an error that says nothing about the pack (HTTP 500) is raised, not named."""
+    from dmipy_replay import ReplayPack
+    gone, old, bare, ok, down = (f"hf://o/n/{n}.rpk" for n in ("gone", "old", "bare", "ok", "down"))
+
+    def opened(uri, workers=8):
+        if uri == gone:
+            raise _not_found()
+        if uri == old:
+            raise ValueError(f"{uri}: This pack carries `walk_params.n_t`, which this revision does not read; re-walk the pack")
+        if uri == down:
+            import urllib.error
+            raise urllib.error.HTTPError(uri, 500, "Server Error", {}, None)
+        return object()
+
+    def table(uri):
+        if uri == bare:
+            raise ValueError(f"{uri}: no declared segment table (walk_params.segments, RPK.md 4.3) in its record")
+        return dict(n=1, n_t=1001, T=0.1)
+
+    monkeypatch.setattr(ReplayPack, "open", staticmethod(opened))
+    monkeypatch.setattr(B, "declared_segments", table)
+    assert B.pack_refusal(gone) == "it is not in its dataset on the Hub"
+    assert "the reader refuses it" in B.pack_refusal(old) and "walk_params.n_t" in B.pack_refusal(old)
+    assert "cannot read its window record" in B.pack_refusal(bare) and "walk_params.segments" in B.pack_refusal(bare)
+    assert B.pack_refusal(ok) is None
+    with pytest.raises(Exception, match="500"):
+        B.pack_refusal(down)
+    local = tmp_path / "p.rpk"
+    assert B.pack_refusal(str(local)) == "there is no such file"
+    local.write_bytes(b"x")
+    assert B.pack_refusal(str(local)) is None
+
+
+def test_the_asset_is_refused_by_name_when_its_revision_or_manifest_is_gone(monkeypatch, tmp_path):
+    """A Hub revision the Hub answers 404 for is named; so is a directory without its manifest; an unreachable Hub is raised."""
+    import types
+    where = "hf://SubstrateCommons/masivar-brain@7556d3900adfe73f2b78b34fffb74902354775e9"
+
+    class Answer(OSError):
+        def __init__(self, status):
+            super().__init__(f"HTTP {status}")
+            self.response = types.SimpleNamespace(status_code=status)
+
+    def gone(w):
+        raise Answer(404)
+
+    monkeypatch.setattr(B, "manifest", gone)
+    assert B.asset_refusal(where) == "the Hub has no manifest.json of SubstrateCommons/masivar-brain at revision 7556d3900adfe73f2b78b34fffb74902354775e9"
+    monkeypatch.setattr(B, "manifest", lambda w: (_ for _ in ()).throw(Answer(503)))
+    with pytest.raises(Answer):
+        B.asset_refusal(where)
+    monkeypatch.undo()
+    assert "has no manifest.json" in B.asset_refusal(str(tmp_path))
+    (tmp_path / "manifest.json").write_text("{}")
+    assert B.asset_refusal(str(tmp_path)) is None
+
+
+def test_the_brain_names_every_input_it_cannot_read_and_says_nothing_when_all_read(monkeypatch):
+    """The page's text lists the asset and each pack of both menus with its reason, under the title and the sentence that
+    the packs are being re-walked for dmipy 3.0; with every input readable there is no text."""
+    cfg = BRAIN_TOML
+    uris = [p["uri"] for t in ("wm", "gm") for p in cfg["packs"][t]]
+    monkeypatch.setattr(B, "asset_refusal", lambda w: "the Hub has no manifest.json of the asset")
+    monkeypatch.setattr(B, "pack_refusal", lambda u: None if u == uris[1] else f"reason of {u.rsplit('/', 1)[-1]}")
+    title, text = B.Brain.unavailable(cfg)
+    assert title == cfg["describe"]["title"] and text.startswith(f"# {title}")
+    assert "its packs are being re-walked for dmipy 3.0" in text and "the Hub has no manifest.json of the asset" in text
+    for u in uris:
+        assert (u in text) == (u != uris[1]), u
+    assert text.count("\n- ") == 4
+    monkeypatch.setattr(B, "asset_refusal", lambda w: None)
+    monkeypatch.setattr(B, "pack_refusal", lambda u: None)
+    assert B.Brain.unavailable(cfg) is None
+
+
+def test_the_brain_page_is_the_text_instead_of_a_crash_and_the_disco_page_is_unaffected(monkeypatch):
+    """With the Hub serving none of the brain's inputs the Space builds a page of the refusal -- the title, the text, nothing
+    to run, and without reading the asset the page would be built from -- while the DiSCo page is built from its own
+    class and never consults the brain."""
+    gr = pytest.importorskip("gradio")
+    from space import app as A
+    from space.sources import disco as D
+    cfg = BRAIN_TOML
+    monkeypatch.setattr(B, "asset_refusal", lambda w: "the Hub has no manifest.json of the asset")
+    monkeypatch.setattr(B, "pack_refusal", lambda u: "it is not in its dataset on the Hub")
+    monkeypatch.setattr(B.Brain, "describe", classmethod(lambda cls, c: (_ for _ in ()).throw(AssertionError("the asset was read"))))
+    demo = A.build(cfg=cfg)
+    texts = [b.value for b in demo.blocks.values() if isinstance(b, gr.Markdown)]
+    assert demo.title == cfg["describe"]["title"] and len(texts) == 1 and "its packs are being re-walked for dmipy 3.0" in texts[0]
+    assert not any(isinstance(b, gr.Button) for b in demo.blocks.values())
+    disco = P.config()
+    assert D.Layout.unavailable(disco) is None and D.Columns.unavailable(disco) is None
+    monkeypatch.setattr(B.Brain, "unavailable", classmethod(lambda cls, c: (_ for _ in ()).throw(AssertionError("the brain was consulted"))))
+    page = A.build(cfg=disco)
+    assert page.title == D.Disco.describe(disco)["title"] and any(isinstance(b, gr.Button) for b in page.blocks.values())

tests/test_deploy_requirements.py

@@ -0,0 +1,40 @@
+"""The ZeroGPU requirements of a configuration are the default's, with the pins of its own file when it has one: one
+family of requirements files."""
+import importlib.util
+import os
+
+ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
+
+
+def _deploy():
+    spec = importlib.util.spec_from_file_location("deploy_tool", os.path.join(ROOT, "tools", "deploy.py"))
+    mod = importlib.util.module_from_spec(spec)
+    spec.loader.exec_module(mod)
+    return mod
+
+
+def test_the_default_is_requirements_zero_as_it_stands():
+    d = _deploy()
+    with open(os.path.join(ROOT, "requirements-zero.txt")) as f:
+        assert d.compose_requirements(ROOT, d.DEFAULT_CONFIG) == f.read()
+
+
+def test_the_brain_space_takes_the_default_pins():
+    """One code runs both pages, so one release of the dmipy packages serves both: the brain Space has no pins of its own."""
+    d = _deploy()
+    assert not os.path.exists(os.path.join(ROOT, d.requirements_for("brain.toml")))
+    assert d.compose_requirements(ROOT, "brain.toml") == d.compose_requirements(ROOT, d.DEFAULT_CONFIG)
+
+
+def test_a_configuration_with_pins_of_its_own_replaces_only_those_lines(tmp_path):
+    d = _deploy()
+    (tmp_path / "requirements-zero.txt").write_text("numpy>=2\n# the stack\njax==1\ndmipy-sim==3.0.0\n")
+    (tmp_path / "requirements-x-zero.txt").write_text("# x was accepted on this one\ndmipy-sim==3.0.1\n")
+    assert d.compose_requirements(str(tmp_path), "x.toml") == "numpy>=2\n# the stack\njax==1\n# x was accepted on this one\ndmipy-sim==3.0.1\n"
+    (tmp_path / "requirements-y-zero.txt").write_text("not-listed==1\n")
+    try:
+        d.compose_requirements(str(tmp_path), "y.toml")
+    except SystemExit as e:
+        assert "not-listed" in str(e)
+    else:
+        raise AssertionError("a pin of a package the default does not list is refused")

tests/test_dmipy_names.py

@@ -0,0 +1,45 @@
+"""Every dmipy name this repository imports resolves in the installed release. A static pass over the imports of the
+package, the tools and the tests -- nothing of them is run -- because most sit inside functions that only a deployment
+calls: a name that moved between releases fails here, in the gate, not on the Space."""
+import ast
+import importlib
+import os
+
+ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
+
+
+def _sources():
+    for top in ("space", "tools", "tests", "app.py"):
+        path = os.path.join(ROOT, top)
+        if os.path.isfile(path):
+            yield path
+            continue
+        for base, dirs, files in os.walk(path):
+            dirs[:] = [d for d in dirs if d != "__pycache__"]
+            yield from (os.path.join(base, f) for f in sorted(files) if f.endswith(".py"))
+
+
+def _imports():
+    """``(file, line, module, name)`` of every ``import dmipy_x.y`` (``name`` None) and ``from dmipy_x.y import name``."""
+    for path in _sources():
+        with open(path, encoding="utf-8") as f:
+            tree = ast.parse(f.read(), path)
+        for node in ast.walk(tree):
+            if isinstance(node, ast.Import):
+                yield from ((path, node.lineno, a.name, None) for a in node.names if a.name.startswith("dmipy"))
+            elif isinstance(node, ast.ImportFrom) and node.module and node.module.startswith("dmipy") and not node.level:
+                yield from ((path, node.lineno, node.module, a.name) for a in node.names)
+
+
+def test_every_dmipy_name_the_repository_imports_resolves():
+    imports = list(_imports())
+    assert imports                                             # the pass reads the package, the tools and the tests
+    unresolved = []
+    for path, line, module, name in imports:
+        try:
+            mod = importlib.import_module(module)
+            if name is not None and not hasattr(mod, name):
+                importlib.import_module(f"{module}.{name}")    # a submodule: `from dmipy_replay.replay import so3`
+        except Exception as e:                                 # an ImportError names the moved name; anything else is a broken module
+            unresolved.append(f"{os.path.relpath(path, ROOT)}:{line}: {module}{'' if name is None else ' import ' + name}: {e}")
+    assert not unresolved, "\n".join(unresolved)

tests/test_gate_job.py

@@ -0,0 +1,23 @@
+"""The acceptance job's script runs the gate: it sets every variable the acceptance tests are gated on, so the paid job
+never ends with every test skipped."""
+import importlib.util
+import os
+import re
+
+ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
+
+
+def _gate_job():
+    spec = importlib.util.spec_from_file_location("gate_job_tool", os.path.join(ROOT, "tools", "gate_job.py"))
+    mod = importlib.util.module_from_spec(spec)
+    spec.loader.exec_module(mod)
+    return mod
+
+
+def test_the_job_script_sets_every_variable_the_acceptance_tests_are_gated_on():
+    with open(os.path.join(ROOT, "tests", "test_acceptance.py")) as f:
+        gates = re.findall(r'skipif\(\s*not os\.environ\.get\("(\w+)"\)', f.read())
+    assert gates, "the acceptance tests are gated on a variable"
+    script = _gate_job().SCRIPT
+    for name in gates:
+        assert re.search(rf"(^|\s){name}=1(\s|$)", script), f"the job script does not set {name}=1"

tests/test_pipeline.py

@@ -13,12 +13,50 @@ def test_the_config_names_every_shell_and_preset_shape():
     shapes = set(CFG["shapes"])
     assert {s["shape"] for s in CFG["disco"]["shells"]} <= shapes
     for name, preset in CFG["presets"].items():
-        assert {s["shape"] for s in preset["shells"]} <= shapes, name
+        assert preset["shape"] in shapes, name
         assert P.preset_protocol(CFG, name).name == name
     for s in CFG["shapes"].values():
         assert 0 < s["delta"] < s["Delta"] < s["TE"]
 
 
+#: the paper's flagship acquisitions (experiments/disco/disco_acq.py): b = 0 rows, shells (b, directions), class
+FLAGSHIP = {"clinical b1000 x 30": (3, [(1000, 30)], "d12-D24"), "research 3-shell x 90": (9, [(1000, 90), (2000, 90), (3000, 90)], "d17-D30"),
+            "Connectome 2.0 b6000 x 60": (6, [(6000, 60)], "d8-D20"), "cosine OGSE 50 Hz b1000 x 60": (6, [(1000, 60)], "cos-f50-T20")}
+
+
+def test_the_presets_are_the_papers_flagship_tables():
+    """Every preset is its stored gradient table, the paper's own (data/protocols/, the same files the paper replays):
+    its b = 0 rows first, then each shell's directions, every row on the preset's class."""
+    assert set(FLAGSHIP) <= set(CFG["presets"])
+    for name, (n_b0, shells, shape) in FLAGSHIP.items():
+        p = P.preset_protocol(CFG, name)
+        assert p.n_b0 == n_b0 and [(s.b, s.n_dirs) for s in p.shells] == shells and {s.shape for s in p.shells} == {shape}
+        m = P.measurements(p, CFG["shapes"])
+        np.testing.assert_allclose(np.linalg.norm(m.dirs, axis=1), 1.0, atol=1e-6)
+
+
+def test_each_class_plays_its_one_waveform():
+    """A class's rows are its one shape scaled to each row's b along each row's direction: the b of every row is its
+    own; a square PGSE class needs the textbook amplitude; the OGSE class keeps the 200 T/m/s edges of its b = 1000
+    build at every b (the shape is fixed), a b = 0 row plays nothing."""
+    from space import classes as C
+    u = np.array([[1.0, 0.0, 0.0], [0.0, 0.6, 0.8], [0.0, 0.0, 1.0]])
+    for name, sh in CFG["shapes"].items():
+        if C.kind(sh) == "tensor":                                     # its rows are level matrices: tests/test_tensor.py
+            continue
+        b = np.array([500.0, 1000.0, 0.0]) if C.kind(sh) == "ogse" else np.array([1000.0, 3000.0, 0.0])
+        seq = C.play(sh, b * 1e6, u)
+        np.testing.assert_allclose(seq.b(), b * 1e6, rtol=1e-5, atol=1e3)
+        G = np.asarray(seq.G, np.float64)
+        assert not G[2].any()
+        if C.kind(sh) == "pgse":
+            assert abs(C.amplitude(sh, 1000) / P.gradient_needed(1000, sh["delta"], sh["Delta"]) - 1) < 0.02, name
+    sh = CFG["shapes"]["cos-f50-T20"]
+    G = np.abs(np.asarray(C.play(sh, [0.5e9, 1e9], u[:2]).G, np.float64)).max(axis=2)
+    np.testing.assert_allclose(G[0] / G[0].max(), G[1] / G[1].max(), atol=1e-6)        # one shape at both b
+    assert 0.2 < C.amplitude(sh, 1000) < 0.3                                             # 0.25 T/m: a 50 Hz cosine is weak in b
+
+
 def test_a_protocols_measurements_are_its_rows_in_order():
     p = P.Protocol((P.Shell("d12-D24", 1000, 30), P.Shell("d8-D20", 6000, 12)), n_b0=2)
     m = P.measurements(p, CFG["shapes"])
@@ -122,9 +160,9 @@ def test_the_volumes_round_trip(tmp_path):
 
 def test_the_backend_and_residency_come_from_the_config_or_the_environment(monkeypatch):
     monkeypatch.delenv("DISCO_BACKEND", raising=False); monkeypatch.delenv("DISCO_RESIDENT", raising=False)
-    assert P.backend(CFG) == CFG["compute"]["backend"] == "jax" and P.resident(CFG) is True
-    monkeypatch.setenv("DISCO_BACKEND", "torch"); monkeypatch.setenv("DISCO_RESIDENT", "0")
-    assert P.backend(CFG) == "torch" and P.resident(CFG) is False
+    assert P.backend(CFG) == CFG["compute"]["backend"] == "jax" and P.resident(CFG) is False
+    monkeypatch.setenv("DISCO_BACKEND", "torch"); monkeypatch.setenv("DISCO_RESIDENT", "1")
+    assert P.backend(CFG) == "torch" and P.resident(CFG) is True
     monkeypatch.setenv("DISCO_BACKEND", "numpy")
     with pytest.raises(ValueError, match="backend"):
         P.backend(CFG)
@@ -268,28 +306,30 @@ def test_floor_stats_are_over_the_voxels_with_signal():
 
 
 def test_the_layout_at_the_pinned_revision_holds_the_configs_classes_and_tiers():
-    """The manifest the Space will download: every [shapes] class with the config's timing, and the tiers when the
-    tissue panel starts on (needs the Hub; skipped offline)."""
+    """The manifest the Space will download: every [shapes] class is a stored partition whose terms are that class's
+    waveform (its template fits them), and the tiers when the tissue panel starts on (needs the Hub; skipped offline)."""
+    import hashlib
     import json
+    from dmipy_replay.replay.time_segments import Partition, Shape
+    from space import classes as C
     hf = pytest.importorskip("huggingface_hub")
     d = CFG["data"]
     import requests
     from huggingface_hub.errors import HfHubHTTPError, LocalEntryNotFoundError, OfflineModeIsEnabled
     try:
-        path = hf.hf_hub_download(d["repo"], f"{d['moments']}/manifest.json", repo_type="dataset", revision=d["revision"])
+        path = hf.hf_hub_download(d["repo"], f"{d['layout']}/manifest.json", repo_type="dataset", revision=d["revision"])
     except (requests.ConnectionError, HfHubHTTPError, LocalEntryNotFoundError, OfflineModeIsEnabled) as e:   # offline, or the Hub unreachable
         pytest.skip(f"the Hub is not reachable: {e!r}"[:120])
-    with open(path) as f:
-        m = json.load(f)
-    assert m["source"]["manifest_sha256"] == d["source_manifest_sha256"]
-    missing = [c for c in CFG["shapes"] if c not in m["columns"]]
+    raw = open(path, "rb").read(); m = json.loads(raw)
+    if d.get("manifest_sha256"):
+        assert hashlib.sha256(raw).hexdigest() == d["manifest_sha256"]
+    missing = [c for c in CFG["shapes"] if c not in m["partitions"]]
     assert not missing, f"the layout at {d['revision'][:12]} lacks {missing}"
-    for name, t in CFG["shapes"].items():
-        enc = m["shapes"][name].get("encoding") or {}
-        for k in ("delta", "Delta", "TE"):
-            v = enc.get(k); v = v[0] if isinstance(v, list) else v
-            if v is not None:
-                assert abs(v - t[k]) < 1e-9, (name, k, v, t[k])
+    for name, sh in CFG["shapes"].items():
+        p = m["partitions"][name]
+        part = Partition(name, tuple(tuple(t["t"]) for t in p["segments"]), tuple(int(t["gate"]) for t in p["segments"]),
+                         tuple(Shape.from_meta(t["shape"], name) for t in p["segments"]))
+        part.levels(C.template(sh))                                       # refused by name unless the waveform is the class's
     if CFG["physics"]["default_on"]:
         assert m.get("tiers"), "the tissue panel starts on but the layout has no tiers"
 
@@ -328,7 +368,7 @@ def test_the_ingredients_are_the_layouts_tier_maps_in_the_pages_units():
             calls.append((shape, tissue is not None, scanner, backend, resident, tuple(b0_direction)))
             return dict(pool={"intra": np.full((2, 1, 1), 0.4), "extra": np.full((2, 1, 1), 0.6)}, contact=np.array([[[-2e-5]], [[np.nan]]]),
                         contact_weight=np.full((2, 1, 1), 0.9), phase=np.zeros((2, 1, 1)), phase_std=np.full((2, 1, 1), 0.1))
-    lay = D.Layout.__new__(D.Layout); lay.moments = _M(); lay.tiers = True; lay.backend = "torch"; lay.unseeded = ("myelin",)
+    lay = D.Layout.__new__(D.Layout); lay.layout = _M(); lay.tiers = True; lay.backend = "torch"; lay.unseeded = ("myelin",); lay.resident = True
     meas = P.measurements(P.Protocol((P.Shell("d12-D24", 1000, 6),), n_b0=1), P.config()["shapes"])
     ph = P.Physics.at(3.0)
     out = D.Layout.ingredients(lay, meas, ph)
@@ -344,27 +384,37 @@ def test_the_ingredients_are_the_layouts_tier_maps_in_the_pages_units():
 
 
 def test_a_machine_reaches_the_layout_as_a_delivery_at_the_placement():
-    """On a stand-in layout: the ideal scanner asks the image for the commanded acquisition (no delivery); a machine
-    asks the layout's delivery for every class at the machine's catalogue entry and the grid centred at the physics'
-    offset, and hands it to that class's image; a machine with its own gradient on a layout without background moments,
-    and with a transmit law on one without RF schedules, is refused by name."""
-    from dmipy_sim.phantom import Grid
+    """On a stand-in layout: the ideal scanner asks the image for the commanded acquisition (the class's waveform, no
+    delivery); a machine asks the layout's delivery for every class with the grid centred at the physics' offset, and
+    hands it to that class's image; the image's normalised signal is scaled by every voxel's walker weight; a machine
+    with its own gradient on a class without background vectors, and any machine on a shaped class (the OGSE), is
+    refused by name."""
+    from dmipy_replay.phantom import Grid
+    from space import classes as C
     seen = []
+
+    class _Part:
+        def __init__(self, constant):
+            self.is_constant = constant
+
     class _M:
         grid = Grid(shape=(2, 1, 1), voxel_size_m=(25e-6,) * 3, origin_m=(0.0, 0.0, 0.0))
-        background = None
-        manifest = {"shapes": {"d12-D24": {}, "d8-D20": {}}}
-        def delivery(self, shape, b, dirs, scanner, grid):
-            seen.append(("delivery", shape, scanner.name, tuple(np.round(np.asarray(grid.isocenter_m) - np.asarray(self.grid.isocenter_m), 6))))
-            return ("delivered", shape)
-        def image(self, shape, b, dirs, *, backend, resident, tissue, scanner, b0_direction, delivered=None):
-            seen.append(("image", shape, delivered, scanner, tuple(b0_direction)))
-            return np.ones((2, 1, 1, len(b))), np.zeros((2, 1, 1))
-    lay = D.Layout.__new__(D.Layout); lay.moments = _M(); lay.tiers = True; lay.backend = "jax"; lay.unseeded = ("myelin",)
-    lay.M0 = np.ones((2, 1, 1)); lay.mask = np.ones((2, 1, 1), bool); lay.pools = ("intra", "extra"); lay.shapes = CFG["shapes"]
+        partitions = {n: _Part(C.kind(sh) != "ogse") for n, sh in CFG["shapes"].items()}
+        manifest = {"partitions": {n: {"group": "g0"} for n in CFG["shapes"]}, "groups": {"g0": {}}}
+
+        def delivery(self, seq, scanner, grid, *, partition):
+            seen.append(("delivery", partition, scanner.name, tuple(np.round(np.asarray(grid.isocenter_m) - np.asarray(self.grid.isocenter_m), 6))))
+            return ("delivered", partition)
+
+        def image(self, seq, *, partition, backend, resident, tissue, scanner, b0_direction, delivered=None):
+            seen.append(("image", partition, delivered, scanner, tuple(b0_direction), np.asarray(seq.G).shape[0]))
+            return np.full((2, 1, 1, np.asarray(seq.G).shape[0]), 0.5), np.zeros((2, 1, 1))
+    lay = D.Layout.__new__(D.Layout); lay.layout = _M(); lay.tiers = True; lay.backend = "jax"; lay.unseeded = ("myelin",); lay.resident = False
+    lay.M0 = np.array([2.0, 4.0]).reshape(2, 1, 1); lay.mask = np.ones((2, 1, 1), bool); lay.pools = ("intra", "extra"); lay.shapes = CFG["shapes"]
     meas = P.measurements(P.Protocol((P.Shell("d12-D24", 1000, 3), P.Shell("d8-D20", 2000, 3)), n_b0=1), CFG["shapes"])
-    D.Layout.replay(lay, meas, P.Physics.at(3.0))
-    assert [s[2] for s in seen if s[0] == "image"] == [None, None]
+    dwi, floor, factor, _ = D.Layout.replay(lay, meas, P.Physics.at(3.0))
+    assert [s[2] for s in seen if s[0] == "image"] == [None, None] and [s[5] for s in seen if s[0] == "image"] == [4, 3]
+    np.testing.assert_allclose(dwi, 1.0); np.testing.assert_allclose(factor.ravel(), [0.25, 0.5])   # S0 = 0.5 W over the fullest voxel's W
     seen.clear()
     ph = P.Physics.at(3.0, scanner="siemens_magnetom_prisma_3T", offset_m=(0.0, 0.0, 0.05))
     D.Layout.validate(lay, P.Protocol((P.Shell("d12-D24", 1000, 3),), n_b0=1), ph)
@@ -372,9 +422,45 @@ def test_a_machine_reaches_the_layout_as_a_delivery_at_the_placement():
     assert [s for s in seen if s[0] == "delivery"] == [("delivery", "d12-D24", "siemens_magnetom_prisma_3T", (0.0, 0.0, -0.05)),
                                                         ("delivery", "d8-D20", "siemens_magnetom_prisma_3T", (0.0, 0.0, -0.05))]
     assert [s[2] for s in seen if s[0] == "image"] == [("delivered", "d12-D24"), ("delivered", "d8-D20")]
+    with pytest.raises(ValueError, match="shaped waveforms"):
+        D.Layout.validate(lay, P.Protocol((P.Shell("cos-f50-T20", 1000, 3),), n_b0=1), ph)
     swoop = P.Physics.at(0.064, scanner="hyperfine_swoop_64mT", b0_direction=(0.0, 1.0, 0.0))
-    with pytest.raises(ValueError, match="background moments"):
-        D.Layout.validate(lay, P.Protocol((P.Shell("d12-D24", 100, 3),), n_b0=1), swoop)
-    lay.moments.background = {"amplitude": 5e-3}
-    with pytest.raises(ValueError, match="RF schedules"):
+    with pytest.raises(ValueError, match="background vectors"):
         D.Layout.validate(lay, P.Protocol((P.Shell("d12-D24", 100, 3),), n_b0=1), swoop)
+    lay.layout.manifest["groups"]["g0"]["background"] = {"column": "bg_g0.npy"}
+    D.Layout.validate(lay, P.Protocol((P.Shell("d12-D24", 100, 3),), n_b0=1), swoop)
+
+
+def test_each_preset_starts_at_its_own_snr():
+    """The b-tensor preset starts at SNR 100 (µFA is noise-limited); every other preset at ``[noise] snr``."""
+    for name, p in CFG["presets"].items():
+        assert P.preset_snr(CFG, name) == (100.0 if "tensor" in p else float(CFG["noise"]["snr"])) and CFG["noise"]["snr"] == 30
+    assert P.preset_snr(CFG, "DiSCo 364") == 30.0
+
+
+def test_a_tractogram_writes_to_tck_and_reads_back(tmp_path):
+    """The .tck download: dmipy-tract's writer read back through dmipy-sim's reader (their conventions agree)."""
+    from dmipy_sim.io.strands import read_tck
+    from dmipy_tract import Tractogram
+    t = Tractogram(np.random.default_rng(0).random((10, 3)).astype(np.float32), np.array([0, 4, 10]), np.array([0, 1]), np.zeros((2, 2), int))
+    t.to_tck(str(tmp_path / "t.tck"))
+    assert [len(s) for s in read_tck(str(tmp_path / "t.tck"), coordinate_unit_m=1e-3)] == [4, 6]
+
+
+def test_the_csd_runs_on_the_configured_backend_alone(monkeypatch):
+    """The FOD field (the response estimate and the fit) runs on the backend it is given: on ``jax`` nothing reaches
+    for PyTorch, which the Docker Space does not install."""
+    import sys
+    monkeypatch.setitem(sys.modules, "torch", None)            # any import of torch now fails
+    p = P.Protocol((P.Shell("d12-D24", 1000, 30),))
+    m = P.measurements(p, CFG["shapes"])
+    rng = np.random.default_rng(3)
+    axes = rng.normal(size=(2, 2, 2, 3))
+    axes /= np.linalg.norm(axes, axis=-1, keepdims=True)
+    g = m.dirs / np.maximum(np.linalg.norm(m.dirs, axis=-1, keepdims=True), 1e-12)
+    cos2 = np.einsum("xyzk,nk->xyzn", axes, g) ** 2
+    b = m.bvals * 1e-3                                         # ms/um^2 x um^2/ms
+    dwi = np.exp(-b * (0.3 + 1.4 * cos2))                     # a prolate tensor along each voxel's axis
+    mask = np.ones(dwi.shape[:3], bool)
+    sh, seconds, _ = P.csd(dwi, m, mask, backend="jax")
+    assert sh.shape == (2, 2, 2, 45) and np.isfinite(sh).all() and seconds > 0

tests/test_tensor.py

@@ -0,0 +1,106 @@
+"""The b-tensor class's encodings (space/tensor.py) without the layout: every row's b-tensor, read off the sequence by
+dmipy-replay, is the encoding asked for; the class's waveform stays inside the scanner it was designed for; the maps
+recover a known microscopic anisotropy from Gaussian compartments played through the same rows."""
+import numpy as np
+import pytest
+
+from space import pipeline as P
+from space import tensor as T
+
+CFG = P.config()
+TC = next(sh for sh in CFG["shapes"].values() if sh.get("kind") == "tensor")
+
+
+def _rows(bvals=(500.0, 1000.0, 2000.0), n=12):
+    return T.rows(TC, list(bvals), n)
+
+
+def test_every_row_is_the_encoding_asked_for():
+    """b to 1e-6 relative and b_delta to 1e-6 (1, -1/2, 0) for every row; the b = 0 row plays nothing."""
+    from dmipy_replay.acquisition.waveforms import btensor_invariants
+    seq, enc, bs = _rows()
+    B = seq.btensor()
+    b, b_delta, _ = btensor_invariants(B)
+    assert enc[0] == "b0" and abs(b[0]) < 1e-6
+    for e in T.ENCODINGS:
+        m = enc == e
+        np.testing.assert_allclose(b[m] / 1e6, bs[m], rtol=1e-6)
+        np.testing.assert_allclose(b_delta[m], T.B_DELTA[e], atol=1e-6)
+    ste = B[enc == "STE"]                                            # spherical: b/3 I at every rotation
+    np.testing.assert_allclose(ste, (bs[enc == "STE"] * 1e6 / 3)[:, None, None] * np.eye(3), rtol=0, atol=1e-6 * 2e9)
+
+
+def test_the_waveform_is_refocused_and_inside_its_scanner():
+    """Each shape's net moment at the echo vanishes and the designed amplitude and slew are within 80 mT/m and
+    200 T/m/s per axis; the spherical row at b = 2000 at the design's own frame needs no more than the design."""
+    W = T.waveform(TC); g = T._grid(); e = np.asarray(g.effective_gate, np.float64)
+    net = np.abs((W * e[:, None]).sum(0) * g.dt)
+    assert (net < 1e-9 * np.abs(W).max() * 0.0535).all(), net
+    assert np.abs(W).max() <= 0.080 * (1 + 1e-6)
+    assert np.abs(np.diff(W, axis=0) / g.dt).max() <= 200.0 * 3 * (1 + 1e-6)   # linear interpolation of a 200 T/m/s design
+
+
+@pytest.mark.parametrize("bvals", [(500.0, 1000.0, 2000.0), (100.0, 700.0, 1400.0, 2000.0)])
+def test_the_maps_recover_a_known_microscopic_anisotropy(bvals):
+    """Powder of randomly oriented axisymmetric Gaussian tensors (D_iso 1 um^2/ms, muFA 0.8) played through the rows'
+    own b-tensors: MD within 3 %, muFA within 0.03 over the b-values (the cumulant fits' own truncation); the powder
+    averages are read by the b-values as given (700 is not 0.7 x 1e3 in floating point)."""
+    seq, enc, bs = _rows(bvals=bvals, n=30)
+    B = seq.btensor() * 1e-9                                           # s/m^2 -> ms/um^2
+    D_iso, mufa = 1.0, 0.8
+    x = mufa ** 2 * D_iso ** 2 / (3 - 2 * mufa ** 2)                  # (D_iso D_Delta)^2 of an axisymmetric tensor of that muFA
+    D_delta = np.sqrt(x) / D_iso
+    rng = np.random.default_rng(0)
+    u = rng.normal(size=(4000, 3)); u /= np.linalg.norm(u, axis=1)[:, None]
+    Dpar, Dperp = D_iso * (1 + 2 * D_delta), D_iso * (1 - D_delta)
+    D = Dperp * np.eye(3)[None] + (Dpar - Dperp) * u[:, :, None] * u[:, None, :]
+    S = np.exp(-np.einsum("mij,cij->mc", B, D)).mean(1)
+    S = S / S[0]
+    pw = T.powder(S[None, None, None, :], enc, bs, list(bvals))
+    md, got = T.md_ufa(pw, list(bvals))
+    assert abs(float(md.ravel()[0]) - D_iso) < 0.03 * D_iso
+    assert abs(float(got.ravel()[0]) - mufa) < 0.03
+
+
+def test_the_tensor_preset_is_linear_and_spherical_rows_on_the_class_and_flags_a_prisma():
+    """The preset: one b = 0 row, then LTE and STE interleaved at b 100, 700, 1400, 2000 with 6, 6, 10, 16 axes, every
+    row's b-tensor the one asked for as dmipy-replay reads it off the sequence the Layout plays; rotated, the rows need
+    more than a Prisma's 80 mT/m on one axis, so the page flags that machine and passes the ideal scanner."""
+    from dmipy_replay.acquisition.waveforms import btensor_invariants
+    from space import classes as C
+    name = next(n for n in CFG["presets"] if "tensor" in CFG["presets"][n])
+    p = P.preset_protocol(CFG, name)
+    assert p.tensor and p.n_b0 == 1 and [(s.encoding, s.b, s.n_dirs) for s in p.shells] == [
+        (e, b, n) for b, n in ((100.0, 6), (700.0, 6), (1400.0, 10), (2000.0, 16)) for e in ("LTE", "STE")]
+    m = P.measurements(p, CFG["shapes"])
+    G = C.tensor_rows(TC, m.encoding, m.bvals, m.dirs)
+    b, b_delta, _ = btensor_invariants(T._sequence(G).btensor())
+    np.testing.assert_allclose(b[1:] / 1e6, m.bvals[1:], rtol=1e-6)
+    np.testing.assert_allclose(b_delta[m.encoding == "LTE"], 1.0, atol=1e-6); np.testing.assert_allclose(b_delta[m.encoding == "STE"], 0.0, atol=1e-6)
+    rows = P.playable(p, CFG["shapes"], "siemens_magnetom_prisma_3T")
+    assert not all(r[-1] for r in rows) and all(r[-1] for r in P.playable(p, CFG["shapes"], None))
+    assert 0.10 < max(r[4] for r in rows) < 0.14                                   # ~124 mT/m at b 2000
+
+
+def test_the_maps_hold_under_noise_on_one_b0_row():
+    """The same voxel as above, 3000 times, with Rician noise at SNR 15 and ONE b = 0 row, every row normalised by it
+    as the page does: the joint fit's shared intercept takes up the b = 0 row's noise, so the median muFA stays within
+    0.08 of the truth and MD within 5 % (a fit through ln S/S0 = 0 at b = 0 reads that noise as curvature)."""
+    bvals = (100.0, 700.0, 1400.0, 2000.0)
+    seq, enc, bs = _rows(bvals=bvals, n=10)
+    B = seq.btensor() * 1e-9
+    D_iso, mufa = 0.8, 0.7
+    x = mufa ** 2 * D_iso ** 2 / (3 - 2 * mufa ** 2); D_delta = np.sqrt(x) / D_iso
+    rng = np.random.default_rng(1)
+    u = rng.normal(size=(4000, 3)); u /= np.linalg.norm(u, axis=1)[:, None]
+    Dpar, Dperp = D_iso * (1 + 2 * D_delta), D_iso * (1 - D_delta)
+    D = Dperp * np.eye(3)[None] + (Dpar - Dperp) * u[:, :, None] * u[:, None, :]
+    S = np.exp(-np.einsum("mij,cij->mc", B, D)).mean(1)
+    sd = 1.0 / 15.0
+    noisy = np.sqrt((S[None] + rng.normal(0, sd, (3000, len(S)))) ** 2 + rng.normal(0, sd, (3000, len(S))) ** 2)
+    noisy = noisy / noisy[:, :1]                                       # the one b = 0 row
+    pw = T.powder(noisy[:, None, None, :], enc, bs, list(bvals))
+    counts = {(e, b): int(((enc == e) & (bs == b)).sum()) for e in ("LTE", "STE") for b in bvals}
+    md, got = T.md_ufa(pw, list(bvals), counts, n_b0=1)
+    assert abs(float(np.nanmedian(got)) - mufa) < 0.08, float(np.nanmedian(got))
+    assert abs(float(np.nanmedian(md)) - D_iso) < 0.05 * D_iso

tests/test_viewers.py

@@ -79,3 +79,15 @@ def test_the_map_slice_draws_plain_and_diverging():
     vol = np.random.default_rng(0).normal(size=(6, 6, 3)); vol[0, 0] = np.nan
     assert V.map_slice(vol, 1, "t").size[0] > 100 and V.map_slice(vol, 1, "d", symmetric=True).size[0] > 100
     assert V.map_slice(np.full((6, 6, 3), np.nan), 1, "empty", symmetric=True).size[0] > 100
+
+
+def test_the_tensor_views_draw_maps_with_the_crossings_and_a_voxels_powder_curve():
+    """The b-tensor results: FA, muFA (crossings outlined) and MD on a slice; a voxel's powder signal against b for
+    A and, dashed, B."""
+    rng = np.random.default_rng(1)
+    bv = [100.0, 700.0, 1400.0, 2000.0]
+    pw = {(e, b): np.exp(-b * 1e-3 * (0.6 if e == "STE" else 0.5)) * np.ones((4, 4, 4)) for e in ("LTE", "STE") for b in bv}
+    m = dict(fa=rng.random((4, 4, 4)), ufa=rng.random((4, 4, 4)), md=rng.random((4, 4, 4)), powder=pw, bvals=bv)
+    cross = np.zeros((4, 4, 4), bool); cross[1:3, 1:3] = True
+    assert V.tensor_maps(m, cross, 2).size[0] > 0
+    assert V.powder_curve(m, (1, 1, 2), other=m).size[0] > 0

tools/build_brain_asset.py

@@ -94,8 +94,7 @@ def fit_msmt_csd(data, bvals, bvecs, mask, sh_order=SH_ORDER, solver="csd_cvxpy"
     those voxels are fit; the rest of ``mask`` is still used to estimate the responses.
 
     Always fits serially (``use_parallel_processing=False``): dmipy-fit's parallel path forks a
-    ``ProcessPoolExecutor`` after the JAX import this fit chain pulls in (``dmipy_sim.replay.so3``, for the SH
-    basis), which JAX warns is unsafe. Thread-level parallelism (BLAS/OSQP) is capped by
+    ``ProcessPoolExecutor``, which JAX warns is unsafe once JAX is imported. Thread-level parallelism (BLAS/OSQP) is capped by
     ``OMP_NUM_THREADS``/``MKL_NUM_THREADS`` in the environment, per this lab's shared-box rule.
 
     Returns ``(sh, fractions, S0_responses, response_models, scheme, seconds)``: ``sh`` is ``(X, Y, Z, n_coef)`` (the

tools/build_brain_fixture.py

@@ -5,9 +5,9 @@ before the MASiVar asset exists. It stays local (BATMAN's licence, dmipy-sim#193
 
 What each file is:
 
-* ``fod_wm.npy`` -- ``wmfod_norm.mif`` (MRtrix3 multi-tissue CSD, lmax 8, the MRtrix3 basis, which is dmipy-sim's
+* ``fod_wm.npy`` -- ``wmfod_norm.mif`` (MRtrix3 multi-tissue CSD, lmax 8, the MRtrix3 basis, which is dmipy-replay's
   required basis and dmipy-fit's) turned from the scanner frame into the image frame by ``R.T`` (``R`` from
-  :meth:`dmipy_sim.phantom.Grid.from_oblique_affine`), zero outside the brain mask.
+  :meth:`dmipy_replay.phantom.Grid.from_oblique_affine`), zero outside the brain mask.
 * ``fractions.npy`` -- WM / GM / CSF from ``5tt_coreg.mif`` (on the 1 mm T1 grid), averaged into every 2.5 mm
   voxel on a 3^3 sub-grid as ``examples/rph/brain_from_csd.py`` does: GM is cortical + sub-cortical GM, WM where
   the FOD has no positive l = 0 coefficient is counted as GM (no FOD, no oriented WM), zero outside the mask.
@@ -19,7 +19,7 @@ What each file is:
 * ``manifest.json`` -- the grid (the image's affine, voxel -> scanner RAS mm), the tutorial's gradient table
   (scanner frame) with the capstone's pulse timing (TE 100 / delta 25 / Delta 55 ms, which the 100 ms packs play).
 
-    python tools/build_brain_fixture.py [--batman ~/dmrai-ws/data/batman] [--out ~/dmrai-ws/data/batman/brain_fixture]
+    python tools/build_brain_fixture.py --batman <dir> --out <dir>        # or DISCO_BATMAN / DISCO_BRAIN_FIXTURE
 """
 import argparse
 import json
@@ -76,13 +76,15 @@ def parcellate(gm, shape):
 
 
 def main(argv=None):
+    from dmipy_replay.phantom import Grid
+    from dmipy_replay.replay.so3 import rotate_sh
     from dmipy_sim.io.mrtrix import read_mif
-    from dmipy_sim.phantom import Grid
-    from dmipy_sim.replay.so3 import rotate_sh
     ap = argparse.ArgumentParser()
-    ap.add_argument("--batman", default=os.path.expanduser("~/dmrai-ws/data/batman"))
-    ap.add_argument("--out", default=os.path.expanduser("~/dmrai-ws/data/batman/brain_fixture"))
+    ap.add_argument("--batman", default=os.environ.get("DISCO_BATMAN"), help="the BATMAN directory (DISCO_BATMAN)")
+    ap.add_argument("--out", default=os.environ.get("DISCO_BRAIN_FIXTURE"), help="the fixture directory (DISCO_BRAIN_FIXTURE)")
     a = ap.parse_args(argv)
+    if not a.batman or not a.out:
+        ap.error("name the BATMAN directory and the output directory (--batman/--out or DISCO_BATMAN/DISCO_BRAIN_FIXTURE)")
     os.makedirs(a.out, exist_ok=True)
     fod = read_mif(os.path.join(a.batman, "wmfod_norm.mif"))
     mask = np.asarray(read_mif(os.path.join(a.batman, "mask_den_unr_preproc_unb.mif")).data, bool)

tools/deploy.py

@@ -15,6 +15,7 @@ bucket is retired), so the deploy also removes a ``DISCO_MOMENTS`` variable and
 """
 import argparse
 import os
+import re
 import shutil
 import subprocess
 import tempfile
@@ -44,6 +45,51 @@ def readme_for(config):
     return "README-zero.md" if config == DEFAULT_CONFIG else f"README-{os.path.splitext(config)[0]}-zero.md"
 
 
+def requirements_for(config):
+    """The file of the pins a configuration replaces in ``requirements-zero.txt`` (:func:`compose_requirements`):
+    ``requirements-<stem>-zero.txt``, when the configuration has one."""
+    return f"requirements-{os.path.splitext(config)[0]}-zero.txt"
+
+
+def _package_of(line):
+    """The distribution a requirements line pins (``name @ url``, ``name==x``, ``name>=x``), or None for a comment."""
+    line = line.strip()
+    if not line or line.startswith("#"):
+        return None
+    return re.split(r"[\s@=<>!~;\[]", line, maxsplit=1)[0].lower().replace("_", "-")
+
+
+def compose_requirements(root, config):
+    """The ZeroGPU ``requirements.txt`` of a configuration: ``requirements-zero.txt`` and, for a configuration that has a
+    file of pins of its own, the lines of ``requirements-<stem>-zero.txt`` replacing the lines of the same package
+    (every other line is the default's, in one place). A configuration without that file takes the default's."""
+    with open(os.path.join(root, "requirements-zero.txt")) as f:
+        lines = f.read().splitlines()
+    own = os.path.join(root, requirements_for(config))
+    if config == DEFAULT_CONFIG or not os.path.exists(own):
+        return "\n".join(lines) + "\n"
+    with open(own) as f:
+        over = f.read().splitlines()
+    pins = {_package_of(l): l for l in over if _package_of(l)}
+    notes = {_package_of(over[i + 1]): over[i] for i in range(len(over) - 1) if _package_of(over[i + 1]) and over[i].lstrip().startswith("#")}
+    out, used = [], set()
+    for l in lines:
+        k = _package_of(l)
+        if k in pins:
+            used.add(k)
+            if k in notes:
+                out.append(notes[k])
+            out.append(pins[k])
+        elif k is None and l.lstrip().startswith("#") and _package_of(lines[min(lines.index(l) + 1, len(lines) - 1)]) in pins:
+            continue                                        # the default's note on a pin the configuration replaces
+        else:
+            out.append(l)
+    unknown = sorted(set(pins) - used)
+    if unknown:
+        raise SystemExit(f"{requirements_for(config)} pins {unknown}, which requirements-zero.txt does not list")
+    return "\n".join(out) + "\n"
+
+
 def set_config(api, space, config):
     """``DISCO_CONFIG`` on the Space: the configuration's name, absent for the default."""
     if config == DEFAULT_CONFIG:
@@ -79,7 +125,8 @@ def main():
         for name in ("app.py", "pyproject.toml"):
             shutil.copy(os.path.join(root, name), staged)
         shutil.copy(os.path.join(root, readme_for(a.config)), os.path.join(staged, "README.md"))
-        shutil.copy(os.path.join(root, "requirements-zero.txt"), os.path.join(staged, "requirements.txt"))
+        with open(os.path.join(staged, "requirements.txt"), "w") as f:
+            f.write(compose_requirements(root, a.config))
         folder = staged
     elif a.base_tag:
         staged = tempfile.mkdtemp(prefix="disco-space-")
@@ -91,7 +138,7 @@ def main():
             g.write(f.read().replace("BASE_TAG", a.base_tag))
         folder = staged
     info = api.upload_folder(folder_path=folder, repo_id=a.space, repo_type="space", ignore_patterns=IGNORE,
-                             delete_patterns=["Dockerfile", "Dockerfile.base", "Dockerfile.space", "entrypoint.sh", "README-*zero.md", "requirements-zero.txt"] if a.zero
+                             delete_patterns=["Dockerfile", "Dockerfile.base", "Dockerfile.space", "entrypoint.sh", "README-*zero.md", "requirements-*zero.txt"] if a.zero
                              else (["Dockerfile.base", "Dockerfile.space"] if a.base_tag else None),
                              commit_message=f"disco-space {sha[:12]}" + (f" on base {a.base_tag[:12]}" if a.base_tag else "") + (f", {a.config}" if a.config != DEFAULT_CONFIG else ""))
     if staged:

tools/gate_job.py

@@ -4,19 +4,43 @@ nothing but pytest. The job runs ``tests/test_acceptance.py`` and prints ``gate.
 prints its log and exits with the job's outcome. The Space is private, so the container gets the caller's Hub token
 as the job secret ``HF_TOKEN`` (the layout download reads it too).
 
-    python tools/gate_job.py [--space rfick/disco] [--flavor l4x1] [--timeout 45m]
+The job starts only once the Space runs the deployed commit (:func:`wait_for_build`): a deploy that changes a variable
+rebuilds the image, and a job started on an image still building errors with no log.
+
+    python tools/gate_job.py [--space rfick/disco] [--flavor l4x1] [--timeout 45m] [--build-timeout 2400]
 """
 import argparse
 import os
 import sys
 import time
 
+POLL_S = 30
+
+
+def wait_for_build(api, space, rev, timeout_s):
+    """Return once ``space`` runs commit ``rev`` (stage RUNNING with that sha, on two polls ``POLL_S`` apart, so a stage
+    read just before a rebuild starts does not count); exit on an error stage or after ``timeout_s``."""
+    t0, seen = time.time(), 0
+    time.sleep(POLL_S)
+    while True:
+        r = api.get_space_runtime(space)
+        stage, sha = r.stage, (r.raw or {}).get("sha")
+        if str(stage).endswith("ERROR"):
+            sys.exit(f"the Space is in {stage} at {sha}: the image to test did not build or start")
+        seen = seen + 1 if (stage == "RUNNING" and sha == rev) else 0
+        if seen >= 2:
+            print(f"{space} runs {rev[:12]} ({time.time() - t0:.0f} s after the deploy)", flush=True)
+            return
+        if time.time() - t0 > timeout_s:
+            sys.exit(f"the Space is still {stage} at {sha} after {timeout_s} s, not RUNNING at {rev}")
+        time.sleep(POLL_S)
+
 SCRIPT = r'''
 set -e
 cd /app
 pip install -q --retries 10 pytest
 python -c "import jax; print('devices', jax.devices())"
-JAX_PLATFORMS=cuda XLA_PYTHON_CLIENT_PREALLOCATE=false XLA_FLAGS=--xla_gpu_deterministic_ops=true HF_HOME=/tmp/hf JAX_COMPILATION_CACHE_DIR=/tmp/jax-cache \
+DISCO_GATE=1 JAX_PLATFORMS=cuda XLA_PYTHON_CLIENT_PREALLOCATE=false XLA_FLAGS=--xla_gpu_deterministic_ops=true HF_HOME=/tmp/hf DMIPY_REPLAY_CACHE=/tmp/replay JAX_COMPILATION_CACHE_DIR=/tmp/jax-cache \
   GATE_OUT=/tmp/gate.json PYTHONPATH=/app pytest tests/test_acceptance.py -s -q -p no:cacheprovider
 echo "=== gate.json"; cat /tmp/gate.json; echo; echo "=== GATE PASSED"
 '''
@@ -27,11 +51,13 @@ def main():
     ap.add_argument("--space", default=os.environ.get("DISCO_SPACE", "rfick/disco"))
     ap.add_argument("--flavor", default="l4x1")
     ap.add_argument("--timeout", default="45m")
+    ap.add_argument("--build-timeout", type=int, default=2400, help="seconds to wait for the Space to run the deployed commit")
     a = ap.parse_args()
     from huggingface_hub import HfApi, get_token
     token = os.environ.get("HF_TOKEN") or get_token()
     api = HfApi(token=token)
     rev = api.repo_info(a.space, repo_type="space").sha
+    wait_for_build(api, a.space, rev, a.build_timeout)
     job = api.run_job(image=f"hf.co/spaces/{a.space}", command=["bash", "-c", SCRIPT],
                       secrets={"HF_TOKEN": token}, flavor=a.flavor, timeout=a.timeout, name=f"disco-space-gate-{rev[:8]}")
     print(f"job {job.id} on {a.flavor}: the image of {a.space} (repo at {rev})", flush=True)

tools/live.py

@@ -1,4 +1,4 @@
-"""The live Space through its Gradio API: the configuration's first preset at SNR 30 with the tissue panel's defaults,
+"""The live Space through its Gradio API: the configuration's first preset at its own starting SNR with the tissue panel's defaults,
 then a custom two-shell protocol on the configuration's first timing class; prints the headline, the timings and the
 seconds the worker spent on the packs' responses the page had not cached and, with ``--stream``, every stage update as it
 arrives; downloads the figures and files to ``$LIVE_OUT``. ``--config`` is the configuration the Space serves
@@ -39,7 +39,7 @@ def main():
     def inputs(preset, *, knob, n_keys=1, ladder=True):
         shells = [True, shapes[0], 1000, 30, 12.0, 24.0, 53.5, True, shapes[0], 3000, 45, 8.0, 20.0, 53.5,
                   False, shapes[0], 3000, 90, 17.0, 30.0, 53.5, False, shapes[0], 6000, 60, 17.0, 30.0, 53.5]
-        return [preset, 2, True, 30, tc["density"][2], tc["max_angle"][2], tc["step"][2], 0, None, knob, A.IDEAL, n_keys, ladder, *physics, *shells]
+        return [preset, 2, True, P.preset_snr(cfg, preset) if preset in cfg.get("presets", {}) else cfg["noise"]["snr"], tc["density"][2], tc["max_angle"][2], tc["step"][2], 0, None, knob, A.IDEAL, n_keys, ladder, *physics, *shells]
 
     c = Client(a.space, token=get_token(), verbose=False, download_files=os.environ.get("LIVE_OUT", "/tmp/disco-live"), httpx_kwargs={"timeout": 900})
     for preset in a.preset or (S.presets(cfg)[0], A.CUSTOM):
@@ -62,7 +62,7 @@ def main():
         rows = named["timings"]["data"] if isinstance(named["timings"], dict) else named["timings"]
         responses = [r[1] for r in rows if r[0] == A.RESPONSES_ROW]
         print(f"the packs' responses on the worker: {responses[0] if responses else 'none (the source prepares nothing)'} s", flush=True)
-        print("tck", named["tck"], "volumes", named["volumes"], flush=True)
+        print("tck", named.get("tck"), "volumes", named.get("volumes"), flush=True)
 
 
 if __name__ == "__main__":

tools/smoke_tensor.py

@@ -0,0 +1,55 @@
+"""The b-tensor preset through the page's own run button without a browser: A at 3 T with B = A with the field off
+(the knob), the explorer's ladder on; every output of :data:`space.app.OUTPUTS` is produced and the b-tensor ones are
+images; then a voxel click, an STE row of the DWI slice and the explorer's B − A of µFA.
+
+    DISCO_LAYOUT=/local/layout python tools/smoke_tensor.py [--snr 30]
+"""
+import argparse
+import json
+import time
+import types
+
+import numpy as np
+
+
+def main():
+    ap = argparse.ArgumentParser(); ap.add_argument("--snr", type=float, default=30.0); ap.add_argument("--knob", default="field tier → off")
+    a = ap.parse_args()
+    from space import app as A
+    from space import pipeline as P
+    st = A._load()
+    if st["error"]:
+        raise SystemExit(st["error"])
+    src, cfg = st["source"], st["cfg"]; S = type(src)
+    preset = next(n for n, v in cfg["presets"].items() if "tensor" in v)
+    panel = S.panel(cfg)
+    vals = {c.name: c.value for row in panel.rows for c in row}
+    values = [vals[f] for f in panel.fields]
+    shells = [False, "d12-D24", 1000, 30, 10.0, 20.0, 60.0] * A.MAX_SHELLS
+    args = (preset, 1, True, a.snr, 2, 30.0, 0.5, 0, None, a.knob, A.IDEAL, 1, True, *values, *shells)
+    t0 = time.perf_counter(); last = None; texts = []
+    for item in A.run_pipeline(A.compute, *args):
+        last = item
+        if isinstance(item, tuple) and isinstance(item[1], str):
+            texts.append(item[1])
+    out = dict(zip(A.OUTPUTS, last)); secs = time.perf_counter() - t0
+    rs = out["result"]
+    kinds = {k: type(v).__name__ for k, v in out.items()}
+    for k in ("tensor_view", "tensor_pick", "tensor_view_b", "dwi_view", "explore_view", "metric_view"):
+        assert out[k] is not None and hasattr(out[k], "size"), (k, kinds[k])
+    assert out["tract_view"] is None and out["mats"] is None
+    evt = types.SimpleNamespace(index=(20 * A.PICK + 3, 19 * A.PICK + 3))
+    curve = A.pick_voxel(rs, 20, evt)
+    ste = A.redraw_slice(rs, 20, 3, True, "A", "STE")
+    ex = A.redraw_explorer(rs, rs["explorer"]["layers"][-1][0], "B − A", A.METRICS[3], 20, 0, "LTE")
+    tm = rs["A"]["tensor"]; sig = np.isfinite(rs["A"]["dwi"][..., 0])
+    report = dict(seconds=secs, headline=out["headline"][:600], stages=texts[-6:], outputs=kinds,
+                  ufa_median=float(np.nanmedian(tm["ufa"][sig])), fa_median=float(np.nanmedian(tm["fa"][sig])),
+                  md_median=float(np.nanmedian(tm["md"][sig])), voxel_curve=curve is not None, ste_slice=ste[0] is not None,
+                  explorer_b_minus_a_mufa=ex[0] is not None, files=[str(f).rsplit("/", 1)[-1] for f in out["volumes"]],
+                  timings=out["timings"])
+    print(json.dumps(report, indent=1, default=str))
+
+
+if __name__ == "__main__":
+    main()