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.

  1. Allocator.hh: There are two types of allocators;
    1. "Stack" comes with preallocated memory that the algorithm can use without needing to malloc or free. Stacks are used in the DRP.
    2. "Heap" behaves as a regular heap. Memory is managed via malloc and free. Heaps are used in the Python mode.
  2. AllocArray.hh: There are two types of allocated arrays designed to use the Stack or the Heap:
    1. The AllocArray() class is a multi-dimensional array that inherits from the XtcData::Array parent class. 
    2. 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:

  1. An allocator is required in the constructor.
  2. The allocator is stored as a private member.
  3. A method to switch to a new allocator is needed.
  4. 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

 

 

  • No labels