{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# Choosing the indexing parameters\n\n`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`).\n\n| Section | Options | What the plots show |\n|---|---|---|\n| 1. Rings | `--rings`, `--tth-tol`, `--cif` | the 2θ profile and the (2θ, η) cake: which rings are clean, how wide they are, what lies between them |\n| 2. Data | `--monitor`, `--gridstep` | intensity per dty row and per ω: flux steps between the rows' scans |\n| 3. 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) |\n| 4. Lit map | `--lit`, `--etacut` | how much of the map is lit against how much intensity it holds; what the η cut removes |\n| 5. Grid | `--grid`, `--max-chance` | chance completeness and size of every grid step, and how much each finer step improves the fit |\n| 6. Pruning | `--prune`, `--min-comp`, `--min-lr`, `--keep` | completeness and likelihood ratio of every orientation, and how many each cut keeps |\n| 7. Occupancy | `--cand`, `--iter`, `--block-gb`, `--beam`, `--censor` | the deviance per iteration, how many voxels still change orientation, and whether the candidates suffice |\n| 8. Model against data | | what the fit leaves unexplained, per row, per ring and spot by spot |\n| 9. Populations | `--occupied`, `--min-frac` | the occupied voxels within the mask, and the populations per voxel |\n\nWhere 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.\n\nThe 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).\n\nWith no dataset given, it renders the 316L [phantom](phantom.ipynb) as in the [indexing tutorial](indexing.ipynb) (about a minute on a laptop CPU), so it runs out of the box.",
   "id": "cell-00"
  },
  {
   "cell_type": "code",
   "metadata": {
    "tags": [
     "parameters"
    ]
   },
   "execution_count": null,
   "outputs": [],
   "source": "# Parameters: the options of python -m anri.index, at its defaults. The sliders below start from these values.\n\n# the dataset, as python -m anri.index <analysisroot> <sample> <dataset>; None renders a phantom to try things on\nanalysisroot = None\nsample = None\ndataset = None\nphase = None  # --phase: name in pars.json (None: the only one)\nparfile = None  # --parfile (None: the DataSet's, else pars/pars.json beside PROCESSED_DATA)\ncif = None  # --cif: CIF of the phase, for structure factors (None: |F|^2 = 1)\nmonitor = None  # --monitor: counter to normalise intensities by, e.g. \"fpico6\" (None: none)\ngridstep = 1  # --gridstep: voxel = gridstep x dty step, rows summed to match\n\n# 1. rings\nn_rings = 6  # --rings\ntth_tol = None  # --tth-tol, deg (None: measured per ring)\n# 3. rotation axis and sample mask\ny0 = None  # --y0 (None: the DataSet's)\nmask = None  # --mask: \"auto\" (threshold of the back-projection), a .npy file, or None (every voxel)\n# 4. lit map\nlit_cut = 1.0  # --lit: x the median non-empty bin\netacut = 0.2  # --etacut\n# 5. grid\ngrid = None  # --grid, deg (None: the coarsest with chance completeness <= max_chance)\nmax_chance = 0.5  # --max-chance\n# 6. pruning\nprune = \"likelihood\"  # --prune: \"likelihood\" or \"completeness\"\nmin_comp = None  # --min-comp (None: the chance level for likelihood, halfway to the maximum for completeness)\nmin_lr = 25.0  # --min-lr\nkeep = 100000  # --keep\n# 7. occupancy\nn_cand = 64  # --cand\nn_iter = 50  # --iter\ncoarse = 1  # --coarse: only passed on to the command line; this notebook fits at full resolution\nblock_gb = 1.0  # --block-gb\nbeam = 0.0  # --beam: FWHM of the beam across dty (0: the 2-row model, close to FWHM = dty step)\ncensor = 0.0  # --censor: counts per bin below which an empty bin agrees with the model (0: off)\n# 9. populations\noccupied_cut = 0.2  # --occupied\nmin_frac = 0.1  # --min-frac\n\nn_cpu = 4  # --n-cpu\noutdir = \".\"  # --outdir, for the command line at the end\n\n# this notebook only\nmax_rings_shown = 10  # rings drawn on the 2theta profile, to choose n_rings from\nrun_fit = True  # False: stop after pruning (sections 7-9 run the full occupancy fit)\ndraw = False  # draw the sample mask by hand in section 3 (needs %matplotlib widget)\ngrid_scan = True  # fit the orientations of every grid step in section 5, to compare their deviance",
   "id": "cell-01"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "import anri.utils\n\nanri.utils.setup(n_cpu=n_cpu)  # before any JAX computation\n\n%matplotlib widget\n\nimport os\nimport shlex\nimport tempfile\nimport time\nfrom functools import partial\n\nimport h5py\nimport ipywidgets as widgets\nimport jax\nimport jax.numpy as jnp\nimport numpy as np\nfrom IPython.display import display\nfrom matplotlib import pyplot as plt\nfrom matplotlib.colors import LogNorm\n\nimport anri.crystal\nimport anri.geom\nimport anri.index as ix\nimport anri.io\n\n# 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\n# histogram in bins R_E x R_O times larger (1 x 1 deg), per dty row\nB_E, B_O = 0.5, 0.25\nR_E, R_O = 2, 4\nN_POP = 4\nMAGMA = plt.colormaps[\"magma\"].with_extremes(bad=\"black\")  # log-scale maps: empty (masked) bins black, not white\nT_START = time.time()\n\n\ndef figure(name, **kwargs):\n    # a figure shown where display(fig.canvas) puts it (under its sliders), not at the end of the cell; re-running the\n    # cell replaces it\n    plt.close(name)\n    with plt.ioff():\n        fig = plt.figure(name, layout=\"constrained\", **kwargs)\n    fig.canvas.header_visible = False\n    return fig\n\n\ndef freeze(fig):\n    # fix every axis's limits as first drawn: the sliders then update the plots without undoing a zoom\n    for a in fig.axes:\n        if getattr(a, \"_colorbar\", None) is None:  # colorbars follow their images' limits\n            a.set_xlim(a.get_xlim())\n            a.set_ylim(a.get_ylim())",
   "id": "cell-02"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## The dataset\n\nWithout a dataset, the 316L phantom of the [indexing tutorial](indexing.ipynb) is rendered into a temporary folder as an ImageD11 dataset (this needs ImageD11). Its rotation axis is at dty = 0.3 µm.",
   "id": "cell-03"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "if analysisroot is None:\n    from ImageD11.sinograms.tensor_map import TensorMap\n\n    truth = TensorMap.from_h5(os.path.join(\"..\", \"..\", \"..\", \"tests\", \"data\", \"phantoms\", \"am316l\", \"am316l_tmap.h5\"))\n    a = truth.phases[0].lattice_parameters[0]\n    wl = 0.2843\n    pars = {\n        \"y_center\": 1023.5, \"y_size\": 75.0, \"tilt_y\": 1e-3, \"z_center\": 1023.5, \"z_size\": 75.0, \"tilt_z\": -2e-3,\n        \"tilt_x\": 0.0, \"distance\": 150e3, \"o11\": -1, \"o12\": 0, \"o21\": 0, \"o22\": -1, \"wavelength\": wl, \"wedge\": 0.0,\n        \"chi\": 0.0, \"omegasign\": 1.0, \"t_x\": 0.0, \"t_y\": 0.0, \"t_z\": 0.0,\n    }\n    y0_true = 0.3\n    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,\n                                    sig_psf=0.5)\n    entries = anri.io.entries_from_tensormap(truth)\n    entries[\"density\"] = np.full(len(entries[\"pos\"]), 30.0)  # counts scale\n    rings8 = ix.ring_table(np.array([a, a, a, 90.0, 90.0, 90.0]), 225, wl, 8)\n    omega_r, dty_r = anri.io.motor_grid((0.0, 180.0), 0.1, (y0_true - 30.0, y0_true + 30.0), 1.0)\n    analysisroot, sample, dataset = tempfile.mkdtemp(), \"phantom\", \"am316l\"\n    sparse = os.path.join(analysisroot, \"phantom_sparse.h5\")\n    t0 = time.time()\n    anri.io.simulate_sparse(sparse, entries, rings8[\"hkls\"], np.ones(len(rings8[\"hkls\"])), geom_r, omega_r, dty_r,\n                            (2048, 2048), max_frames=61)\n    cell_r = {\"cell__a\": a, \"cell__b\": a, \"cell__c\": a, \"cell_alpha\": 90.0, \"cell_beta\": 90.0, \"cell_gamma\": 90.0,\n              \"cell_lattice_[P,A,B,C,I,F,R]\": 225}\n    pars_r = anri.io.write_pars(os.path.join(analysisroot, \"pars\"), pars, {\"316L\": cell_r})\n    anri.io.write_dataset(sparse, analysisroot, sample, dataset, y0=y0_true, parfile=pars_r)\n    outdir = analysisroot  # the results go beside the phantom, not into the current folder\n    print(f\"phantom rendered in {time.time() - t0:.0f} s into {analysisroot}\")",
   "id": "cell-04"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# as python -m anri.index: the DataSet, its sparse pixels, and the phase from pars.json\ndsname = f\"{sample}_{dataset}\"\ndsfile = os.path.join(analysisroot, sample, dsname, f\"{dsname}_dataset.h5\")\nsparsefile = os.path.join(analysisroot, sample, dsname, f\"{dsname}_sparse.h5\")\nds = anri.io.read_dataset(dsfile)\nif ds[\"sparsefile\"] and os.path.exists(ds[\"sparsefile\"]):\n    sparsefile = ds[\"sparsefile\"]\nY0_DS = ds[\"y0\"] if y0 is None else y0\nif Y0_DS is None:\n    raise ValueError(f\"{dsfile} has no y0: set y0 in the parameters (section 3 then helps to check it)\")\npfile = parfile or ds[\"parfile\"]\nif not pfile or not os.path.exists(pfile):  # e.g. processed elsewhere: pars/ beside PROCESSED_DATA\n    root = os.path.abspath(analysisroot)\n    while os.path.basename(root) != \"PROCESSED_DATA\" and root != os.path.dirname(root):\n        root = os.path.dirname(root)\n    pfile = os.path.join(os.path.dirname(root), \"pars\", \"pars.json\")\ngeo, phase_name, cell = anri.io.read_pars_json(pfile, phase)\nlpars = np.array([cell[k] for k in (\"cell__a\", \"cell__b\", \"cell__c\", \"cell_alpha\", \"cell_beta\", \"cell_gamma\")])\nsg = int(cell[\"cell_lattice_[P,A,B,C,I,F,R]\"])  # a space-group number is needed, not a centring letter\n\nybin, yedge, oedge = ds[\"ybincens\"], ds[\"ybinedges\"], ds[\"obinedges\"]\nystep0, nk0 = float(np.median(np.diff(ybin))), len(ybin)\nYSTEP, DTY0, NK = gridstep * ystep0, float(ybin[0]) + 0.5 * (gridstep - 1) * ystep0, -(-nk0 // gridstep)\nOM0, OSTEP = float(oedge[0]), float(np.median(np.diff(oedge)))\nN_E, N_O = round(360 / B_E), round(float(oedge[-1] - oedge[0]) / B_O)\nN_O -= N_O % R_O\nbins = (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\nn_e, n_o = bins[2], bins[3]\ndty_c = DTY0 + np.arange(NK) * YSTEP  # row centres\nom_c = OM0 + (np.arange(n_o) + 0.5) * bins[1]  # omega bin centres of the fit's histogram\n\nWL = geo[\"wavelength\"]\ngeom = 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)\ngeom = {k: jnp.asarray(v, jnp.float32) if np.issubdtype(np.asarray(v).dtype, np.floating) else v\n        for k, v in geom.items()}\nB64 = anri.crystal.B_matrix(lpars)\nB = B64.astype(np.float32)\nops = anri.crystal.laue_rotations(anri.crystal.symmetry_matrices(sg), B64)\nif cif is None:\n    structure = None\nelse:\n    import Dans_Diffraction\n\n    structure = Dans_Diffraction.Crystal(cif)\n\nwith h5py.File(sparsefile, \"r\") as h:\n    groups = list(h.keys())\n    n_pix = {g: int(h[g][\"nnz\"][()].sum()) for g in groups}\nchunk = int(min(1 << 24, 1 << max(10, int(np.ceil(np.log2(max(max(n_pix.values()), 1)))))))\nspatial = anri.io.read_spatial(ds)  # the DataSet's detector distortion maps (e2dx/e2dy or detector file), if any\nspatial_src = ds[\"detectorh5\"] or (f\"{ds['e2dxfile']}, {ds['e2dyfile']}\" if spatial is not None else \"\")\nmonitor_ref = None\nif monitor:  # one reference for every scan, the mean, as python -m anri.index\n    mons = anri.io.read_monitor(sparsefile, groups, monitor, ds[\"masterfile\"])\n    monitor_ref = float(np.mean(np.concatenate(list(mons.values()))))\n\n\ndef stream(groups_):\n    # sparse pixels of some scans, a chunk at a time, read ahead in a thread, as python -m anri.index\n    return anri.io.prefetch(anri.io.stream_sparse(sparsefile, yedge, ds[\"omegamotor\"], ds[\"dtymotor\"], chunk, groups_,\n                                                  gridstep, ds[\"dty\"], ds[\"scans\"], monitor, monitor_ref,\n                                                  ds[\"masterfile\"], ds[\"omega\"], spatial))\n\n\n\n\n# the sinogram's rows, as the DataSet has binned the frames (ds.dty, ds.omega: [rows, frames]); stream_sparse takes\n# each frame's dty and omega from there. How far dty moves within a row, in dty steps:\ndty_rows = None if ds[\"dty\"] is None else np.asarray(ds[\"dty\"], float)\ndrift = 0.0 if dty_rows is None else float(np.ptp(dty_rows, axis=1).max()) / ystep0\nprint(f\"dataset {dsfile}\\nsparse pixels {sparsefile}: {len(groups)} scans, {sum(n_pix.values()) / 1e6:.1f}M pixels\")\nadvice = (\": its frames alias between neighbouring rows (gaps and doubled rows in the sinogram); replace ds.dty by \"\n          \"its mean over each row in ImageD11, and save the DataSet\")\nprint(f\"DataSet: {len(ds['scans'] or [])} rows ({', '.join((ds['scans'] or [])[:2])}, ...); its dty moves by up to \"\n      f\"{drift:.2f} dty steps within a row\" + (advice if drift > 0.2 else \"\"))\n# where the beam really was: each row's DataSet dty plus dty's motion within the row, read from the sparse file (fly\n# and helical scans, where dty moves while omega turns); the model follows it (scan[\"ddty\"]), as python -m anri.index\nddty = None\nif dty_rows is not None and ds[\"scans\"] is not None and ds[\"omega\"] is not None:\n    dd_ = ix.dty_offsets(anri.io.read_frame_dty(ds, sparsefile), ds[\"omega\"], yedge, gridstep, OM0, bins[1], n_o, DTY0,\n                         YSTEP)\n    if np.abs(dd_).max() > 1e-3 * YSTEP:\n        ddty = jnp.asarray(dd_, jnp.float32)\nprint(\"dty within the rows: \" + (\"constant (a step scan)\" if ddty is None else\n      f\"moves up to {float(jnp.abs(ddty).max()) / YSTEP:.2f} dty steps from each row's dty; the model follows it\"))\nprint(f\"spatial correction: {spatial_src or 'none (the DataSet names no e2dx/e2dy or detector file)'}\")\nprint(f\"pars {pfile}: phase {phase_name}, lattice {', '.join(f'{v:g}' for v in lpars)}, space group {sg}, \"\n      f\"{len(ops)} Laue-group rotations\")\nprint(f\"wavelength {WL:.5f}, distance {geo['distance']:g}; y0 {Y0_DS:.6g}; dty {DTY0:.6g} + {NK} x {YSTEP:.6g}; \"\n      f\"omega {OM0:.4g} .. {oedge[-1]:.4g} in {len(oedge) - 1} frames of {OSTEP:.4g}\")",
   "id": "cell-05"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 1. Rings: `--rings`, `--tth-tol`, `--cif`\n\nPixels 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.\n\nWhat to look for:\n\n- **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.\n- **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.\n- **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.\n- **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.\n- **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.",
   "id": "cell-06"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "rings_all = ix.ring_table(lpars, sg, WL, max_rings_shown, structure)\nn_all = len(rings_all[\"tth\"])\nsample_groups = [groups[i] for i in np.unique(np.linspace(0, len(groups) - 1, min(len(groups), 9)).round().astype(int))]\n\n\n@partial(jax.jit, static_argnames=(\"n_t\", \"n_eta\"))\ndef cake(x, val, lo, step, n_t, n_eta):\n    # sum pixel intensities in (2theta, eta) bins: n_t of width step from lo, and n_eta over 360 deg from -180\n    it = jnp.floor((x[:, 0] - lo) / step).astype(jnp.int32)\n    ie = jnp.floor((x[:, 1] + 180.0) * n_eta / 360.0).astype(jnp.int32) % n_eta\n    ok = (it >= 0) & (it < n_t)\n    idx = jnp.where(ok, it * n_eta + ie, n_t * n_eta)\n    return jax.ops.segment_sum(jnp.where(ok, val, 0.0), idx, n_t * n_eta + 1)[:-1].reshape(n_t, n_eta)\n\n\nTTH_STEP, N_CAKE = 0.002, 180  # 2theta bins as ix.ring_profile; eta in 2 deg bins\ntth_lo = float(rings_all[\"tth\"][0]) - 0.5\nn_t = int(np.ceil((float(rings_all[\"tth\"][-1]) + 0.5 - tth_lo) / TTH_STEP))\nC = jnp.zeros((n_t, N_CAKE), jnp.float32)\nt0 = time.time()\nfor slow, fast, om, _, val in stream(sample_groups):\n    m = len(val)\n\n    def pad(a, m=m):\n        return jnp.asarray(np.pad(np.asarray(a, np.float32), (0, chunk - m)))\n\n    C = C + cake(ix.pixel_angles(pad(slow), pad(fast), pad(om), geom), pad(val), tth_lo, TTH_STEP, n_t, N_CAKE)\nC = np.asarray(C)\nprof = C.sum(1)  # the 2theta profile, as ix.ring_profile measures it\ntth_x = tth_lo + (np.arange(n_t) + 0.5) * TTH_STEP\nprint(f\"profile of {len(sample_groups)} scans: {time.time() - t0:.0f} s\")\n\n# each ring's reflections: a representative hkl, the multiplicity and the mean |F|^2\nring_of_h = rings_all[\"ring_j\"][::2]\nring_info = []\nfor r in range(n_all):\n    h = rings_all[\"hkls\"][ring_of_h == r]\n    hkl = tuple(int(v) for v in sorted(np.abs(h[0]), reverse=True))\n    ring_info.append((hkl, len(h), float(rings_all[\"F2\"][ring_of_h == r].mean())))\ngap = np.diff(rings_all[\"tth\"])\n\nfig = figure(\"rings\", figsize=(14, 9))\nsf = fig.subfigures(2, 1, height_ratios=[2.3, 0.9])\nax_p, ax_c = sf[0].subplots(2, 1, sharex=True)\nax_p.semilogy(tth_x, np.maximum(prof, 0.5), lw=0.7, color=\"k\")\nfor t in rings_all[\"tth\"]:\n    ax_p.axvline(t, color=\"0.6\", ls=\"--\", lw=0.6)\nax_p.set(ylabel=\"counts\", title=\"2θ profile: every ring (grey), windows of the rings used (orange)\")\nk = 5  # 0.01 deg in 2theta\nCb = C[: n_t // k * k].reshape(-1, k, N_CAKE).sum(1)\nax_c.imshow(np.ma.masked_less_equal(Cb.T, 0), aspect=\"auto\", origin=\"lower\", cmap=MAGMA, norm=LogNorm(),\n            extent=(tth_lo, tth_lo + len(Cb) * k * TTH_STEP, -180, 180), interpolation=\"nearest\")\nax_c.set(xlabel=\"2θ (°)\", ylabel=\"η (°)\", title=\"cake: rings should be straight; windows in cyan\")\naxz = np.atleast_1d(sf[1].subplots(1, n_all, sharex=True, sharey=True))\nfor r, a in enumerate(axz):\n    z = np.abs(tth_x - rings_all[\"tth\"][r]) < 0.15\n    a.semilogy(tth_x[z] - rings_all[\"tth\"][r], np.maximum(prof[z], 0.5), lw=0.8, color=\"k\")\n    a.set_title(f\"{r}: {ring_info[r][0]}\", fontsize=9)\n    a.tick_params(labelsize=7)\naxz[0].set_ylabel(\"counts\")\nsf[1].supxlabel(\"2θ − ring (°); orange: window, blue: measured centre; grey: ring not used\", fontsize=9)\ndrawn = []  # the windows, redrawn when the sliders move\n\n\ndef show_rings(n_rings, tth_tol):\n    tth = rings_all[\"tth\"][:n_rings]\n    off, hw = ix.ring_widths(prof, tth_lo, TTH_STEP, tth)\n    tol = np.abs(off) + hw if tth_tol == 0 else np.full(n_rings, tth_tol)\n    inside = np.zeros(n_t, bool)\n    for t, w in zip(tth, tol):\n        inside |= np.abs(tth_x - t) < w\n    print(\" ring  hkl        mult  |F|^2  2theta   offset  half-width  tolerance  gap to next  counts in window\")\n    for r in range(n_rings):\n        win = np.abs(tth_x - tth[r]) < tol[r]\n        g = f\"{gap[r]:11.3f}\" if r < len(gap) else \"          -\"\n        warn = \"  <- windows overlap\" if r < len(gap) and r + 1 < n_rings and tol[r] + tol[r + 1] > gap[r] else \"\"\n        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}  \"\n              f\"{off[r]:+7.3f}  {hw[r]:10.3f}  {tol[r]:9.3f}  {g}  {prof[win].sum():.3g}{warn}\")\n    used = tth_x < tth[-1] + 0.5\n    print(f\"{prof[inside].sum() / prof[used].sum() * 100:.1f}% of the intensity up to the last ring used is in a \"\n          \"window\")\n    for art in drawn:\n        art.remove()\n    drawn.clear()\n    for r in range(n_rings):\n        drawn.append(ax_p.axvspan(tth[r] - tol[r], tth[r] + tol[r], color=\"tab:orange\", alpha=0.4))\n        drawn.extend(ax_c.axvline(tth[r] + sg_ * tol[r], color=\"c\", lw=0.6) for sg_ in (-1, 1))\n        drawn.append(axz[r].axvspan(-tol[r], tol[r], color=\"tab:orange\", alpha=0.4))\n        drawn.append(axz[r].axvline(off[r], color=\"tab:blue\", lw=0.8))\n    for r, a in enumerate(axz):\n        a.set_facecolor(\"white\" if r < n_rings else \"0.9\")\n    fig.canvas.draw_idle()\n\n\nw_rings = widgets.interactive(\n    show_rings,\n    n_rings=widgets.IntSlider(value=n_rings, min=1, max=n_all, description=\"n_rings\", continuous_update=False),\n    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,\n                                readout_format=\".3f\", description=\"tth_tol\", continuous_update=False),\n)\nfreeze(fig)\nprint(\"tth_tol 0: measured per ring\")\ndisplay(w_rings, fig.canvas)",
   "id": "cell-07"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 2. The data: `--monitor`, `--gridstep`\n\nAll 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).\n\nThen 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.\n\n`--gridstep` sums rows in groups to make voxels larger: fewer voxels and rows, a much cheaper fit, coarser maps.\n\nEach 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.",
   "id": "cell-08"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "n_rings_used = w_rings.kwargs[\"n_rings\"]\nrings = ix.ring_table(lpars, sg, WL, n_rings_used, structure)\nring_off, ring_hw = ix.ring_widths(prof, tth_lo, TTH_STEP, rings[\"tth\"])\nrings[\"hw\"] = ring_hw\ntth_tol_used = w_rings.kwargs[\"tth_tol\"] or None\ntth_tol_r = np.abs(ring_off) + ring_hw if tth_tol_used is None else np.full(n_rings_used, tth_tol_used)\nprint(f\"{len(rings['hkls'])} hkls in {n_rings_used} rings; tolerances (deg) {np.round(tth_tol_r, 3)}\")\nn_cells = n_rings_used * n_e * n_o * NK\nprint(f\"histogram: {n_cells / 1e6:.0f}M bins, {n_cells * 4 / 1e9:.2f} GB (x ~4 in the fit)\")\n\nt0 = time.time()\nH_lit, H = ix.histogram_pixels(stream(groups), geom, rings[\"tth\"], tth_tol_r, OM0,\n                               [((B_E, B_O, N_E, N_O), 1), (bins, NK)], chunk)\nHs = H_lit.reshape(n_rings_used, N_E, N_O)\nH4 = np.asarray(H).reshape(n_rings_used, n_e, n_o, NK)  # [ring, eta, omega, row] on the host, for the plots\nprint(f\"histograms: {time.time() - t0:.0f} s; {float(jnp.mean(H > 0)) * 100:.1f}% of the bins non-empty\")",
   "id": "cell-09"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "row_tot = H4.sum((0, 1, 2))\nom_tot = H4.sum((0, 1, 3))\nnb = np.convolve(row_tot, [0.5, 0.0, 0.5], \"same\")  # the mean of each row's neighbours\nlit_row = row_tot > 0.05 * row_tot.max()\nlit_row[[0, -1]] = False\njump = np.where(lit_row, row_tot / np.maximum(nb, 1e-30), np.nan)\n\nn_ax = 4 if monitor else 3\nfig, ax = plt.subplots(1, n_ax, figsize=(5 * n_ax, 3.6), layout=\"constrained\")\nif dty_rows is not None:\n    for v in dty_rows[:: max(1, len(dty_rows) // 10)]:\n        ax[-1].plot(np.linspace(0, 1, len(v)), (v - v.mean()) / ystep0, lw=0.7)\nax[-1].set(xlabel=\"fraction of the row\", ylabel=\"dty − row mean (dty steps)\", title=\"the DataSet's dty within a row\")\nax[0].plot(dty_c, row_tot, \".-\", ms=3)\nax0 = ax[0].twinx()\nax0.plot(dty_c, jump, \".\", color=\"tab:red\", ms=3)\nax0.axhline(1, color=\"tab:red\", lw=0.5)\nax0.set_ylabel(\"row / mean of its neighbours\", color=\"tab:red\")\nax[0].set(xlabel=\"dty\", ylabel=\"intensity\", title=\"per dty row\")\nax[1].plot(om_c, om_tot, lw=0.8)\nax[1].set(xlabel=\"ω (°)\", ylabel=\"intensity\", title=\"per ω bin, all rows\")\nif monitor:\n    for g in groups[:: max(1, len(groups) // 20)]:\n        ax[2].plot(mons[g] / monitor_ref, lw=0.5)\n    ax[2].set(xlabel=\"frame\", ylabel=f\"{monitor} / mean\", title=f\"{monitor}, some scans\")\nplt.show()\nprint(f\"row-to-row scatter (rows with intensity): {np.nanstd(jump) * 100:.1f}%; sample edges excluded\")",
   "id": "cell-10"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 3. Rotation axis: `--y0`\n\ny0 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.\n\nHere 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.\n\nCheck 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.",
   "id": "cell-11"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "sino = H4.sum((0, 1))  # [omega, row], every ring and eta\nw_om = sino.sum(1)\ngood = w_om > 0.05 * w_om.max()\ncom = sino @ dty_c / np.maximum(w_om, 1e-30)  # intensity-weighted mean dty\n\n\ndef mid(p, frac=0.05):\n    # the middle of one projection: halfway between its edges, where it crosses frac x its maximum\n    thr = frac * p.max()\n    k = np.flatnonzero(p > thr)\n    i0, i1 = k[0], k[-1]\n    lo = dty_c[i0] - YSTEP * (p[i0] - thr) / max(p[i0] - (p[i0 - 1] if i0 > 0 else 0.0), 1e-30)\n    hi = dty_c[i1] + YSTEP * (p[i1] - thr) / max(p[i1] - (p[i1 + 1] if i1 < NK - 1 else 0.0), 1e-30)\n    return 0.5 * (lo + hi)\n\n\nmids = np.array([mid(p) if g else np.nan for p, g in zip(sino, good)])\nA = np.stack([np.ones(n_o), np.sin(np.radians(om_c)), np.cos(np.radians(om_c))], 1)\ncoef = np.linalg.lstsq(A[good], mids[good], rcond=None)[0]\ncoef_com = np.linalg.lstsq(A[good], com[good], rcond=None)[0]\ny0_fit = float(coef[0])\nedge = max(row_tot[0], row_tot[-1]) / row_tot.max()\nprint(f\"y0 from the sample's edges: {y0_fit:.4f} (DataSet / parameters: {Y0_DS:.4f}; \"\n      f\"{(y0_fit - Y0_DS) / YSTEP:+.2f} dty steps); from the mean dty: {coef_com[0]:.4f}\")\nif edge > 0.1:\n    print(f\"the first or last row holds {edge * 100:.0f}% of the brightest row's intensity: the sample may extend past \"\n          \"the scan, which biases the fit\")\n\n\n_, pad_ds = anri.geom.sino_shift_and_pad(Y0_DS, NK, DTY0, YSTEP)\nnr_ds = NK + pad_ds\npos_ds = anri.geom.recon_positions(nr_ds, YSTEP).astype(jnp.float32)\nsino_log = np.log1p(H4).sum((0, 1))  # [omega, row]: every spot counts about the same, as in tomo_2_map\nsino_f = jnp.asarray(ix.ramp_filter(sino_log), jnp.float32)\nom_j = jnp.asarray(om_c, jnp.float32)\nrec_ds = np.asarray(ix.backproject(sino_f, om_j, pos_ds, Y0_DS, DTY0, YSTEP, ddty)).reshape(nr_ds, nr_ds)\nvmax = np.percentile(rec_ds, 99.5)\next = (-nr_ds / 2 * YSTEP, nr_ds / 2 * YSTEP) * 2\next_sino = (OM0, OM0 + n_o * bins[1], DTY0 - YSTEP / 2, DTY0 + (NK - 0.5) * YSTEP)  # omega, dty\n\nfig = figure(\"y0\", figsize=(15, 4.6))\nax = fig.subplots(1, 3)\nax[2].sharex(ax[1])\nax[2].sharey(ax[1])\nax[0].imshow(sino_log.T, aspect=\"auto\", origin=\"lower\", cmap=\"magma\", extent=ext_sino, interpolation=\"nearest\")\nax[0].plot(om_c[good], com[good], \".\", color=\"0.6\", ms=2, label=\"mean dty\")\nax[0].plot(om_c[good], mids[good], \".\", color=\"c\", ms=2, label=\"middle of the edges\")\nax[0].plot(om_c, A @ coef, color=\"w\", lw=1, label=f\"fit: y0 = {y0_fit:.3f}\")\ny0_line = ax[0].axhline(Y0_DS, color=\"tab:green\", lw=1, ls=\"--\", label=\"y0 (slider)\")\nax[0].legend(fontsize=8, loc=\"upper right\")\nax[0].set(xlabel=\"ω (°)\", ylabel=\"dty\", title=\"sinogram of log(1 + counts), all rings\")\nax[1].imshow(rec_ds, origin=\"lower\", cmap=\"gray\", vmin=0, vmax=vmax, extent=ext)\nax[1].set(title=f\"back-projection, y0 = {Y0_DS:.6g} (DataSet)\", xlabel=\"voxels (dty units)\")\nim_y0 = ax[2].imshow(rec_ds, origin=\"lower\", cmap=\"gray\", vmin=0, vmax=vmax, extent=ext)\nax[2].set(xlabel=\"voxels (dty units)\")\n\n\ndef show_y0(y0_try):\n    im_y0.set_data(np.asarray(ix.backproject(sino_f, om_j, pos_ds, y0_try, DTY0, YSTEP, ddty)).reshape(nr_ds, nr_ds))\n    y0_line.set_ydata([y0_try, y0_try])\n    ax[2].set_title(f\"back-projection, y0 = {y0_try:.6g} (slider)\")\n    fig.canvas.draw_idle()\n\n\nlo_, hi_ = min(Y0_DS, y0_fit) - 10 * YSTEP, max(Y0_DS, y0_fit) + 10 * YSTEP\nw_y0 = widgets.interactive(show_y0, y0_try=widgets.FloatSlider(value=Y0_DS, min=lo_, max=hi_, step=YSTEP / 20,\n                                                               readout_format=\".6g\", description=\"y0\",\n                                                               continuous_update=False))\nfreeze(fig)\ndisplay(w_y0, fig.canvas)",
   "id": "cell-12"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "Y0 = w_y0.kwargs[\"y0_try\"]  # the rotation axis used from here on\n_, PAD = anri.geom.sino_shift_and_pad(Y0, NK, DTY0, YSTEP)\nNR = NK + PAD  # the voxel grid: NR x NR, centred on the rotation axis, as ImageD11 pads its reconstructions\nNV = NR * NR\nprint(f\"y0 = {Y0:.4f}; voxels {NR} x {NR} of {YSTEP:g} ({NK} dty rows + pad {PAD})\")",
   "id": "cell-13"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "### The sample mask: `--mask`\n\nThe 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`).\n\nWith 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.",
   "id": "cell-13a"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "scan0 = {\"y0\": Y0, \"dty0\": DTY0, \"ystep\": YSTEP, \"n_rows\": NK, \"om0\": OM0}\nif ddty is not None:\n    scan0[\"ddty\"] = ddty\nrec = ix.reconstruct(H, n_rings_used, n_e, n_o, scan0, bins[1], NR)  # [NR, NR], reconstruction order\nt_otsu = ix.otsu(rec)\nvmax_r = np.percentile(rec, 99.5)\nfig = figure(\"mask\", figsize=(14, 4.4))\nax = fig.subplots(1, 3)\nax[1].sharex(ax[0])\nax[1].sharey(ax[0])\nax[0].imshow(rec, origin=\"lower\", cmap=\"gray\", vmin=0, vmax=vmax_r)\nax[0].set(title=f\"back-projection, y0 = {Y0:.6g}\", xticks=[], yticks=[])\nim_m = ax[1].imshow(rec, origin=\"lower\", cmap=\"gray\", vmin=0, vmax=vmax_r)\nax[1].set(xticks=[], yticks=[])\nax[2].hist(rec.ravel(), bins=200, log=True, color=\"0.5\")\nax[2].axvline(t_otsu, color=\"k\", ls=\":\", label=\"Otsu\")\nv_t = ax[2].axvline(t_otsu, color=\"tab:red\", label=\"threshold\")\nax[2].set(xlabel=\"back-projection\", ylabel=\"voxels\")\nax[2].legend(fontsize=8)\n\n\ndef show_mask(threshold):\n    m = ix.threshold_mask(rec, threshold)\n    im_m.set_data(np.where(m, rec, np.nan))\n    v_t.set_xdata([threshold, threshold])\n    ax[1].set_title(f\"mask: {m.sum()} of {NV} voxels\")\n    fig.canvas.draw_idle()\n\n\nw_mask = widgets.interactive(show_mask, threshold=widgets.FloatSlider(\n    value=t_otsu, min=float(np.percentile(rec, 1)), max=float(np.percentile(rec, 99.9)),\n    step=float(np.ptp(rec)) / 500, readout_format=\".3g\", description=\"threshold\", continuous_update=False))\nfreeze(fig)\ndisplay(w_mask, fig.canvas)\nif draw:\n    get_mask = ix.draw_mask(rec)  # draw the outline on this figure, then run the next cell",
   "id": "cell-13b"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# the mask used from here on: drawn, thresholded at the slider (mask = \"auto\"), a file, or every voxel\nmask_used = mask\nif draw:\n    mask_v = get_mask()\nelif mask == \"auto\":\n    mask_v = ix.threshold_mask(rec, w_mask.kwargs[\"threshold\"])\nelif mask is not None:\n    mask_v = np.load(mask).astype(bool)\nelse:\n    mask_v = np.ones((NR, NR), bool)\nif draw or (mask == \"auto\" and not np.isclose(w_mask.kwargs[\"threshold\"], t_otsu)):\n    mask_used = os.path.join(outdir, f\"{dsname}_mask.npy\")  # for --mask\n    np.save(mask_used, mask_v)\nmask_v = mask_v.ravel()\nprint(f\"mask: {mask_v.sum()} of {NV} voxels fitted\" + (f\" (--mask {mask_used})\" if mask_used else \" (no mask)\"))",
   "id": "cell-13c"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 4. Lit map: `--lit` and `--etacut`\n\nThe 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.",
   "id": "cell-14"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "med = float(jnp.median(Hs[Hs > 0]))\nvals = np.sort(np.asarray(Hs).ravel())\ncum = np.concatenate([[0.0], np.cumsum(vals, dtype=np.float64)])\nmults = np.logspace(-1, 2, 61)\nidx = np.searchsorted(vals, mults * med, side=\"right\")\narea = 1 - idx / vals.size  # fraction of the bins above the cut\nheld = 1 - cum[idx] / cum[-1]  # fraction of the intensity in them\nHs_np = np.asarray(Hs)\nv_pos = vals[np.searchsorted(vals, 0.0, side=\"right\")]  # the smallest non-zero bin\n\nfig = figure(\"lit\", figsize=(15, 4.6))\nax = fig.subplots(1, 3, width_ratios=[1, 1.3, 1.3])\nax[2].sharex(ax[1])\nax[2].sharey(ax[1])\nax[0].loglog(mults, area, label=\"lit area\")\nax[0].loglog(mults, held, label=\"intensity held\")\nv_lit = ax[0].axvline(lit_cut, color=\"tab:red\")\nax[0].set(xlabel=\"lit cut (x median non-empty bin)\", ylabel=\"fraction\", title=f\"median non-empty bin: {med:.3g}\")\nax[0].legend()\ne = (OM0, OM0 + N_O * B_O, -180, 180)\nim_int = ax[1].imshow(np.ma.masked_less_equal(Hs_np[0], 0), aspect=\"auto\", origin=\"lower\", cmap=MAGMA,\n                      norm=LogNorm(v_pos, vals[-1]), extent=e, interpolation=\"nearest\")\nax[1].set(xlabel=\"ω (°)\", ylabel=\"η (°)\")\nim_lit = ax[2].imshow(np.zeros(Hs_np[0].shape), aspect=\"auto\", origin=\"lower\", cmap=\"gray_r\", vmin=0, vmax=1,\n                      extent=e, interpolation=\"nearest\")\nax[2].set(xlabel=\"ω (°)\")\nec = np.degrees(np.arcsin(etacut))\nfor a in ax[1:]:\n    for el in (ec, -ec, 180 - ec, ec - 180):  # where |sin eta| = etacut (the parameter)\n        a.axhline(el, color=\"tab:cyan\", lw=0.6, ls=\":\")\n\n\ndef show_lit(lit_cut, ring):\n    i = np.searchsorted(vals, lit_cut * med, side=\"right\")\n    lit_np = Hs_np[ring] > lit_cut * med\n    v_lit.set_xdata([lit_cut, lit_cut])\n    im_int.set_data(np.ma.masked_less_equal(Hs_np[ring], 0))\n    im_lit.set_data(lit_np.astype(np.float32))\n    ax[1].set_title(f\"ring {ring}: intensity, rows summed\")\n    ax[2].set_title(f\"ring {ring}: lit, {lit_np.mean() * 100:.2f}% of its bins\")\n    fig.canvas.draw_idle()\n    print(f\"lit: {(1 - i / vals.size) * 100:.2f}% of the bins, holding {(1 - cum[i] / cum[-1]) * 100:.1f}% of the \"\n          f\"intensity; per ring {np.round((Hs_np > lit_cut * med).mean((1, 2)) * 100, 2)}%\")\n\n\nw_lit = widgets.interactive(\n    show_lit,\n    lit_cut=widgets.FloatLogSlider(value=lit_cut, base=10, min=-1, max=2, step=0.05, description=\"lit\",\n                                   continuous_update=False),\n    ring=widgets.IntSlider(value=0, min=0, max=n_rings_used - 1, description=\"ring\", continuous_update=False),\n)\nfreeze(fig)\ndisplay(w_lit, fig.canvas)",
   "id": "cell-15"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "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).",
   "id": "cell-16"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "eta_x = -180.0 + (np.arange(N_E) + 0.5) * B_E\nprof_eta = Hs_np.sum((0, 2))\n\nfig = figure(\"etacut\", figsize=(12, 3.6))\nax = fig.subplots()\nax.plot(eta_x, prof_eta, lw=0.8, color=\"k\")\nax.set(xlabel=\"η (°)\", ylabel=\"intensity, all rings\", xlim=(-180, 180), title=\"red: not used\")\nax2 = ax.twinx()\nax2.semilogy(eta_x, 1 / np.maximum(np.abs(np.sin(np.radians(eta_x))), 1e-3), color=\"tab:blue\", lw=0.8)\nax2.set_ylim(1, 50)\nax2.set_ylabel(\"ω smearing, 1 / |sin η|\", color=\"tab:blue\")\nshade = []\n\n\ndef show_etacut(etacut):\n    cut = np.abs(np.sin(np.radians(eta_x))) <= etacut\n    for art in shade:\n        art.remove()\n    shade[:] = [ax.fill_between(eta_x, 0, prof_eta.max(), where=cut, color=\"tab:red\", alpha=0.2)]\n    fig.canvas.draw_idle()\n    print(f\"etacut {etacut:.2f}: {2 * np.arcsin(etacut) / np.pi * 100:.1f}% of the reflections and \"\n          f\"{prof_eta[cut].sum() / prof_eta.sum() * 100:.1f}% of the measured intensity not used; spots smeared up to \"\n          f\"{1 / max(etacut, 1e-3):.1f} x in ω are kept\")\n\n\nw_eta = widgets.interactive(show_etacut, etacut=widgets.FloatSlider(value=etacut, min=0.0, max=0.7, step=0.01,\n                                                                    description=\"etacut\", continuous_update=False))\nfreeze(fig)\ndisplay(w_eta, fig.canvas)",
   "id": "cell-17"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 5. Orientation grid: `--grid`, `--max-chance`\n\nThe 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).\n\nThe 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.",
   "id": "cell-18"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "lit_used, etacut_used = w_lit.kwargs[\"lit_cut\"], w_eta.kwargs[\"etacut\"]\nlit = {\"table\": ix.lit_table(Hs > lit_used * med), \"om0\": OM0, \"bins\": (B_E, B_O, N_E, N_O), \"frame_step\": OSTEP,\n       \"etacut\": etacut_used}\nrng = np.random.default_rng(0)  # as choose_grid\nchance, n_grid = {}, {}\nt0 = time.time()\nfor s in ix.GRID_STEPS:\n    if n_grid and n_grid[ix.GRID_STEPS[0]] * (ix.GRID_STEPS[0] / s) ** 3 > 4e6:\n        print(f\"grid {s} deg: skipped, ~{n_grid[ix.GRID_STEPS[0]] * (ix.GRID_STEPS[0] / s) ** 3 / 1e6:.0f}M \"\n              \"orientations\")\n        continue\n    U_s, d_s = anri.crystal.orientation_grid(s, ops)\n    smp = U_s[rng.choice(len(U_s), min(len(U_s), 1 << 14), replace=False)]\n    chance[s], n_grid[s] = float(np.median(ix.completeness_of(smp, d_s, B, rings, geom, lit))), len(U_s)\n    print(f\"grid {s} deg: {len(U_s)} orientations, up to {d_s:.2f} deg from the truth; chance completeness \"\n          f\"{chance[s]:.2f}\")\ndel U_s\nprint(f\"{time.time() - t0:.0f} s\")\n\n\ndef auto_step(max_chance):\n    # the step choose_grid picks: the coarsest at or below max_chance, else the finest\n    return next((s for s in chance if chance[s] <= max_chance), ix.GRID_STEPS[-1])\n\n\nfig = figure(\"grid\", figsize=(7, 3.6))\nax = fig.subplots()\nxs = [str(s) for s in chance]\nbars = ax.bar(xs, list(chance.values()))\nh_chance = ax.axhline(max_chance, color=\"k\", ls=\":\")\nax.set(xlabel=\"grid step (°)\", ylabel=\"chance completeness\", ylim=(0, 1))\nax2 = ax.twinx()\nax2.semilogy(xs, list(n_grid.values()), \"o-\", color=\"0.4\")\nax2.set_ylabel(\"orientations\", color=\"0.4\")\n\n\ndef show_grid(max_chance):\n    st = grid if grid is not None else auto_step(max_chance)\n    for b, s in zip(bars, chance):\n        b.set_color(\"tab:red\" if s == st else \"tab:blue\")\n    h_chance.set_ydata([max_chance, max_chance])\n    ax.set_title(f\"grid step {st}°\" + (\" (fixed by grid in the parameters)\" if grid is not None else \" (automatic)\"))\n    fig.canvas.draw_idle()\n\n\nw_grid = widgets.interactive(show_grid, max_chance=widgets.FloatSlider(\n    value=max_chance, min=0.05, max=0.95, step=0.05, description=\"max_chance\", continuous_update=False))\nfreeze(fig)\ndisplay(w_grid, fig.canvas)",
   "id": "cell-19"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**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).\n\nBelow, 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.",
   "id": "cell-19a"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "if grid_scan:\n    d_rows = H.reshape(-1, NK).sum(1)  # the histogram with its dty rows summed\n    dev_s, n_lr_s = {}, {}\n    for s in [*chance, 0.75]:\n        t0 = time.time()\n        U_s, d_s = anri.crystal.orientation_grid(s, ops)\n        _, comp_s, info_s = ix.prune(U_s, d_s, B, rings, geom, lit, None, keep)\n        pre_s = np.flatnonzero(comp_s > info_s[\"chance\"])\n        _, lr_s, mu_s = ix.orientation_mlem(d_rows, U_s[pre_s], B, rings, geom, (*bins, OM0), etacut_used,\n                                            log=lambda m: None, censor=censor, return_model=True)\n        dev_s[s], n_lr_s[s] = float(ix.deviance(d_rows, jnp.asarray(mu_s), censor)), int(np.sum(lr_s > min_lr))\n        print(f\"grid {s} deg: {len(pre_s)} orientations above chance; deviance {dev_s[s]:.5g}; {n_lr_s[s]} with \"\n              f\"likelihood ratio > {min_lr:g} ({time.time() - t0:.0f} s)\")\n    del U_s, comp_s\n    steps_s = list(dev_s)\n    fig, ax = plt.subplots(1, 2, figsize=(11, 3.6), layout=\"constrained\")\n    ax[0].plot(steps_s, [dev_s[s] / dev_s[min(steps_s)] for s in steps_s], \"o-\")\n    ax[0].set(xlabel=\"grid step (°)\", ylabel=f\"deviance / deviance at {min(steps_s)}°\",\n              title=\"orientation fit to the row-summed data\")\n    ax[0].invert_xaxis()\n    ax[1].plot(steps_s, [n_lr_s[s] for s in steps_s], \"o-\", color=\"0.4\")\n    ax[1].set(xlabel=\"grid step (°)\", ylabel=f\"orientations with likelihood ratio > {min_lr:g}\", yscale=\"log\")\n    ax[1].invert_xaxis()\n    plt.show()",
   "id": "cell-19b"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 6. Pruning: `--prune`, `--min-comp`, `--min-lr`, `--keep`\n\nEvery grid orientation gets its completeness. Real orientations sit well above the chance level; the bulk of the grid sits at it.\n\n- `--prune completeness` keeps everything at or above `min_comp` (default halfway from chance to the maximum), the most complete first, at most `keep`.\n- `--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.",
   "id": "cell-20"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "step = grid if grid is not None else auto_step(w_grid.kwargs[\"max_chance\"])\nU_grid, delta = anri.crystal.orientation_grid(step, ops)\nt0 = time.time()\n_, comp, info = ix.prune(U_grid, delta, B, rings, geom, lit, None, keep)\nhalfway = info[\"min_comp\"]  # the default cut of --prune completeness\nprint(f\"grid {step} deg: {len(U_grid)} orientations; completeness ({time.time() - t0:.0f} s): \"\n      f\"chance {info['chance']:.2f}, 99th percentile {np.percentile(comp, 99):.2f}, max {comp.max():.2f}; \"\n      f\"halfway {halfway:.2f}\")\ncomp_sorted = np.sort(comp)\n\n\ndef default_cut(prune):\n    return info[\"chance\"] if prune == \"likelihood\" else halfway\n\n\nfig = figure(\"completeness\", figsize=(13, 3.6))\nax = fig.subplots(1, 2, sharex=True)\nax[0].hist(comp, bins=100, log=True, color=\"tab:blue\")\nax[0].axvline(info[\"chance\"], color=\"k\", ls=\":\", label=\"chance\")\nax[0].axvline(halfway, color=\"0.5\", ls=\"--\", label=\"halfway\")\nv_comp = [a.axvline(0.0, color=\"tab:red\", label=\"min_comp\") for a in ax]\nax[0].set(xlabel=\"completeness\", ylabel=\"orientations\")\nax[0].legend()\ncs = np.linspace(0, 1, 201)\nax[1].semilogy(cs, np.maximum(comp.size - np.searchsorted(comp_sorted, cs), 0.5), color=\"tab:blue\")\nax[1].set(xlabel=\"min_comp\", ylabel=\"orientations above\")\n\n\ndef show_comp(prune, min_comp):\n    n_above = comp.size - np.searchsorted(comp_sorted, min_comp, side=\"left\" if prune == \"completeness\" else \"right\")\n    for v in v_comp:\n        v.set_xdata([min_comp, min_comp])\n    fig.canvas.draw_idle()\n    if prune == \"likelihood\":\n        print(f\"{n_above} orientations above {min_comp:.2f} go to the orientation fit\")\n    else:\n        print(f\"{n_above} orientations at or above {min_comp:.2f}; {min(n_above, keep)} kept (keep = {keep})\")\n\n\nw_comp = widgets.interactive(\n    show_comp,\n    prune=widgets.Dropdown(options=[\"likelihood\", \"completeness\"], value=prune, description=\"prune\"),\n    min_comp=widgets.FloatSlider(value=min_comp if min_comp is not None else default_cut(prune), min=0.0, max=1.0,\n                                 step=0.01, description=\"min_comp\", continuous_update=False),\n)\n\n\ndef reset_cut(change):\n    w_comp.children[1].value = default_cut(change[\"new\"])  # a new prune mode starts from its default min_comp\n\n\nw_comp.children[0].observe(reset_cut, names=\"value\")\nfreeze(fig)\nprint(f\"min_comp defaults: likelihood {info['chance']:.2f} (chance), completeness {halfway:.2f} (halfway)\")\ndisplay(w_comp, fig.canvas)",
   "id": "cell-21"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "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.\n\nIn 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.",
   "id": "cell-22"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "prune_used, min_comp_used = w_comp.kwargs[\"prune\"], w_comp.kwargs[\"min_comp\"]\nif min_comp is None and np.isclose(min_comp_used, default_cut(prune_used)):\n    min_comp_used = None  # the command line's own default\ncut = default_cut(prune_used) if min_comp_used is None else min_comp_used\nif prune_used == \"likelihood\":\n    pre = np.flatnonzero(comp > cut)\n    d_rows = H.reshape(-1, NK).sum(1)  # the histogram with its dty rows summed\n    t0 = time.time()\n    g_occ, lr = ix.orientation_mlem(d_rows, U_grid[pre], B, rings, geom, (*bins, OM0), etacut_used, censor=censor)\n    print(f\"orientation fit of {len(pre)} orientations: {time.time() - t0:.0f} s\")\n    lr_sorted = np.sort(lr)\n\n    # each orientation's misorientation to the nearest one with a higher likelihood ratio (the strongest 5000)\n    top = np.argsort(lr)[::-1][: min(len(pre), 5000)]\n    q = anri.crystal.mat_to_quat(U_grid[pre[top]]).astype(np.float32)  # [M, 4]\n    qs = anri.crystal.mat_to_quat(U_grid[pre[top]][:, None] @ ops[None]).astype(np.float32)  # [M, S, 4]\n    dist = np.full(len(top), 180.0)\n    for s0 in range(0, len(top), 256):\n        c = np.abs(np.einsum(\"ik,jsk->ijs\", q[s0 : s0 + 256], qs)).max(2)  # cos(half the misorientation)\n        c[np.arange(c.shape[0])[:, None] <= np.arange(len(top))[None] - s0] = -1.0  # only stronger ones\n        dist[s0 : s0 + 256] = np.degrees(2 * np.arccos(np.clip(c.max(1), -1.0, 1.0)))\n    dist[0] = 180.0\n    lr_top = np.maximum(lr[top], 0.1)\n\n    fig = figure(\"likelihood\", figsize=(16, 3.9))\n    ax = fig.subplots(1, 3)\n    ax[2].sharey(ax[0])\n    ax[0].scatter(comp[pre], np.maximum(lr, 0.1), s=1, alpha=0.4)\n    h_lr = [ax[0].axhline(min_lr, color=\"tab:red\", label=\"min_lr\"), ax[2].axhline(min_lr, color=\"tab:red\")]\n    ax[0].axvline(cut, color=\"tab:orange\", ls=\"--\", label=\"min_comp\")\n    ax[0].set(yscale=\"log\", xlabel=\"completeness\", ylabel=\"likelihood ratio\")\n    ax[0].legend(loc=\"upper left\", fontsize=8)\n    ts = np.logspace(-1, max(1, np.log10(max(lr.max(), 10))), 100)\n    ax[1].loglog(ts, np.maximum(lr.size - np.searchsorted(lr_sorted, ts, side=\"right\"), 0.5))\n    v_lr = ax[1].axvline(min_lr, color=\"tab:red\")\n    h_keep = ax[1].axhline(keep, color=\"0.5\", ls=\":\", label=\"keep\")\n    ax[1].set(xlabel=\"min_lr\", ylabel=\"orientations kept\", ylim=(0.5, max(2 * len(pre), 2 * keep)))\n    ax[1].legend(fontsize=8)\n    sc_out = ax[2].scatter(dist, lr_top, s=2, color=\"0.75\", label=\"not kept\")\n    sc_in = ax[2].scatter(dist, lr_top, s=3, color=\"tab:blue\", label=\"kept\")\n    ax[2].axvline(1.8 * step, color=\"0.4\", ls=\"--\", label=\"population radius\")\n    ax[2].set(xscale=\"log\", xlabel=\"misorientation to the nearest stronger orientation (°)\",\n              title=f\"the {len(top)} strongest\")\n    ax[2].legend(fontsize=8, loc=\"lower right\")\n\n    def show_lr(min_lr, keep):\n        ok = lr > min_lr\n        n_ok = int(ok.sum())\n        k = pre[ok][np.argsort(lr[ok])[::-1][:keep]] if n_ok else pre[:0]\n        sel = np.isin(pre[top], k)\n        for h_ in h_lr:\n            h_.set_ydata([min_lr, min_lr])\n        v_lr.set_xdata([min_lr, min_lr])\n        h_keep.set_ydata([keep, keep])\n        sc_in.set_offsets(np.c_[dist[sel], lr_top[sel]])\n        sc_out.set_offsets(np.c_[dist[~sel], lr_top[~sel]])\n        fig.canvas.draw_idle()\n        capped = f\" (capped by keep = {keep})\" if n_ok > keep else \"\"\n        print(f\"{n_ok} with likelihood ratio > {min_lr:.3g}; {len(k)} kept{capped}; \"\n              f\"{int(np.sum(sel & (dist > 1.8 * step)))} kept peaks (beyond 1.8 grid steps from anything stronger)\")\n\n    w_lr = widgets.interactive(\n        show_lr,\n        min_lr=widgets.FloatLogSlider(value=min_lr, base=10, min=0, max=4, step=0.05, description=\"min_lr\",\n                                      continuous_update=False),\n        keep=widgets.BoundedIntText(value=keep, min=1, max=10**8, description=\"keep\"),\n    )\n    freeze(fig)\n    display(w_lr, fig.canvas)\nelse:\n    print(\"prune = completeness: no orientation fit\")",
   "id": "cell-23"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "if prune_used == \"likelihood\":\n    min_lr_used, keep_used = w_lr.kwargs[\"min_lr\"], w_lr.kwargs[\"keep\"]\n    above = np.flatnonzero(lr > min_lr_used)\n    kept = pre[above[np.argsort(lr[above])[::-1]][:keep_used]]\nelse:\n    min_lr_used, keep_used = min_lr, keep\n    above = np.flatnonzero(comp >= cut)\n    kept = above[np.argsort(comp[above])[::-1]][:keep_used]\nU_kept = U_grid[kept]\nprint(f\"{len(U_kept)} orientations kept for the occupancy fit\")",
   "id": "cell-24"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 7. Occupancy: `--cand`, `--iter`, `--block-gb`, `--beam`, `--censor`\n\nEach 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.\n\n`--block-gb` only sets how many voxels are handled at once (memory, not results).\n\n- `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.\n- `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.",
   "id": "cell-25"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "if run_fit:\n    pred = ix.predictions(U_kept, B, rings, geom, etacut_used)\n    pos = np.asarray(anri.geom.recon_positions(NR, YSTEP), np.float32)  # every voxel\n    pos_fit = pos[mask_v]  # the voxels fitted\n    n_fit = len(pos_fit)\n    scan = {\"y0\": Y0, \"dty0\": DTY0, \"ystep\": YSTEP, \"n_rows\": NK, \"om0\": OM0}\n    if ddty is not None:  # fly, helical scans\n        scan[\"ddty\"] = ddty\n    beam_used = beam\n    bins_fit = bins\n    if beam_used > 0:  # the beam's profile across dty, as python -m anri.index\n        scan.update({\"sig_beam\": beam_used / (2 * np.sqrt(2 * np.log(2))), \"width_beam\": (gridstep - 1) * ystep0,\n                     \"voxel\": YSTEP})\n        bins_fit = (*bins, ix.beam_rows(scan))\n    n_corners = 4 * bins_fit[4] if len(bins_fit) > 4 else 8\n    ring_j = jnp.asarray(rings[\"ring_j\"])\n    n_j = pred[0].shape[1]\n    vb_all = ix.block_voxels(16, n_j, block_gb * 1e9, n_corners)\n    vb = ix.block_voxels(n_cand, n_j, block_gb * 1e9, n_corners)\n    pos_p = ix.pad_voxels(pos_fit, max(vb_all, vb))\n    print(f\"{n_fit} voxels x {len(U_kept)} orientations x {n_j} reflections; occupancies {n_fit * n_cand * 8 / 1e9:.2f} \"\n          f\"GB, blocks of {vb} voxels\" + (f\"; beam FWHM {beam_used:g}: {bins_fit[4]} rows per voxel\" if beam_used > 0\n                                          else \"; the 2-row model\"))\n    if beam_used > 0:  # the weights per row for a voxel at a few places between two rows, mid-scan\n        fig, ax = plt.subplots(figsize=(7, 3), layout=\"constrained\")\n        one = jnp.ones((1, 1, 1))\n        om_mid = OM0 + n_o * bins[1] / 2\n        for k_, off in enumerate((0.0, 0.25, 0.5)):\n            # 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\n            ylab = Y0 - DTY0 - (NK // 2 + off) * YSTEP\n            p_ = jnp.asarray([[ylab * np.sin(np.radians(om_mid)), ylab * np.cos(np.radians(om_mid)), 0.0]], jnp.float32)\n            for nb, ls in ((bins_fit[4], \"-\"), (0, \":\")):\n                idx_, wt_ = ix.system(one * 10.0, one * om_mid, one > 0, one, jnp.asarray([0]), p_, scan, *bins, nb)\n                idx_, wt_ = np.asarray(idx_).ravel(), np.asarray(wt_).ravel()\n                w_row = np.bincount(idx_[idx_ >= 0] % NK, weights=wt_[idx_ >= 0], minlength=NK)\n                ax.plot(dty_c - (DTY0 + NK // 2 * YSTEP), w_row, ls, marker=\"o\", ms=3, color=f\"C{k_}\",\n                        label=f\"voxel {off:g} step off a row\" + (\"\" if nb else \", 2-row model\"))\n        ax.set(xlim=(-3 * YSTEP, 3 * YSTEP), xlabel=\"dty - dty of the middle row\", ylabel=\"weight\",\n               title=f\"weights per row, at omega {om_mid:.4g}\")\n        ax.legend(fontsize=7)\n        plt.show()\n    t0 = time.time()\n    f0, cand = ix.candidates(H, pred, ring_j, pos_p, scan, bins_fit, n_cand, vb_all)\n    print(f\"candidates: {time.time() - t0:.0f} s\")",
   "id": "cell-26"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "if run_fit:  # anri.index.mlem, recording the deviance at every iteration and the map at some\n    norm = jnp.maximum(ix.backward(jnp.ones(H.shape[0], H.dtype), cand, pred, ring_j, pos_p, scan, bins_fit, vb), 1e-30)\n    fj = f0\n    dev, moved, f_prev, it_prev = [], {}, f0, 0\n    snap_at = sorted({i for i in (1, 2, 5, 10, 20, 50, 100, 200, 500) if i < n_iter} | {n_iter})\n    t0 = time.time()\n    for it in range(n_iter):\n        Af = ix.forward(fj, cand, pred, ring_j, pos_p, scan, bins_fit, H.shape[0], vb)\n        dev.append(float(ix.deviance(H, Af, censor)))\n        fj = fj * ix.backward(ix.censored_ratio(H, Af, censor), cand, pred, ring_j, pos_p, scan, bins_fit, vb) / norm\n        if it + 1 in snap_at:  # the share of each voxel's occupancy that moved since the last snapshot, per iteration\n            pa = f_prev / jnp.maximum(f_prev.sum(1, keepdims=True), 1e-30)\n            pb = fj / jnp.maximum(fj.sum(1, keepdims=True), 1e-30)\n            moved[it + 1] = np.full(NV, np.nan)\n            moved[it + 1][mask_v] = np.asarray(0.5 * jnp.abs(pb - pa).sum(1))[:n_fit] / (it + 1 - it_prev)\n            f_prev, it_prev = fj, it + 1\n            print(f\"  iteration {it + 1}: deviance {dev[-1]:.5g} ({time.time() - t0:.0f} s)\")\n    model = ix.forward(fj, cand, pred, ring_j, pos_p, scan, bins_fit, H.shape[0], vb)\n    f_v, cand_v = np.zeros((NV, n_cand), np.float32), np.zeros((NV, n_cand), np.int32)  # per voxel\n    f_v[mask_v], cand_v[mask_v] = np.asarray(fj)[:n_fit], np.asarray(cand)[:n_fit]  # outside the mask: empty",
   "id": "cell-27"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "**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.\n\n**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.",
   "id": "cell-28"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "if run_fit:\n    tot = f_v.sum(1)\n    occ = mask_v & (tot > occupied_cut * np.percentile(tot[mask_v], 99))\n\n    its = list(moved)[1:]  # the first snapshot is from the candidates' starting point\n    fig, ax = plt.subplots(1, 4, figsize=(18, 4), layout=\"constrained\")\n    ax[0].plot(np.arange(1, n_iter + 1), dev)\n    ax[0].set(xscale=\"log\", yscale=\"log\", xlabel=\"iteration\", ylabel=\"deviance\")\n    ax0 = ax[0].twinx()\n    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)\n    ax0.set_ylabel(\"relative decrease per iteration\", color=\"tab:red\")\n    if its:\n        for q, ls in ((50, \"-\"), (90, \"--\"), (99, \":\")):\n            ax[1].semilogy(its, [np.percentile(moved[i][occ], q) for i in its], \"o\" + ls, color=\"k\", ms=3,\n                           label=f\"{q}th percentile\")\n        ax[1].set(xscale=\"log\", xticks=its, xticklabels=[str(i) for i in its], xlabel=\"iteration\",\n                  ylabel=\"share moved per iteration\", title=\"occupancy moving between orientations\")\n        ax[1].minorticks_off()\n        ax[1].legend(fontsize=8)\n        im = ax[2].imshow(np.where(occ, moved[its[-1]], np.nan).reshape(NR, NR), origin=\"lower\", cmap=\"Reds\",\n                          norm=LogNorm(1e-5, 1e-1))\n        ax[2].set(title=f\"share moved per iteration, at {its[-1]}\", xticks=[], yticks=[])\n        fig.colorbar(im, ax=ax[2], shrink=0.8)\n    iv = np.flatnonzero(occ)\n    iv = iv[:: max(1, len(iv) // 200000)]  # a sample of the occupied voxels is enough\n    share = np.sort(f_v[iv], 1)[:, ::-1] / np.maximum(tot[iv, None], 1e-30)\n    for q, ls in ((50, \"-\"), (90, \"--\"), (99, \":\")):\n        ax[3].semilogy(np.arange(1, n_cand + 1), np.maximum(np.percentile(share, q, 0), 1e-6), ls, color=\"k\",\n                       label=f\"{q}th percentile\")\n    ax[3].set(xlabel=\"candidate rank\", ylabel=\"share of the voxel's occupancy\", title=f\"cand = {n_cand}\")\n    ax[3].legend(fontsize=8)\n    plt.show()\n    tail = share[:, -max(1, n_cand // 8) :].sum(1)\n    print(f\"the last {max(1, n_cand // 8)} candidates hold {np.median(tail) * 100:.2f}% of an occupied voxel's \"\n          f\"occupancy (median), {np.percentile(tail, 99) * 100:.2f}% (99th percentile)\")",
   "id": "cell-29"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 8. Model against data\n\nThe fitted histogram `A f` against the data.\n\n- **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.\n- **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.\n\nOnly η bins the fit models are used in the ratios (|sin η| above `etacut`).",
   "id": "cell-30"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "if run_fit:\n    M4 = np.asarray(model).reshape(n_rings_used, n_e, n_o, NK)\n    eta_c = -180.0 + (np.arange(n_e) + 0.5) * bins[0]\n    use_e = np.abs(np.sin(np.radians(eta_c))) > etacut_used\n    d_row, m_row = H4[:, use_e].sum((0, 1, 2)), M4[:, use_e].sum((0, 1, 2))\n    lit_rows = m_row > 0.05 * m_row.max()\n    row_ratio = np.where(lit_rows, d_row / np.maximum(m_row, 1e-30), np.nan)\n    d_ring, m_ring = H4[:, use_e].sum((1, 2, 3)), M4[:, use_e].sum((1, 2, 3))\n    fig, ax = plt.subplots(1, 2, figsize=(13, 3.6), layout=\"constrained\", width_ratios=[2, 1])\n    ax[0].plot(dty_c, row_ratio, \".-\", ms=3)\n    ax[0].axhline(1, color=\"k\", lw=0.5)\n    ax[0].set(xlabel=\"dty\", ylabel=\"measured / fitted\", title=\"per dty row\")\n    ax[1].bar(np.arange(n_rings_used), m_ring / d_ring)\n    ax[1].set(xlabel=\"ring\", ylabel=\"fitted / measured\", title=\"per ring\")\n    plt.show()\n    dv = row_ratio[lit_rows] - 1\n    print(f\"per row: rms {np.sqrt(np.mean(dv**2)):.3f}, row-to-row rms \"\n          f\"{np.sqrt(np.nanmean(np.diff(row_ratio) ** 2)):.3f}\")\n\n    # data, model and residual: one ring with the rows summed (eta against omega, top), and the sinogram of an eta\n    # window (bottom); every panel shares omega, each row its other axis\n    fig = figure(\"residual\", figsize=(16, 8))\n    ax = fig.subplots(2, 3, sharex=True)\n    for row in range(2):\n        for j in (1, 2):\n            ax[row, j].sharey(ax[row, 0])\n    om_ext = (OM0, OM0 + n_o * bins[1])\n    exts = ((*om_ext, -180, 180), (*om_ext, DTY0 - YSTEP / 2, DTY0 + (NK - 0.5) * YSTEP))\n    shapes = ((n_e, n_o), (NK, n_o))\n    ims = []\n    for row in range(2):\n        z = np.zeros(shapes[row])\n        ims.append([ax[row, j].imshow(z, aspect=\"auto\", origin=\"lower\", extent=exts[row], interpolation=\"nearest\",\n                                      cmap=MAGMA if j < 2 else \"RdBu_r\") for j in range(3)])\n        fig.colorbar(ims[row][2], ax=ax[row, 2], shrink=0.8)\n        ax[row, 0].set_ylabel((\"η (°)\", \"dty\")[row])\n        for j in range(3):\n            ax[row, j].set_xlabel(\"ω (°)\")\n    eta_lines = [[a.axhline(0.0, color=\"c\", lw=0.6) for _ in range(2)] for a in ax[0]]\n\n    def show_residual(ring, eta):\n        ie = (eta_c >= eta[0]) & (eta_c < eta[1])\n        pairs = ((H4[ring].sum(2), M4[ring].sum(2)), (H4[ring, ie].sum(0).T, M4[ring, ie].sum(0).T))\n        for row, (dd, mm) in enumerate(pairs):\n            top = max(dd.max(), 1.0)\n            for j, im in enumerate((dd, mm)):\n                ims[row][j].set_data(np.ma.masked_less_equal(im, 0))\n                ims[row][j].set_norm(LogNorm(top * 1e-2, top))\n            r = dd - mm\n            v = max(np.percentile(np.abs(r), 99.9), 1e-9)\n            ims[row][2].set_data(r)\n            ims[row][2].set_clim(-v, v)\n            for j, t in enumerate((\"data\", \"model\", \"data − model\")):\n                ax[row, j].set_title(f\"ring {ring}, {t}\" + (\"\" if row == 0 else f\", η {eta[0]}..{eta[1]}°\"))\n        for pair in eta_lines:\n            for ln, e_ in zip(pair, eta):\n                ln.set_ydata([e_, e_])\n        fig.canvas.draw_idle()\n\n    w_res = widgets.interactive(\n        show_residual,\n        ring=widgets.IntSlider(value=0, min=0, max=n_rings_used - 1, description=\"ring\", continuous_update=False),\n        eta=widgets.IntRangeSlider(value=(80, 100), min=-180, max=180, step=1, description=\"η window\",\n                                   continuous_update=False),\n    )\n    freeze(fig)\n    display(w_res, fig.canvas)",
   "id": "cell-31"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 9. Populations: `--occupied`, `--min-frac`\n\nEach 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).\n\n- `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.)\n- `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.",
   "id": "cell-32"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "if run_fit:\n    t0 = time.time()\n    frac, U_pop, spread, n_members = ix.populations(f_v, cand_v, U_kept, ops, 1.8 * step, p=N_POP)\n    print(f\"populations: {time.time() - t0:.0f} s\")\n    p99 = np.percentile(tot[mask_v], 99)  # within the mask\n    f2 = frac[mask_v & (tot > occupied_cut * p99), 1:].ravel()  # in the voxels occupied at the parameters' cut\n\n    fig = figure(\"populations\", figsize=(15, 9))\n    ax = fig.subplots(2, 3)\n    maps = (ax[0, 0], ax[0, 2], ax[1, 0], ax[1, 2])  # they share their axes; the histograms do not\n    for a in maps[1:]:\n        a.sharex(maps[0])\n        a.sharey(maps[0])\n    im = ax[0, 0].imshow(tot.reshape(NR, NR), origin=\"lower\", cmap=\"viridis\")\n    fig.colorbar(im, ax=ax[0, 0], shrink=0.8)\n    im_out = ax[0, 0].imshow(np.zeros((NR, NR, 4)), origin=\"lower\")  # the voxels not occupied, greyed out\n    ax[0, 0].set(title=\"total occupancy; grey: not occupied\")\n    ax[0, 1].hist(tot / p99, bins=np.linspace(0, max(1.5, (tot / p99).max()), 150), log=True)\n    v_occ = ax[0, 1].axvline(occupied_cut, color=\"tab:red\")\n    ax[0, 1].set(xlabel=\"total occupancy / 99th percentile\", ylabel=\"voxels\")\n    im_n = ax[0, 2].imshow(np.zeros((NR, NR)), origin=\"lower\", cmap=\"Blues\", vmin=0, vmax=N_POP)\n    ax[0, 2].set(title=\"populations per voxel\")\n    fig.colorbar(im_n, ax=ax[0, 2], shrink=0.8)\n    im_f2 = ax[1, 0].imshow(np.zeros((NR, NR)), origin=\"lower\", cmap=\"magma\", vmin=0, vmax=0.5)\n    ax[1, 0].set(title=\"fraction of the second population\")\n    fig.colorbar(im_f2, ax=ax[1, 0], shrink=0.8)\n    ax[1, 1].hist(f2[f2 > 0], bins=np.linspace(0.0, 0.5, 101), log=True)\n    v_frac = ax[1, 1].axvline(min_frac, color=\"tab:red\")\n    ax[1, 1].set(xlabel=\"fraction of populations 2 to 4\", ylabel=\"populations\")\n    im_sp = ax[1, 2].imshow(np.zeros((NR, NR)), origin=\"lower\", cmap=\"viridis\", vmin=0, vmax=max(2 * step, 1e-3))\n    ax[1, 2].set(title=\"spread of the main population (°)\")\n    fig.colorbar(im_sp, ax=ax[1, 2], shrink=0.8)\n    for a in maps:\n        a.set(xticks=[], yticks=[])\n\n    def show_pop(occupied_cut, min_frac):\n        occ = mask_v & (tot > occupied_cut * p99)\n        present = (frac >= min_frac) & occ[:, None]\n        present[:, 0] = occ\n        n_pop = present.sum(1)\n        grey = np.zeros((NV, 4))\n        grey[~occ] = (0.85, 0.85, 0.85, 0.85)\n        im_out.set_data(grey.reshape(NR, NR, 4))\n        im_n.set_data(np.where(occ, n_pop, np.nan).reshape(NR, NR))\n        im_f2.set_data(np.where(occ, np.where(present[:, 1], frac[:, 1], 0), np.nan).reshape(NR, NR))\n        im_sp.set_data(np.where(occ, spread[:, 0], np.nan).reshape(NR, NR))\n        v_occ.set_xdata([occupied_cut, occupied_cut])\n        v_frac.set_xdata([min_frac, min_frac])\n        fig.canvas.draw_idle()\n        print(f\"{occ.sum()} of {NV} voxels occupied; populations per occupied voxel: \"\n              + \", \".join(f\"{k}: {np.mean(n_pop[occ] == k) * 100:.1f}%\" for k in range(1, N_POP + 1)))\n\n    w_pop = widgets.interactive(\n        show_pop,\n        occupied_cut=widgets.FloatSlider(value=occupied_cut, min=0.0, max=1.0, step=0.01, description=\"occupied\",\n                                         continuous_update=False),\n        min_frac=widgets.FloatSlider(value=min_frac, min=0.0, max=0.5, step=0.01, description=\"min_frac\",\n                                     continuous_update=False),\n    )\n    freeze(fig)\n    display(w_pop, fig.canvas)",
   "id": "cell-33"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "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.",
   "id": "cell-34"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "if run_fit:\n    occupied_used, min_frac_used = w_pop.kwargs[\"occupied_cut\"], w_pop.kwargs[\"min_frac\"]\n    occ = mask_v & (tot > occupied_used * p99)\n    comp_main = np.full(NV, np.nan)\n    comp_main[occ] = ix.completeness_of(U_pop[occ, 0], delta, B, rings, geom, lit)\n    fig, ax = plt.subplots(1, 2, figsize=(11, 5), layout=\"constrained\", sharex=True, sharey=True)\n    im = ax[1].imshow(comp_main.reshape(NR, NR), origin=\"lower\", cmap=\"viridis\", vmin=0, vmax=1)\n    ax[1].set(title=\"completeness of the main population\", xticks=[], yticks=[])\n    fig.colorbar(im, ax=ax[1], shrink=0.8)\n    try:\n        maps = {\"UBI\": np.where(occ[:, None, None], np.linalg.inv(U_pop[:, 0] @ B), np.nan).reshape(NR, NR, 3, 3),\n                \"phase_ids\": np.where(occ, 0, -1).reshape(NR, NR)}\n        tmap = anri.io.tensormap_from_recon(maps, lpars, sg, phase_name, YSTEP)\n        tmap.get_ipf_maps()\n        ax[0].imshow(np.where((tmap.phase_ids[0] >= 0)[..., None], tmap.ipf_z[0], 1.0), origin=\"lower\")\n        ax[0].set(title=\"main population, IPF-Z\")\n    except ImportError as e:\n        ax[0].set(title=f\"no IPF map ({e.name} is not installed)\")\n    ax[0].set(xticks=[], yticks=[])\n    plt.show()",
   "id": "cell-35"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## The command line, and the results\n\nEvery 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.\n\nThen 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.",
   "id": "cell-36"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "# every value chosen above (sliders included), named as python -m anri.index names its options\nif not run_fit:\n    occupied_used, min_frac_used = occupied_cut, min_frac\noptions = {\n    \"analysisroot\": analysisroot, \"sample\": sample, \"dataset\": dataset, \"phase\": phase, \"parfile\": parfile,\n    \"rings\": n_rings_used, \"grid\": grid, \"max_chance\": w_grid.kwargs[\"max_chance\"], \"keep\": keep_used,\n    \"min_comp\": min_comp_used, \"prune\": prune_used, \"min_lr\": min_lr_used, \"cif\": cif, \"lit\": lit_used,\n    \"etacut\": etacut_used, \"tth_tol\": tth_tol_used, \"iter\": n_iter, \"cand\": n_cand, \"coarse\": coarse,\n    \"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,\n    \"y0\": None if ds[\"y0\"] is not None and np.isclose(Y0, ds[\"y0\"]) else Y0, \"gridstep\": gridstep,\n    \"outdir\": outdir, \"n_cpu\": n_cpu, \"check\": False,\n}\ncmd = [\"python\", \"-m\", \"anri.index\", analysisroot, sample, dataset]\nfor k, v in options.items():\n    if v is None or v is False or k in (\"analysisroot\", \"sample\", \"dataset\"):\n        continue\n    cmd += [\"--\" + k.replace(\"_\", \"-\")] + ([] if v is True else [f\"{v:.10g}\" if isinstance(v, float) else str(v)])\ncommand = shlex.join(cmd)\nprint(command)",
   "id": "cell-37"
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": "if run_fit:  # the results, written as python -m anri.index writes them\n    import datetime\n\n    from anri.index.__main__ import _git_commit, save_results, write_toml\n\n    occ = mask_v & (tot > occupied_used * p99)\n    present = (frac >= min_frac_used) & occ[:, None]\n    present[:, 0] = occ  # the main population always\n    comp_pop = np.zeros(frac.shape, np.float32)\n    comp_pop[present] = ix.completeness_of(U_pop[present], delta, B, rings, geom, lit)\n    tag = os.path.join(outdir, f\"{dsname}_index\")\n    results = {\"f\": f_v, \"cand\": cand_v, \"U\": U_kept, \"comp\": comp[kept], \"frac\": frac, \"U_pop\": U_pop,\n               \"spread\": spread, \"n\": n_members, \"comp_pop\": comp_pop, \"occupied\": occ, \"present\": present,\n               \"pos\": pos, \"grid_step\": step, \"delta\": delta, \"row_ratio\": row_ratio, \"row_data\": d_row,\n               \"row_model\": m_row, \"y0\": Y0, \"mask\": mask_v}\n    save_results(tag, results, B, lpars, sg, phase_name, YSTEP)\n    run_info = {\n        \"run\": {\"command\": command, \"notebook\": \"indexing_parameters.ipynb\", \"anri_version\": anri.VERSION,\n                **_git_commit(), \"cwd\": os.getcwd(),\n                \"date\": datetime.datetime.now(datetime.timezone.utc).astimezone().isoformat(timespec=\"seconds\")},\n        \"options\": options,\n        \"resolved\": {\"dataset\": dsfile, \"sparsefile\": sparsefile, \"parfile\": pfile, \"phase\": phase_name,\n                     \"lattice\": list(lpars), \"space_group\": sg, \"wavelength\": WL, \"y0\": Y0, \"voxel_size\": YSTEP,\n                     \"voxels\": NR, \"omega_step\": OSTEP, \"grid_step\": step, \"grid_auto\": grid is None,\n                     \"grid_worst_case_deg\": delta, \"rings_tth_deg\": list(rings[\"tth\"]), \"tth_tol_deg\": list(tth_tol_r),\n                     \"n_hkls\": len(rings[\"hkls\"]), \"spatial_correction\": spatial_src, \"beam_fwhm\": beam_used},\n        \"results\": {\"chance_completeness\": info[\"chance\"], \"min_comp_used\": cut, \"orientations_kept\": len(U_kept),\n                    \"voxels_occupied\": int(occ.sum()), \"seconds\": time.time() - T_START},\n    }\n    write_toml(f\"{tag}_params.toml\", run_info)\n    print(f\"-> {tag}_params.toml\")\nelse:\n    print(\"run_fit = False: nothing indexed here, so nothing written; run the command above\")",
   "id": "cell-38"
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  },
  "nbsphinx": {
   "execute": "never"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}