Implementation of PeakFinder8 on GPU¶
The peakfinder8 is the core algorithm for assessing the quality of a single frame in serial crystallography and was initially implemented in C++ within the cheetah [1]
This algorithm is called peakfinder8 because it consits of 8 subsequent steps perfromed on evry single frame:
perfrom the azimuthal integration with uncertainety propagation
discard pixels which differ by more than N-sigma from the mean and cycle to 1 about 3 to 5 times
pick all pixels with I > mean + min(N*sigma, noise)
such pixel is a peak if it is the maximum of the 3x3 or 5x5 patch and there are connected pixels in the patch with their intensity above the previous threshold.
subtract background and sum the signal over the patch
return the index of the peak, the integrated signal and the center of mass of the peak
exclude neighboring peaks (un-implemented)
Validate the frame if there are enough peaks found.
There is a attempt to implement peakfinder8 on GPU within the pyFAI. The steps 1+2 correspond to the sigma-clipping algorithm and enforce an azimuthal, normal distribution for the background.
This tutorial demontrates how peak-finding can be called from Jupyter notebooks and what are the performances expected. Finally, the performances will be compared with the reference implementation.
[1] A. Barty, R. A. Kirian, F. R. N. C. Maia, M. Hantke, C. H. Yoon, T. A. White, and H. N. Chapman, “Cheetah: software for high-throughput reduction and analysis of serial femtosecond x-ray diffraction data”, J Appl Crystallogr, vol. 47, pp. 1118-1131 (2014)
[1]:
%matplotlib inline
# use `widget` for better user experience; `inline` is for documentation generation
[2]:
import os
import sys
import shutil
import posixpath
import numpy
import glob
from matplotlib.pylab import subplots
import fabio
import pyFAI, pyFAI.azimuthalIntegrator
from pyFAI.gui import jupyter
from pyFAI import units
import pyopencl
from pyFAI.opencl.peak_finder import OCL_PeakFinder
from pyFAI.test.utilstest import UtilsTest
from pyFAI.containers import ErrorModel
import time
start_time = time.perf_counter()
os.environ["PYOPENCL_COMPILER_OUTPUT"] = "1"
[3]:
fimg = fabio.open(UtilsTest.getimage("Pilatus6M.cbf"))
mask = numpy.logical_or(fimg.data>65000, fimg.data<0)
print(f"Number of masked pixels: {mask.sum()}")
Number of masked pixels: 527055
[4]:
det = pyFAI.detector_factory("Pilatus6M")
det.mask = mask
[5]:
dimg = fimg.data.copy()
dimg[mask] = 0
fig,ax = subplots(1, 2, figsize=(10,5))
jupyter.display(dimg, ax=ax[0])
jupyter.display(dimg, ax=ax[1])
ax[1].set_xlim(1500, 1800)
ax[1].set_ylim(850, 1020)
pass
[6]:
ponifile = UtilsTest.getimage("Pilatus6M.poni")
ai = pyFAI.load(ponifile)
ai.detector = det
print(ai)
Detector Pilatus 6M PixelSize= 1.720e-04, 1.720e-04 m BottomRight (3)
Wavelength= 1.033200e-10 m
SampleDetDist= 3.000000e-01 m PONI= 2.254060e-01, 2.285880e-01 m rot1=0.000000 rot2=0.000000 rot3=0.000000 rad
DirectBeamDist= 300.000 mm Center: x=1329.000, y=1310.500 pix Tilt= 0.000° tiltPlanRotation= 0.000° 𝛌= 1.033Å
[7]:
kwargs = {"data": fimg.data,
"npt":1000,
"method": ("no", "csr", "opencl"),
"polarization_factor": 0.99,
"unit":"r_mm", }
ax = jupyter.plot1d(ai.integrate1d(**kwargs))
ax.errorbar(*ai.sigma_clip_ng(error_model="azimuthal", **kwargs), label="sigma-clip")
_=ax.legend()
[8]:
ctx = pyopencl.create_some_context(interactive=False)
print(f"Using {ctx.devices[0].name}")
Using NVIDIA RTX A2000 12GB
[9]:
unit = units.to_unit("r_mm")
image_size = det.shape[0] * det.shape[1]
integrator = ai.setup_sparse_integrator(ai.detector.shape, 1000, mask=mask, unit=unit, split="no", algo="CSR", scale=False)
polarization = ai._cached_array["last_polarization"]
pf = OCL_PeakFinder(integrator.lut,
image_size=image_size,
bin_centers=integrator.bin_centers,
radius=ai._cached_array[unit.name.split("_")[0] + "_center"],
mask=mask,
ctx=ctx,
# block_size=512,
unit=unit)
kwargs = {"data":fimg.data,
"error_model": ErrorModel.parse("azimuthal"),
"polarization":polarization.array,
"polarization_checksum": polarization.checksum}
print(f"Number of high intensity pixels at stage #3:\t{pf.count_intense(**kwargs ,cycle=5, cutoff_pick=3.0)}\n\
Number of peaks identified at stage #6:\t\t{pf._count_peak(**kwargs, cycle=5, cutoff_peak=3.0)}")
Number of high intensity pixels at stage #3: 19464
Number of peaks identified at stage #6: 340
[10]:
from silx.opencl.processing import ProfileDescription, EventDescription
def average_opencl_runtime(events):
stats = {}
total_time = 0.0
for e in events:
if isinstance(e, ProfileDescription):
name = e[0]
t0 = e[1]
t1 = e[2]
elif isinstance(e, EventDescription) or "__len__" in dir(e) and len(e) == 2:
name = e[0]
pr = e[1].profile
t0 = pr.start
t1 = pr.end
else:
name = "?"
t0 = e.profile.start
t1 = e.profile.end
et = 1e-6 * (t1 - t0)
total_time += et
if name in stats:
stats[name].append(et)
else:
stats[name] = [et]
return total_time/max(len(stats[i]) for i in stats)
[11]:
# Performance measurement of the pixel recording (stage 1->3)
pf.reset_log()
pf.set_profiling(True)
timeit_count_intense = %timeit -o pf.count(**kwargs, cycle=3, cutoff_pick=3, noise=1)
print("\n".join(pf.log_profile(True)))
print(f"Overhead due to Python: {(0.001*average_opencl_runtime(pf.events)/timeit_count_intense.average-1.0)*-100:.1f}%")
pf.set_profiling(False)
9.51 ms ± 40.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
OpenCL kernel profiling statistics in milliseconds for: OCL_PeakFinder
Kernel name (count): min median max mean std
copy H->D image_raw ( 811): 2.778 2.864 4.162 2.919 0.163
memset_ng ( 811): 0.003 0.003 0.328 0.007 0.013
corrections4a ( 811): 0.601 0.604 0.975 0.609 0.020
csr_sigma_clip4 ( 811): 3.993 4.317 6.967 4.440 0.441
memset counter ( 811): 0.003 0.003 0.011 0.003 0.001
find_intense ( 811): 0.835 0.845 0.905 0.850 0.015
copy D->H counter ( 811): 0.001 0.001 0.003 0.001 0.000
________________________________________________________________________________
Total OpenCL execution time : 7160.324ms
Overhead due to Python: 7.2%
The overhead from calling OpenCL from Python is as low as 8%
[12]:
# Performance measurement of the pixel recording (stage 1->6)
pf.reset_log()
pf.set_profiling(True)
timeit_gpu = %timeit -o pf.peakfinder8(**kwargs, cycle=3, cutoff_peak=3, noise=1, connected=3, patch_size=3)
print("\n".join(pf.log_profile(True)))
print(f"Overhead due to Python: {(0.001*average_opencl_runtime(pf.events)/timeit_gpu.average-1.0)*-100:.1f}%")
pf.set_profiling(False)
9.71 ms ± 84 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
OpenCL kernel profiling statistics in milliseconds for: OCL_PeakFinder
Kernel name (count): min median max mean std
copy H->D image_raw ( 811): 2.775 2.879 4.333 2.942 0.189
memset_ng ( 811): 0.003 0.003 0.017 0.005 0.004
corrections4a ( 811): 0.597 0.604 1.018 0.609 0.019
csr_sigma_clip4 ( 811): 4.002 4.392 5.757 4.473 0.380
memset counter ( 811): 0.003 0.003 0.004 0.003 0.000
peakfinder ( 811): 0.843 0.860 0.925 0.862 0.011
copy D->H counter ( 811): 0.001 0.002 0.002 0.002 0.000
copy D->H peak positions ( 811): 0.001 0.002 0.008 0.002 0.000
copy D->H peak descriptor ( 811): 0.001 0.002 0.010 0.002 0.000
________________________________________________________________________________
Total OpenCL execution time : 7216.835ms
Overhead due to Python: 8.4%
The overhead from calling OpenCL from Python is as low as 10% (lower performances due to memory allocation)
[13]:
# Visualization of the performances:
res8 = pf.peakfinder8(fimg.data, error_model="azimuthal",
cycle=3, cutoff_peak=3, noise=2, connected=9, patch_size=5)
print(len(res8), res8.dtype)
128 [('index', '<i4'), ('intensity', '<f4'), ('sigma', '<f4'), ('pos0', '<f4'), ('pos1', '<f4')]
[14]:
width = fimg.shape[-1]
y = res8["index"] // width
x = res8["index"] % width
fig, ax = subplots(1, 2, figsize=(10,5))
jupyter.display(dimg, ax=ax[0])
jupyter.display(dimg, ax=ax[1])
ax[0].plot(x, y, ".", label="maxi")
ax[1].plot(x, y, ".")
ax[0].plot(res8["pos1"], res8["pos0"], ".g", label="peak")
ax[1].plot(res8["pos1"], res8["pos0"], ".g")
ax[1].set_xlim(1500, 1800)
ax[1].set_ylim(850, 1020)
_=ax[0].legend()
[15]:
fig,ax = subplots()
res=ax.hist(res8["intensity"], 100, )
Comparison with the original “peakfinder8”¶
This algorithm has a python wrapper available from https://github.com/tjlane/peakfinder8
The next cells installs a local version of the Cython-binded peakfinder8 from github.
Nota: This is a quick & dirty solution.
[16]:
targeturl = "https://github.com/kif/peakfinder8"
targetdir = posixpath.split(targeturl)[-1]
if os.path.exists(targetdir):
shutil.rmtree(targetdir, ignore_errors=True)
pwd = os.getcwd()
try:
os.system("git clone " + targeturl)
os.chdir(targetdir)
os.system(f"'{sys.executable}' setup.py build")
finally:
os.chdir(pwd)
sys.path.append(pwd+"/"+glob.glob(f"{targetdir}/build/lib*")[0])
from ssc.peakfinder8_extension import peakfinder_8
Cloning into 'peakfinder8'...
Compiling ext/peakfinder8/peakfinder8_extension.pyx because it changed.
[1/1] Cythonizing ext/peakfinder8/peakfinder8_extension.pyx
/media/edgar1993a/TOSHIBA EXT/venvs/py310_/lib/python3.10/site-packages/Cython/Compiler/Main.py:381: FutureWarning: Cython directive 'language_level' not set, using '3str' for now (Py3). This has changed from earlier releases! File: /home/edgar1993a/work/pyFAI_2024/doc/source/usage/tutorial/Separation/peakfinder8/ext/peakfinder8/peakfinder8_extension.pyx
tree = Parsing.p_module(s, pxd, full_module_name)
In file included from /media/edgar1993a/TOSHIBA EXT/venvs/py310_/lib/python3.10/site-packages/numpy/core/include/numpy/ndarraytypes.h:1929,
from /media/edgar1993a/TOSHIBA EXT/venvs/py310_/lib/python3.10/site-packages/numpy/core/include/numpy/ndarrayobject.h:12,
from /media/edgar1993a/TOSHIBA EXT/venvs/py310_/lib/python3.10/site-packages/numpy/core/include/numpy/arrayobject.h:5,
from ext/peakfinder8/peakfinder8_extension.cpp:1235:
/media/edgar1993a/TOSHIBA EXT/venvs/py310_/lib/python3.10/site-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:17:2: warning: #warning "Using deprecated NumPy API, disable it with " "#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-Wcpp]
17 | #warning "Using deprecated NumPy API, disable it with " \
| ^~~~~~~
ext/peakfinder8/peakfinders.cpp: In function ‘int peakfinder3(tPeakList*, float*, char*, long int, long int, long int, long int, float, float, long int, long int, long int)’:
ext/peakfinder8/peakfinders.cpp:319:15: warning: variable ‘thisr’ set but not used [-Wunused-but-set-variable]
319 | float thisr;
| ^~~~~
ext/peakfinder8/peakfinders.cpp:176:8: warning: variable ‘total’ set but not used [-Wunused-but-set-variable]
176 | float total;
| ^~~~~
ext/peakfinder8/peakfinders.cpp: In function ‘int peakfinder8(tPeakList*, float*, char*, float*, long int, long int, long int, long int, float, float, long int, long int, long int)’:
ext/peakfinder8/peakfinders.cpp:463:8: warning: variable ‘total’ set but not used [-Wunused-but-set-variable]
463 | float total;
| ^~~~~
ext/peakfinder8/peakfinders.cpp:507:7: warning: variable ‘lminr’ set but not used [-Wunused-but-set-variable]
507 | long lminr, lmaxr;
| ^~~~~
ext/peakfinder8/peakfinders.cpp:522:33: warning: argument 1 value ‘18446744072709551617’ exceeds maximum object size 9223372036854775807 [-Walloc-size-larger-than=]
522 | float *rsigma = (float*) calloc(lmaxr, sizeof(float));
| ~~~~~~^~~~~~~~~~~~~~~~~~~~~~
In file included from /usr/include/c++/9/bits/std_abs.h:38,
from /usr/include/c++/9/cmath:47,
from /usr/include/c++/9/math.h:36,
from ext/peakfinder8/peakfinders.cpp:5:
/usr/include/stdlib.h:542:14: note: in a call to allocation function ‘void* calloc(size_t, size_t)’ declared here
542 | extern void *calloc (size_t __nmemb, size_t __size)
| ^~~~~~
ext/peakfinder8/peakfinders.cpp:523:34: warning: argument 1 value ‘18446744072709551617’ exceeds maximum object size 9223372036854775807 [-Walloc-size-larger-than=]
523 | float *roffset = (float*) calloc(lmaxr, sizeof(float));
| ~~~~~~^~~~~~~~~~~~~~~~~~~~~~
In file included from /usr/include/c++/9/bits/std_abs.h:38,
from /usr/include/c++/9/cmath:47,
from /usr/include/c++/9/math.h:36,
from ext/peakfinder8/peakfinders.cpp:5:
/usr/include/stdlib.h:542:14: note: in a call to allocation function ‘void* calloc(size_t, size_t)’ declared here
542 | extern void *calloc (size_t __nmemb, size_t __size)
| ^~~~~~
ext/peakfinder8/peakfinders.cpp:524:31: warning: argument 1 value ‘18446744072709551617’ exceeds maximum object size 9223372036854775807 [-Walloc-size-larger-than=]
524 | long *rcount = (long*) calloc(lmaxr, sizeof(long));
| ~~~~~~^~~~~~~~~~~~~~~~~~~~~
In file included from /usr/include/c++/9/bits/std_abs.h:38,
from /usr/include/c++/9/cmath:47,
from /usr/include/c++/9/math.h:36,
from ext/peakfinder8/peakfinders.cpp:5:
/usr/include/stdlib.h:542:14: note: in a call to allocation function ‘void* calloc(size_t, size_t)’ declared here
542 | extern void *calloc (size_t __nmemb, size_t __size)
| ^~~~~~
ext/peakfinder8/peakfinders.cpp:525:37: warning: argument 1 value ‘18446744072709551617’ exceeds maximum object size 9223372036854775807 [-Walloc-size-larger-than=]
525 | float *rthreshold = (float*) calloc(lmaxr, sizeof(float));
| ~~~~~~^~~~~~~~~~~~~~~~~~~~~~
In file included from /usr/include/c++/9/bits/std_abs.h:38,
from /usr/include/c++/9/cmath:47,
from /usr/include/c++/9/math.h:36,
from ext/peakfinder8/peakfinders.cpp:5:
/usr/include/stdlib.h:542:14: note: in a call to allocation function ‘void* calloc(size_t, size_t)’ declared here
542 | extern void *calloc (size_t __nmemb, size_t __size)
| ^~~~~~
[17]:
%%time
#Create some compatibility layer:
img = fimg.data.astype("float32")
r = ai._cached_array['r_center'].astype("float32")
# r = numpy.ones_like(img)
imask = (1-mask).astype("int8")
max_num_peaks = 1000
asic_nx = img.shape[-1]
asic_ny = img.shape[0]
nasics_x = 1
nasics_y = 1
adc_threshold = 2.0
minimum_snr = 3.0
min_pixel_count = 9
max_pixel_count = 999
local_bg_radius = 3
accumulated_shots = 1
min_res = 0
max_res = 3000
CPU times: user 14.5 ms, sys: 16.5 ms, total: 31.1 ms
Wall time: 28.4 ms
[18]:
%%time
ref = peakfinder_8(max_num_peaks,
img, imask, r,
asic_nx, asic_ny, nasics_x, nasics_y,
adc_threshold, minimum_snr,
min_pixel_count, max_pixel_count, local_bg_radius)
CPU times: user 185 ms, sys: 12.1 ms, total: 197 ms
Wall time: 196 ms
[19]:
timeit_cpu = %timeit -o peakfinder_8(max_num_peaks, img, imask, r, asic_nx, asic_ny, nasics_x, nasics_y, adc_threshold, minimum_snr,min_pixel_count, max_pixel_count, local_bg_radius)
176 ms ± 3.49 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
[20]:
print("Number of peak found: ", len(ref[0]), len(ref[1]), len(ref[2]))
print("Speed up of GPU vs CPU: ", timeit_cpu.best/timeit_gpu.best)
Number of peak found: 995 995 995
Speed up of GPU vs CPU: 17.747103983352165
[21]:
#Display the peaks
fig, ax = subplots(1, 2, figsize=(10,5))
jupyter.display(dimg, ax=ax[0])
jupyter.display(dimg, ax=ax[1])
ax[0].plot(ref[0], ref[1], ".g")
ax[1].plot(ref[0], ref[1], ".g")
ax[1].set_xlim(1500, 1800)
ax[1].set_ylim(850, 1020)
pass
Conclusion¶
The re-implementation of peakfinder8 in pyFAI takes advantage of the many parallel threads available on GPU which makes it 20 times faster than the original implementation in C++. Despite this algorithm has been re-designed for GPU, it can also run on CPU but it would not be optimized there thus it is likely to be slower.
The results obtained with the Python/OpenCL implementation looks better, this is probably due to a slightly different threshold \(I > mean + max(N*sigma, noise)\) instead of \(I > max(noise, mean + N*sigma)\).
[22]:
print(f"Total execution time: {time.perf_counter()-start_time:.3f}s")
Total execution time: 47.771s