<jupyter_start><jupyter_text>Analysing simulated serial crystallography data<jupyter_code>%matplotlib inline
import io
import h5py
from h5glance import H5Glance
import matplotlib.pyplot as plt
import numpy as np
import os
import os.path as osp
import pandas as pd
from scipy.stats import norm

from CrystFEL_Jupyter_utilities.stream_read import search_crystals_parameters
from cfel_geom import CrystfelGeom<jupyter_output><empty_output><jupyter_text>Input data<jupyter_code>data_dir = "../../../data/simulation/"
data_file = osp.join(data_dir, "xstal/xstal_diffr200.h5")
data_f = h5py.File(data_file, 'r')
H5Glance(data_f['data/0000001'])<jupyter_output><empty_output><jupyter_text>Here's one frame from the simulation:<jupyter_code>geom = CrystfelGeom.from_crystfel_geom('1panel.geom')
frame = data_f['data/0000001/data'][:]
fig, ax = plt.subplots(figsize=(12, 10))
geom.plot_data_fast(frame[np.newaxis], axis_units='m', ax=ax, vmin=0, vmax=10000);<jupyter_output><empty_output><jupyter_text>Running CrystFEL<jupyter_code>with open('test.lst', 'w') as f:
    f.write(f'{data_file}\n')
import modules
modules.load('exfel', 'exfel_crystfel')
!time indexamajig -i test.lst -g 1panel.geom -o out.stream --peaks=zaef -j80<jupyter_output>Unrecognised field 'px'
Unrecognised field 'py'
Unrecognised field 'pix_width'
Unrecognised field 'd'
WARNING: You did not specify --int-radius.
WARNING: I will use the default values, which are probably not appropriate for your patterns.
No indexing methods specified.  I will try to automatically detect the available methods.
To disable auto-detection of indexing methods, specify which methods to use with --indexing=<methods>.
Use --indexing=none to disable indexing and integration.
No reference unit cell provided.
WARNING: Forcing --no-check-cell because reference unit cell parameters were not given.
WARNING: Forcing all indexing methods to use "-nocell", because reference cell parameters were not given.
List of indexing methods:
   0: asdf-nolatt-nocell        (asdf - no prior information)
Indexing parameters:
                  Check unit cell parameters: off
                        Check peak alignment: on
                   Refine indexing solutions: on
 Multi-lattice indexing ("d[...]<jupyter_text>Identified peaks<jupyter_code>def read_stream_table(fh, end_marker):
    buf = io.StringIO()
    for line in fh:
        if line.startswith(end_marker):
            buf.seek(0)
            return pd.read_fwf(buf)
        
        buf.write(line)
def read_peak_lists(stream_file):
    stream_chunks = {}

    with open(stream_file, 'r') as f:
        in_peaks_table = False
        chunk_data = {}
        for line in f:
            if '-- Begin chunk --' in line:
                chunk_data = {}
            elif '-- End chunk --' in line:
                if 'event' in chunk_data:
                    stream_chunks[chunk_data['event']] = chunk_data

            elif line.startswith('Event:'):
                chunk_data['event'] = line.split(':', 1)[1].strip()

            elif line.startswith('Peaks from peak search'):
                chunk_data['peaks'] = read_stream_table(f, 'End of peak list')
    
    return stream_chunks
peaks_found = read_peak_lists('out.stream')
peak_tbl = peaks_found['0000001//']['peaks']
peak_tbl
module_no = np.zeros(len(peak_tbl), dtype=np.uint32)
peak_coords = geom.data_coords_to_positions(module_no, peak_tbl['ss/px'], peak_tbl['fs/px'])
peak_x = peak_coords[:, 0]
peak_y = peak_coords[:, 1]
fig, ax = plt.subplots(figsize=(12,10))
geom.plot_data_fast(frame[np.newaxis], axis_units='m', ax=ax, vmin=0, vmax=10000)
ax.scatter(peak_x, peak_y, s=80, marker='o', facecolor='none', edgecolor='0.7');<jupyter_output><empty_output><jupyter_text>Examine unit cell parameters<jupyter_code>params_lst = search_crystals_parameters('out.stream')
params = pd.DataFrame(params_lst).drop(columns='name').rename(columns={'alfa': 'alpha'})
params
fig, axes = plt.subplots(2, 3, figsize=(12, 8))

for name, ax in zip('abc', axes[0]):
    ax.set_title(name)
    ax.set_xlabel("Å")
    ax.hist(params[name])

for name, ax in zip(('alpha', 'beta', 'gamma'), axes[1]):
    ax.set_title(f"$\\{name}$")
    ax.set_xlabel('degrees')
    ax.hist(params[name])

fig.tight_layout()<jupyter_output><empty_output><jupyter_text>Estimate unit cell parametersWe can estimate the value of each parameter by fitting a normal distribution to it:<jupyter_code>def plot_param(name, ax, units):
    data = params[name]
    ax.set_xlabel(units)
    ax.hist(data, density=True, color='0.8')
    
    μ, σ = norm.fit(data)
    ax.text(0.05, 0.95, f"μ = {μ:.04}, σ = {σ:.03}",
            transform=ax.transAxes, fontsize=12)
    x = np.linspace(data.min(), data.max())
    ax.plot(x, norm.pdf(x, μ, σ), linestyle='--')
fig, axes = plt.subplots(2, 3, figsize=(12, 8))

for name, ax in zip('abc', axes[0]):
    plot_param(name, ax, "Å")
    ax.set_title(name)
    

for name, ax in zip(('alpha', 'beta', 'gamma'), axes[1]):
    plot_param(name, ax, 'degrees')
    ax.set_title(f"$\\{name}$")

fig.tight_layout()<jupyter_output><empty_output><jupyter_text>Simulated extreme noiseLet's repeat it with a simulation with a small number of extreme values randomly added to the results.<jupyter_code>data_file = osp.join(data_dir, "superlarge_data_reduced_2.h5")
with open('extreme_noise.lst', 'w') as f:
    f.write(f'{data_file}\n')
!time indexamajig -i extreme_noise.lst -g 1panel.geom -o out-xn.stream --peaks=zaef -j80
peaks_found = read_peak_lists('out-xn.stream')
peak_tbl = peaks_found['0000001//']['peaks']
module_no = np.zeros(len(peak_tbl), dtype=np.uint32)
peak_coords = geom.data_coords_to_positions(module_no, peak_tbl['ss/px'], peak_tbl['fs/px'])
peak_x = peak_coords[:, 0]
peak_y = peak_coords[:, 1]
with h5py.File(data_file, 'r') as data_f:
    frame = data_f['data/0000001/data'][:]

fig, ax = plt.subplots(figsize=(12,10))
geom.plot_data_fast(frame[np.newaxis], axis_units='m', ax=ax, vmin=0, vmax=10000)
ax.scatter(peak_x, peak_y, s=80, marker='o', facecolor='none', edgecolor='0.7');
params_lst = search_crystals_parameters('out-xn.stream')
params = pd.DataFrame(params_lst).drop(columns='name').rename(columns={'alfa': 'alpha'})
fig, axes = plt.subplots(2, 3, figsize=(12, 8))

for name, ax in zip('abc', axes[0]):
    ax.set_title(name)
    ax.set_xlabel("Å")
    ax.hist(params[name])

for name, ax in zip(('alpha', 'beta', 'gamma'), axes[1]):
    ax.set_title(f"$\\{name}$")
    ax.set_xlabel('degrees')
    ax.hist(params[name])

fig.tight_layout()<jupyter_output><empty_output><jupyter_text>AGIPD multi-module layout<jupyter_code>data_file = osp.join(data_dir, "panoscDEMO/diffr200.h5")
geom_file = "agipd_simple_2d.geom"<jupyter_output><empty_output><jupyter_text>Make a virtual dataset to present data as a 3D array.<jupyter_code>vlayout = h5py.VirtualLayout((200, 8192, 128), dtype='f4')

frames_dir = osp.join(data_dir, "panoscDEMO/diffr200")

for i, filename in enumerate(sorted(os.listdir(frames_dir))):
    with h5py.File(os.path.join(frames_dir, filename), 'r') as f:
        vsrc = h5py.VirtualSource(f['data/data'])
    
    vlayout[i] = vsrc


with h5py.File('vdata.h5', 'w') as f:
    f.create_virtual_dataset('data/data', vlayout)
data_file = 'vdata.h5'
data_f = h5py.File(data_file, "r")<jupyter_output><empty_output><jupyter_text>Here's one frame from the simulation:<jupyter_code>from extra_geom import AGIPD_1MGeometry
geom = AGIPD_1MGeometry.from_crystfel_geom(geom_file)

frame = data_f['data/data'][0].reshape(16, 512, 128)
fig, ax = plt.subplots(figsize=(12, 10))
geom.plot_data_fast(frame, axis_units='m', ax=ax, vmax=50);
with open('multi-panel.lst', 'w') as f:
    f.write(f'{data_file}\n')
!time indexamajig -i multi-panel.lst -g agipd_simple_2d.geom -o multi-panel.stream \
    --peaks=peakfinder8 --threshold 20 -j80
peaks_found = read_peak_lists('multi-panel.stream')
peak_tbl = peaks_found['//0']['peaks']
module_no = np.zeros(len(peak_tbl), dtype=np.uint32)
peak_coords = geom.data_coords_to_positions(module_no, peak_tbl['ss/px'], peak_tbl['fs/px'])
peak_x = peak_coords[:, 0]
peak_y = peak_coords[:, 1]
fig, ax = plt.subplots(figsize=(12,10))
geom.plot_data_fast(frame, axis_units='m', ax=ax, vmin=0, vmax=100)
ax.scatter(peak_x, peak_y, s=80, marker='o', facecolor='none', edgecolor='0.7');
params_lst = search_crystals_parameters('multi-panel.stream')
params = pd.DataFrame(params_lst).drop(columns='name').rename(columns={'alfa': 'alpha'})
fig, axes = plt.subplots(2, 3, figsize=(12, 8))

for name, ax in zip('abc', axes[0]):
    ax.set_title(name)
    ax.set_xlabel("Å")
    ax.hist(params[name])

for name, ax in zip(('alpha', 'beta', 'gamma'), axes[1]):
    ax.set_title(f"$\\{name}$")
    ax.set_xlabel('degrees')
    ax.hist(params[name])

fig.tight_layout()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15279672578444747, 'fraction_numerical': 0.024215552523874487, 'mean_word_length': 4.802770448548813, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 45, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '652907', 'license': 'MIT', 'n_tokens_mistral': 3474, 'n_tokens_neox': 3274, 'n_words': 680, 'path': 'src/controller/panoscDEMO/CrystFEL indexing.ipynb', 'repo_name': 'ejcjason/crystalProject', 'size': '8796'}