disco-space 8c9afd1593ec

README.md

@@ -11,12 +11,16 @@ license: mit
 short_description: DiSCo replay to connectome on the shared GPU pool
 ---
 
-# The DiSCo Space on ZeroGPU
+# The DiSCo replay Space (`rfick/disco-replay`)
 
-The same application as `rfick/disco` (dmrai-lab/disco-space) on Hugging Face's shared GPU pool: the compute backend
-is PyTorch (`DISCO_BACKEND=torch`: dmipy-sim's shape-moment image, dmipy-fit's `csd_tournier07_torch`, dmipy-tract's
-torch tracker), the GPU exists for the duration of one run (`spaces.GPU`), and the layout's tiles are rebuilt on the
-device inside every call from the layout downloaded from the Hub when the container started. Deterministic algorithms are on and TF32 off inside the
+The DiSCo page (dmrai-lab/disco-space) on Hugging Face's shared GPU pool, serving the v2 DiSCo walk through its
+time-segment layout (`SubstrateCommons/disco-time-segments`): the compute backend is PyTorch (`DISCO_BACKEND=torch`:
+dmipy-sim's time-segment image, dmipy-fit's `csd_tournier07_torch`, dmipy-tract's torch tracker), the GPU exists for
+the duration of one run (`spaces.GPU`), and the layout's packed words are streamed to the device a chunk at a time
+inside every call (decoded in the kernel) from the layout downloaded from the Hub when the container started. Its
+acquisitions are the paper's DiSCo flagship set: DiSCo's own 364 measurements, clinical, research, Connectome 2.0 and
+a cosine OGSE at 50 Hz, each one stored waveform of the layout. (Formerly `rfick/disco-zero`; the old address
+redirects.) Deterministic algorithms are on and TF32 off inside the
 call. Source and pins: https://github.com/dmrai-lab/disco-space (`requirements-zero.txt`, `app.py`).
 
 ## GPU quota per visitor

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

@@ -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

@@ -17,7 +17,7 @@ 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-sim: main at the merge of #828 (the columnar field window; after #823, the streamed time-segment image)
+dmipy-sim @ git+https://github.com/dmrai-lab/dmipy-sim@8d7b97f982e0c10e744d0065dd2cafbcc1c6d781
 dmipy-fit @ git+https://github.com/dmrai-lab/dmipy-fit@c3dcc546595cbde8e8e7c7bb2ff4e87f099d5dab
 dmipy-tract @ git+https://github.com/dmrai-lab/dmipy-tract@7da22c3a9473c50a08cc7c6246ec16b68a64d877

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
@@ -523,19 +588,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():
@@ -624,7 +733,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])
@@ -657,6 +768,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 +783,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 +819,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 +846,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 +859,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

@@ -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,114 @@
+"""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-sim's own 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):
+    import dmipy_sim as d
+    sh = dict(key)
+    k = kind(sh)
+    if k == "pgse":
+        return d.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 d.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 d.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."""
+    from dmipy_sim.acquisition.waveforms import calc_b
+    t = template(sh)
+    G = np.asarray(t.G, np.float64)[0, :, 2]
+    peak = float(np.abs(G).max())
+    return G / peak, float(calc_b(t)[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_sim.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"

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 = "afe5589d09a85b70119773cf2948dcc001fa8433"   # the dataset revision the layout is read at (with the background moments, disco-space#29)
-source_manifest_sha256 = "90cbd8f238b3946c224b765e9afeb748817f735410cded75e719bf600911d236"   # 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_sim.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 = "8277ff05ea8237a98f1f8bdfe6a86ce37c4fcdc3"
+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-sim's own
+# 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
@@ -88,7 +134,7 @@ pools = ["intra", "extra"]          # the pack's seeded pools (its spec also nam
 # The scanner menu after "ideal": catalogue machines of dmipy-sim (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/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
@@ -1015,13 +1076,18 @@ def limits(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
 
@@ -1120,6 +1186,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 +1200,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_sim.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_sim.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/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-sim'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
@@ -299,7 +303,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 +347,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 +413,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-sim's time-segment layout (``dmipy_sim.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_sim.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 None, 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 +479,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_sim.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 +517,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_sim.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 +556,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,17 +598,17 @@ 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."""
+    :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
 
@@ -534,15 +621,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 +640,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."""
+        """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."""
         import dmipy_sim as d
+        from dmipy_sim.acquisition.scanner_sequence import ScannerSequence
+        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 = d.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``
@@ -583,7 +684,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():
+    import dmipy_sim as d
+    return d.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_sim.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_sim.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-sim's own b-tensor of the sequence, and the peak gradient per axis (T/m) of every row."""
+    from dmipy_sim.acquisition.waveforms import btensor_invariants, calc_btensor
+    B = calc_btensor(seq)
+    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,66 @@ 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 tensor_waveform(G, dt, gate, label):
+    """The class's three shapes as played (physical gradient, mT/m) over the echo time, the 180 marked."""
+    plt = _mpl()
+    t = np.arange(G.shape[0]) * dt * 1e3
+    fig, ax = plt.subplots(figsize=(7.5, 2.8))
+    for k, c in enumerate(AXIS_COLOURS):
+        ax.plot(t, G[:, k] * 1e3, color=c, lw=1.2, label=f"ψ{k + 1}")
+    flip = np.nonzero(np.diff(np.asarray(gate)) != 0)[0]
+    for f in flip:
+        ax.axvline(t[f] + 0.5 * (t[1] - t[0]), color="0.5", ls=":", lw=1)
+    ax.set_xlabel("t (ms)"); ax.set_ylabel("G (mT/m)"); ax.set_title(label, fontsize=9); ax.legend(fontsize=8, ncol=3)
+    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
@@ -219,7 +221,7 @@ def test_the_estimated_seconds_grow_with_the_run_and_stay_in_the_pools_window():
     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 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 +232,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 +286,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_pipeline.py

@@ -13,12 +13,51 @@ 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 dmipy_sim.acquisition.waveforms import calc_b
+    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(calc_b(seq), 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 +161,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 +307,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_sim.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 +369,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 +385,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."""
+    """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_sim.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 +423,17 @@ 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

tests/test_tensor.py

@@ -0,0 +1,108 @@
+"""The b-tensor class's encodings (space/tensor.py) without the layout: every row's b-tensor, read off the sequence by
+dmipy-sim, 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_sim.acquisition.waveforms import btensor_invariants, calc_btensor
+    seq, enc, bs = _rows()
+    B = calc_btensor(seq)
+    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)
+    from dmipy_sim.acquisition.waveforms import calc_btensor
+    B = calc_btensor(seq) * 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-sim 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_sim.acquisition.waveforms import btensor_invariants, calc_btensor
+    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(calc_btensor(T._sequence(G)))
+    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)
+    from dmipy_sim.acquisition.waveforms import calc_btensor
+    B = calc_btensor(seq) * 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/deploy.py

@@ -44,6 +44,12 @@ 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 ZeroGPU requirements of a configuration: ``requirements-zero.txt`` for the default (the DiSCo page's pins),
+    ``requirements-<stem>-zero.txt`` else (the brain Space keeps the pins its packs were accepted on)."""
+    return "requirements-zero.txt" if config == DEFAULT_CONFIG else f"requirements-{os.path.splitext(config)[0]}-zero.txt"
+
+
 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 +85,7 @@ 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"))
+        shutil.copy(os.path.join(root, requirements_for(a.config)), os.path.join(staged, "requirements.txt"))
         folder = staged
     elif a.base_tag:
         staged = tempfile.mkdtemp(prefix="disco-space-")
@@ -91,7 +97,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/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()