Measured frames against a map’s forward model#

Give this notebook an ImageD11 dataset and a map of it (a TensorMap from the indexer, ImageD11 or a refinement), and it shows the measured frames next to the frames the map predicts, for any row (dty) and frame (ω), with sliders. It is the most direct check of a map: a map that is right reproduces the spots where they are, with their shapes and intensities, frame by frame.

The model simulates the dataset’s own scan (python -m anri.simulate prepare --dataset): every frame’s ω and dty from the DataSet, its detector and geometry, and the map’s phases. Rows are rendered when first shown (a few seconds on a GPU, after the first compile) and kept. The measured side is the dataset’s sparse file (what the indexer sees, above the segmentation cut) or the raw detector frames from the master file; both are moved by the DataSet’s detector distortion maps to where an ideal detector would see them, as the model is. Each row’s model is scaled onto the data by one factor (least squares over the pixels the data hold or the scaled model puts above the cut), so the overall brightness of the map does not matter.

It needs the dataset (not distributed with anri), so it is stored without outputs. With interactive = True (the default) the sliders redraw the plots; it needs ipywidgets, ipympl and matplotlib.

[ ]:

# Parameters

dataset_file = None  # an ImageD11 <sample>_<dataset>_dataset.h5 (required)
map_file = None  # a TensorMap (.h5) or python -m anri.index's _entries.npz (required)
remap = None  # {old prefix: new prefix} for the paths the DataSet records, e.g. data mounted elsewhere
pars_file = None  # geometry (.par or pars.json; None: the DataSet's)
y0 = None  # where the rotation axis is in the beam (None: the DataSet's); use the y0 the map was made with
cif = None  # --cif: one path, or {phase name: path}
optics = "al-crl"  # ID11's focusing optics: "al-crl" or "si-nanolens"
window, max_frames = (3, 7, 7), 31  # as python -m anri.simulate render
cut = 1.0  # the segmentation's threshold: model pixels below it do not count against the data in the scale fit
outdir = None  # where the prepare stage file goes (None: a temporary folder)
n_cpu = 4

# this notebook only
interactive = True  # sliders (ipympl); False: static figures at the starting values
row = None  # the row shown first (None: the middle row)
frame = None  # the frame shown first (None: the model's brightest in that row)
source = "sparse"  # "sparse" (the segmented pixels) or "raw" (the detector frames from the master file)

[ ]:

import anri.utils

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

import os
import tempfile
import time
from types import SimpleNamespace

import ipywidgets as widgets
import numpy as np
from IPython import get_ipython
from IPython.display import display
from matplotlib import pyplot as plt

import anri.io
import anri.simulate as sim
from anri.simulate import diagnostics as dg
from anri.simulate import plots as pl

get_ipython().run_line_magic("matplotlib", "widget" if interactive else "inline")
T0 = time.time()


def log(msg):
    print(f"[{time.time() - T0:7.1f} s] {msg}", flush=True)


def figure(name, **kwargs):
    plt.close(name)
    with plt.ioff():
        fig = plt.figure(name, layout="constrained", **kwargs)
    if interactive:
        fig.canvas.header_visible = False
    return fig


def interact(fig, update, **controls):
    # sliders that call a plot's update, which redraws it and returns a summary; their values are in .kwargs
    if not interactive:
        values = {k: c.value for k, c in controls.items()}
        print(update(**values))
        display(fig)
        return SimpleNamespace(kwargs=values)
    w = widgets.interactive(lambda **kw: print(update(**kw)), **controls)
    display(w, fig.canvas)
    return w

The dataset and the model#

The DataSet gives the scan, the sparse and raw files and the detector distortion; stage_prepare(..., dataset=...) prepares the map to simulate exactly that scan. A map made by python -m anri.index is centred on the y0 its fit used (python -m anri.index show <file> prints it): give that y0, or the model’s spots move with the rotation.

[ ]:

if dataset_file is None or map_file is None:
    raise ValueError("set dataset_file and map_file in the parameters")
ds = anri.io.imaged11.read_dataset(dataset_file, remap)
spatial = anri.io.imaged11.read_spatial(ds)
outdir = outdir or tempfile.mkdtemp()
prepare_file = sim.stages.stage_prepare(map_file, pars_file, os.path.join(outdir, "compare_prepare.h5"), dataset=dataset_file,
                                remap=remap, y0=y0, cif=cif, optics=optics, log=log)
p = sim.stages.load_prepare(prepare_file)
n_rows, n_frames = p["omega"].shape
print(f"{n_rows} rows x {n_frames} frames; detector distortion: {'yes' if spatial is not None else 'none'}")

Frames#

Data, model and both overlaid (data red, model cyan: white where they agree) for frames frames around frame of row, summed. The panels share their zoom: use the toolbar’s zoom on one, and the others follow. What to look for:

  • Positions: spots in the same place at the same frame. A shift that grows with the distance from the rotation axis points to y0; one that rotates with ω, to the map’s position or orientation.

  • Red spots (data only): grains or phases the map lacks, or orientations it has wrong.

  • Cyan spots (model only): orientations or voxels the data do not support.

  • Shapes: a measured spot wider or longer than the model’s is spread the model lacks (orientation spread within the voxels, the beam’s convergence, the detector’s point spread); smeared along the ring it is orientation, across it strain.

  • Intensities: the row’s scale factor is printed; spot-to-spot differences within a row point to structure factors, absorption or the voxels’ densities.

[ ]:

rows_seen = {}  # row: (observed, model, scale), filled as rows are shown


def row_data(r):
    if r not in rows_seen:
        t = time.time()
        obs, mod = dg.observed_row(ds, r, spatial), dg.model_row(p, r, window, max_frames)
        rows_seen[r] = (obs, mod, dg.row_scale(obs, mod, cut))
        log(f"row {r}: {len(obs['value'])} measured and {len(mod['value'])} model pixels ({time.time() - t:.1f} s)")
    return rows_seen[r]


state = {"source": source}


def images(r, frames):
    obs, mod, s = row_data(r)
    data = dg.raw_image(ds, r, frames, spatial) if state["source"] == "raw" else dg.frames_image(obs, frames)
    return data, s * dg.frames_image(mod, frames), f"{state['source']}; model scaled by {s:.3g}"


row = n_rows // 2 if row is None else row
if frame is None:  # the model's brightest pixel in that row
    _, mod0, _ = row_data(row)
    frame = int(mod0["frame"][np.argmax(mod0["value"])])
fig = figure("frames", figsize=(18, 6.5))
update = pl.plot_frame_compare(fig, images, p["det_shape"])


def show_frames(row, frame, frames, source):
    state["source"] = source
    return update(row, frame, frames)


w = interact(fig, show_frames,
             row=widgets.IntSlider(value=row, min=0, max=n_rows - 1, description="row", continuous_update=False),
             frame=widgets.IntSlider(value=frame, min=0, max=n_frames - 1, description="frame",
                                     continuous_update=False),
             frames=widgets.IntSlider(value=3, min=1, max=31, step=2, description="frames", continuous_update=False),
             source=widgets.Dropdown(options=["sparse", "raw"], value=source, description="data"))

Zoom on the model’s brightest spot in the row shown first (static: the figure above, cropped).

[ ]:

obs, mod, s = row_data(row)
k = int(np.argmax(mod["value"]))
r0, c0 = divmod(int(mod["pixel"][k]), p["det_shape"][1])
fig = figure("zoom", figsize=(18, 6.2))
zoom = pl.plot_frame_compare(fig, images, p["det_shape"])
print(zoom(row, int(mod["frame"][k]), 3))
for a in fig.axes[:3]:
    a.set_xlim(c0 - 60, c0 + 60)
    a.set_ylim(r0 + 60, r0 - 60)
display(fig.canvas if interactive else fig)