We describe a pattern for writing an algorithm that can operate in both the DRP mode and Python mode. In the DRP mode, the algorithm is written in C++ and the memory allocation is done on the PEBBLE (a ring buffer of preallocated memory) with the idea of removing malloc and free commands that hurt performance. In Python mode, the algorithm uses the standard heap and exposes the output arrays to python (without copying its contents). In psalg, we propose to expose on arrays of fundamental types and the ownership of the output arrays gets passed on to python.
Here is an example of converting peak finder v3r3 algorithm into a DRP/python compatible code. Let's start with two classes that the algorithm needs to include in the C++ header.
- Allocator.hh: There are two types of allocators;
- "Stack" comes with preallocated memory that the algorithm can use without needing to malloc or free. Stacks are used in the DRP.
- "Heap" behaves as a regular heap. Memory is managed via malloc and free. Heaps are used in the Python mode.
- AllocArray.hh: There are two types of allocated arrays designed to use the Stack or the Heap:
- The AllocArray() class is a multi-dimensional array that inherits from the XtcData::Array parent class.
- The AllocArray1D() class is a 1-dimensional array that has methods that mimic the std::vector.
// lcls2/psalg/psalg/include/Allocator.hh class Stack:public Allocator class Heap:public Allocator // lcls2/psalg/psalg/include/AllocArray.hh template <typename T> class AllocArray{...} template <typename T> class AllocArray1D{...}
PeakFinderAlgos class was modified in the following way:
- An allocator is required in the constructor.
- The allocator is stored as a private member.
- A method to switch to a new allocator is needed.
- Replace all std::vector<T> with AllocArray1D<T>
// lcls2/psalg/psalg/include/peakFinderAlgos.h PeakFinderAlgos(Allocator *allocator, ...); // constructor private: Allocator *m_allocator; // a private member of PeakFinderAlgos public: void setAllocator(Allocator *allocator); // a method to switch to a new allocator per event AllocArray1D<float> rows; // Replaces std::vector<float> rows;
Cython code for exposing C++ AllocArray1D to python without copying. You will need PyAllocArray1D class:
# lcls2/psana/psana/peakFinder/peakFinder.pyx cdef class PyAllocArray1D: cdef void* shape_ptr cdef void* data_ptr cdef int size cdef int typenum cdef init(self, void* buffer_ptr, void* data_ptr, int size, int typenum): self.shape_ptr = buffer_ptr self.data_ptr = data_ptr self.size = size self.typenum = typenum def __array__(self): """ this is called when numpy converts this object to np.array.""" cdef cnp.npy_intp shape[1] shape[0] = <cnp.npy_intp> self.size # Create a 1D array, of length 'size' ndarray = cnp.PyArray_SimpleNewFromData(1, shape, self.typenum, self.data_ptr) # Increment the reference count, as the above assignement was done in # C, and Python does not know that there is this additional reference Py_INCREF(self) return ndarray def __dealloc__(self): """ Frees the array. This is called by Python when all the references to the object are gone. """ free(<void*>self.shape_ptr)
Here's how to use PyAllocArray1D in Cython. In this code snippet, we expose C++ peakFinderAlgos.rows AllocArray1D to Python rows_cgrav as a numpy array. NOTE: When a Python user calls getPeaks(), rows_cgrav is returned that persists till the user no longer needs it:
# lcls2/psana/psana/peakFinder/peakFinder.pyx def getPeaks(self): cdef cnp.ndarray rows_cgrav arr0 = PyAllocArray1D() # self.cptr is a pointer to peakFinderAlgos object arr0.init(<void*> self.cptr.rows.shape(), <void*> self.cptr.rows.data(), self.cptr.numPeaksSelected, cnp.NPY_FLOAT) rows_cgrav = np.array(arr0, copy=False) rows_cgrav.base = <PyObject*> arr0 self.cptr.rows.incRefCnt() # increment ref count so that C++ doesn't delete array return rows_cgrav
Here's setup.py for the peakFinder algorithm:
# lcls2/psana/setup.py from Cython.Build import cythonize ext = Extension("peakFinder", packages=['psana.peakfinder',], sources=["psana/peakFinder/peakFinder.pyx", "../psalg/psalg/src/PeakFinderAlgos.cpp", "../psalg/psalg/src/LocalExtrema.cpp"], language="c++", extra_compile_args=['-std=c++11'], include_dirs=[np.get_include(), "../install/include"] ) setup(name="peakFinder", ext_modules=cythonize(ext))
For wrapping C++ objects of non-fundamental types, here's the pattern in Cython for exposing it to python (courtesy of Mikhail Dubrovin). This makes a copy of the object and hands it off to python.
# lcls2/psana/psana/peakFinder/peakFinder.pyx cdef extern from "../../../psalg/psalg/include/PeakFinderAlgos.h" namespace "psalgos": cdef cppclass Peak : Peak() except + Peak(const Peak& o) except + Peak operator=(const Peak& rhs) except + # https://groups.google.com/forum/#!topic/cython-users/39Nwqsksdto @staticmethod cdef factory(Peak cpp_obj): py_obj = py_peak.__new__(py_peak, _make_obj=False) (<py_peak> py_obj).cptr = new Peak(cpp_obj) # C++ copy constructor return py_obj