Choosing the indexing parameters

python -m anri.index runs from a handful of options. This notebook runs the same steps one at a time on your dataset, with a diagnostic plot at each step, so that you can see what each option does to your data before a long run, in the spirit of ImageD11’s S3DXRD notebooks (tomo_1_index, tomo_2_map, pbp_1_indexing, pbp_3_refinement).

Section	Options	What the plots show
Rings	--rings, --tth-tol, --cif	the 2θ profile and the (2θ, η) cake: which rings are clean, how wide they are, what lies between them
Data	--monitor, --gridstep	intensity per dty row and per ω: flux steps between the rows’ scans
Rotation axis, mask	--y0, --mask	the sinogram’s centre of mass, a back-projection for any y0, and the sample mask from it (thresholded or drawn)
Lit map	--lit, --etacut	how much of the map is lit against how much intensity it holds; what the η cut removes
Grid	--grid, --max-chance	chance completeness and size of every grid step, and how much each finer step improves the fit
Pruning	--prune, --min-comp, --min-lr, --keep	completeness and likelihood ratio of every orientation, and how many each cut keeps
Occupancy	--cand, --iter, --block-gb, --beam, --censor	the deviance per iteration, how many voxels still change orientation, and whether the candidates suffice
Model against data		what the fit leaves unexplained, per row, per ring and spot by spot
Populations	--occupied, --min-frac	the occupied voxels within the mask, and the populations per voxel

Where a choice is cheap to redo, a slider redraws the plot, and the value on the slider is what the later cells use. The last cell prints the python -m anri.index command with every value chosen.

The notebook mirrors anri/index/__main__.py: the same binning, geometry, and calls into anri.index, so a choice made here behaves the same in the command-line run. Run it top to bottom; after moving a slider, re-run the cells below it. Figures are interactive (ipympl, %matplotlib widget): zoom and pan with each figure’s toolbar, and the zoom stays while the sliders move; panels that show the same quantity share their axes, so they zoom together. It needs ipywidgets, ipympl and matplotlib. Sections 7 to 9 run the full occupancy fit, which costs as much here as in python -m anri.index (set run_fit = False to stop after pruning on big datasets).

With no dataset given, it renders the 316L phantom as in the indexing tutorial (about a minute on a laptop CPU), so it runs out of the box.

# Parameters: the options of python -m anri.index, at its defaults. The sliders below start from these values.

# the dataset, as python -m anri.index <analysisroot> <sample> <dataset>; None renders a phantom to try things on
analysisroot = None
sample = None
dataset = None
phase = None  # --phase: name in pars.json (None: the only one)
parfile = None  # --parfile (None: the DataSet's, else pars/pars.json beside PROCESSED_DATA)
cif = None  # --cif: CIF of the phase, for structure factors (None: |F|^2 = 1)
monitor = None  # --monitor: counter to normalise intensities by, e.g. "fpico6" (None: none)
gridstep = 1  # --gridstep: voxel = gridstep x dty step, rows summed to match

# 1. rings
n_rings = 6  # --rings
tth_tol = None  # --tth-tol, deg (None: measured per ring)
# 3. rotation axis and sample mask
y0 = None  # --y0 (None: the DataSet's)
mask = None  # --mask: "auto" (threshold of the back-projection), a .npy file, or None (every voxel)
# 4. lit map
lit_cut = 1.0  # --lit: x the median non-empty bin
etacut = 0.2  # --etacut
# 5. grid
grid = None  # --grid, deg (None: the coarsest with chance completeness <= max_chance)
max_chance = 0.5  # --max-chance
# 6. pruning
prune = "likelihood"  # --prune: "likelihood" or "completeness"
min_comp = None  # --min-comp (None: the chance level for likelihood, halfway to the maximum for completeness)
min_lr = 25.0  # --min-lr
keep = 100000  # --keep
# 7. occupancy
n_cand = 64  # --cand
n_iter = 50  # --iter
coarse = 1  # --coarse: only passed on to the command line; this notebook fits at full resolution
block_gb = 1.0  # --block-gb
beam = 0.0  # --beam: FWHM of the beam across dty (0: the 2-row model, close to FWHM = dty step)
censor = 0.0  # --censor: counts per bin below which an empty bin agrees with the model (0: off)
# 9. populations
occupied_cut = 0.2  # --occupied
min_frac = 0.1  # --min-frac

n_cpu = 4  # --n-cpu
outdir = "."  # --outdir, for the command line at the end

# this notebook only
max_rings_shown = 10  # rings drawn on the 2theta profile, to choose n_rings from
run_fit = True  # False: stop after pruning (sections 7-9 run the full occupancy fit)
draw = False  # draw the sample mask by hand in section 3 (needs %matplotlib widget)
grid_scan = True  # fit the orientations of every grid step in section 5, to compare their deviance

import anri.utils

anri.utils.setup(n_cpu=n_cpu)  # before any JAX computation

%matplotlib widget

import os
import shlex
import tempfile
import time
from functools import partial

import h5py
import ipywidgets as widgets
import jax
import jax.numpy as jnp
import numpy as np
from IPython.display import display
from matplotlib import pyplot as plt
from matplotlib.colors import LogNorm

import anri.crystal
import anri.geom
import anri.index as ix
import anri.io

# the binning of python -m anri.index: the lit map in (eta, omega) bins of 0.5 x 0.25 deg, all rows summed; the fit's
# histogram in bins R_E x R_O times larger (1 x 1 deg), per dty row
B_E, B_O = 0.5, 0.25
R_E, R_O = 2, 4
N_POP = 4
MAGMA = plt.colormaps["magma"].with_extremes(bad="black")  # log-scale maps: empty (masked) bins black, not white
T_START = time.time()


def figure(name, **kwargs):
    # a figure shown where display(fig.canvas) puts it (under its sliders), not at the end of the cell; re-running the
    # cell replaces it
    plt.close(name)
    with plt.ioff():
        fig = plt.figure(name, layout="constrained", **kwargs)
    fig.canvas.header_visible = False
    return fig


def freeze(fig):
    # fix every axis's limits as first drawn: the sliders then update the plots without undoing a zoom
    for a in fig.axes:
        if getattr(a, "_colorbar", None) is None:  # colorbars follow their images' limits
            a.set_xlim(a.get_xlim())
            a.set_ylim(a.get_ylim())

The dataset

Without a dataset, the 316L phantom of the indexing tutorial is rendered into a temporary folder as an ImageD11 dataset (this needs ImageD11). Its rotation axis is at dty = 0.3 µm.

if analysisroot is None:
    from ImageD11.sinograms.tensor_map import TensorMap

    truth = TensorMap.from_h5(os.path.join("..", "..", "..", "tests", "data", "phantoms", "am316l", "am316l_tmap.h5"))
    a = truth.phases[0].lattice_parameters[0]
    wl = 0.2843
    pars = {
        "y_center": 1023.5, "y_size": 75.0, "tilt_y": 1e-3, "z_center": 1023.5, "z_size": 75.0, "tilt_z": -2e-3,
        "tilt_x": 0.0, "distance": 150e3, "o11": -1, "o12": 0, "o21": 0, "o22": -1, "wavelength": wl, "wedge": 0.0,
        "chi": 0.0, "omegasign": 1.0, "t_x": 0.0, "t_y": 0.0, "t_z": 0.0,
    }
    y0_true = 0.3
    geom_r = anri.io.geom_from_pars(pars, y0_true, wl * 2e-4 / 2.355, 5e-5, 5e-5, sig_beam=1.4 / 2.355, voxel_size=0.5,
                                    sig_psf=0.5)
    entries = anri.io.entries_from_tensormap(truth)
    entries["density"] = np.full(len(entries["pos"]), 30.0)  # counts scale
    rings8 = ix.ring_table(np.array([a, a, a, 90.0, 90.0, 90.0]), 225, wl, 8)
    omega_r, dty_r = anri.io.motor_grid((0.0, 180.0), 0.1, (y0_true - 30.0, y0_true + 30.0), 1.0)
    analysisroot, sample, dataset = tempfile.mkdtemp(), "phantom", "am316l"
    sparse = os.path.join(analysisroot, "phantom_sparse.h5")
    t0 = time.time()
    anri.io.simulate_sparse(sparse, entries, rings8["hkls"], np.ones(len(rings8["hkls"])), geom_r, omega_r, dty_r,
                            (2048, 2048), max_frames=61)
    cell_r = {"cell__a": a, "cell__b": a, "cell__c": a, "cell_alpha": 90.0, "cell_beta": 90.0, "cell_gamma": 90.0,
              "cell_lattice_[P,A,B,C,I,F,R]": 225}
    pars_r = anri.io.write_pars(os.path.join(analysisroot, "pars"), pars, {"316L": cell_r})
    anri.io.write_dataset(sparse, analysisroot, sample, dataset, y0=y0_true, parfile=pars_r)
    outdir = analysisroot  # the results go beside the phantom, not into the current folder
    print(f"phantom rendered in {time.time() - t0:.0f} s into {analysisroot}")

# as python -m anri.index: the DataSet, its sparse pixels, and the phase from pars.json
dsname = f"{sample}_{dataset}"
dsfile = os.path.join(analysisroot, sample, dsname, f"{dsname}_dataset.h5")
sparsefile = os.path.join(analysisroot, sample, dsname, f"{dsname}_sparse.h5")
ds = anri.io.read_dataset(dsfile)
if ds["sparsefile"] and os.path.exists(ds["sparsefile"]):
    sparsefile = ds["sparsefile"]
Y0_DS = ds["y0"] if y0 is None else y0
if Y0_DS is None:
    raise ValueError(f"{dsfile} has no y0: set y0 in the parameters (section 3 then helps to check it)")
pfile = parfile or ds["parfile"]
if not pfile or not os.path.exists(pfile):  # e.g. processed elsewhere: pars/ beside PROCESSED_DATA
    root = os.path.abspath(analysisroot)
    while os.path.basename(root) != "PROCESSED_DATA" and root != os.path.dirname(root):
        root = os.path.dirname(root)
    pfile = os.path.join(os.path.dirname(root), "pars", "pars.json")
geo, phase_name, cell = anri.io.read_pars_json(pfile, phase)
lpars = np.array([cell[k] for k in ("cell__a", "cell__b", "cell__c", "cell_alpha", "cell_beta", "cell_gamma")])
sg = int(cell["cell_lattice_[P,A,B,C,I,F,R]"])  # a space-group number is needed, not a centring letter

ybin, yedge, oedge = ds["ybincens"], ds["ybinedges"], ds["obinedges"]
ystep0, nk0 = float(np.median(np.diff(ybin))), len(ybin)
YSTEP, DTY0, NK = gridstep * ystep0, float(ybin[0]) + 0.5 * (gridstep - 1) * ystep0, -(-nk0 // gridstep)
OM0, OSTEP = float(oedge[0]), float(np.median(np.diff(oedge)))
N_E, N_O = round(360 / B_E), round(float(oedge[-1] - oedge[0]) / B_O)
N_O -= N_O % R_O
bins = (B_E * R_E, B_O * R_O, N_E // R_E, N_O // R_O)  # the fit's histogram: bin widths and counts in eta, omega
n_e, n_o = bins[2], bins[3]
dty_c = DTY0 + np.arange(NK) * YSTEP  # row centres
om_c = OM0 + (np.arange(n_o) + 0.5) * bins[1]  # omega bin centres of the fit's histogram

WL = geo["wavelength"]
geom = anri.io.geom_from_pars(geo, Y0_DS, WL * 2e-3 / 2.355, 1.5e-4, 1.5e-4, sig_beam=YSTEP / 2.355, voxel_size=YSTEP)
geom = {k: jnp.asarray(v, jnp.float32) if np.issubdtype(np.asarray(v).dtype, np.floating) else v
        for k, v in geom.items()}
B64 = anri.crystal.B_matrix(lpars)
B = B64.astype(np.float32)
ops = anri.crystal.laue_rotations(anri.crystal.symmetry_matrices(sg), B64)
if cif is None:
    structure = None
else:
    import Dans_Diffraction

    structure = Dans_Diffraction.Crystal(cif)

with h5py.File(sparsefile, "r") as h:
    groups = list(h.keys())
    n_pix = {g: int(h[g]["nnz"][()].sum()) for g in groups}
chunk = int(min(1 << 24, 1 << max(10, int(np.ceil(np.log2(max(max(n_pix.values()), 1)))))))
spatial = anri.io.read_spatial(ds)  # the DataSet's detector distortion maps (e2dx/e2dy or detector file), if any
spatial_src = ds["detectorh5"] or (f"{ds['e2dxfile']}, {ds['e2dyfile']}" if spatial is not None else "")
monitor_ref = None
if monitor:  # one reference for every scan, the mean, as python -m anri.index
    mons = anri.io.read_monitor(sparsefile, groups, monitor, ds["masterfile"])
    monitor_ref = float(np.mean(np.concatenate(list(mons.values()))))


def stream(groups_):
    # sparse pixels of some scans, a chunk at a time, read ahead in a thread, as python -m anri.index
    return anri.io.prefetch(anri.io.stream_sparse(sparsefile, yedge, ds["omegamotor"], ds["dtymotor"], chunk, groups_,
                                                  gridstep, ds["dty"], ds["scans"], monitor, monitor_ref,
                                                  ds["masterfile"], ds["omega"], spatial))




# the sinogram's rows, as the DataSet has binned the frames (ds.dty, ds.omega: [rows, frames]); stream_sparse takes
# each frame's dty and omega from there. How far dty moves within a row, in dty steps:
dty_rows = None if ds["dty"] is None else np.asarray(ds["dty"], float)
drift = 0.0 if dty_rows is None else float(np.ptp(dty_rows, axis=1).max()) / ystep0
print(f"dataset {dsfile}\nsparse pixels {sparsefile}: {len(groups)} scans, {sum(n_pix.values()) / 1e6:.1f}M pixels")
advice = (": its frames alias between neighbouring rows (gaps and doubled rows in the sinogram); replace ds.dty by "
          "its mean over each row in ImageD11, and save the DataSet")
print(f"DataSet: {len(ds['scans'] or [])} rows ({', '.join((ds['scans'] or [])[:2])}, ...); its dty moves by up to "
      f"{drift:.2f} dty steps within a row" + (advice if drift > 0.2 else ""))
# where the beam really was: each row's DataSet dty plus dty's motion within the row, read from the sparse file (fly
# and helical scans, where dty moves while omega turns); the model follows it (scan["ddty"]), as python -m anri.index
ddty = None
if dty_rows is not None and ds["scans"] is not None and ds["omega"] is not None:
    dd_ = ix.dty_offsets(anri.io.read_frame_dty(ds, sparsefile), ds["omega"], yedge, gridstep, OM0, bins[1], n_o, DTY0,
                         YSTEP)
    if np.abs(dd_).max() > 1e-3 * YSTEP:
        ddty = jnp.asarray(dd_, jnp.float32)
print("dty within the rows: " + ("constant (a step scan)" if ddty is None else
      f"moves up to {float(jnp.abs(ddty).max()) / YSTEP:.2f} dty steps from each row's dty; the model follows it"))
print(f"spatial correction: {spatial_src or 'none (the DataSet names no e2dx/e2dy or detector file)'}")
print(f"pars {pfile}: phase {phase_name}, lattice {', '.join(f'{v:g}' for v in lpars)}, space group {sg}, "
      f"{len(ops)} Laue-group rotations")
print(f"wavelength {WL:.5f}, distance {geo['distance']:g}; y0 {Y0_DS:.6g}; dty {DTY0:.6g} + {NK} x {YSTEP:.6g}; "
      f"omega {OM0:.4g} .. {oedge[-1]:.4g} in {len(oedge) - 1} frames of {OSTEP:.4g}")

1. Rings: --rings, --tth-tol, --cif

Pixels count only within each ring’s 2θ window, so the windows decide what the indexer sees. As python -m anri.index does, the profile is measured on 9 scans spread over the dty range, and each ring’s offset and half-width (holding 95% of its intensity above the local background) are measured from it; the window is |offset| + half-width, unless tth_tol is set.

What to look for:

• Rings that overlap or touch (the gap column, or windows that meet in the plots): their spots cannot be told apart, so a reflection is predicted in both. Use fewer rings, or a tolerance below the gap.

• Intensity between rings, or a ring that is not this phase’s: another phase, or a wrong lattice. It is ignored by the indexer, but if it falls in a window it lights the map for nothing.

• Detector distortion widens every ring: a pixel or two of distortion is about a ring’s width. The pixels are corrected with the DataSet’s e2dx/e2dy or detector file, as ImageD11 corrects its peaks (the dataset cell says which); if it names none, the widths include the distortion.

• A ring that wobbles with η in the cake (a sinusoid): the geometry (beam centre, tilts) is off. Fix the calibration rather than widening tth_tol, which would also take in more of the neighbours.

• Weak rings add predictions but little intensity: each ring costs time in every step. Six to eight rings are usually enough; more help small grains and lower symmetry.

rings_all = ix.ring_table(lpars, sg, WL, max_rings_shown, structure)
n_all = len(rings_all["tth"])
sample_groups = [groups[i] for i in np.unique(np.linspace(0, len(groups) - 1, min(len(groups), 9)).round().astype(int))]


@partial(jax.jit, static_argnames=("n_t", "n_eta"))
def cake(x, val, lo, step, n_t, n_eta):
    # sum pixel intensities in (2theta, eta) bins: n_t of width step from lo, and n_eta over 360 deg from -180
    it = jnp.floor((x[:, 0] - lo) / step).astype(jnp.int32)
    ie = jnp.floor((x[:, 1] + 180.0) * n_eta / 360.0).astype(jnp.int32) % n_eta
    ok = (it >= 0) & (it < n_t)
    idx = jnp.where(ok, it * n_eta + ie, n_t * n_eta)
    return jax.ops.segment_sum(jnp.where(ok, val, 0.0), idx, n_t * n_eta + 1)[:-1].reshape(n_t, n_eta)


TTH_STEP, N_CAKE = 0.002, 180  # 2theta bins as ix.ring_profile; eta in 2 deg bins
tth_lo = float(rings_all["tth"][0]) - 0.5
n_t = int(np.ceil((float(rings_all["tth"][-1]) + 0.5 - tth_lo) / TTH_STEP))
C = jnp.zeros((n_t, N_CAKE), jnp.float32)
t0 = time.time()
for slow, fast, om, _, val in stream(sample_groups):
    m = len(val)

    def pad(a, m=m):
        return jnp.asarray(np.pad(np.asarray(a, np.float32), (0, chunk - m)))

    C = C + cake(ix.pixel_angles(pad(slow), pad(fast), pad(om), geom), pad(val), tth_lo, TTH_STEP, n_t, N_CAKE)
C = np.asarray(C)
prof = C.sum(1)  # the 2theta profile, as ix.ring_profile measures it
tth_x = tth_lo + (np.arange(n_t) + 0.5) * TTH_STEP
print(f"profile of {len(sample_groups)} scans: {time.time() - t0:.0f} s")

# each ring's reflections: a representative hkl, the multiplicity and the mean |F|^2
ring_of_h = rings_all["ring_j"][::2]
ring_info = []
for r in range(n_all):
    h = rings_all["hkls"][ring_of_h == r]
    hkl = tuple(int(v) for v in sorted(np.abs(h[0]), reverse=True))
    ring_info.append((hkl, len(h), float(rings_all["F2"][ring_of_h == r].mean())))
gap = np.diff(rings_all["tth"])

fig = figure("rings", figsize=(14, 9))
sf = fig.subfigures(2, 1, height_ratios=[2.3, 0.9])
ax_p, ax_c = sf[0].subplots(2, 1, sharex=True)
ax_p.semilogy(tth_x, np.maximum(prof, 0.5), lw=0.7, color="k")
for t in rings_all["tth"]:
    ax_p.axvline(t, color="0.6", ls="--", lw=0.6)
ax_p.set(ylabel="counts", title="2θ profile: every ring (grey), windows of the rings used (orange)")
k = 5  # 0.01 deg in 2theta
Cb = C[: n_t // k * k].reshape(-1, k, N_CAKE).sum(1)
ax_c.imshow(np.ma.masked_less_equal(Cb.T, 0), aspect="auto", origin="lower", cmap=MAGMA, norm=LogNorm(),
            extent=(tth_lo, tth_lo + len(Cb) * k * TTH_STEP, -180, 180), interpolation="nearest")
ax_c.set(xlabel="2θ (°)", ylabel="η (°)", title="cake: rings should be straight; windows in cyan")
axz = np.atleast_1d(sf[1].subplots(1, n_all, sharex=True, sharey=True))
for r, a in enumerate(axz):
    z = np.abs(tth_x - rings_all["tth"][r]) < 0.15
    a.semilogy(tth_x[z] - rings_all["tth"][r], np.maximum(prof[z], 0.5), lw=0.8, color="k")
    a.set_title(f"{r}: {ring_info[r][0]}", fontsize=9)
    a.tick_params(labelsize=7)
axz[0].set_ylabel("counts")
sf[1].supxlabel("2θ − ring (°); orange: window, blue: measured centre; grey: ring not used", fontsize=9)
drawn = []  # the windows, redrawn when the sliders move


def show_rings(n_rings, tth_tol):
    tth = rings_all["tth"][:n_rings]
    off, hw = ix.ring_widths(prof, tth_lo, TTH_STEP, tth)
    tol = np.abs(off) + hw if tth_tol == 0 else np.full(n_rings, tth_tol)
    inside = np.zeros(n_t, bool)
    for t, w in zip(tth, tol):
        inside |= np.abs(tth_x - t) < w
    print(" ring  hkl        mult  |F|^2  2theta   offset  half-width  tolerance  gap to next  counts in window")
    for r in range(n_rings):
        win = np.abs(tth_x - tth[r]) < tol[r]
        g = f"{gap[r]:11.3f}" if r < len(gap) else "          -"
        warn = "  <- windows overlap" if r < len(gap) and r + 1 < n_rings and tol[r] + tol[r + 1] > gap[r] else ""
        print(f"{r:5d}  {str(ring_info[r][0]):10s} {ring_info[r][1]:4d}  {ring_info[r][2]:5.2f}  {tth[r]:6.3f}  "
              f"{off[r]:+7.3f}  {hw[r]:10.3f}  {tol[r]:9.3f}  {g}  {prof[win].sum():.3g}{warn}")
    used = tth_x < tth[-1] + 0.5
    print(f"{prof[inside].sum() / prof[used].sum() * 100:.1f}% of the intensity up to the last ring used is in a "
          "window")
    for art in drawn:
        art.remove()
    drawn.clear()
    for r in range(n_rings):
        drawn.append(ax_p.axvspan(tth[r] - tol[r], tth[r] + tol[r], color="tab:orange", alpha=0.4))
        drawn.extend(ax_c.axvline(tth[r] + sg_ * tol[r], color="c", lw=0.6) for sg_ in (-1, 1))
        drawn.append(axz[r].axvspan(-tol[r], tol[r], color="tab:orange", alpha=0.4))
        drawn.append(axz[r].axvline(off[r], color="tab:blue", lw=0.8))
    for r, a in enumerate(axz):
        a.set_facecolor("white" if r < n_rings else "0.9")
    fig.canvas.draw_idle()


w_rings = widgets.interactive(
    show_rings,
    n_rings=widgets.IntSlider(value=n_rings, min=1, max=n_all, description="n_rings", continuous_update=False),
    tth_tol=widgets.FloatSlider(value=tth_tol or 0.0, min=0.0, max=max(0.3, 2 * (tth_tol or 0.0)), step=0.002,
                                readout_format=".3f", description="tth_tol", continuous_update=False),
)
freeze(fig)
print("tth_tol 0: measured per ring")
display(w_rings, fig.canvas)

2. The data: --monitor, --gridstep

All the pixels go once into two histograms, as in python -m anri.index: the lit map (0.5° in η, 0.25° in ω, all rows summed) for pruning, and the fit’s data (1° × 1°, per dty row). This is the one pass over the whole dataset (re-run it after changing the rings).

Then the total intensity of each dty row and of each ω bin. Across dty, a sample gives a smooth profile; steps from one row to the next mean that the flux changed between the rows’ scans (e.g. a refill), which makes ring artefacts centred on the rotation axis: normalise with --monitor (e.g. fpico6). Along ω, a slow drift is beam decay within a rotation. Sharp spikes in either are single strong spots, not flux.

--gridstep sums rows in groups to make voxels larger: fewer voxels and rows, a much cheaper fit, coarser maps.

Each frame’s ω and dty come from the DataSet (ds.omega, ds.dty, reshaped by ImageD11 to the sinogram’s rows; a fly scan in one scan is split into its rotations, 1.1::[0:1440], …), not from the sparse file’s raw readings. The last plot is how the DataSet’s dty moves within each row. In a fly scan over ω and dty together, the readings drift slowly through each rotation; binned by them, each frame falls in the row of the moment, so the rows alias: gaps in some rows of the sinogram and doubled intensity in others, and a scalloped sample edge. ImageD11 only reshapes the readings: if they drift here, replace ds.dty by its mean over each row in ImageD11, save the DataSet, and re-run from the top.

n_rings_used = w_rings.kwargs["n_rings"]
rings = ix.ring_table(lpars, sg, WL, n_rings_used, structure)
ring_off, ring_hw = ix.ring_widths(prof, tth_lo, TTH_STEP, rings["tth"])
rings["hw"] = ring_hw
tth_tol_used = w_rings.kwargs["tth_tol"] or None
tth_tol_r = np.abs(ring_off) + ring_hw if tth_tol_used is None else np.full(n_rings_used, tth_tol_used)
print(f"{len(rings['hkls'])} hkls in {n_rings_used} rings; tolerances (deg) {np.round(tth_tol_r, 3)}")
n_cells = n_rings_used * n_e * n_o * NK
print(f"histogram: {n_cells / 1e6:.0f}M bins, {n_cells * 4 / 1e9:.2f} GB (x ~4 in the fit)")

t0 = time.time()
H_lit, H = ix.histogram_pixels(stream(groups), geom, rings["tth"], tth_tol_r, OM0,
                               [((B_E, B_O, N_E, N_O), 1), (bins, NK)], chunk)
Hs = H_lit.reshape(n_rings_used, N_E, N_O)
H4 = np.asarray(H).reshape(n_rings_used, n_e, n_o, NK)  # [ring, eta, omega, row] on the host, for the plots
print(f"histograms: {time.time() - t0:.0f} s; {float(jnp.mean(H > 0)) * 100:.1f}% of the bins non-empty")

row_tot = H4.sum((0, 1, 2))
om_tot = H4.sum((0, 1, 3))
nb = np.convolve(row_tot, [0.5, 0.0, 0.5], "same")  # the mean of each row's neighbours
lit_row = row_tot > 0.05 * row_tot.max()
lit_row[[0, -1]] = False
jump = np.where(lit_row, row_tot / np.maximum(nb, 1e-30), np.nan)

n_ax = 4 if monitor else 3
fig, ax = plt.subplots(1, n_ax, figsize=(5 * n_ax, 3.6), layout="constrained")
if dty_rows is not None:
    for v in dty_rows[:: max(1, len(dty_rows) // 10)]:
        ax[-1].plot(np.linspace(0, 1, len(v)), (v - v.mean()) / ystep0, lw=0.7)
ax[-1].set(xlabel="fraction of the row", ylabel="dty − row mean (dty steps)", title="the DataSet's dty within a row")
ax[0].plot(dty_c, row_tot, ".-", ms=3)
ax0 = ax[0].twinx()
ax0.plot(dty_c, jump, ".", color="tab:red", ms=3)
ax0.axhline(1, color="tab:red", lw=0.5)
ax0.set_ylabel("row / mean of its neighbours", color="tab:red")
ax[0].set(xlabel="dty", ylabel="intensity", title="per dty row")
ax[1].plot(om_c, om_tot, lw=0.8)
ax[1].set(xlabel="ω (°)", ylabel="intensity", title="per ω bin, all rows")
if monitor:
    for g in groups[:: max(1, len(groups) // 20)]:
        ax[2].plot(mons[g] / monitor_ref, lw=0.5)
    ax[2].set(xlabel="frame", ylabel=f"{monitor} / mean", title=f"{monitor}, some scans")
plt.show()
print(f"row-to-row scatter (rows with intensity): {np.nanstd(jump) * 100:.1f}%; sample edges excluded")

3. Rotation axis: --y0

y0 is the dty at which the rotation axis is in the beam. A voxel at (x, y) in the sample is in the beam at dty = y0 − (x sin ω + y cos ω), so any point fixed in the sample traces y0 + b sin ω + c cos ω in the sinogram, and a fit of that curve gives y0. ImageD11’s tomo_2_map fits each grain’s sinogram this way.

Here the point is the middle of the sample’s projection: at each ω, halfway between the first and last dty rows above 5% of that projection’s maximum (the sample’s edges, interpolated between rows). The edges are sharp whichever grains diffract, so this is much steadier than the intensity-weighted mean dty (also shown), which jumps as grains light up and go dark: on the phantom (true y0 0.3, 1 µm steps) the edges give 0.19, the mean −0.03. Both need the whole sample within the scanned dty range.

Check it on the reconstruction, made as ImageD11’s tomo_2_map makes its whole-sample mask: a sinogram of log intensities, so that every spot counts about the same however bright (ImageD11 sums log(intensity) of its 2D peaks; here, log(1 + counts) of each (ring, η) bin of the fit’s histogram, where a spot fills about one bin), and filtered back-projection with a Hamming-windowed ramp, which tames the noise. With raw intensities the few brightest Bragg peaks dominate and smear into streaks. With the wrong y0, the sample’s edge blurs and the network of grain boundaries inside it (voxels where two grains’ spots add up) doubles and smears; on the phantom that shows clearly at 3 dty steps, faintly at 1. The slider moves y0 by up to 10 steps either way; the reconstruction is on the indexer’s voxel grid, with its model of where a voxel is in the beam.

sino = H4.sum((0, 1))  # [omega, row], every ring and eta
w_om = sino.sum(1)
good = w_om > 0.05 * w_om.max()
com = sino @ dty_c / np.maximum(w_om, 1e-30)  # intensity-weighted mean dty


def mid(p, frac=0.05):
    # the middle of one projection: halfway between its edges, where it crosses frac x its maximum
    thr = frac * p.max()
    k = np.flatnonzero(p > thr)
    i0, i1 = k[0], k[-1]
    lo = dty_c[i0] - YSTEP * (p[i0] - thr) / max(p[i0] - (p[i0 - 1] if i0 > 0 else 0.0), 1e-30)
    hi = dty_c[i1] + YSTEP * (p[i1] - thr) / max(p[i1] - (p[i1 + 1] if i1 < NK - 1 else 0.0), 1e-30)
    return 0.5 * (lo + hi)


mids = np.array([mid(p) if g else np.nan for p, g in zip(sino, good)])
A = np.stack([np.ones(n_o), np.sin(np.radians(om_c)), np.cos(np.radians(om_c))], 1)
coef = np.linalg.lstsq(A[good], mids[good], rcond=None)[0]
coef_com = np.linalg.lstsq(A[good], com[good], rcond=None)[0]
y0_fit = float(coef[0])
edge = max(row_tot[0], row_tot[-1]) / row_tot.max()
print(f"y0 from the sample's edges: {y0_fit:.4f} (DataSet / parameters: {Y0_DS:.4f}; "
      f"{(y0_fit - Y0_DS) / YSTEP:+.2f} dty steps); from the mean dty: {coef_com[0]:.4f}")
if edge > 0.1:
    print(f"the first or last row holds {edge * 100:.0f}% of the brightest row's intensity: the sample may extend past "
          "the scan, which biases the fit")


_, pad_ds = anri.geom.sino_shift_and_pad(Y0_DS, NK, DTY0, YSTEP)
nr_ds = NK + pad_ds
pos_ds = anri.geom.recon_positions(nr_ds, YSTEP).astype(jnp.float32)
sino_log = np.log1p(H4).sum((0, 1))  # [omega, row]: every spot counts about the same, as in tomo_2_map
sino_f = jnp.asarray(ix.ramp_filter(sino_log), jnp.float32)
om_j = jnp.asarray(om_c, jnp.float32)
rec_ds = np.asarray(ix.backproject(sino_f, om_j, pos_ds, Y0_DS, DTY0, YSTEP, ddty)).reshape(nr_ds, nr_ds)
vmax = np.percentile(rec_ds, 99.5)
ext = (-nr_ds / 2 * YSTEP, nr_ds / 2 * YSTEP) * 2
ext_sino = (OM0, OM0 + n_o * bins[1], DTY0 - YSTEP / 2, DTY0 + (NK - 0.5) * YSTEP)  # omega, dty

fig = figure("y0", figsize=(15, 4.6))
ax = fig.subplots(1, 3)
ax[2].sharex(ax[1])
ax[2].sharey(ax[1])
ax[0].imshow(sino_log.T, aspect="auto", origin="lower", cmap="magma", extent=ext_sino, interpolation="nearest")
ax[0].plot(om_c[good], com[good], ".", color="0.6", ms=2, label="mean dty")
ax[0].plot(om_c[good], mids[good], ".", color="c", ms=2, label="middle of the edges")
ax[0].plot(om_c, A @ coef, color="w", lw=1, label=f"fit: y0 = {y0_fit:.3f}")
y0_line = ax[0].axhline(Y0_DS, color="tab:green", lw=1, ls="--", label="y0 (slider)")
ax[0].legend(fontsize=8, loc="upper right")
ax[0].set(xlabel="ω (°)", ylabel="dty", title="sinogram of log(1 + counts), all rings")
ax[1].imshow(rec_ds, origin="lower", cmap="gray", vmin=0, vmax=vmax, extent=ext)
ax[1].set(title=f"back-projection, y0 = {Y0_DS:.6g} (DataSet)", xlabel="voxels (dty units)")
im_y0 = ax[2].imshow(rec_ds, origin="lower", cmap="gray", vmin=0, vmax=vmax, extent=ext)
ax[2].set(xlabel="voxels (dty units)")


def show_y0(y0_try):
    im_y0.set_data(np.asarray(ix.backproject(sino_f, om_j, pos_ds, y0_try, DTY0, YSTEP, ddty)).reshape(nr_ds, nr_ds))
    y0_line.set_ydata([y0_try, y0_try])
    ax[2].set_title(f"back-projection, y0 = {y0_try:.6g} (slider)")
    fig.canvas.draw_idle()


lo_, hi_ = min(Y0_DS, y0_fit) - 10 * YSTEP, max(Y0_DS, y0_fit) + 10 * YSTEP
w_y0 = widgets.interactive(show_y0, y0_try=widgets.FloatSlider(value=Y0_DS, min=lo_, max=hi_, step=YSTEP / 20,
                                                               readout_format=".6g", description="y0",
                                                               continuous_update=False))
freeze(fig)
display(w_y0, fig.canvas)

Y0 = w_y0.kwargs["y0_try"]  # the rotation axis used from here on
_, PAD = anri.geom.sino_shift_and_pad(Y0, NK, DTY0, YSTEP)
NR = NK + PAD  # the voxel grid: NR x NR, centred on the rotation axis, as ImageD11 pads its reconstructions
NV = NR * NR
print(f"y0 = {Y0:.4f}; voxels {NR} x {NR} of {YSTEP:g} ({NK} dty rows + pad {PAD})")

The sample mask: --mask

The fit only needs the voxels the sample fills. Here is the back-projection at the y0 just chosen, and a mask from it as ImageD11’s tomo_2_map makes its whole-sample mask: the pixels above a threshold (Otsu’s by default; the slider moves it), their largest connected region, and its convex hull. Or draw it by hand (draw = True in the parameters): click the sample’s outline on the image that appears, close it on its first corner, then run the next cell (anri.index.draw_mask, as ImageD11’s InteractiveMask).

With a mask, the voxels outside it are not fitted (fewer voxels, no noise orientations in the air around the sample), and the --occupied cut of section 9 works within it. mask = "auto" in the parameters uses the threshold at the slider; a drawn mask, or a threshold moved off Otsu’s, is saved as a .npy file for --mask. Check that the mask covers the whole sample: a sample that fills the scan leaves no air to threshold against, and needs no mask.

scan0 = {"y0": Y0, "dty0": DTY0, "ystep": YSTEP, "n_rows": NK, "om0": OM0}
if ddty is not None:
    scan0["ddty"] = ddty
rec = ix.reconstruct(H, n_rings_used, n_e, n_o, scan0, bins[1], NR)  # [NR, NR], reconstruction order
t_otsu = ix.otsu(rec)
vmax_r = np.percentile(rec, 99.5)
fig = figure("mask", figsize=(14, 4.4))
ax = fig.subplots(1, 3)
ax[1].sharex(ax[0])
ax[1].sharey(ax[0])
ax[0].imshow(rec, origin="lower", cmap="gray", vmin=0, vmax=vmax_r)
ax[0].set(title=f"back-projection, y0 = {Y0:.6g}", xticks=[], yticks=[])
im_m = ax[1].imshow(rec, origin="lower", cmap="gray", vmin=0, vmax=vmax_r)
ax[1].set(xticks=[], yticks=[])
ax[2].hist(rec.ravel(), bins=200, log=True, color="0.5")
ax[2].axvline(t_otsu, color="k", ls=":", label="Otsu")
v_t = ax[2].axvline(t_otsu, color="tab:red", label="threshold")
ax[2].set(xlabel="back-projection", ylabel="voxels")
ax[2].legend(fontsize=8)


def show_mask(threshold):
    m = ix.threshold_mask(rec, threshold)
    im_m.set_data(np.where(m, rec, np.nan))
    v_t.set_xdata([threshold, threshold])
    ax[1].set_title(f"mask: {m.sum()} of {NV} voxels")
    fig.canvas.draw_idle()


w_mask = widgets.interactive(show_mask, threshold=widgets.FloatSlider(
    value=t_otsu, min=float(np.percentile(rec, 1)), max=float(np.percentile(rec, 99.9)),
    step=float(np.ptp(rec)) / 500, readout_format=".3g", description="threshold", continuous_update=False))
freeze(fig)
display(w_mask, fig.canvas)
if draw:
    get_mask = ix.draw_mask(rec)  # draw the outline on this figure, then run the next cell

# the mask used from here on: drawn, thresholded at the slider (mask = "auto"), a file, or every voxel
mask_used = mask
if draw:
    mask_v = get_mask()
elif mask == "auto":
    mask_v = ix.threshold_mask(rec, w_mask.kwargs["threshold"])
elif mask is not None:
    mask_v = np.load(mask).astype(bool)
else:
    mask_v = np.ones((NR, NR), bool)
if draw or (mask == "auto" and not np.isclose(w_mask.kwargs["threshold"], t_otsu)):
    mask_used = os.path.join(outdir, f"{dsname}_mask.npy")  # for --mask
    np.save(mask_used, mask_v)
mask_v = mask_v.ravel()
print(f"mask: {mask_v.sum()} of {NV} voxels fitted" + (f" (--mask {mask_used})" if mask_used else " (no mask)"))

4. Lit map: --lit and --etacut

The lit map says where there is intensity: a bin is lit above lit × the median non-empty bin. Completeness, the fraction of an orientation’s predicted spots that land on lit bins, is the first filter on the orientation grid, so the lit map should be sparse: the more of it is lit, the more wrong orientations match by chance (section 5). But a high cut drops weak spots, i.e. small grains. The left plot is the trade-off: the lit area against the intensity it holds, as the cut rises. A good cut is near the elbow, where the area falls fast and the intensity slowly.

med = float(jnp.median(Hs[Hs > 0]))
vals = np.sort(np.asarray(Hs).ravel())
cum = np.concatenate([[0.0], np.cumsum(vals, dtype=np.float64)])
mults = np.logspace(-1, 2, 61)
idx = np.searchsorted(vals, mults * med, side="right")
area = 1 - idx / vals.size  # fraction of the bins above the cut
held = 1 - cum[idx] / cum[-1]  # fraction of the intensity in them
Hs_np = np.asarray(Hs)
v_pos = vals[np.searchsorted(vals, 0.0, side="right")]  # the smallest non-zero bin

fig = figure("lit", figsize=(15, 4.6))
ax = fig.subplots(1, 3, width_ratios=[1, 1.3, 1.3])
ax[2].sharex(ax[1])
ax[2].sharey(ax[1])
ax[0].loglog(mults, area, label="lit area")
ax[0].loglog(mults, held, label="intensity held")
v_lit = ax[0].axvline(lit_cut, color="tab:red")
ax[0].set(xlabel="lit cut (x median non-empty bin)", ylabel="fraction", title=f"median non-empty bin: {med:.3g}")
ax[0].legend()
e = (OM0, OM0 + N_O * B_O, -180, 180)
im_int = ax[1].imshow(np.ma.masked_less_equal(Hs_np[0], 0), aspect="auto", origin="lower", cmap=MAGMA,
                      norm=LogNorm(v_pos, vals[-1]), extent=e, interpolation="nearest")
ax[1].set(xlabel="ω (°)", ylabel="η (°)")
im_lit = ax[2].imshow(np.zeros(Hs_np[0].shape), aspect="auto", origin="lower", cmap="gray_r", vmin=0, vmax=1,
                      extent=e, interpolation="nearest")
ax[2].set(xlabel="ω (°)")
ec = np.degrees(np.arcsin(etacut))
for a in ax[1:]:
    for el in (ec, -ec, 180 - ec, ec - 180):  # where |sin eta| = etacut (the parameter)
        a.axhline(el, color="tab:cyan", lw=0.6, ls=":")


def show_lit(lit_cut, ring):
    i = np.searchsorted(vals, lit_cut * med, side="right")
    lit_np = Hs_np[ring] > lit_cut * med
    v_lit.set_xdata([lit_cut, lit_cut])
    im_int.set_data(np.ma.masked_less_equal(Hs_np[ring], 0))
    im_lit.set_data(lit_np.astype(np.float32))
    ax[1].set_title(f"ring {ring}: intensity, rows summed")
    ax[2].set_title(f"ring {ring}: lit, {lit_np.mean() * 100:.2f}% of its bins")
    fig.canvas.draw_idle()
    print(f"lit: {(1 - i / vals.size) * 100:.2f}% of the bins, holding {(1 - cum[i] / cum[-1]) * 100:.1f}% of the "
          f"intensity; per ring {np.round((Hs_np > lit_cut * med).mean((1, 2)) * 100, 2)}%")


w_lit = widgets.interactive(
    show_lit,
    lit_cut=widgets.FloatLogSlider(value=lit_cut, base=10, min=-1, max=2, step=0.05, description="lit",
                                   continuous_update=False),
    ring=widgets.IntSlider(value=0, min=0, max=n_rings_used - 1, description="ring", continuous_update=False),
)
freeze(fig)
display(w_lit, fig.canvas)

Reflections near η = 0 and 180° (12 and 6 o’clock) are poorly placed in ω: a small rotation δ of the crystal moves a spot by up to δ / |sin η| in ω, so near the rotation axis spots smear over many frames and the grid’s tolerances blow up. Reflections with |sin η| ≤ etacut are not used anywhere (dotted lines above). Raising it costs reflections (2 arcsin(etacut) / π of every ring) and the intensity there (right).

eta_x = -180.0 + (np.arange(N_E) + 0.5) * B_E
prof_eta = Hs_np.sum((0, 2))

fig = figure("etacut", figsize=(12, 3.6))
ax = fig.subplots()
ax.plot(eta_x, prof_eta, lw=0.8, color="k")
ax.set(xlabel="η (°)", ylabel="intensity, all rings", xlim=(-180, 180), title="red: not used")
ax2 = ax.twinx()
ax2.semilogy(eta_x, 1 / np.maximum(np.abs(np.sin(np.radians(eta_x))), 1e-3), color="tab:blue", lw=0.8)
ax2.set_ylim(1, 50)
ax2.set_ylabel("ω smearing, 1 / |sin η|", color="tab:blue")
shade = []


def show_etacut(etacut):
    cut = np.abs(np.sin(np.radians(eta_x))) <= etacut
    for art in shade:
        art.remove()
    shade[:] = [ax.fill_between(eta_x, 0, prof_eta.max(), where=cut, color="tab:red", alpha=0.2)]
    fig.canvas.draw_idle()
    print(f"etacut {etacut:.2f}: {2 * np.arcsin(etacut) / np.pi * 100:.1f}% of the reflections and "
          f"{prof_eta[cut].sum() / prof_eta.sum() * 100:.1f}% of the measured intensity not used; spots smeared up to "
          f"{1 / max(etacut, 1e-3):.1f} x in ω are kept")


w_eta = widgets.interactive(show_etacut, etacut=widgets.FloatSlider(value=etacut, min=0.0, max=0.7, step=0.01,
                                                                    description="etacut", continuous_update=False))
freeze(fig)
display(w_eta, fig.canvas)

5. Orientation grid: --grid, --max-chance

The grid covers one fundamental zone. Its step sets how far the truth can be from the nearest grid point, so the matching tolerances, and with them the completeness that wrong orientations reach by chance: the chance completeness, the median over a sample of the grid (most grid orientations are wrong). python -m anri.index takes the coarsest step whose chance completeness is at most max_chance. A finer grid separates real orientations from decoys better and gives more precise orientations, but has many more orientations (the cost of pruning grows with it).

The samples are drawn as anri.index.choose_grid draws them, so the values match the command line’s log. The lit map and η cut chosen above are used; re-run this cell after changing them.

lit_used, etacut_used = w_lit.kwargs["lit_cut"], w_eta.kwargs["etacut"]
lit = {"table": ix.lit_table(Hs > lit_used * med), "om0": OM0, "bins": (B_E, B_O, N_E, N_O), "frame_step": OSTEP,
       "etacut": etacut_used}
rng = np.random.default_rng(0)  # as choose_grid
chance, n_grid = {}, {}
t0 = time.time()
for s in ix.GRID_STEPS:
    if n_grid and n_grid[ix.GRID_STEPS[0]] * (ix.GRID_STEPS[0] / s) ** 3 > 4e6:
        print(f"grid {s} deg: skipped, ~{n_grid[ix.GRID_STEPS[0]] * (ix.GRID_STEPS[0] / s) ** 3 / 1e6:.0f}M "
              "orientations")
        continue
    U_s, d_s = anri.crystal.orientation_grid(s, ops)
    smp = U_s[rng.choice(len(U_s), min(len(U_s), 1 << 14), replace=False)]
    chance[s], n_grid[s] = float(np.median(ix.completeness_of(smp, d_s, B, rings, geom, lit))), len(U_s)
    print(f"grid {s} deg: {len(U_s)} orientations, up to {d_s:.2f} deg from the truth; chance completeness "
          f"{chance[s]:.2f}")
del U_s
print(f"{time.time() - t0:.0f} s")


def auto_step(max_chance):
    # the step choose_grid picks: the coarsest at or below max_chance, else the finest
    return next((s for s in chance if chance[s] <= max_chance), ix.GRID_STEPS[-1])


fig = figure("grid", figsize=(7, 3.6))
ax = fig.subplots()
xs = [str(s) for s in chance]
bars = ax.bar(xs, list(chance.values()))
h_chance = ax.axhline(max_chance, color="k", ls=":")
ax.set(xlabel="grid step (°)", ylabel="chance completeness", ylim=(0, 1))
ax2 = ax.twinx()
ax2.semilogy(xs, list(n_grid.values()), "o-", color="0.4")
ax2.set_ylabel("orientations", color="0.4")


def show_grid(max_chance):
    st = grid if grid is not None else auto_step(max_chance)
    for b, s in zip(bars, chance):
        b.set_color("tab:red" if s == st else "tab:blue")
    h_chance.set_ydata([max_chance, max_chance])
    ax.set_title(f"grid step {st}°" + (" (fixed by grid in the parameters)" if grid is not None else " (automatic)"))
    fig.canvas.draw_idle()


w_grid = widgets.interactive(show_grid, max_chance=widgets.FloatSlider(
    value=max_chance, min=0.05, max=0.95, step=0.05, description="max_chance", continuous_update=False))
freeze(fig)
display(w_grid, fig.canvas)

How fine a grid the fit needs. Chance completeness says how coarse a grid the pruning can work with. The fit itself predicts point spots from grid orientations into 1° × 1° bins: a grain between grid points has its spots up to the grid’s worst case away, and a grid much coarser than the bins cannot represent it (MLEM then spreads it over several grid points, differently from voxel to voxel, and grain boundaries go noisy).

Below, for each grid step, the orientations above chance completeness are fitted to the row-summed data (the orientation fit of section 6). The deviance falls as the grid gets finer while the grid limits the fit, and flattens where the data no longer tell finer grids apart. Choose the coarsest step on the flat part (grid in the parameters; 0.75° is shown as a check that 1° is fine enough). Set grid_scan = False to skip this on big datasets: it costs about one orientation fit per step.

if grid_scan:
    d_rows = H.reshape(-1, NK).sum(1)  # the histogram with its dty rows summed
    dev_s, n_lr_s = {}, {}
    for s in [*chance, 0.75]:
        t0 = time.time()
        U_s, d_s = anri.crystal.orientation_grid(s, ops)
        _, comp_s, info_s = ix.prune(U_s, d_s, B, rings, geom, lit, None, keep)
        pre_s = np.flatnonzero(comp_s > info_s["chance"])
        _, lr_s, mu_s = ix.orientation_mlem(d_rows, U_s[pre_s], B, rings, geom, (*bins, OM0), etacut_used,
                                            log=lambda m: None, censor=censor, return_model=True)
        dev_s[s], n_lr_s[s] = float(ix.deviance(d_rows, jnp.asarray(mu_s), censor)), int(np.sum(lr_s > min_lr))
        print(f"grid {s} deg: {len(pre_s)} orientations above chance; deviance {dev_s[s]:.5g}; {n_lr_s[s]} with "
              f"likelihood ratio > {min_lr:g} ({time.time() - t0:.0f} s)")
    del U_s, comp_s
    steps_s = list(dev_s)
    fig, ax = plt.subplots(1, 2, figsize=(11, 3.6), layout="constrained")
    ax[0].plot(steps_s, [dev_s[s] / dev_s[min(steps_s)] for s in steps_s], "o-")
    ax[0].set(xlabel="grid step (°)", ylabel=f"deviance / deviance at {min(steps_s)}°",
              title="orientation fit to the row-summed data")
    ax[0].invert_xaxis()
    ax[1].plot(steps_s, [n_lr_s[s] for s in steps_s], "o-", color="0.4")
    ax[1].set(xlabel="grid step (°)", ylabel=f"orientations with likelihood ratio > {min_lr:g}", yscale="log")
    ax[1].invert_xaxis()
    plt.show()

6. Pruning: --prune, --min-comp, --min-lr, --keep

Every grid orientation gets its completeness. Real orientations sit well above the chance level; the bulk of the grid sits at it.

• --prune completeness keeps everything at or above min_comp (default halfway from chance to the maximum), the most complete first, at most keep.

• --prune likelihood (default) fits one occupancy per orientation above min_comp (default: the chance level) to the row-summed data, and keeps those whose likelihood ratio is above min_lr: how much the fit worsens without that orientation (about χ² with one degree of freedom for one that is not there, so 25 is about 5σ). Raising min_comp here drops decoys, but also small grains, whose completeness is lower.

step = grid if grid is not None else auto_step(w_grid.kwargs["max_chance"])
U_grid, delta = anri.crystal.orientation_grid(step, ops)
t0 = time.time()
_, comp, info = ix.prune(U_grid, delta, B, rings, geom, lit, None, keep)
halfway = info["min_comp"]  # the default cut of --prune completeness
print(f"grid {step} deg: {len(U_grid)} orientations; completeness ({time.time() - t0:.0f} s): "
      f"chance {info['chance']:.2f}, 99th percentile {np.percentile(comp, 99):.2f}, max {comp.max():.2f}; "
      f"halfway {halfway:.2f}")
comp_sorted = np.sort(comp)


def default_cut(prune):
    return info["chance"] if prune == "likelihood" else halfway


fig = figure("completeness", figsize=(13, 3.6))
ax = fig.subplots(1, 2, sharex=True)
ax[0].hist(comp, bins=100, log=True, color="tab:blue")
ax[0].axvline(info["chance"], color="k", ls=":", label="chance")
ax[0].axvline(halfway, color="0.5", ls="--", label="halfway")
v_comp = [a.axvline(0.0, color="tab:red", label="min_comp") for a in ax]
ax[0].set(xlabel="completeness", ylabel="orientations")
ax[0].legend()
cs = np.linspace(0, 1, 201)
ax[1].semilogy(cs, np.maximum(comp.size - np.searchsorted(comp_sorted, cs), 0.5), color="tab:blue")
ax[1].set(xlabel="min_comp", ylabel="orientations above")


def show_comp(prune, min_comp):
    n_above = comp.size - np.searchsorted(comp_sorted, min_comp, side="left" if prune == "completeness" else "right")
    for v in v_comp:
        v.set_xdata([min_comp, min_comp])
    fig.canvas.draw_idle()
    if prune == "likelihood":
        print(f"{n_above} orientations above {min_comp:.2f} go to the orientation fit")
    else:
        print(f"{n_above} orientations at or above {min_comp:.2f}; {min(n_above, keep)} kept (keep = {keep})")


w_comp = widgets.interactive(
    show_comp,
    prune=widgets.Dropdown(options=["likelihood", "completeness"], value=prune, description="prune"),
    min_comp=widgets.FloatSlider(value=min_comp if min_comp is not None else default_cut(prune), min=0.0, max=1.0,
                                 step=0.01, description="min_comp", continuous_update=False),
)


def reset_cut(change):
    w_comp.children[1].value = default_cut(change["new"])  # a new prune mode starts from its default min_comp


w_comp.children[0].observe(reset_cut, names="value")
freeze(fig)
print(f"min_comp defaults: likelihood {info['chance']:.2f} (chance), completeness {halfway:.2f} (halfway)")
display(w_comp, fig.canvas)

The orientation fit, for --prune likelihood: every orientation above min_comp gets one occupancy, fitted by MLEM to the histogram with its dty rows summed (it explains overlapping spots jointly, with Lorentz, polarisation and |F|²), then its likelihood ratio.

In the plots: real orientations reach high likelihood ratios. On the right, each kept orientation’s likelihood ratio against its misorientation to the nearest kept orientation with a higher ratio. A grain shows as one strong point far from anything stronger (top right: a peak), with a halo of weaker points one or two grid steps from it, its neighbours on the grid soaking up what the grid cannot fit exactly. The number of peaks (beyond the population radius, 1.8 grid steps) is roughly the number of distinct orientations: grains and twins. Weak peaks far from everything are small grains or decoys; min_lr decides which are kept.

prune_used, min_comp_used = w_comp.kwargs["prune"], w_comp.kwargs["min_comp"]
if min_comp is None and np.isclose(min_comp_used, default_cut(prune_used)):
    min_comp_used = None  # the command line's own default
cut = default_cut(prune_used) if min_comp_used is None else min_comp_used
if prune_used == "likelihood":
    pre = np.flatnonzero(comp > cut)
    d_rows = H.reshape(-1, NK).sum(1)  # the histogram with its dty rows summed
    t0 = time.time()
    g_occ, lr = ix.orientation_mlem(d_rows, U_grid[pre], B, rings, geom, (*bins, OM0), etacut_used, censor=censor)
    print(f"orientation fit of {len(pre)} orientations: {time.time() - t0:.0f} s")
    lr_sorted = np.sort(lr)

    # each orientation's misorientation to the nearest one with a higher likelihood ratio (the strongest 5000)
    top = np.argsort(lr)[::-1][: min(len(pre), 5000)]
    q = anri.crystal.mat_to_quat(U_grid[pre[top]]).astype(np.float32)  # [M, 4]
    qs = anri.crystal.mat_to_quat(U_grid[pre[top]][:, None] @ ops[None]).astype(np.float32)  # [M, S, 4]
    dist = np.full(len(top), 180.0)
    for s0 in range(0, len(top), 256):
        c = np.abs(np.einsum("ik,jsk->ijs", q[s0 : s0 + 256], qs)).max(2)  # cos(half the misorientation)
        c[np.arange(c.shape[0])[:, None] <= np.arange(len(top))[None] - s0] = -1.0  # only stronger ones
        dist[s0 : s0 + 256] = np.degrees(2 * np.arccos(np.clip(c.max(1), -1.0, 1.0)))
    dist[0] = 180.0
    lr_top = np.maximum(lr[top], 0.1)

    fig = figure("likelihood", figsize=(16, 3.9))
    ax = fig.subplots(1, 3)
    ax[2].sharey(ax[0])
    ax[0].scatter(comp[pre], np.maximum(lr, 0.1), s=1, alpha=0.4)
    h_lr = [ax[0].axhline(min_lr, color="tab:red", label="min_lr"), ax[2].axhline(min_lr, color="tab:red")]
    ax[0].axvline(cut, color="tab:orange", ls="--", label="min_comp")
    ax[0].set(yscale="log", xlabel="completeness", ylabel="likelihood ratio")
    ax[0].legend(loc="upper left", fontsize=8)
    ts = np.logspace(-1, max(1, np.log10(max(lr.max(), 10))), 100)
    ax[1].loglog(ts, np.maximum(lr.size - np.searchsorted(lr_sorted, ts, side="right"), 0.5))
    v_lr = ax[1].axvline(min_lr, color="tab:red")
    h_keep = ax[1].axhline(keep, color="0.5", ls=":", label="keep")
    ax[1].set(xlabel="min_lr", ylabel="orientations kept", ylim=(0.5, max(2 * len(pre), 2 * keep)))
    ax[1].legend(fontsize=8)
    sc_out = ax[2].scatter(dist, lr_top, s=2, color="0.75", label="not kept")
    sc_in = ax[2].scatter(dist, lr_top, s=3, color="tab:blue", label="kept")
    ax[2].axvline(1.8 * step, color="0.4", ls="--", label="population radius")
    ax[2].set(xscale="log", xlabel="misorientation to the nearest stronger orientation (°)",
              title=f"the {len(top)} strongest")
    ax[2].legend(fontsize=8, loc="lower right")

    def show_lr(min_lr, keep):
        ok = lr > min_lr
        n_ok = int(ok.sum())
        k = pre[ok][np.argsort(lr[ok])[::-1][:keep]] if n_ok else pre[:0]
        sel = np.isin(pre[top], k)
        for h_ in h_lr:
            h_.set_ydata([min_lr, min_lr])
        v_lr.set_xdata([min_lr, min_lr])
        h_keep.set_ydata([keep, keep])
        sc_in.set_offsets(np.c_[dist[sel], lr_top[sel]])
        sc_out.set_offsets(np.c_[dist[~sel], lr_top[~sel]])
        fig.canvas.draw_idle()
        capped = f" (capped by keep = {keep})" if n_ok > keep else ""
        print(f"{n_ok} with likelihood ratio > {min_lr:.3g}; {len(k)} kept{capped}; "
              f"{int(np.sum(sel & (dist > 1.8 * step)))} kept peaks (beyond 1.8 grid steps from anything stronger)")

    w_lr = widgets.interactive(
        show_lr,
        min_lr=widgets.FloatLogSlider(value=min_lr, base=10, min=0, max=4, step=0.05, description="min_lr",
                                      continuous_update=False),
        keep=widgets.BoundedIntText(value=keep, min=1, max=10**8, description="keep"),
    )
    freeze(fig)
    display(w_lr, fig.canvas)
else:
    print("prune = completeness: no orientation fit")

if prune_used == "likelihood":
    min_lr_used, keep_used = w_lr.kwargs["min_lr"], w_lr.kwargs["keep"]
    above = np.flatnonzero(lr > min_lr_used)
    kept = pre[above[np.argsort(lr[above])[::-1]][:keep_used]]
else:
    min_lr_used, keep_used = min_lr, keep
    above = np.flatnonzero(comp >= cut)
    kept = above[np.argsort(comp[above])[::-1]][:keep_used]
U_kept = U_grid[kept]
print(f"{len(U_kept)} orientations kept for the occupancy fit")

7. Occupancy: --cand, --iter, --block-gb, --beam, --censor

Each voxel keeps its cand best orientations (by the first MLEM update from flat occupancy), and MLEM fits their occupancies, every voxel jointly. This is anri.index.fit_occupancy written out, so that the deviance can be recorded at every iteration and the map at a few. It costs the same as in python -m anri.index: the estimate below is for memory; time grows with voxels × candidates × reflections × iterations.

--block-gb only sets how many voxels are handled at once (memory, not results).

• beam: the beam’s FWHM across dty. Each voxel is spread over the dty rows by the beam’s profile (a Gaussian of that FWHM) integrated over the voxel at each ω, so a beam wider than the step (overfocusing) puts a voxel into more rows. With FWHM = step it is close to the default, a linear spread over the 2 nearest rows (beam = 0), and maps come out as good (AP1_1, by the figure of merit), at ~2.4x the cost of the candidate pass. The plot shows the weights per row for a voxel at a few places between two rows.

• censor: ImageD11 keeps only pixels above its segmentation cut, so an empty bin means every pixel in it was below the cut, not zero. Where the model predicts fewer than censor counts in an empty bin, the bin counts as agreeing with the model (in both orientation fits and the pruning’s likelihood ratio). That stops weak orientations, whose faint spots fell below the cut, from being driven to zero. Roughly the cut times the pixel-frames a faint spot covers; 0 turns it off.

if run_fit:
    pred = ix.predictions(U_kept, B, rings, geom, etacut_used)
    pos = np.asarray(anri.geom.recon_positions(NR, YSTEP), np.float32)  # every voxel
    pos_fit = pos[mask_v]  # the voxels fitted
    n_fit = len(pos_fit)
    scan = {"y0": Y0, "dty0": DTY0, "ystep": YSTEP, "n_rows": NK, "om0": OM0}
    if ddty is not None:  # fly, helical scans
        scan["ddty"] = ddty
    beam_used = beam
    bins_fit = bins
    if beam_used > 0:  # the beam's profile across dty, as python -m anri.index
        scan.update({"sig_beam": beam_used / (2 * np.sqrt(2 * np.log(2))), "width_beam": (gridstep - 1) * ystep0,
                     "voxel": YSTEP})
        bins_fit = (*bins, ix.beam_rows(scan))
    n_corners = 4 * bins_fit[4] if len(bins_fit) > 4 else 8
    ring_j = jnp.asarray(rings["ring_j"])
    n_j = pred[0].shape[1]
    vb_all = ix.block_voxels(16, n_j, block_gb * 1e9, n_corners)
    vb = ix.block_voxels(n_cand, n_j, block_gb * 1e9, n_corners)
    pos_p = ix.pad_voxels(pos_fit, max(vb_all, vb))
    print(f"{n_fit} voxels x {len(U_kept)} orientations x {n_j} reflections; occupancies {n_fit * n_cand * 8 / 1e9:.2f} "
          f"GB, blocks of {vb} voxels" + (f"; beam FWHM {beam_used:g}: {bins_fit[4]} rows per voxel" if beam_used > 0
                                          else "; the 2-row model"))
    if beam_used > 0:  # the weights per row for a voxel at a few places between two rows, mid-scan
        fig, ax = plt.subplots(figsize=(7, 3), layout="constrained")
        one = jnp.ones((1, 1, 1))
        om_mid = OM0 + n_o * bins[1] / 2
        for k_, off in enumerate((0.0, 0.25, 0.5)):
            # a voxel at lab y = x sin(omega) + y cos(omega) is in the beam at dty = y0 - lab y: off rows past the middle
            ylab = Y0 - DTY0 - (NK // 2 + off) * YSTEP
            p_ = jnp.asarray([[ylab * np.sin(np.radians(om_mid)), ylab * np.cos(np.radians(om_mid)), 0.0]], jnp.float32)
            for nb, ls in ((bins_fit[4], "-"), (0, ":")):
                idx_, wt_ = ix.system(one * 10.0, one * om_mid, one > 0, one, jnp.asarray([0]), p_, scan, *bins, nb)
                idx_, wt_ = np.asarray(idx_).ravel(), np.asarray(wt_).ravel()
                w_row = np.bincount(idx_[idx_ >= 0] % NK, weights=wt_[idx_ >= 0], minlength=NK)
                ax.plot(dty_c - (DTY0 + NK // 2 * YSTEP), w_row, ls, marker="o", ms=3, color=f"C{k_}",
                        label=f"voxel {off:g} step off a row" + ("" if nb else ", 2-row model"))
        ax.set(xlim=(-3 * YSTEP, 3 * YSTEP), xlabel="dty - dty of the middle row", ylabel="weight",
               title=f"weights per row, at omega {om_mid:.4g}")
        ax.legend(fontsize=7)
        plt.show()
    t0 = time.time()
    f0, cand = ix.candidates(H, pred, ring_j, pos_p, scan, bins_fit, n_cand, vb_all)
    print(f"candidates: {time.time() - t0:.0f} s")

if run_fit:  # anri.index.mlem, recording the deviance at every iteration and the map at some
    norm = jnp.maximum(ix.backward(jnp.ones(H.shape[0], H.dtype), cand, pred, ring_j, pos_p, scan, bins_fit, vb), 1e-30)
    fj = f0
    dev, moved, f_prev, it_prev = [], {}, f0, 0
    snap_at = sorted({i for i in (1, 2, 5, 10, 20, 50, 100, 200, 500) if i < n_iter} | {n_iter})
    t0 = time.time()
    for it in range(n_iter):
        Af = ix.forward(fj, cand, pred, ring_j, pos_p, scan, bins_fit, H.shape[0], vb)
        dev.append(float(ix.deviance(H, Af, censor)))
        fj = fj * ix.backward(ix.censored_ratio(H, Af, censor), cand, pred, ring_j, pos_p, scan, bins_fit, vb) / norm
        if it + 1 in snap_at:  # the share of each voxel's occupancy that moved since the last snapshot, per iteration
            pa = f_prev / jnp.maximum(f_prev.sum(1, keepdims=True), 1e-30)
            pb = fj / jnp.maximum(fj.sum(1, keepdims=True), 1e-30)
            moved[it + 1] = np.full(NV, np.nan)
            moved[it + 1][mask_v] = np.asarray(0.5 * jnp.abs(pb - pa).sum(1))[:n_fit] / (it + 1 - it_prev)
            f_prev, it_prev = fj, it + 1
            print(f"  iteration {it + 1}: deviance {dev[-1]:.5g} ({time.time() - t0:.0f} s)")
    model = ix.forward(fj, cand, pred, ring_j, pos_p, scan, bins_fit, H.shape[0], vb)
    f_v, cand_v = np.zeros((NV, n_cand), np.float32), np.zeros((NV, n_cand), np.int32)  # per voxel
    f_v[mask_v], cand_v[mask_v] = np.asarray(fj)[:n_fit], np.asarray(cand)[:n_fit]  # outside the mask: empty

Iterations. The deviance flattens early, but thin features (twins one or two voxels thick, small grains) keep changing long after. The middle plots show how much of each occupied voxel’s occupancy moved between orientations from one snapshot to the next (half the summed change of its occupancies, as a share of its total, per iteration), and where it still moves at the end. If that is still falling steeply at the last iteration, more iterations will change the map.

Candidates. Each voxel’s occupancy, sorted, as a share of the voxel’s total: if the last candidates still hold a noticeable share in many voxels, the voxel wanted more than cand orientations.

if run_fit:
    tot = f_v.sum(1)
    occ = mask_v & (tot > occupied_cut * np.percentile(tot[mask_v], 99))

    its = list(moved)[1:]  # the first snapshot is from the candidates' starting point
    fig, ax = plt.subplots(1, 4, figsize=(18, 4), layout="constrained")
    ax[0].plot(np.arange(1, n_iter + 1), dev)
    ax[0].set(xscale="log", yscale="log", xlabel="iteration", ylabel="deviance")
    ax0 = ax[0].twinx()
    ax0.semilogy(np.arange(2, n_iter + 1), np.maximum(-np.diff(dev) / np.array(dev[1:]), 1e-9), color="tab:red", lw=0.8)
    ax0.set_ylabel("relative decrease per iteration", color="tab:red")
    if its:
        for q, ls in ((50, "-"), (90, "--"), (99, ":")):
            ax[1].semilogy(its, [np.percentile(moved[i][occ], q) for i in its], "o" + ls, color="k", ms=3,
                           label=f"{q}th percentile")
        ax[1].set(xscale="log", xticks=its, xticklabels=[str(i) for i in its], xlabel="iteration",
                  ylabel="share moved per iteration", title="occupancy moving between orientations")
        ax[1].minorticks_off()
        ax[1].legend(fontsize=8)
        im = ax[2].imshow(np.where(occ, moved[its[-1]], np.nan).reshape(NR, NR), origin="lower", cmap="Reds",
                          norm=LogNorm(1e-5, 1e-1))
        ax[2].set(title=f"share moved per iteration, at {its[-1]}", xticks=[], yticks=[])
        fig.colorbar(im, ax=ax[2], shrink=0.8)
    iv = np.flatnonzero(occ)
    iv = iv[:: max(1, len(iv) // 200000)]  # a sample of the occupied voxels is enough
    share = np.sort(f_v[iv], 1)[:, ::-1] / np.maximum(tot[iv, None], 1e-30)
    for q, ls in ((50, "-"), (90, "--"), (99, ":")):
        ax[3].semilogy(np.arange(1, n_cand + 1), np.maximum(np.percentile(share, q, 0), 1e-6), ls, color="k",
                       label=f"{q}th percentile")
    ax[3].set(xlabel="candidate rank", ylabel="share of the voxel's occupancy", title=f"cand = {n_cand}")
    ax[3].legend(fontsize=8)
    plt.show()
    tail = share[:, -max(1, n_cand // 8) :].sum(1)
    print(f"the last {max(1, n_cand // 8)} candidates hold {np.median(tail) * 100:.2f}% of an occupied voxel's "
          f"occupancy (median), {np.percentile(tail, 99) * 100:.2f}% (99th percentile)")

8. Model against data

The fitted histogram A f against the data.

• Per dty row (as python -m anri.index logs it): measured / fitted intensity. Steps from row to row point to the flux varying between the rows’ scans (--monitor); a smooth trend with radius, to the model.

• Per ring and spot by spot: the data, the model and the difference for one ring with every row summed (η against ω), and the sinogram of an η window. Compact positive spots in the difference are intensity that no kept orientation explains: a missed grain (lower min_lr, lit or min_comp, or a finer grid) or another phase. Negative spots are predicted where there is nothing: decoys, or wrong orientations. A smooth residual all along a ring is background, which the model does not have.

Only η bins the fit models are used in the ratios (|sin η| above etacut).

if run_fit:
    M4 = np.asarray(model).reshape(n_rings_used, n_e, n_o, NK)
    eta_c = -180.0 + (np.arange(n_e) + 0.5) * bins[0]
    use_e = np.abs(np.sin(np.radians(eta_c))) > etacut_used
    d_row, m_row = H4[:, use_e].sum((0, 1, 2)), M4[:, use_e].sum((0, 1, 2))
    lit_rows = m_row > 0.05 * m_row.max()
    row_ratio = np.where(lit_rows, d_row / np.maximum(m_row, 1e-30), np.nan)
    d_ring, m_ring = H4[:, use_e].sum((1, 2, 3)), M4[:, use_e].sum((1, 2, 3))
    fig, ax = plt.subplots(1, 2, figsize=(13, 3.6), layout="constrained", width_ratios=[2, 1])
    ax[0].plot(dty_c, row_ratio, ".-", ms=3)
    ax[0].axhline(1, color="k", lw=0.5)
    ax[0].set(xlabel="dty", ylabel="measured / fitted", title="per dty row")
    ax[1].bar(np.arange(n_rings_used), m_ring / d_ring)
    ax[1].set(xlabel="ring", ylabel="fitted / measured", title="per ring")
    plt.show()
    dv = row_ratio[lit_rows] - 1
    print(f"per row: rms {np.sqrt(np.mean(dv**2)):.3f}, row-to-row rms "
          f"{np.sqrt(np.nanmean(np.diff(row_ratio) ** 2)):.3f}")

    # data, model and residual: one ring with the rows summed (eta against omega, top), and the sinogram of an eta
    # window (bottom); every panel shares omega, each row its other axis
    fig = figure("residual", figsize=(16, 8))
    ax = fig.subplots(2, 3, sharex=True)
    for row in range(2):
        for j in (1, 2):
            ax[row, j].sharey(ax[row, 0])
    om_ext = (OM0, OM0 + n_o * bins[1])
    exts = ((*om_ext, -180, 180), (*om_ext, DTY0 - YSTEP / 2, DTY0 + (NK - 0.5) * YSTEP))
    shapes = ((n_e, n_o), (NK, n_o))
    ims = []
    for row in range(2):
        z = np.zeros(shapes[row])
        ims.append([ax[row, j].imshow(z, aspect="auto", origin="lower", extent=exts[row], interpolation="nearest",
                                      cmap=MAGMA if j < 2 else "RdBu_r") for j in range(3)])
        fig.colorbar(ims[row][2], ax=ax[row, 2], shrink=0.8)
        ax[row, 0].set_ylabel(("η (°)", "dty")[row])
        for j in range(3):
            ax[row, j].set_xlabel("ω (°)")
    eta_lines = [[a.axhline(0.0, color="c", lw=0.6) for _ in range(2)] for a in ax[0]]

    def show_residual(ring, eta):
        ie = (eta_c >= eta[0]) & (eta_c < eta[1])
        pairs = ((H4[ring].sum(2), M4[ring].sum(2)), (H4[ring, ie].sum(0).T, M4[ring, ie].sum(0).T))
        for row, (dd, mm) in enumerate(pairs):
            top = max(dd.max(), 1.0)
            for j, im in enumerate((dd, mm)):
                ims[row][j].set_data(np.ma.masked_less_equal(im, 0))
                ims[row][j].set_norm(LogNorm(top * 1e-2, top))
            r = dd - mm
            v = max(np.percentile(np.abs(r), 99.9), 1e-9)
            ims[row][2].set_data(r)
            ims[row][2].set_clim(-v, v)
            for j, t in enumerate(("data", "model", "data − model")):
                ax[row, j].set_title(f"ring {ring}, {t}" + ("" if row == 0 else f", η {eta[0]}..{eta[1]}°"))
        for pair in eta_lines:
            for ln, e_ in zip(pair, eta):
                ln.set_ydata([e_, e_])
        fig.canvas.draw_idle()

    w_res = widgets.interactive(
        show_residual,
        ring=widgets.IntSlider(value=0, min=0, max=n_rings_used - 1, description="ring", continuous_update=False),
        eta=widgets.IntRangeSlider(value=(80, 100), min=-180, max=180, step=1, description="η window",
                                   continuous_update=False),
    )
    freeze(fig)
    display(w_res, fig.canvas)

9. Populations: --occupied, --min-frac

Each voxel’s occupied candidates are grouped into populations (candidates within 1.8 grid steps of the most occupied one join it, and so on), each with a fraction of the voxel, a mean orientation and a spread (which includes the grid’s own spacing, so it is an upper bound).

• occupied: voxels in the mask (section 3) count as sample above this × the 99th percentile of the total occupancy in the mask; voxels outside it are never occupied. The histogram of the total occupancy should show the sample well apart from the air around it; put the cut in the gap. (The raw occupancy is always saved.)

• min_frac: populations holding at least this fraction of a voxel are reported. Second populations are real where two grains (or a grain and its twin) share a voxel, at boundaries and in lamellae; a scatter of weak ones all over the grains is noise from the coarse model (the beam is wider than a voxel, the grid coarser than the truth). Put the cut above that floor.

if run_fit:
    t0 = time.time()
    frac, U_pop, spread, n_members = ix.populations(f_v, cand_v, U_kept, ops, 1.8 * step, p=N_POP)
    print(f"populations: {time.time() - t0:.0f} s")
    p99 = np.percentile(tot[mask_v], 99)  # within the mask
    f2 = frac[mask_v & (tot > occupied_cut * p99), 1:].ravel()  # in the voxels occupied at the parameters' cut

    fig = figure("populations", figsize=(15, 9))
    ax = fig.subplots(2, 3)
    maps = (ax[0, 0], ax[0, 2], ax[1, 0], ax[1, 2])  # they share their axes; the histograms do not
    for a in maps[1:]:
        a.sharex(maps[0])
        a.sharey(maps[0])
    im = ax[0, 0].imshow(tot.reshape(NR, NR), origin="lower", cmap="viridis")
    fig.colorbar(im, ax=ax[0, 0], shrink=0.8)
    im_out = ax[0, 0].imshow(np.zeros((NR, NR, 4)), origin="lower")  # the voxels not occupied, greyed out
    ax[0, 0].set(title="total occupancy; grey: not occupied")
    ax[0, 1].hist(tot / p99, bins=np.linspace(0, max(1.5, (tot / p99).max()), 150), log=True)
    v_occ = ax[0, 1].axvline(occupied_cut, color="tab:red")
    ax[0, 1].set(xlabel="total occupancy / 99th percentile", ylabel="voxels")
    im_n = ax[0, 2].imshow(np.zeros((NR, NR)), origin="lower", cmap="Blues", vmin=0, vmax=N_POP)
    ax[0, 2].set(title="populations per voxel")
    fig.colorbar(im_n, ax=ax[0, 2], shrink=0.8)
    im_f2 = ax[1, 0].imshow(np.zeros((NR, NR)), origin="lower", cmap="magma", vmin=0, vmax=0.5)
    ax[1, 0].set(title="fraction of the second population")
    fig.colorbar(im_f2, ax=ax[1, 0], shrink=0.8)
    ax[1, 1].hist(f2[f2 > 0], bins=np.linspace(0.0, 0.5, 101), log=True)
    v_frac = ax[1, 1].axvline(min_frac, color="tab:red")
    ax[1, 1].set(xlabel="fraction of populations 2 to 4", ylabel="populations")
    im_sp = ax[1, 2].imshow(np.zeros((NR, NR)), origin="lower", cmap="viridis", vmin=0, vmax=max(2 * step, 1e-3))
    ax[1, 2].set(title="spread of the main population (°)")
    fig.colorbar(im_sp, ax=ax[1, 2], shrink=0.8)
    for a in maps:
        a.set(xticks=[], yticks=[])

    def show_pop(occupied_cut, min_frac):
        occ = mask_v & (tot > occupied_cut * p99)
        present = (frac >= min_frac) & occ[:, None]
        present[:, 0] = occ
        n_pop = present.sum(1)
        grey = np.zeros((NV, 4))
        grey[~occ] = (0.85, 0.85, 0.85, 0.85)
        im_out.set_data(grey.reshape(NR, NR, 4))
        im_n.set_data(np.where(occ, n_pop, np.nan).reshape(NR, NR))
        im_f2.set_data(np.where(occ, np.where(present[:, 1], frac[:, 1], 0), np.nan).reshape(NR, NR))
        im_sp.set_data(np.where(occ, spread[:, 0], np.nan).reshape(NR, NR))
        v_occ.set_xdata([occupied_cut, occupied_cut])
        v_frac.set_xdata([min_frac, min_frac])
        fig.canvas.draw_idle()
        print(f"{occ.sum()} of {NV} voxels occupied; populations per occupied voxel: "
              + ", ".join(f"{k}: {np.mean(n_pop[occ] == k) * 100:.1f}%" for k in range(1, N_POP + 1)))

    w_pop = widgets.interactive(
        show_pop,
        occupied_cut=widgets.FloatSlider(value=occupied_cut, min=0.0, max=1.0, step=0.01, description="occupied",
                                         continuous_update=False),
        min_frac=widgets.FloatSlider(value=min_frac, min=0.0, max=0.5, step=0.01, description="min_frac",
                                     continuous_update=False),
    )
    freeze(fig)
    display(w_pop, fig.canvas)

The main population’s IPF colours (with ImageD11 and orix installed) and its completeness: the fraction of its predicted reflections that land on lit bins. Low completeness inside the sample points to a wrong orientation, or a voxel the data barely see.

if run_fit:
    occupied_used, min_frac_used = w_pop.kwargs["occupied_cut"], w_pop.kwargs["min_frac"]
    occ = mask_v & (tot > occupied_used * p99)
    comp_main = np.full(NV, np.nan)
    comp_main[occ] = ix.completeness_of(U_pop[occ, 0], delta, B, rings, geom, lit)
    fig, ax = plt.subplots(1, 2, figsize=(11, 5), layout="constrained", sharex=True, sharey=True)
    im = ax[1].imshow(comp_main.reshape(NR, NR), origin="lower", cmap="viridis", vmin=0, vmax=1)
    ax[1].set(title="completeness of the main population", xticks=[], yticks=[])
    fig.colorbar(im, ax=ax[1], shrink=0.8)
    try:
        maps = {"UBI": np.where(occ[:, None, None], np.linalg.inv(U_pop[:, 0] @ B), np.nan).reshape(NR, NR, 3, 3),
                "phase_ids": np.where(occ, 0, -1).reshape(NR, NR)}
        tmap = anri.io.tensormap_from_recon(maps, lpars, sg, phase_name, YSTEP)
        tmap.get_ipf_maps()
        ax[0].imshow(np.where((tmap.phase_ids[0] >= 0)[..., None], tmap.ipf_z[0], 1.0), origin="lower")
        ax[0].set(title="main population, IPF-Z")
    except ImportError as e:
        ax[0].set(title=f"no IPF map ({e.name} is not installed)")
    ax[0].set(xticks=[], yticks=[])
    plt.show()

The command line, and the results

Every value chosen above (sliders included), as options of python -m anri.index: the same command indexes another dataset (or this one again) with them. Options left to the data (--tth-tol, --min-comp, --grid when automatic) are only given if set; --y0 only if it differs from the DataSet’s.

Then the results of this notebook’s run are written as python -m anri.index writes them (with the same function), in outdir: <sample>_<dataset>_index_params.toml (the command above, every option and the values resolved from the data), _index.npz (occupancies and populations), _index_entries.npz (every population as map entries, for the renderer) and, with ImageD11, _index_tmap.h5 (a TensorMap of the main population, and .xdmf for ParaView). With the phantom, outdir is the phantom’s temporary folder.

# every value chosen above (sliders included), named as python -m anri.index names its options
if not run_fit:
    occupied_used, min_frac_used = occupied_cut, min_frac
options = {
    "analysisroot": analysisroot, "sample": sample, "dataset": dataset, "phase": phase, "parfile": parfile,
    "rings": n_rings_used, "grid": grid, "max_chance": w_grid.kwargs["max_chance"], "keep": keep_used,
    "min_comp": min_comp_used, "prune": prune_used, "min_lr": min_lr_used, "cif": cif, "lit": lit_used,
    "etacut": etacut_used, "tth_tol": tth_tol_used, "iter": n_iter, "cand": n_cand, "coarse": coarse,
    "monitor": monitor, "mask": mask_used, "beam": beam or None, "censor": censor or None, "occupied": occupied_used, "min_frac": min_frac_used, "block_gb": block_gb,
    "y0": None if ds["y0"] is not None and np.isclose(Y0, ds["y0"]) else Y0, "gridstep": gridstep,
    "outdir": outdir, "n_cpu": n_cpu, "check": False,
}
cmd = ["python", "-m", "anri.index", analysisroot, sample, dataset]
for k, v in options.items():
    if v is None or v is False or k in ("analysisroot", "sample", "dataset"):
        continue
    cmd += ["--" + k.replace("_", "-")] + ([] if v is True else [f"{v:.10g}" if isinstance(v, float) else str(v)])
command = shlex.join(cmd)
print(command)

if run_fit:  # the results, written as python -m anri.index writes them
    import datetime

    from anri.index.__main__ import _git_commit, save_results, write_toml

    occ = mask_v & (tot > occupied_used * p99)
    present = (frac >= min_frac_used) & occ[:, None]
    present[:, 0] = occ  # the main population always
    comp_pop = np.zeros(frac.shape, np.float32)
    comp_pop[present] = ix.completeness_of(U_pop[present], delta, B, rings, geom, lit)
    tag = os.path.join(outdir, f"{dsname}_index")
    results = {"f": f_v, "cand": cand_v, "U": U_kept, "comp": comp[kept], "frac": frac, "U_pop": U_pop,
               "spread": spread, "n": n_members, "comp_pop": comp_pop, "occupied": occ, "present": present,
               "pos": pos, "grid_step": step, "delta": delta, "row_ratio": row_ratio, "row_data": d_row,
               "row_model": m_row, "y0": Y0, "mask": mask_v}
    save_results(tag, results, B, lpars, sg, phase_name, YSTEP)
    run_info = {
        "run": {"command": command, "notebook": "indexing_parameters.ipynb", "anri_version": anri.VERSION,
                **_git_commit(), "cwd": os.getcwd(),
                "date": datetime.datetime.now(datetime.timezone.utc).astimezone().isoformat(timespec="seconds")},
        "options": options,
        "resolved": {"dataset": dsfile, "sparsefile": sparsefile, "parfile": pfile, "phase": phase_name,
                     "lattice": list(lpars), "space_group": sg, "wavelength": WL, "y0": Y0, "voxel_size": YSTEP,
                     "voxels": NR, "omega_step": OSTEP, "grid_step": step, "grid_auto": grid is None,
                     "grid_worst_case_deg": delta, "rings_tth_deg": list(rings["tth"]), "tth_tol_deg": list(tth_tol_r),
                     "n_hkls": len(rings["hkls"]), "spatial_correction": spatial_src, "beam_fwhm": beam_used},
        "results": {"chance_completeness": info["chance"], "min_comp_used": cut, "orientations_kept": len(U_kept),
                    "voxels_occupied": int(occ.sum()), "seconds": time.time() - T_START},
    }
    write_toml(f"{tag}_params.toml", run_info)
    print(f"-> {tag}_params.toml")
else:
    print("run_fit = False: nothing indexed here, so nothing written; run the command above")