Wilson plots generated from sparse datasets

The input data comes from the sparsification of the zenodo dataset available at: https://zenodo.org/record/3834335

filename="/mnt/data/jungfrau/sparse_512_3_0_2_fit_q2.h5"
import os, time
os.chdir(os.path.dirname(filename))
start_time = time.perf_counter()

---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
Cell In[1], line 3
      1 filename="/mnt/data/jungfrau/sparse_512_3_0_2_fit_q2.h5"
      2 import os, time
----> 3 os.chdir(os.path.dirname(filename))
      4 start_time = time.perf_counter()

FileNotFoundError: [Errno 2] No such file or directory: '/mnt/data/jungfrau'

%matplotlib nbagg
import json, pyFAI
import hdf5plugin
import h5py
import numpy
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
from matplotlib.pyplot import subplots

root = h5py.File(filename, "r")
de = root.attrs["default"]
entry = root[de]
nx_data = entry[entry.attrs["default"]]
dict(nx_data)

radius = nx_data["radius"]
I_bg = nx_data["background_avg"]
resolution = 10/numpy.sqrt(radius)
resolution[:4],resolution[-4:]

fig, ax = subplots()
ax.plot(radius, numpy.log(I_bg[0]))
ax.set_ylabel("$log(\overline{I_{bg}})$")
ax.set_xlabel("$d^*² (nm^{-2})$")
ax.set_title("Wilson-like plot on frame #0")
pass

frame_ptr = nx_data["frame_ptr"]
npix = frame_ptr[1:]-frame_ptr[:-1]
fig,ax = subplots()
ax.plot(npix)
ax.set_xlabel("Frame number")
ax.set_ylabel("Number of recorded pixels")
ax.set_title("Number of peaks per frame")
frame = numpy.argmax(npix)
print(f"Frame index with largest number of pixel recorded: {frame}")

mask = nx_data["mask"]
shape = mask.shape
size = numpy.prod(shape)

nticks = 7
fig, ax = subplots()
ax.plot(radius, numpy.log(I_bg[frame]), label="Background")
tick_value = numpy.linspace(radius[0], radius[-1], nticks)
ax.set_xticks(tick_value)
tick_label = [f"{10/numpy.sqrt(i) if i>0 else numpy.nan:6.4f}" for i in tick_value]
ax.set_xticklabels(tick_label)
ax.set_xlabel("Resolution ($\mathrm{\AA}$)")
ax.set_ylabel("$log(\overline{I_{bg}})$")
ax.set_title(f"Wilson-like plot on frame #{frame}")
pass

indexes = nx_data["index"][frame_ptr[frame]:frame_ptr[frame+1]]
intensities= nx_data["intensity"][frame_ptr[frame]:frame_ptr[frame+1]]
signal = numpy.zeros(size)
norm = numpy.zeros(size)
signal[indexes] = intensities
norm[indexes] = 1

%%time
sparsify = nx_data.parent["sparsify"]
config = json.loads(sparsify["configuration/data"][()])
ai = pyFAI.load(config["geometry"])
engine = ai.setup_sparse_integrator(shape=ai.detector.shape,
                      npt=resolution.size,
                      mask=numpy.logical_not(numpy.isfinite(mask)),
                      unit="d*2_nm^-2",
                      split='no',
                      algo='CSR',
                      scale=True)

from scipy.sparse import csr as csr_module
csr = csr_module.csr_matrix(engine.lut, shape=(resolution.size, size))
csr

nticks = 7
fig, ax = subplots()
ax.plot(radius, numpy.log(I_bg[frame]), label="Background")
tick_value = numpy.linspace(radius[0], radius[-1], nticks)
ax.set_xticks(tick_value)
tick_label = [f"{10/numpy.sqrt(i) if i>0 else numpy.nan:6.4f}" for i in tick_value]
ax.set_xticklabels(tick_label)
ax.set_xlabel("Resolution ($\mathrm{\AA}$)")
#ax.plot(radius, numpy.log(csr.dot(signal)/csr.dot(norm)-I_bg[frame]), label="Peaks-Background")
ax.plot(radius, numpy.log(csr.dot(signal)/csr.dot(norm)), label="Peaks")
ax.set_ylabel("$log(\overline{I})$")
ax.legend()
ax.set_title(f"Wilson-like plot for frame #{frame}")
pass

def update_plot(frame):
    ax.set_title(f"Wilson-like plot for frame #{frame}")
    ax.lines[0].set_data(radius, numpy.log(I_bg[frame]))
    signal = numpy.zeros(size)
    norm = numpy.zeros(size)
    indexes = nx_data["index"][frame_ptr[frame]:frame_ptr[frame+1]]
    intensities= nx_data["intensity"][frame_ptr[frame]:frame_ptr[frame+1]]
    signal[indexes] = intensities
    norm[indexes] = 1
    #ax.lines[0].set_data(radius, numpy.log(csr.dot(signal)/csr.dot(norm)-I_bg[frame]))
    ax.lines[1].set_data(radius, numpy.log(csr.dot(signal)/csr.dot(norm)))
    return frame
interact(update_plot, frame=widgets.IntSlider(min=0, max=len(npix)-1, step=1, value=1));

root.close()
print(f"Total execution time: {time.perf_counter()-start_time:6.3f} s")