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= 172µm, 172µ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)
10.5 ms ± 502 μ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.896 3.747 17.709 4.068 1.260
memset_ng ( 811): 0.003 0.015 0.072 0.015 0.008
corrections4a ( 811): 0.600 0.605 0.716 0.610 0.015
csr_sigma_clip4 ( 811): 4.009 4.057 5.168 4.149 0.191
memset counter ( 811): 0.002 0.003 0.009 0.003 0.001
find_intense ( 811): 0.835 0.844 0.947 0.848 0.015
copy D->H counter ( 811): 0.001 0.001 0.002 0.001 0.000
________________________________________________________________________________
Total OpenCL execution time : 7862.400ms
Overhead due to Python: 7.8%
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)
12.6 ms ± 635 μ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.859 3.530 8.728 3.834 0.964
memset_ng ( 811): 0.003 0.014 0.248 0.012 0.012
corrections4a ( 811): 0.599 0.604 0.715 0.608 0.014
csr_sigma_clip4 ( 811): 4.028 4.090 8.674 4.224 0.304
memset counter ( 811): 0.003 0.003 0.012 0.003 0.001
peakfinder ( 811): 3.033 3.197 5.532 3.249 0.165
copy D->H counter ( 811): 0.001 0.002 0.004 0.002 0.000
copy D->H peak positions ( 811): 0.002 0.002 0.005 0.002 0.000
copy D->H peak descriptor ( 811): 0.001 0.001 0.006 0.002 0.000
________________________________________________________________________________
Total OpenCL execution time : 9678.920ms
Overhead due to Python: 5.5%
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'...
/home/edgar1993a/miniforge3/envs/ewoks/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/doc/source/usage/tutorial/Separation/peakfinder8/ext/peakfinder8/peakfinder8_extension.pyx
tree = Parsing.p_module(s, pxd, full_module_name)
Compiling ext/peakfinder8/peakfinder8_extension.pyx because it changed.
[1/1] Cythonizing ext/peakfinder8/peakfinder8_extension.pyx
In file included from /home/edgar1993a/miniforge3/envs/ewoks/lib/python3.10/site-packages/numpy/_core/include/numpy/ndarraytypes.h:1909,
from /home/edgar1993a/miniforge3/envs/ewoks/lib/python3.10/site-packages/numpy/_core/include/numpy/ndarrayobject.h:12,
from /home/edgar1993a/miniforge3/envs/ewoks/lib/python3.10/site-packages/numpy/_core/include/numpy/arrayobject.h:5,
from ext/peakfinder8/peakfinder8_extension.cpp:1270:
/home/edgar1993a/miniforge3/envs/ewoks/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 16.9 ms, sys: 12 ms, total: 28.9 ms
Wall time: 26.8 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 178 ms, sys: 15.9 ms, total: 194 ms
Wall time: 193 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)
184 ms ± 1.37 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: 15.258822549471374
[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: 52.966s