diff --git a/bootstrap.py b/bootstrap.py index 57adc14..0059409 100755 --- a/bootstrap.py +++ b/bootstrap.py @@ -1,4 +1,4 @@ -#!/usr/bin/env python +#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Bootstrap helps you to test scripts without installing them diff --git a/dynamix/correlator/common.py b/dynamix/correlator/common.py index e60ecf0..82084d4 100644 --- a/dynamix/correlator/common.py +++ b/dynamix/correlator/common.py @@ -1,11 +1,11 @@ import numpy as np from os import linesep - -import pyopencl.array as parray -from pyopencl.tools import dtype_to_ctype +from .. import resources +resources.silx_integration() from silx.opencl.common import pyopencl as cl from silx.opencl.processing import OpenclProcessing, KernelContainer - +import pyopencl.array as parray +from pyopencl.tools import dtype_to_ctype class BaseCorrelator(object): "Abstract base class for all Correlators" @@ -193,23 +193,3 @@ def _reset_arrays(self, arrays_names): if old_array is not None: setattr(self, array_name, old_array) setattr(self, old_array_name, None) - - # Overwrite OpenclProcessing.compile_kernel, as it does not support - # kernels outside silx/opencl/resources - def compile_kernels(self, kernel_files=None, compile_options=None): - kernel_files = kernel_files or self.kernel_files - - allkernels_src = [] - for kernel_file in kernel_files: - with open(kernel_file) as fid: - kernel_src = fid.read() - allkernels_src.append(kernel_src) - allkernels_src = linesep.join(allkernels_src) - - compile_options = compile_options or self.get_compiler_options() - try: - self.program = cl.Program(self.ctx, allkernels_src).build(options=compile_options) - except (cl.MemoryError, cl.LogicError) as error: - raise MemoryError(error) - else: - self.kernels = KernelContainer(self.program) diff --git a/dynamix/correlator/dense.py b/dynamix/correlator/dense.py index 2bb2d53..2a7b7c9 100644 --- a/dynamix/correlator/dense.py +++ b/dynamix/correlator/dense.py @@ -4,7 +4,7 @@ import pyopencl.array as parray from os import path from multiprocessing import cpu_count -from ..utils import nextpow2, updiv, get_opencl_srcfile, get_next_power +from ..utils import nextpow2, get_opencl_srcfile, get_next_power, CorrelationResult from .common import OpenclCorrelator, BaseCorrelator from silx.math.fft.fftw import FFTW @@ -26,7 +26,6 @@ NCPU = cpu_count() -CorrelationResult = namedtuple("CorrelationResult", "res dev") def py_dense_correlator(xpcs_data, mask, calc_std=False): @@ -144,7 +143,7 @@ def correlate(self, frames, calc_std=False): class DenseCorrelator(OpenclCorrelator): - kernel_files = ["densecorrelator.cl"] + kernel_files = ["dynamix:opencl/densecorrelator.cl"] def __init__( self, shape, nframes, @@ -167,9 +166,8 @@ def __init__( self._allocate_arrays() def _setup_kernels(self): - kernel_files = list(map(get_opencl_srcfile, self.kernel_files)) self.compile_kernels( - kernel_files=kernel_files, + kernel_files=self.kernel_files, compile_options=[ "-DIMAGE_WIDTH=%d" % self.shape[1], "-DNUM_BINS=%d" % self.n_bins, diff --git a/dynamix/correlator/event.py b/dynamix/correlator/event.py index f254884..1f73998 100644 --- a/dynamix/correlator/event.py +++ b/dynamix/correlator/event.py @@ -1,12 +1,25 @@ +__authors__ = {"Pierre Paleo", "Jerome Kieffer"} +__contact__ = "Jerome.Kieffer@ESRF.eu" +__license__ = "MIT" +__copyright__ = "European Synchrotron Radiation Facility, Grenoble, France" +__date__ = "05/07/2021" +__status__ = "stable" + import numpy as np -# from silx.opencl.common import pyopencl as cl -import pyopencl.array as parray -from ..utils import get_opencl_srcfile -from .common import OpenclCorrelator +import logging +from silx.opencl.processing import BufferDescription, EventDescription +from silx.opencl.common import pyopencl +from pyopencl import array as parray +from pyopencl.tools import dtype_to_ctype +from ..utils import get_opencl_srcfile, Compacted +from ..tools.decorators import timeit +from .common import OpenclCorrelator, OpenclProcessing +logger = logging.getLogger(__name__) + class EventCorrelator(OpenclCorrelator): - kernel_files = ["evtcorrelator.cl"] + kernel_files = ["dynamix:opencl/evtcorrelator.cl"] """ A class to compute the correlation function for sparse XPCS data. @@ -57,20 +70,17 @@ def __init__( self._allocate_events_arrays() self._setup_kernels() - def _set_events_size(self, max_events_count, total_events_count): self.max_events_count = max_events_count self.total_events_count = total_events_count or max_events_count * np.prod(self.shape) - def _setup_kernels(self): - kernel_files = list(map(get_opencl_srcfile, self.kernel_files)) self.compile_kernels( - kernel_files=kernel_files, + kernel_files=self.kernel_files, compile_options=[ "-DIMAGE_WIDTH=%d" % self.shape[1], "-DDTYPE=%s" % self.c_dtype, - "-DSUM_WG_SIZE=%d" % 1024, # TODO tune ? + "-DSUM_WG_SIZE=%d" % 1024, # TODO tune ? "-DMAX_EVT_COUNT=%d" % self.max_events_count, "-DSCALE_FACTOR=%f" % self.scale_factors[1], "-DNUM_BINS=%d" % self.n_bins, @@ -80,15 +90,14 @@ def _setup_kernels(self): self.normalization_kernel = self.kernels.get_kernel("normalize_correlation") self.grid = self.shape[::-1] - self.wg = None # tune ? - + self.wg = None # tune ? def _allocate_events_arrays(self, is_reallocating=False): tot_nnz = self.total_events_count self.d_vol_times = parray.zeros(self.queue, tot_nnz, dtype=np.int32) self.d_vol_data = parray.zeros(self.queue, tot_nnz, dtype=self.dtype) - self.d_offsets = parray.zeros(self.queue, np.prod(self.shape)+1, dtype=np.uint32) + self.d_offsets = parray.zeros(self.queue, np.prod(self.shape) + 1, dtype=np.uint32) self._old_d_vol_times = None self._old_d_vol_data = None @@ -100,7 +109,6 @@ def _allocate_events_arrays(self, is_reallocating=False): self.d_res = parray.zeros(self.queue, self.output_shape, np.float32) self.d_scale_factors = parray.to_device(self.queue, np.array(list(self.scale_factors.values()), dtype=np.float32)) - def _check_event_arrays(self, vol_times, vol_data, offsets): for arr in [vol_times, vol_data, offsets]: assert arr.ndim == 1 @@ -114,7 +122,6 @@ def _check_event_arrays(self, vol_times, vol_data, offsets): else: raise ValueError("Too many events and allow_reallocate was set to False") - def correlate(self, vol_times, vol_data, offsets): self._check_event_arrays(vol_times, vol_data, offsets) self._set_data({ @@ -145,7 +152,7 @@ def correlate(self, vol_times, vol_data, offsets): evt = self.normalization_kernel( self.queue, (self.nframes, self.n_bins), - None, # tune wg ? + None, # tune wg ? self.d_res_int.data, self.d_res.data, self.d_sums.data, @@ -160,7 +167,6 @@ def correlate(self, vol_times, vol_data, offsets): return self.d_res.get() - class FramesCompressor(object): """ A class for compressing frames on-the-fly. @@ -213,14 +219,12 @@ def __init__(self, shape, nframes, max_nnz, dtype=np.int8): self.npix = np.prod(self.shape) self._init_events_datastructure() - def _init_events_datastructure(self): self.events_counter = np.zeros(self.npix, dtype=np.int32) self.events = np.zeros((self.npix, self.max_nnz), dtype=self.dtype) self.times = np.zeros((self.npix, self.max_nnz), dtype=np.int32) self.frames_counter = 0 - @staticmethod def compress_all_stack(frames): """ @@ -235,11 +239,10 @@ def compress_all_stack(frames): res_times = times[nnz_indices[-1]] offsets = np.cumsum((frames > 0).sum(axis=0).ravel()) - res_offsets = np.zeros(np.prod(frames.shape[1:])+1, dtype=np.uint32) + res_offsets = np.zeros(np.prod(frames.shape[1:]) + 1, dtype=np.uint32) res_offsets[1:] = offsets[:] - return res_data, res_times, res_offsets - + return Compacted(res_data, res_times, res_offsets) def process_frame(self, frame): """ @@ -253,7 +256,6 @@ def process_frame(self, frame): self.events_counter[mask] += 1 self.frames_counter += 1 - def get_compacted_events(self, wait_for_all_frames=True): """ Compact all compressed frames into three 1D structures: @@ -275,7 +277,135 @@ def get_compacted_events(self, wait_for_all_frames=True): events = self.events.ravel()[m] times = self.times.ravel()[m] - return events, times, offsets + return Compacted(events, times, offsets) +class OpenclCompressor(OpenclProcessing): + """See doc of FramesCompressor""" + + kernel_files = ["dynamix:opencl/compactor.cl"] + + def __init__(self, shape, nframes, max_nnz, dtype=np.int8, slab_size=100, + ctx=None, devicetype="all", platformid=None, deviceid=None, + block_size=None, memory=None, profile=False): + OpenclProcessing.__init__(self, ctx=ctx, devicetype=devicetype, platformid=platformid, deviceid=deviceid, + block_size=block_size, memory=memory, profile=profile) + self.shape = shape + self.nframes = nframes + self.dtype = dtype + self.max_nnz = max_nnz + self.npix = np.prod(self.shape) + self.slab_size = slab_size + self.frames_counter = 0 + self.grid = None + self._allocate_buffers() + self._setup_kernel() + self.reset() + + def _allocate_buffers(self): + buffers = [BufferDescription("counter", (self.npix,), np.uint32, None), + BufferDescription("slab", (self.slab_size, self.npix), self.dtype, None), + BufferDescription("events", (self.npix, self.max_nnz), self.dtype, None), + BufferDescription("times", (self.npix, self.max_nnz), np.int32, None), + ] + + self.allocate_buffers(buffers, use_array=True) + + def _setup_kernel(self): + self.compile_kernels( + kernel_files=self.kernel_files, + compile_options=[ + f"-DDTYPE={dtype_to_ctype(self.dtype)}", + ] + ) + self.grid = self.shape[::-1] + def reset(self): + """Reset the compressor""" + self.frames_counter = 0 + self.cl_mem["slab"].fill(0) + self.cl_mem["times"].fill(0), + self.cl_mem["events"].fill(0), + self.cl_mem["counter"].fill(0) + + # @timeit + def process_stack(self, stack): + self.reset() + for start in range(0, stack.shape[0], self.slab_size): + end = min(stack.shape[0], start+self.slab_size) + self.process_slab(stack[start:end], start, end) + + def process_slab(self, frames, start, end): + """ + Compress a single frame, and update the events stack state. + """ + nframes = end-start + if frames.ndim == 3: + assert nframes == frames.shape[0] + frames = frames.reshape((nframes, -1)) + elif nframes == 1: + frames = frames.reshape((1, -1)) + assert start >= self.frames_counter, "Process slabs in incremental order !" + if nframesself.slab_size: + logger.warning("Too many frames in slab: %s slab_size: %s, dropping some of them", + frames.shape[0], self.slab_size) + self.cl_mem["slab"].set(frames[:self.slab_size]) + else: + self.cl_mem["slab"].set(frames) + + evt = self.program.compactor1(self.queue, self.grid, None, + self.cl_mem["slab"].data, + np.int32(start), + np.int32(end), + np.int32(self.shape[1]), + np.int32(self.shape[0]), + np.int32(self.max_nnz), + self.cl_mem["events"].data, + self.cl_mem["times"].data, + self.cl_mem["counter"].data, + ) + + self.frames_counter = end + if self.profile: + self.events.append(EventDescription("compactor1", evt)) + + def process_frame(self, frame): + """ + Compress a single frame, and update the events stack state. + """ + self.process_slab(frame, self.frames_counter, self.frames_counter+1) + + # @timeit + def get_compacted_events(self, wait_for_all_frames=True): + """ + Compact all compressed frames into three 1D structures: + - values: 1D array containing all the nonzero data points + - times: time indices corresponding to nonzero data points + - offsets: values[offsets[k]:offsets[k+1]] corresponds to all values>0 for pixel index k + """ + if self.frames_counter < self.nframes - 1 and wait_for_all_frames: + raise RuntimeError( + "Not all frames were compressed yet (%d/%d)" % + (self.frames_counter, self.nframes) + ) + counter = self.cl_mem["counter"].get() + + #TODO: implement in OpenCL as well ! + offsets = np.zeros(self.npix + 1, dtype=np.uint32) + offsets[1:] = np.cumsum(counter) + + values = self.cl_mem["events"].get().ravel() + times = self.cl_mem["times"].get().ravel() + msk = values > 0 + values = values[msk] + times = times[msk] + return Compacted(values, times, offsets) + + def compress_all_stack(self, frames): + """ + Build the whole event structure assuming that all the frames are given. + """ + self.process_stack(frames) + return self.get_compacted_events() \ No newline at end of file diff --git a/dynamix/correlator/event_y.py b/dynamix/correlator/event_y.py index c802398..9756556 100644 --- a/dynamix/correlator/event_y.py +++ b/dynamix/correlator/event_y.py @@ -1,107 +1,103 @@ #! /usr/bin/env python3 -#wxpcs code that works with ini file -import sys -#sys.path.append("/data/id10/inhouse/Programs/PyXPCS_project/wxpcs") -#sys.path.append("/users/chushkin/Documents/Analysis/Glass_school_2019/wxpcs") -#sys.path.append("/users/chushkin/Documents/Programs/PyXPCS_project/wxpcs") +# wxpcs code that works with ini file import numpy as np import numba as nb from .WXPCS import fecorrt import time from collections import namedtuple - -import os import psutil -nproc = psutil.cpu_count(logical=True)#len(os.sched_getaffinity(0)) * psutil.cpu_count(logical=False) // psutil.cpu_count(logical=True) +nproc = psutil.cpu_count(logical=True) # len(os.sched_getaffinity(0)) * psutil.cpu_count(logical=False) // psutil.cpu_count(logical=True) CorrelationResult = namedtuple("CorrelationResult", "res dev trc") + ########## Numba implementation ################################## @nb.jit(nopython=True, parallel=True, fastmath=True) -def ncorrelate(evs,tms,cnt,q,n_frames): +def ncorrelate(evs, tms, cnt, q, n_frames): qm = q.max() - cc = np.zeros((qm,n_frames,n_frames),np.float32) - mint = np.zeros((qm,n_frames),np.float32) + cc = np.zeros((qm, n_frames, n_frames), np.float32) + mint = np.zeros((qm, n_frames), np.float32) ll = cnt.size k = 0 for i in range(ll): - qp = q[i]-1 - if qp>=0: - for j in range(cnt[i]): + qp = q[i] - 1 + if qp >= 0: + for j in range(cnt[i]): t0 = tms[k] - n = k+1 - mint[qp,t0] += evs[k] - for f in range(j+1,cnt[i],1): - cc[qp,t0,tms[n]] += evs[k]*evs[n] - n += 1 - k += 1 + n = k + 1 + mint[qp, t0] += evs[k] + for f in range(j + 1, cnt[i], 1): + cc[qp, t0, tms[n]] += evs[k] * evs[n] + n += 1 + k += 1 else: k += cnt[i] - return cc,mint + return cc, mint + @nb.jit(nopython=True, parallel=True, fastmath=True) -def ncorrelatep(evs,tms,cnt,q,n_frames,nproc): +def ncorrelatep(evs, tms, cnt, q, n_frames, nproc): qm = q.max() - cc = np.zeros((nproc,qm,n_frames,n_frames),np.float32) - mint = np.zeros((nproc,qm,n_frames),np.float32) - vk = np.zeros((cnt.size+1),np.uint32) + cc = np.zeros((nproc, qm, n_frames, n_frames), np.float32) + mint = np.zeros((nproc, qm, n_frames), np.float32) + vk = np.zeros((cnt.size + 1), np.uint32) for i in range(cnt.size): - vk[i+1] = vk[i]+cnt[i] + vk[i + 1] = vk[i] + cnt[i] for proc in nb.prange(nproc): for i in range(proc, cnt.size, nproc): k = vk[i] - qp = q[i]-1 - if qp>=0: + qp = q[i] - 1 + if qp >= 0: for j in range(cnt[i]): t0 = tms[k] - n = k+1 - mint[proc, qp,t0] += evs[k] - for f in range(j+1,cnt[i],1): - cc[proc, qp,t0,tms[n]] += evs[k]*evs[n] - n += 1 + n = k + 1 + mint[proc, qp, t0] += evs[k] + for f in range(j + 1, cnt[i], 1): + cc[proc, qp, t0, tms[n]] += evs[k] * evs[n] + n += 1 k += 1 else: k += cnt[i] - for i in range(1,nproc,1): + for i in range(1, nproc, 1): cc[0,:,:,:] += cc[i,:,:,:] mint[0,:,:] += mint[i,:,:] - return cc[0,:,:,:],mint[0,:,:] + return cc[0,:,:,:], mint[0,:,:] + #### The fastes numba implementation and does not use extra memory ### @nb.jit(nopython=True, parallel=True, fastmath=True) -def ncorrelatepm(evs,tms,cnt,q,n_frames,nproc): +def ncorrelatepm(evs, tms, cnt, q, n_frames, nproc): qm = q.max() - cc = np.zeros((qm,n_frames,n_frames),np.float32) - mint = np.zeros((qm,n_frames),np.float32) - vk = np.zeros((cnt.size+1),np.uint32) + cc = np.zeros((qm, n_frames, n_frames), np.float32) + mint = np.zeros((qm, n_frames), np.float32) + vk = np.zeros((cnt.size + 1), np.uint32) for i in range(cnt.size): - vk[i+1] = vk[i]+cnt[i] - tstep = int(np.ceil(tms.max()/nproc)) + vk[i + 1] = vk[i] + cnt[i] + tstep = int(np.ceil(tms.max() / nproc)) print(tstep) for proc in nb.prange(nproc): for i in range(cnt.size): k = vk[i] - qp = q[i]-1 - if qp>=0: + qp = q[i] - 1 + if qp >= 0: for j in range(cnt[i]): t0 = tms[k] - if t0 >= tstep * proc and t0 < tstep * (proc+1): - #if t0 % nproc == proc: - n = k+1 - mint[qp,t0] += evs[k] - for f in range(j+1,cnt[i],1): - cc[qp,t0,tms[n]] += evs[k]*evs[n] - n += 1 + if t0 >= tstep * proc and t0 < tstep * (proc + 1): + # if t0 % nproc == proc: + n = k + 1 + mint[qp, t0] += evs[k] + for f in range(j + 1, cnt[i], 1): + cc[qp, t0, tms[n]] += evs[k] * evs[n] + n += 1 k += 1 else: k += cnt[i] - return cc,mint - + return cc, mint ##### Event_correlator standard several qs using Numba ######### -def nbecorrts_q(events,times,cnt,qqmask, n_frames, calc_std=False, ttcf_par=0): +def nbecorrts_q(events, times, cnt, qqmask, n_frames, calc_std=False, ttcf_par=0): """ Calculation of the event correlation function using Numba :param events: 1D array of events values @@ -114,221 +110,213 @@ def nbecorrts_q(events,times,cnt,qqmask, n_frames, calc_std=False, ttcf_par=0): :return: CorrelationResult structure with correlation functions """ print("start calculation of the correlation functions") - t0 = time.time() + t0 = time.perf_counter() + + qqmask = np.ravel(qqmask)[cnt > 0] + cnt = cnt[cnt > 0] - qqmask = np.ravel(qqmask)[cnt>0] - cnt = cnt[cnt>0] - - #num, mint = ncorrelate(events,times,cnt,qqmask,n_frames) print("Number of used processors: ", nproc) - num, mint = ncorrelatepm(events,times,cnt,qqmask,n_frames,nproc)# parallel - #num, mint = ncorrelate(events,times,cnt,qqmask,n_frames)# paralle - res = [] + num, mint = ncorrelatepm(events, times, cnt, qqmask, n_frames, nproc) # parallel qm = qqmask.max() - res = np.zeros((qm,n_frames-1), dtype=np.float32) + res = np.zeros((qm, n_frames - 1), dtype=np.float32) if calc_std: dev = np.zeros_like(res) - for q in range(1,qm+1,1): - for_norm = qqmask[qqmask==q].size - cor = num[q-1,:,:]#/q[q==2].size - s = np.reshape(mint[q-1,:],(n_frames,1))#/for_norm - print("Number of pixels for q %d is %d" % (q,for_norm)) - print("Average photons per frame for q %d is %f" % (q,np.mean(s))) - print("Average photons per frame per pixel for q %d is %f" % (q,np.mean(s)/for_norm)) - - norm = np.dot(s,s.T)/n_frames - cor = cor*for_norm/n_frames - for i in range(1,n_frames): - dia = np.diag(cor,k=i) - sdia = np.diag(norm,k=i) + for q in range(1, qm + 1, 1): + for_norm = qqmask[qqmask == q].size + cor = num[q - 1,:,:] # /q[q==2].size + s = np.reshape(mint[q - 1,:], (n_frames, 1)) # /for_norm + print("Number of pixels for q %d is %d" % (q, for_norm)) + print("Average photons per frame for q %d is %f" % (q, np.mean(s))) + print("Average photons per frame per pixel for q %d is %f" % (q, np.mean(s) / for_norm)) + + norm = np.dot(s, s.T) / n_frames + cor = cor * for_norm / n_frames + for i in range(1, n_frames): + dia = np.diag(cor, k=i) + sdia = np.diag(norm, k=i) ind = np.where(np.isfinite(dia)) - res[q-1,i-1] = np.mean(dia[ind])/np.mean(sdia[ind]) + res[q - 1, i - 1] = np.mean(dia[ind]) / np.mean(sdia[ind]) if calc_std: - dev[q-1,i-1] = np.std(dia[ind]/sdia[ind])/len(sdia[ind])**0.5 + dev[q - 1, i - 1] = np.std(dia[ind] / sdia[ind]) / len(sdia[ind]) ** 0.5 if ttcf_par == q: - trc = cor/norm - tmp = np.diag(trc,k=-5) - if np.sum(tmp)<1: + trc = cor / norm + tmp = np.diag(trc, k=-5) + if np.sum(tmp) < 1: trc += np.rot90(np.fliplr(trc)) - tmp = np.diag(trc,k=1) - tmp = np.mean(tmp[tmp>0]) - for j in range(n_frames-1): - trc[j,j]=tmp + tmp = np.diag(trc, k=1) + tmp = np.mean(tmp[tmp > 0]) + for j in range(n_frames - 1): + trc[j, j] = tmp del tmp - print("Total correlation time %2.2f sec" % (time.time()-t0)) - #return wcf,mtcf,cor + print("Total correlation time %2.2f sec" % (time.perf_counter() - t0)) + # return wcf,mtcf,cor if ttcf_par == 0: - trc = 0 + trc = 0 if not(calc_std): - dev = 0 + dev = 0 return CorrelationResult(res, dev, trc) ########### Event_correlator standard several qs ################ -def ecorrts_q(pixels,s,qqmask, calc_std=False, ttcf_par=0): +def ecorrts_q(pixels, s, qqmask, calc_std=False, ttcf_par=0): print("start calculation of the correlation functions") - t0 = time.time() - lpixels = len(pixels)#number of frames + t0 = time.perf_counter() + lpixels = len(pixels) # number of frames rpixels = range(lpixels) - t=[] + t = [] for t1 in rpixels: - t += [t1]*s[t1] - pix = np.concatenate(pixels).ravel()-1#to fit with in1d function result - pix = np.array(pix,dtype=np.uint32) - t = np.array(t,dtype=np.uint32) - indpi = np.lexsort((t,pix)) + t += [t1] * s[t1] + pix = np.concatenate(pixels).ravel() - 1 # to fit with in1d function result + pix = np.array(pix, dtype=np.uint32) + t = np.array(t, dtype=np.uint32) + indpi = np.lexsort((t, pix)) t = t[indpi] - pix = pix[indpi] - res = [] + pix = pix[indpi] qqmask = np.ravel(qqmask) qm = qqmask.max() - res = np.zeros((qm,lpixels-1), dtype=np.float32) + res = np.zeros((qm, lpixels - 1), dtype=np.float32) if calc_std: dev = np.zeros_like(res) - for q in range(1,qm+1,1): - for_norm = len(qqmask[qqmask==q]) - indx = np.where(qqmask==q)[0]#here I can make +1 or pix-1 - indpi = np.in1d(pix,indx) + for q in range(1, qm + 1, 1): + for_norm = len(qqmask[qqmask == q]) + indx = np.where(qqmask == q)[0] # here I can make +1 or pix-1 + indpi = np.in1d(pix, indx) tq = t[indpi] - pixq = pix[indpi] - s0 = np.array(s)+0 - s, tt = np.histogram(tq, bins=lpixels) - print("Number of pixels for q %d is %d" % (q,for_norm)) - print("Average photons per frame for q %d is %f" % (q,np.mean(s))) - print("Average photons per frame per pixel for q %d is %f" % (q,np.mean(s)/for_norm)) + pixq = pix[indpi] + s, _ = np.histogram(tq, bins=lpixels) + print("Number of pixels for q %d is %d" % (q, for_norm)) + print("Average photons per frame for q %d is %f" % (q, np.mean(s))) + print("Average photons per frame per pixel for q %d is %f" % (q, np.mean(s) / for_norm)) lenpi = len(pixq) - cor = np.zeros((lpixels,lpixels),dtype=np.uint32) - cor = fecorrt(pixq,tq,cor,lenpi,lpixels)#fortran module + cor = np.zeros((lpixels, lpixels), dtype=np.uint32) + cor = fecorrt(pixq, tq, cor, lenpi, lpixels) # fortran module lens = len(s) - s = np.array(s,dtype=np.float32) - cor = np.array(cor,dtype=np.float32) - s.shape = lens,1 - #norm = dot(s,flipud(s.T))/lpixels - norm = np.dot(s,s.T)/lpixels - cor = cor*for_norm/lpixels - for i in range(1,lpixels): - dia = np.diag(cor,k=i) - sdia = np.diag(norm,k=i) + s = np.array(s, dtype=np.float32) + cor = np.array(cor, dtype=np.float32) + s.shape = lens, 1 + # norm = dot(s,flipud(s.T))/lpixels + norm = np.dot(s, s.T) / lpixels + cor = cor * for_norm / lpixels + for i in range(1, lpixels): + dia = np.diag(cor, k=i) + sdia = np.diag(norm, k=i) ind = np.where(np.isfinite(dia)) - res[q-1,i-1] = np.mean(dia[ind])/np.mean(sdia[ind]) + res[q - 1, i - 1] = np.mean(dia[ind]) / np.mean(sdia[ind]) if calc_std: - dev[q-1,i-1] = np.std(dia[ind]/sdia[ind])/len(sdia[ind])**0.5 - #wcf = x+0 - #wcf[np.isnan(wcf)] = 0.01 - #mtcf = cftomt_testing(x) - #mtcf = tools.cftomt(x) - #res.append(wcf) + dev[q - 1, i - 1] = np.std(dia[ind] / sdia[ind]) / len(sdia[ind]) ** 0.5 + # wcf = x+0 + # wcf[np.isnan(wcf)] = 0.01 + # mtcf = cftomt_testing(x) + # mtcf = tools.cftomt(x) + # res.append(wcf) if ttcf_par == q: - trc = cor/norm - tmp = np.diag(trc,k=1) - tmp = np.mean(tmp[tmp>0]) - for j in range(lpixels-1): - trc[j,j]=tmp + trc = cor / norm + tmp = np.diag(trc, k=1) + tmp = np.mean(tmp[tmp > 0]) + for j in range(lpixels - 1): + trc[j, j] = tmp del tmp - print("Total correlation time %2.2f sec" % (time.time()-t0)) - #return wcf,mtcf,cor + print("Total correlation time %2.2f sec" % (time.perf_counter() - t0)) + # return wcf,mtcf,cor if ttcf_par == 0: - trc = 0 + trc = 0 if not(calc_std): - dev = 0 + dev = 0 return CorrelationResult(res, dev, trc) ########### Event_correlator standard ################ -def ecorrts(pixels,s,for_norm): +def ecorrts(pixels, s, for_norm): print("start calculation of the correlation functions") - t0 = time.time() + t0 = time.perf_counter() lpixels = len(pixels) rpixels = range(lpixels) - t=[] + t = [] for t1 in rpixels: - t += [t1]*s[t1] - #print("time for pre loop "+str(time()-timee)) + t += [t1] * s[t1] + # print("time for pre loop "+str(time()-timee)) pix = np.concatenate(pixels).ravel() - #print('start sorting') - times = time() - pix = np.array(pix,dtype=np.int32) + # print('start sorting') + pix = np.array(pix, dtype=np.int32) t = np.array(t) - indpi = np.lexsort((t,pix)) + indpi = np.lexsort((t, pix)) t = t[indpi] - pix = pix[indpi] - #print('sorting took '+ str(time()-times)) - #print('start main loop') - timem = time.time() + pix = pix[indpi] + # print('sorting took '+ str(time()-times)) + # print('start main loop') lenpi = len(pix) - cor = np.zeros((lpixels,lpixels),dtype=np.uint32) - timef = time.time() - cor = fecorrt(pix,t,cor,lenpi,lpixels)#fortran module + cor = np.zeros((lpixels, lpixels), dtype=np.uint32) + cor = fecorrt(pix, t, cor, lenpi, lpixels) # fortran module # to split for multicores - #ind = np.where(pix<133128)[0][-1] - #cor = fecorrt(pix[:ind],t[:ind],cor,len(pix[:ind]),lpixels)#fortran module - #cor1 = np.zeros((lpixels,lpixels),dtype=np.uint32) - #cor1 = fecorrt(pix[ind:],t[ind:],cor1,len(pix[ind:]),lpixels)#fortran module - #cor = cor+cor1 - #print("time for correlating "+str(time()-timem)) - print("average photons per frame "+str(np.mean(s))) - #print(cor.min(),cor.max()) + # ind = np.where(pix<133128)[0][-1] + # cor = fecorrt(pix[:ind],t[:ind],cor,len(pix[:ind]),lpixels)#fortran module + # cor1 = np.zeros((lpixels,lpixels),dtype=np.uint32) + # cor1 = fecorrt(pix[ind:],t[ind:],cor1,len(pix[ind:]),lpixels)#fortran module + # cor = cor+cor1 + # print("time for correlating "+str(time()-timem)) + print("average photons per frame " + str(np.mean(s))) + # print(cor.min(),cor.max()) lens = len(s) - s = np.array(s,dtype=np.float32) - cor = np.array(cor,dtype=np.float32) - s.shape = lens,1 - #norm = dot(s,flipud(s.T))/lpixels - norm = np.dot(s,s.T)/lpixels - cor = cor*for_norm/lpixels - x = np.ones((lpixels-1,3)) - x[:,0] = np.arange(1,lpixels) - for i in range(1,lpixels): - dia = np.diag(cor,k=i) - sdia = np.diag(norm,k=i) - ind = np.where(np.isfinite(dia)) - x[i-1,1] = np.mean(dia[ind])/np.mean(sdia[ind]) - x[i-1,2] = np.std(dia[ind]/sdia[ind])/len(sdia[ind])**0.5 - wcf = x+0 + s = np.array(s, dtype=np.float32) + cor = np.array(cor, dtype=np.float32) + s.shape = lens, 1 + # norm = dot(s,flipud(s.T))/lpixels + norm = np.dot(s, s.T) / lpixels + cor = cor * for_norm / lpixels + x = np.ones((lpixels - 1, 3)) + x[:, 0] = np.arange(1, lpixels) + for i in range(1, lpixels): + dia = np.diag(cor, k=i) + sdia = np.diag(norm, k=i) + ind = np.where(np.isfinite(dia)) + x[i - 1, 1] = np.mean(dia[ind]) / np.mean(sdia[ind]) + x[i - 1, 2] = np.std(dia[ind] / sdia[ind]) / len(sdia[ind]) ** 0.5 + wcf = x + 0 mtcf = cftomt_testing(x) cor /= norm - print("Total time for correlating %2.2f sec" % (time.time()-t0)) - return wcf,mtcf,cor + print("Total time for correlating %2.2f sec" % (time.perf_counter() - t0)) + return wcf, mtcf, cor ############################################################################## + def cftomt_testing(d): - par = 16 - tmp = d[par:,:3] - nt = [] - nd = [] - nse = [] - for i in range(par): - nt.append(d[i,0]) - nd.append(d[i,1]) - nse.append(d[i,2]) - while len(tmp[:,0])>=par: - ntmp = (tmp[:-1,:]+tmp[1:,:])/2 - for i in range(0,par,2): - nt.append(ntmp[i,0]) - nd.append(ntmp[i,1]) - nse.append(ntmp[i,2]) - tmp = ntmp[par:-1:2,:3] - x = np.array([nt,nd,nse]).T - return x + par = 16 + tmp = d[par:,:3] + nt = [] + nd = [] + nse = [] + for i in range(par): + nt.append(d[i, 0]) + nd.append(d[i, 1]) + nse.append(d[i, 2]) + while len(tmp[:, 0]) >= par: + ntmp = (tmp[:-1,:] + tmp[1:,:]) / 2 + for i in range(0, par, 2): + nt.append(ntmp[i, 0]) + nd.append(ntmp[i, 1]) + nse.append(ntmp[i, 2]) + tmp = ntmp[par:-1:2,:3] + x = np.array([nt, nd, nse]).T + return x def gpu_ecorr(events, times, offsets, shape, nframes, dtype, qqmask, F): print("start calculation of the correlation functions") from dynamix.correlator.event import EventCorrelator - t0 = time.time() - E = EventCorrelator(shape, nframes, F.max_nnz, dtype=dtype, total_events_count=events.size, scale_factor=qqmask[qqmask==1].size) + t0 = time.perf_counter() + E = EventCorrelator(shape, nframes, F.max_nnz, dtype=dtype, total_events_count=events.size, scale_factor=qqmask[qqmask == 1].size) result = E.correlate(times, events, offsets) - print("Total correlation time %2.2f sec" % (time.time()-t0)) - res = np.zeros((1,nframes-1), dtype=np.float32) - dev = np.zeros((1,nframes-1), dtype=np.float32) + print("Total correlation time %2.2f sec" % (time.perf_counter() - t0)) + res = np.zeros((1, nframes - 1), dtype=np.float32) + dev = np.zeros((1, nframes - 1), dtype=np.float32) res[0,:] = result[0][1:] - dev[0,:] = result[0][1:]*1e-3 - trc = 0 + dev[0,:] = result[0][1:] * 1e-3 + trc = 0 return CorrelationResult(res, dev, trc) - def gpu_ecorr_q(events, times, offsets, shape, nframes, dtype, qqmask, max_nnz): """ Calculation of the event correlation function using GPU @@ -346,18 +334,17 @@ def gpu_ecorr_q(events, times, offsets, shape, nframes, dtype, qqmask, max_nnz): print("start calculation of the correlation functions") from dynamix.correlator.event import EventCorrelator - t0 = time.time() + t0 = time.perf_counter() qm = qqmask.max() - res = np.zeros((qm,nframes-1),np.float32) + res = np.zeros((qm, nframes - 1), np.float32) dev = np.zeros_like(res) E = EventCorrelator(shape, nframes, max_nnz, qmask=qqmask, dtype=dtype, total_events_count=events.size) result = E.correlate(times, events, offsets) - for q in range(1,qm+1,1): - print("Number of pixels for q %d is %d" % (q,E.scale_factors[q])) - res[q-1,:] = result[q-1][1:]#1 time correlation function - dev[q-1,:] = result[q-1][1:]*1e-3#error of the 1 time correlation function - print("Total correlation time %2.2f sec" % (time.time()-t0)) + for q in range(1, qm + 1, 1): + print("Number of pixels for q %d is %d" % (q, E.scale_factors[q])) + res[q - 1,:] = result[q - 1][1:] # 1 time correlation function + dev[q - 1,:] = result[q - 1][1:] * 1e-3 # error of the 1 time correlation function + print("Total correlation time %2.2f sec" % (time.perf_counter() - t0)) trc = 0 return CorrelationResult(res, dev, trc) - diff --git a/dynamix/correlator/intensity.py b/dynamix/correlator/intensity.py index 5446311..3eb0a2c 100644 --- a/dynamix/correlator/intensity.py +++ b/dynamix/correlator/intensity.py @@ -55,7 +55,7 @@ def _determine_kernel(self): def _setup_kernels(self): self.kernel_name = self._determine_kernel() self.compile_kernels( - get_opencl_srcfile("correlator.cl"), + "dynamix:opencl/correlator.cl", compile_options=[ "-DIMAGE_WIDTH=%d" % self.shape[1], "-DDTYPE=%s" % self.c_dtype, diff --git a/dynamix/correlator/test/test_event.py b/dynamix/correlator/test/test_event.py index a6e6714..853ccda 100644 --- a/dynamix/correlator/test/test_event.py +++ b/dynamix/correlator/test/test_event.py @@ -28,50 +28,120 @@ __authors__ = ["P. Paleo"] __license__ = "MIT" -__date__ = "04/09/2019" +__date__ = "02/07/2021" -from time import time +import time import logging import unittest -import numpy as np -from dynamix.test.utils import XPCSDataset -from dynamix.correlator.event import EventCorrelator, FramesCompressor - -logging.basicConfig(level=logging.DEBUG) +import numpy +from ...test.utils import XPCSDataset, DatasetDescription +from ..event import EventCorrelator, FramesCompressor, OpenclCompressor +from ...tools.decorators import timeit +np = numpy logger = logging.getLogger(__name__) -logger.setLevel(logging.DEBUG) + +@timeit +def build_random_stack(shape=(10,11), nframes=111, dtype="uint8", nnz=10): + dtype = numpy.dtype(dtype) + npix = numpy.prod(shape) + maxi = numpy.iinfo(dtype).max + stack = np.zeros((nframes,)+(shape), dtype=dtype) + nnz2d = numpy.random.randint(1,nnz, size=shape) + pix_ptr = numpy.cumsum(nnz2d) + last = pix_ptr[-1] + pix_ptr = numpy.concatenate(([0],pix_ptr)) + values = numpy.random.randint(0, maxi, last) + times = numpy.random.randint(0,nframes,last) + for i in range(npix): + line = stack[:, i//shape[-1], i%shape[-1]] + line[times[pix_ptr[i]:pix_ptr[i+1]]] = values[pix_ptr[i]:pix_ptr[i+1]] + return stack + +def build_fake_dataset(shape=(233,231), nframes=513, dtype="uint8", nnz=17): + "Build a completely fake dataset to test sparsification" + data = build_random_stack(shape, nframes, dtype, nnz) + + self = XPCSDataset() + self.name = "fake" + descr = {"name":"fake", + "data":data, + "result":None, "description":f"Fake dataset with nnz={nnz}", + "nframes":nframes, + "dtype":dtype, + "frame_shape":shape, + "bins":None} + self.dataset_desc = DatasetDescription(**descr) + self.data = data + self.data_file = None + self.result_file = None + self.result = None + return self class TestEventDataStructure(unittest.TestCase): def setUp(self): - logger.debug("Loading data") - self.dataset = XPCSDataset("andor_1024_3k") - self.shape = self.dataset.dataset_desc.frame_shape - self.nframes = self.dataset.dataset_desc.nframes - self.compressor = FramesCompressor(self.shape, self.nframes, 330, np.int8) + logger.debug("TestEventDataStructure.setUp") + #self.dataset = XPCSDataset("andor_1024_3k") + #self.nnz = 330 + self.nnz = 10 + self.dtype = np.dtype("int8") + self.shape = (101, 103) + self.nframes = 257 + self.dataset = build_fake_dataset(dtype=self.dtype, shape=self.shape, nnz=self.nnz, nframes=self.nframes) def tearDown(self): - pass - + self.dataset = self.shape = self.nframes = self.dtype = self.nnz = None + + @property + def numpy_compressor(self): + "Single use compressor" + return FramesCompressor(self.shape, self.nframes, self.nnz, self.dtype) + + @property + def opencl_compressor(self): + "Single use compressor" + return OpenclCompressor(self.shape, self.nframes, self.nnz, self.dtype) + + + @timeit def compute_reference_datastructure(self): - logger.debug("Computing reference data compaction to events") - return self.compressor.compress_all_stack(self.dataset.data) - + return self.numpy_compressor.compress_all_stack(self.dataset.data) def test_progressive_compression(self): # Simulate progressive acquisition + compaction of frames - logger.debug("Computing progressive data compaction") - for frame in self.dataset.data: - self.compressor.process_frame(frame) + logger.debug("test_progressive_compression") + # Compare with reference implementation (compact all frames in a single pass) + t0 = time.perf_counter() ref_data, ref_times, ref_offsets = self.compute_reference_datastructure() - data, times, offsets = self.compressor.get_compacted_events() - - self.assertTrue(np.allclose(data, ref_data)) - self.assertTrue(np.allclose(times, ref_times)) - self.assertTrue(np.allclose(offsets, ref_offsets)) - + t1 = time.perf_counter() + + for compressor in (self.numpy_compressor, self.opencl_compressor): + name = compressor.__class__.__name__ + logger.debug("Working with %s", name) + for frame in self.dataset.data: + compressor.process_frame(frame) + data, times, offsets = compressor.get_compacted_events() + # print(f"data: {data.shape}, {data.dtype}") + # print(f"times: {times.shape}, {times.dtype}") + # print(f"offsets: {offsets.shape}, {offsets.dtype}") + + self.assertTrue(np.allclose(data, ref_data), msg=name) + self.assertTrue(np.allclose(times, ref_times), msg=name) + self.assertTrue(np.allclose(offsets, ref_offsets), msg=name) + + # test full-compression with OpenCL + compressor.reset() + t2 = time.perf_counter() + compressor.process_stack(self.dataset.data) + data, times, offsets = compressor.get_compacted_events() + t3 = time.perf_counter() + logger.info("OpenCL provides a speed-up of a factor %.3fx", (t1-t0)/(t3-t2)) + self.assertTrue(np.allclose(data, ref_data), msg="process_stack") + self.assertTrue(np.allclose(times, ref_times), msg="process_stack") + self.assertTrue(np.allclose(offsets, ref_offsets), msg="process_stack") + class TestEvent(unittest.TestCase): @@ -86,15 +156,17 @@ def setUp(self): self.shape = self.dataset.dataset_desc.frame_shape self.nframes = self.dataset.dataset_desc.nframes self.ref = self.dataset.result + self.events_struct = None self.compact_frames() self.tol = 5e-3 def tearDown(self): - pass + self.dataset = None + self.events_struct = None def compact_frames(self): logger.debug("Compacting frames") - self.frames_compressor = FramesCompressor( + self.frames_compressor = OpenclCompressor( self.shape, self.nframes, self.max_nnz, @@ -126,23 +198,21 @@ def test_event_correlator(self): profile=True ) # Correlate - t0 = time() + t0 = time.perf_counter() res = self.correlator.correlate( vol_times, vol_data, offsets ) - logger.info("OpenCL event correlator took %.1f ms" % ((time() - t0)*1e3)) + logger.info("OpenCL event correlator took %.1f ms", ((time.perf_counter() - t0)*1e3)) self.compare(res, "OpenCL event correlator") def suite(): testsuite = unittest.TestSuite() - testsuite.addTest( - # unittest.defaultTestLoader.loadTestsFromTestCase(TestEventDataStructure) - unittest.defaultTestLoader.loadTestsFromTestCase(TestEvent) - ) + testsuite.addTest(unittest.defaultTestLoader.loadTestsFromTestCase(TestEventDataStructure)) + testsuite.addTest(unittest.defaultTestLoader.loadTestsFromTestCase(TestEvent)) return testsuite diff --git a/dynamix/resources/__init__.py b/dynamix/resources/__init__.py index e69de29..fbeeb73 100644 --- a/dynamix/resources/__init__.py +++ b/dynamix/resources/__init__.py @@ -0,0 +1,137 @@ +# coding: utf-8 +# /*########################################################################## +# +# Copyright (C) 2016-2018 European Synchrotron Radiation Facility +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +# THE SOFTWARE. +# +# ###########################################################################*/ +"""Access project's data and documentation files. + +All access to data and documentation files MUST be made through the functions +of this modules to ensure access across different distribution schemes: + +- Installing from source or from wheel +- Installing package as a zip (through the use of pkg_resources) +- Linux packaging willing to install data files (and doc files) in + alternative folders. In this case, this file must be patched. +- Frozen fat binary application using pyFAI (frozen with cx_Freeze or py2app). + This needs special care for the resource files in the setup: + +- With cx_Freeze, add pyFAI/resources to include_files: + +.. code-block:: python + + import pyFAI.resources + pyFAI_include_files = (os.path.dirname(pyFAI.resources.__file__), + os.path.join('pyFAI', 'resources')) + setup(..., + options={'build_exe': {'include_files': [pyFAI_include_files]}} + ) + +- With py2app, add pyFAI in the packages list of the py2app options: + +.. code-block:: python + + setup(..., + options={'py2app': {'packages': ['pyFAI']}} + ) +""" + +__authors__ = ["V.A. Sole", "Thomas Vincent"] +__license__ = "MIT" +__date__ = "02/07/2021" + + +import os +import sys +import logging + +logger = logging.getLogger(__name__) + +# pkg_resources is useful when this package is stored in a zip +# When pkg_resources is not available, the resources dir defaults to the +# directory containing this module. +try: + import pkg_resources +except ImportError: + logger.debug("Backtrace", exc_info=True) + pkg_resources = None + + +# For packaging purpose, patch this variable to use an alternative directory +# E.g., replace with _RESOURCES_DIR = '/usr/share/pyFAI/data' +_RESOURCES_DIR = None + +# For packaging purpose, patch this variable to use an alternative directory +# E.g., replace with _RESOURCES_DIR = '/usr/share/pyFAI/doc' +# Not in use, uncomment when functionnality is needed +# _RESOURCES_DOC_DIR = None + +# cx_Freeze forzen support +# See http://cx-freeze.readthedocs.io/en/latest/faq.html#using-data-files +if getattr(sys, 'frozen', False): + # Running in a frozen application: + # We expect resources to be located either in a pyFAI/resources/ dir + # relative to the executable or within this package. + _dir = os.path.join(os.path.dirname(sys.executable), 'dynamix', 'resources') + if os.path.isdir(_dir): + _RESOURCES_DIR = _dir + + +def resource_filename(resource): + """Return filename corresponding to resource. + + resource can be the name of either a file or a directory. + The existence of the resource is not checked. + + :param str resource: Resource path relative to resource directory + using '/' path separator. + :return: Absolute resource path in the file system + """ + # Not in use, uncomment when functionnality is needed + # If _RESOURCES_DOC_DIR is set, use it to get resources in doc/ subflodler + # from an alternative directory. + # if _RESOURCES_DOC_DIR is not None and (resource is 'doc' or + # resource.startswith('doc/')): + # # Remove doc folder from resource relative path + # return os.path.join(_RESOURCES_DOC_DIR, *resource.split('/')[1:]) + + if _RESOURCES_DIR is not None: # if set, use this directory + return os.path.join(_RESOURCES_DIR, *resource.split('/')) + elif pkg_resources is None: # Fallback if pkg_resources is not available + return os.path.join(os.path.abspath(os.path.dirname(__file__)), + *resource.split('/')) + else: # Preferred way to get resources as it supports zipfile package + return pkg_resources.resource_filename(__name__, resource) + + +_integrated = False + + +def silx_integration(): + """Provide `dynamix` resources accessible throug silx using a prefix.""" + global _integrated + if _integrated: + return + import silx.resources + silx.resources.register_resource_directory("dynamix", + __name__, + _RESOURCES_DIR) + _integrated = True diff --git a/dynamix/resources/opencl/compactor.cl b/dynamix/resources/opencl/compactor.cl new file mode 100644 index 0000000..bba8ca5 --- /dev/null +++ b/dynamix/resources/opencl/compactor.cl @@ -0,0 +1,46 @@ +/* Compact some frames into the sparse data-structure. + * + * + */ + +/* Launched with (image_width, image_height) grid of threads. + * One thread handle one line of events, so threads are indexes by frame pixel indices. + * + * Parameters: + * slab: contains a sub-stack of input images + * timestamp_first and timestamp_last are the indices of begin and end of the sub-stack. + * image_width and image_height contain the size of the images + * nnz is the maximum expected number of events in a pixel + * times_array is am (y,x,nnz) array with the timestamps of pixels with signal + * data_array is am (y,x,nnz) array with the values associated with times_array + * counter is an (y,x) array with the number of active points found in the stack + * + */ +kernel void compactor1( const global DTYPE* slab, + const int timestamp_first, + const int timestamp_last, + const int image_width, + const int image_height, + const int nnz, + global DTYPE* data_array, + global int* times_array, + global uint* counter){ + uint x = get_global_id(0); + uint y = get_global_id(1); + if ((x >= image_width) || (y >= image_height)){ + return; + } + uint pos = y*image_width + x; + uint cnt = counter[pos]; + for (int k=0; k0){ + if (cnt0] = 0 - saxs = saxs[mask<1] - qh = qh[mask<1] - qmax = np.arange(int(qh.min()),int(qh.max())+1,1)#this is correct + q = np.float32(np.sqrt(X ** 2 + Y ** 2)) + qh = np.int16(q + 0.5) # better match with data + # qh = np.int16(q)#better match with pyfai + q[mask > 0] = 0 + saxs = saxs[mask < 1] + qh = qh[mask < 1] + qmax = np.arange(int(qh.min()), int(qh.max()) + 1, 1) # this is correct ring_brightness, radius = np.histogram(qh, weights=saxs, bins=qmax) rings, radius = np.histogram(qh, bins=qmax) - radi = np.zeros((len(radius)-1,2)) - radi[:,0] = radius[:-1]#(radius[:-1]+radius[1:])/2.0 - radi[:,1] = ring_brightness/rings - new_saxs = q*0 - f1 = q-np.array(q,np.uint16) - ind = np.array(q,np.uint16)-int(radius[0]) - ind[mask>0] = 0 - val = radi[:,1] - val = np.append(val,val[-2:]) - ind[ind>radius[-1]]=0 - #print(len(val),ind.max()) - new_saxs[mask<1] = val[ind[mask<1]+1]*f1[mask<1] + val[ind[mask<1]]*(1-f1[mask<1]) + radi = np.zeros((len(radius) - 1, 2)) + radi[:, 0] = radius[:-1] # (radius[:-1]+radius[1:])/2.0 + radi[:, 1] = ring_brightness / rings + new_saxs = q * 0 + f1 = q - np.array(q, np.uint16) + ind = np.array(q, np.uint16) - int(radius[0]) + ind[mask > 0] = 0 + val = radi[:, 1] + val = np.append(val, val[-2:]) + ind[ind > radius[-1]] = 0 + # print(len(val),ind.max()) + new_saxs[mask < 1] = val[ind[mask < 1] + 1] * f1[mask < 1] + val[ind[mask < 1]] * (1 - f1[mask < 1]) return radi, q, new_saxs + ######### Beam center ####### -def beam_center(data,mask,cx,cy,lc=60,lw=10): +def beam_center(data, mask, cx, cy, lc=60, lw=10): """ Finding the beam center :param data: 2D array of saxs pattern @@ -62,86 +64,87 @@ def beam_center(data,mask,cx,cy,lc=60,lw=10): :return: cx, cy refined beam center """ - data = np.ma.array(data,mask=mask) - ind = np.where((data>0)&(mask<1)) + data = np.ma.array(data, mask=mask) + ind = np.where((data > 0) & (mask < 1)) - rad, r_q, new_saxs = radi(data,mask,cx,cy)#radial averaging - err = np.abs(data-new_saxs)[ind].mean()/np.mean(data[ind])#error - print("Initial center cx=%.2f, cy=%.2f, err=%1.5f" % (cx,cy,err)) + rad, r_q, new_saxs = radi(data, mask, cx, cy) # radial averaging + err = np.abs(data - new_saxs)[ind].mean() / np.mean(data[ind]) # error + print("Initial center cx=%.2f, cy=%.2f, err=%1.5f" % (cx, cy, err)) ### show the stripes #### - sdata = np.zeros(data.shape,np.uint32) - sdata[int(cy-lc-lw/2):int(cy-lc+lw/2+1),:] += np.uint32(1) - sdata[int(cy+lc-lw/2):int(cy+lc+lw/2+1),:] += np.uint32(1) - sdata[:,int(cx-lc-lw/2):int(cx-lc+lw/2+1)] += np.uint32(1) - sdata[:,int(cx+lc-lw/2):int(cx+lc+lw/2+1)] += np.uint32(1) + sdata = np.zeros(data.shape, np.uint32) + sdata[int(cy - lc - lw / 2):int(cy - lc + lw / 2 + 1),:] += np.uint32(1) + sdata[int(cy + lc - lw / 2):int(cy + lc + lw / 2 + 1),:] += np.uint32(1) + sdata[:, int(cx - lc - lw / 2):int(cx - lc + lw / 2 + 1)] += np.uint32(1) + sdata[:, int(cx + lc - lw / 2):int(cx + lc + lw / 2 + 1)] += np.uint32(1) plt.figure() with np.errstate(divide="ignore", invalid="ignore"): - plt.imshow(np.log10(data),cmap='jet') - plt.imshow(sdata,cmap='gray_r',alpha=0.3) - plt.plot(cx,cy,'r+') + plt.imshow(np.log10(data), cmap='jet') + plt.imshow(sdata, cmap='gray_r', alpha=0.3) + plt.plot(cx, cy, 'r+') ##### Find horizontal center x ########################## - vl1 = np.sum(data[:,int(cx-lc-lw/2):int(cx-lc+lw/2+1)],1)/(lw+1)#vertical line 1 + vl1 = np.sum(data[:, int(cx - lc - lw / 2):int(cx - lc + lw / 2 + 1)], 1) / (lw + 1) # vertical line 1 verr0 = 1e+6 ### pixel precision #### - for llc in range(lc-10,lc+10,1): - vl2 = np.sum(data[:,int(cx+llc-lw/2):int(cx+llc+lw/2+1)],1)/(lw+1)#vertical line - verr = np.mean(np.abs(vl1-vl2)) - if verr=0: - vl2 = (vl20-data[:,int(cx+nlc-lw/2)]*f+data[:,int(cx+nlc+lw/2+1)+1]*f)/(lw+1) + # for f in np.arange(-0.99,1.0,0.01): + for f in np.arange(-0.51, 0.52, 0.01): + if f >= 0: + vl2 = (vl20 - data[:, int(cx + nlc - lw / 2)] * f + data[:, int(cx + nlc + lw / 2 + 1) + 1] * f) / (lw + 1) else: - vl2 = (vl20-data[:,int(cx+nlc-lw/2-1)]*f+data[:,int(cx+nlc+lw/2+1)]*f)/(lw+1) - verr = np.mean(np.abs(vl1-vl2)) - if verr=0: - vl2 = (vl20-data[int(cy+nlc-lw/2),:]*f+data[int(cy+nlc+lw/2+1)+1,:]*f)/(lw+1) + if f >= 0: + vl2 = (vl20 - data[int(cy + nlc - lw / 2),:] * f + data[int(cy + nlc + lw / 2 + 1) + 1,:] * f) / (lw + 1) else: - vl2 = (vl20-data[int(cy+nlc-lw/2-1),:]*f+data[int(cy+nlc+lw/2+1),:]*f)/(lw+1) - verr = np.mean(np.abs(vl1-vl2)) - if verr=par: - ntmp = (tmp[:-1,:]+tmp[1:,:])/2 - for i in range(0,par,2): - nt.append(ntmp[i,0]) - nd.append(ntmp[i,1]) - nse.append(ntmp[i,2]) + nt.append(d[i, 0]) + nd.append(d[i, 1]) + nse.append(d[i, 2]) + while len(tmp[:, 0]) >= par: + ntmp = (tmp[:-1,:] + tmp[1:,:]) / 2 + for i in range(0, par, 2): + nt.append(ntmp[i, 0]) + nd.append(ntmp[i, 1]) + nse.append(ntmp[i, 2]) tmp = ntmp[par:-1:2,:3] - x = np.array([nt,nd,nse]).T + x = np.array([nt, nd, nse]).T return x @@ -177,7 +180,6 @@ def make_q(config): sname = config["sample_description"]["name"] - #### Data location ###################################################################################### datdir = config["data_location"]["data_dir"] @@ -191,7 +193,6 @@ def make_q(config): df2 = int(config["data_location"]["last_dark"]) savdir = config["data_location"]["result_dir"] - #### Experimental setup ###################################################################################### geometry = config["exp_setup"]["geometry"] @@ -208,7 +209,6 @@ def make_q(config): #### Correlator info ###################################################################################### - correlator = config["correlator"]["method"] lth = int(config["correlator"]["low_threshold"]) bADU = int(config["correlator"]["bottom_ADU"]) @@ -224,14 +224,14 @@ def make_q(config): flatfield_file = config["detector"]["flatfield"] ################################################################################################################### print("Making qs") - + try: - data = readdata.readnpz(savdir+sname+"_2D.npz") + data = readdata.readnpz(savdir + sname + "_2D.npz") except: - print("Cannot read "+savdir+sname+"_2D.npz") + print("Cannot read " + savdir + sname + "_2D.npz") sys.exit() - #if beamstop_mask != 'none': + # if beamstop_mask != 'none': # try: # bmask = np.abs(EdfMethods.loadedf(beamstop_mask))#reads edf and npy # bmask[bmask>1] = 1 @@ -240,59 +240,55 @@ def make_q(config): ### read beamstop mask #### bmask = read_beamstop_mask(beamstop_mask) - + ### read detector mask ### - mask = read_det_mask(mask_file,detector) - - mask[bmask>0] = 1 - - data[mask>0] = 0 + mask = read_det_mask(mask_file, detector) + mask[bmask > 0] = 1 - ind = np.where((data>0)&(mask<1)) + data[mask > 0] = 0 - data = np.ma.array(data,mask=mask) - qmask = mask*0 + ind = np.where((data > 0) & (mask < 1)) + data = np.ma.array(data, mask=mask) + qmask = mask * 0 - #t0 = time.time() - rad, r_q, new_saxs = radi(data,mask,cx,cy)#radial averaging - #print("Calculation time %3.4f sec" % (time.time()-t0)) + # t0 = time.time() + rad, r_q, new_saxs = radi(data, mask, cx, cy) # radial averaging + # print("Calculation time %3.4f sec" % (time.time()-t0)) new_saxs[np.isnan(new_saxs)] = 0 - np.save(savdir+sname+"_gaus.npy",np.array(new_saxs,np.float32)) - + np.save(savdir + sname + "_gaus.npy", np.array(new_saxs, np.float32)) - rad[:,0] = 4*np.pi/lambdaw*np.sin(np.arctan(pix_size*rad[:,0]/distance)/2.0)#q vector calculation - r_q = 4*np.pi/lambdaw*np.sin(np.arctan(pix_size*r_q/distance)/2.0)#q vector calculation - width_p = np.tan(np.arcsin(width_q*lambdaw/4/np.pi)*2)*distance/pix_size - qmask = np.array((r_q-first_q+width_q/2)/width_q+1,np.uint16) + rad[:, 0] = 4 * np.pi / lambdaw * np.sin(np.arctan(pix_size * rad[:, 0] / distance) / 2.0) # q vector calculation + r_q = 4 * np.pi / lambdaw * np.sin(np.arctan(pix_size * r_q / distance) / 2.0) # q vector calculation + width_p = np.tan(np.arcsin(width_q * lambdaw / 4 / np.pi) * 2) * distance / pix_size + qmask = np.array((r_q - first_q + width_q / 2) / width_q + 1, np.uint16) print("Number of Qs %d" % number_q) - #print("Width of ROI is %1.1f pixels" % width_p) - np.savetxt(savdir+sname+"_1D.dat",rad) - - #qmask = np.array((r_q-first_q+width_q/2)/width_q+1,np.uint16) - #qmask[mask>0] = 0 - #np.save(savdir+sname+"_qmask.npy",np.array(qmask,np.uint16)) - #qmask[qmask>number_q] = 0 - - - #qp = np.linspace(first_q,first_q+(number_q-1)*width_q,number_q) - qp = np.linspace(first_q,first_q+(number_q-1)*step_q,number_q) - qmask = mask*0 - i_qp = qp*0 + # print("Width of ROI is %1.1f pixels" % width_p) + np.savetxt(savdir + sname + "_1D.dat", rad) + + # qmask = np.array((r_q-first_q+width_q/2)/width_q+1,np.uint16) + # qmask[mask>0] = 0 + # np.save(savdir+sname+"_qmask.npy",np.array(qmask,np.uint16)) + # qmask[qmask>number_q] = 0 + + # qp = np.linspace(first_q,first_q+(number_q-1)*width_q,number_q) + qp = np.linspace(first_q, first_q + (number_q - 1) * step_q, number_q) + qmask = mask * 0 + i_qp = qp * 0 i_bands = [] n = 0 for i in qp: - i_qp[n] = rad[(np.abs(rad[:,0]-i)==(np.abs(rad[:,0]-i).min())),1] - ind1 = np.where((rad[:,0]>=i-width_q/2)&(rad[:,0]<=i+width_q/2))[0] + i_qp[n] = rad[(np.abs(rad[:, 0] - i) == (np.abs(rad[:, 0] - i).min())), 1] + ind1 = np.where((rad[:, 0] >= i - width_q / 2) & (rad[:, 0] <= i + width_q / 2))[0] i_bands.append(ind1) n += 1 - indq = np.where(np.abs(r_q-i)<=width_q/2) + indq = np.where(np.abs(r_q - i) <= width_q / 2) qmask[indq] = n - qmask[mask>0] = 0 - np.save(savdir+sname+"_qmask.npy",np.array(qmask,np.uint16)) + qmask[mask > 0] = 0 + np.save(savdir + sname + "_qmask.npy", np.array(qmask, np.uint16)) return @@ -306,11 +302,11 @@ def test_dir(dir_name): test_savdir = os.path.dirname(dir_name) if not os.path.exists(test_savdir): os.makedirs(test_savdir) - print("Create",test_savdir) + print("Create", test_savdir) return -def format_result(CorrelationResult,qqmask,flatfield,cdata,dt,ttcf_par): +def format_result(CorrelationResult, qqmask, flatfield, cdata, dt, ttcf_par): """ format the obained result to usable form :param CorrelationResult: nametuple containing the results @@ -325,38 +321,39 @@ def format_result(CorrelationResult,qqmask,flatfield,cdata,dt,ttcf_par): res = [] n = 0 nframes = CorrelationResult.res[0].size - cf = np.zeros((nframes,3),np.float32) - cf[:,0] = np.arange(1,nframes+1,1)*dt + cf = np.zeros((nframes, 3), np.float32) + cf[:, 0] = np.arange(1, nframes + 1, 1) * dt number_q = qqmask.max() - for q in range(1,number_q+1,1): - fcorrection = (flatfield[qqmask==q]**2).mean()/flatfield[qqmask==q].mean()**2 - correction = (cdata[qqmask==q]**2).mean()/cdata[qqmask==q].mean()**2 * fcorrection - cf[:,1] = CorrelationResult.res[n,:]/correction # make correct baseline#CPU - cf[:,2] = CorrelationResult.dev[n,:]/correction # make correct baseline#CPU + for q in range(1, number_q + 1, 1): + fcorrection = (flatfield[qqmask == q] ** 2).mean() / flatfield[qqmask == q].mean() ** 2 + correction = (cdata[qqmask == q] ** 2).mean() / cdata[qqmask == q].mean() ** 2 * fcorrection + cf[:, 1] = CorrelationResult.res[n,:] / correction # make correct baseline#CPU + cf[:, 2] = CorrelationResult.dev[n,:] / correction # make correct baseline#CPU if ttcf_par == q: - trc = CorrelationResult.trc/correction + trc = CorrelationResult.trc / correction cfm = cftomt(cf) res.append(cfm) if q == 1: save_cf = cfm else: - save_cf = np.append(save_cf,cfm[:,1:],axis=1) + save_cf = np.append(save_cf, cfm[:, 1:], axis=1) n += 1 trc = CorrelationResult.trc return res, save_cf, trc - -def save_cf(file_name,save_cf,qp): - q_title='#q values 1/A:' + +def save_cf(file_name, save_cf, qp): + q_title = '#q values 1/A:' for q in qp: - q_title = q_title+" %.5f " % q - q_title=q_title+'\n' - f=open(file_name,'w') + q_title = q_title + " %.5f " % q + q_title = q_title + '\n' + f = open(file_name, 'w') f.write(q_title) np.savetxt(f, save_cf) f.close() return + def events(data, mNp): """ convert 3D data matrix to events list @@ -365,36 +362,36 @@ def events(data, mNp): :return: pixels, s list of pixels and intgrated intensity per frame """ - #t0 = time.time() + # t0 = time.time() s = [] pixels = [] sshape = np.shape(data) - if len(sshape) == 3 : + if len(sshape) == 3: nframes, nx, ny = np.shape(data) - nx = nx*ny + nx = nx * ny for i in range(nframes): matr = np.ravel(data[i,:,:]) - msumpix,mpix = eigerpix(matr,mNp,nx) + msumpix, mpix = eigerpix(matr, mNp, nx) mpix = mpix[:msumpix] pixels.append(mpix) s.append(msumpix) - if len(sshape) == 2: + if len(sshape) == 2: nx, ny = np.shape(data) - nx = nx*ny + nx = nx * ny matr = np.ravel(data) - msumpix,mpix = eigerpix(matr,mNp,nx) + msumpix, mpix = eigerpix(matr, mNp, nx) mpix = mpix[:msumpix] pixels.append(mpix) s.append(msumpix) - #print("Compaction time %f" % (time.time()-t0)) + # print("Compaction time %f" % (time.time()-t0)) - return pixels, s + return pixels, s -def make_cdata(file_name,config): +def make_cdata(file_name, config): """ make smooth data :param file_name: filename that can contain cdata @@ -402,15 +399,16 @@ def make_cdata(file_name,config): :return: cdata, smooth 2D image for correction """ if os.path.isfile(file_name): - cdata = np.array(np.load(file_name),np.float32) + cdata = np.array(np.load(file_name), np.float32) else: make_q(config) - cdata = np.array(np.load(file_name),np.float32) + cdata = np.array(np.load(file_name), np.float32) cdata[np.isnan(cdata)] = 0 return cdata -def read_det_mask(det_mask,detector): + +def read_det_mask(det_mask, detector): """ read detector mask file :param det_mask: filename that contains detector mask @@ -419,26 +417,27 @@ def read_det_mask(det_mask,detector): """ if det_mask != 'none': try: - dmask = np.abs(EdfMethods.loadedf(det_mask))#reads edf and npy - dmask[dmask>1] = 1 + dmask = np.abs(EdfMethods.loadedf(det_mask)) # reads edf and npy + dmask[dmask > 1] = 1 except: print("Cannot read detector mask %s, skip" % det_mask) else: if detector == 'Andor': - dshape=(1024,1024) + dshape = (1024, 1024) elif detector == 'maxipix': - dshape=(516,516) + dshape = (516, 516) elif detector == 'eiger500K': - dshape=(1024,514) + dshape = (1024, 514) elif detector == 'eiger4m': - dshape=(2162,2068) - else: + dshape = (2162, 2068) + else: print("Needs detector mask file") - sys.exit() - dmask = np.zeros(dshape,dtype='uint8') + sys.exit() + dmask = np.zeros(dshape, dtype='uint8') return dmask -def read_qmask(qmask_file,mask, number_q): + +def read_qmask(qmask_file, mask, number_q): """ read qmask file :param qmask_file: filename that contains qmask @@ -450,18 +449,19 @@ def read_qmask(qmask_file,mask, number_q): if qmask_file != 'none': try: qqmask = EdfMethods.loadedf(qmask_file) - qqmask[mask>0] = 0 - qqmask[qqmask>number_q] = 0 + qqmask[mask > 0] = 0 + qqmask[qqmask > number_q] = 0 except: print("Cannot read qmask %s, skip" % qmask_file) - qqmask = mask*0 - qqmask[mask<1] = 1 + qqmask = mask * 0 + qqmask[mask < 1] = 1 else: - qqmask = mask*0 - qqmask[mask<1] = 1 - + qqmask = mask * 0 + qqmask[mask < 1] = 1 + return qqmask + def read_beamstop_mask(beamstop_mask): """ read beamstop mask file @@ -471,15 +471,16 @@ def read_beamstop_mask(beamstop_mask): """ if beamstop_mask != 'none': try: - bmask = np.abs(EdfMethods.loadedf(beamstop_mask))#reads edf and npy - bmask[bmask>1] = 1 + bmask = np.abs(EdfMethods.loadedf(beamstop_mask)) # reads edf and npy + bmask[bmask > 1] = 1 except: print("Cannot read beamstop mask %s, skip" % beamstop_mask) else: bmask = 0 return bmask - -def reduce_matrix(data,qqmask,cdata,flatfield): + + +def reduce_matrix(data, qqmask, cdata, flatfield): """ read beamstop mask file :param data: 3D array of frames @@ -489,19 +490,20 @@ def reduce_matrix(data,qqmask,cdata,flatfield): :return: data,qqmask,cdata,flatfield reduced array of original data """ - - ind = np.where(qqmask>0) + + ind = np.where(qqmask > 0) sindy = ind[0].min() eindy = ind[0].max() sindx = ind[1].min() eindx = ind[1].max() - data = data[:,sindy:eindy,sindx:eindx] - qqmask = qqmask[sindy:eindy,sindx:eindx] - cdata = cdata[sindy:eindy,sindx:eindx] - flatfield = flatfield[sindy:eindy,sindx:eindx] - print("Reduced data size is %2.2f Gigabytes" % (data.size*data.itemsize/1024**3)) + data = data[:, sindy:eindy, sindx:eindx] + qqmask = qqmask[sindy:eindy, sindx:eindx] + cdata = cdata[sindy:eindy, sindx:eindx] + flatfield = flatfield[sindy:eindy, sindx:eindx] + print("Reduced data size is %2.2f Gigabytes" % (data.size * data.itemsize / 1024 ** 3)) + + return data, qqmask, cdata, flatfield - return data,qqmask,cdata,flatfield def data_compaction(data): """ compact data for event correlation @@ -512,22 +514,22 @@ def data_compaction(data): """ t0 = time.time() events = [] - times = [] - offsets = [0] - t,x,y = data.shape - y = x*y - data = np.reshape(data,(t,y)) + times = [] + offsets = [0] + t, x, y = data.shape + y = x * y + data = np.reshape(data, (t, y)) for pix in range(y): - ind = np.where(data[:,pix]>0) - events.append(data[ind[0],pix]) + ind = np.where(data[:, pix] > 0) + events.append(data[ind[0], pix]) times.append(ind[0]) offsets.append(ind[0].size) - events = np.array(np.concatenate(events).ravel(),np.int8) - times = np.array(np.concatenate(times).ravel(),np.int32) - offsets = np.array(np.cumsum(offsets),np.uint32) - print("Compaction time %f" % (time.time()-t0)) + events = np.array(np.concatenate(events).ravel(), np.int8) + times = np.array(np.concatenate(times).ravel(), np.int32) + offsets = np.array(np.cumsum(offsets), np.uint32) + print("Compaction time %f" % (time.time() - t0)) return events, times, offsets diff --git a/dynamix/utils.py b/dynamix/utils.py index 47ad8d2..60a7159 100644 --- a/dynamix/utils.py +++ b/dynamix/utils.py @@ -2,6 +2,16 @@ from bisect import bisect import os import numpy as np +from collections import namedtuple + +CorrelationResult = namedtuple("CorrelationResult", "res dev") +Compacted = namedtuple("Compacted", "values times pixel_ptr") +""" Compacted form of a XPCS stack containing 3 1D structure: + - values: pixel values for nonzero data points. dtype: same as reference stack + - times: time indices corresponding to nonzero data points. dtype: uint32 + - pixel_ptr: Pixel k's information is stored in values[pixel_ptr[k]:pixel_ptr[k+1]] and times[pixel_ptr[k]:pixel_ptr[k+1]] +""" + def get_folder_path(foldername=""): _file_dir = os.path.dirname(os.path.realpath(__file__)) diff --git a/run_tests.py b/run_tests.py index 6007344..0636b17 100755 --- a/run_tests.py +++ b/run_tests.py @@ -1,4 +1,4 @@ -#!/usr/bin/env python +#!/usr/bin/env python3 # coding: utf-8 # /*########################################################################## # @@ -32,7 +32,7 @@ """ __authors__ = ["Jérôme Kieffer", "Thomas Vincent"] -__date__ = "02/03/2018" +__date__ = "02/07/2021" __license__ = "MIT" import distutils.util diff --git a/setup.py b/setup.py index 0328317..6ede548 100644 --- a/setup.py +++ b/setup.py @@ -1,4 +1,4 @@ -#!/usr/bin/env python +#!/usr/bin/env python3 #from setuptools import setup, Extension from numpy.distutils.core import setup, Extension