mirror of https://github.com/silx-kit/pyFAI.git
Add sparse builder cython class
This commit is contained in:
Valentin Valls
parent
74277e48d6
commit
66c63663b4
|
|
@ -42,3 +42,4 @@
|
|||
|
||||
#discard generated C-files
|
||||
/pyFAI/ext/*.c
|
||||
/pyFAI/ext/*.cpp
|
||||
|
|
|
|||
|
|
@ -24,7 +24,7 @@
|
|||
|
||||
__authors__ = ["V. Valls"]
|
||||
__license__ = "MIT"
|
||||
__date__ = "31/08/2017"
|
||||
__date__ = "10/07/2018"
|
||||
|
||||
from numpy.distutils.misc_util import Configuration
|
||||
import platform
|
||||
|
|
@ -100,6 +100,13 @@ def configuration(parent_package='', top_path=None):
|
|||
for ext_config in ext_modules:
|
||||
config.add_extension(**ext_config)
|
||||
|
||||
config.add_extension('sparse_builder',
|
||||
sources=['sparse_builder.pyx'],
|
||||
include_dirs=[numpy.get_include()],
|
||||
language='c++',
|
||||
extra_link_args=['-fopenmp'],
|
||||
extra_compile_args=['-fopenmp'])
|
||||
|
||||
return config
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -0,0 +1,342 @@
|
|||
import numpy
|
||||
cimport numpy as cnumpy
|
||||
|
||||
from libcpp.vector cimport vector
|
||||
from libcpp.list cimport list as clist
|
||||
from libcpp cimport bool
|
||||
from libc.math cimport fabs
|
||||
cimport libc.stdlib
|
||||
cimport libc.string
|
||||
|
||||
from cython.parallel import prange
|
||||
from cython.operator cimport dereference
|
||||
from cython.operator cimport preincrement
|
||||
cimport cython
|
||||
from cython cimport floating
|
||||
|
||||
cdef double EPS32 = (1.0 + numpy.finfo(numpy.float32).eps)
|
||||
|
||||
|
||||
cdef struct pixel_t:
|
||||
cnumpy.int32_t index
|
||||
cnumpy.float32_t coef
|
||||
|
||||
cdef cppclass PixelElementaryBlock:
|
||||
vector[cnumpy.int32_t] _indexes
|
||||
vector[cnumpy.float32_t] _coefs
|
||||
int _size
|
||||
int _max_size
|
||||
|
||||
PixelElementaryBlock(int size) nogil:
|
||||
this._indexes.reserve(size)
|
||||
this._coefs.reserve(size)
|
||||
this._size = 0
|
||||
this._max_size = size
|
||||
|
||||
void push(pixel_t &pixel) nogil:
|
||||
this._indexes.push_back(pixel.index)
|
||||
this._coefs.push_back(pixel.coef)
|
||||
this._size += 1
|
||||
|
||||
int size() nogil:
|
||||
return this._size
|
||||
|
||||
bool is_full() nogil:
|
||||
return this._size == this._max_size
|
||||
|
||||
cdef cppclass PixelBlock:
|
||||
clist[PixelElementaryBlock*] _blocks
|
||||
int _block_size
|
||||
|
||||
PixelBlock(int block_size) nogil:
|
||||
this._block_size = block_size
|
||||
|
||||
__dealloc__() nogil:
|
||||
cdef:
|
||||
PixelElementaryBlock* element
|
||||
int i = 0
|
||||
clist[PixelElementaryBlock*].iterator it
|
||||
it = this._blocks.begin()
|
||||
while it != this._blocks.end():
|
||||
element = dereference(it)
|
||||
del element
|
||||
preincrement(it)
|
||||
this._blocks.clear()
|
||||
|
||||
void push(pixel_t &pixel) nogil:
|
||||
cdef:
|
||||
PixelElementaryBlock *block
|
||||
if _blocks.size() == 0 or this._blocks.back().is_full():
|
||||
block = new PixelElementaryBlock(size=this._block_size)
|
||||
this._blocks.push_back(block)
|
||||
block = this._blocks.back()
|
||||
block.push(pixel)
|
||||
|
||||
int size() nogil:
|
||||
cdef:
|
||||
int size = 0
|
||||
clist[PixelElementaryBlock*].iterator it
|
||||
it = this._blocks.begin()
|
||||
while it != this._blocks.end():
|
||||
size += dereference(it).size()
|
||||
preincrement(it)
|
||||
return size
|
||||
|
||||
cnumpy.int32_t[:] index_array():
|
||||
cdef:
|
||||
clist[PixelElementaryBlock*].iterator it
|
||||
PixelElementaryBlock* block
|
||||
int size
|
||||
int begin
|
||||
cnumpy.int32_t[:] data
|
||||
cnumpy.int32_t[:] data2
|
||||
|
||||
size = this.size()
|
||||
data = numpy.empty(size, dtype=numpy.int32)
|
||||
|
||||
begin = 0
|
||||
it = this._blocks.begin()
|
||||
while it != this._blocks.end():
|
||||
block = dereference(it)
|
||||
if block.size() != 0:
|
||||
data2 = numpy.array(block._indexes, dtype=numpy.int32)
|
||||
data[begin:begin + block.size()] = data2
|
||||
begin += block.size()
|
||||
preincrement(it)
|
||||
return data
|
||||
|
||||
cnumpy.float32_t[:] coef_array():
|
||||
cdef:
|
||||
clist[PixelElementaryBlock*].iterator it
|
||||
PixelElementaryBlock* block
|
||||
int size
|
||||
int begin
|
||||
cnumpy.float32_t[:] data
|
||||
cnumpy.float32_t[:] data2
|
||||
|
||||
size = this.size()
|
||||
data = numpy.empty(size, dtype=numpy.float32)
|
||||
|
||||
begin = 0
|
||||
it = this._blocks.begin()
|
||||
while it != this._blocks.end():
|
||||
block = dereference(it)
|
||||
if block.size() != 0:
|
||||
data2 = numpy.array(block._coefs, dtype=numpy.float32)
|
||||
data[begin:begin + block.size()] = data2
|
||||
begin += block.size()
|
||||
preincrement(it)
|
||||
return data
|
||||
|
||||
|
||||
cdef cppclass PixelBin:
|
||||
clist[pixel_t] _pixels
|
||||
PixelBlock *_pixels_in_block
|
||||
|
||||
PixelBin(int block_size) nogil:
|
||||
if block_size > 0:
|
||||
this._pixels_in_block = new PixelBlock(block_size)
|
||||
else:
|
||||
this._pixels_in_block = NULL
|
||||
|
||||
__dealloc__() nogil:
|
||||
if this._pixels_in_block != NULL:
|
||||
del this._pixels_in_block
|
||||
this._pixels_in_block = NULL
|
||||
else:
|
||||
this._pixels.clear()
|
||||
|
||||
void push(pixel_t &pixel) nogil:
|
||||
if this._pixels_in_block != NULL:
|
||||
this._pixels_in_block.push(pixel)
|
||||
else:
|
||||
this._pixels.push_back(pixel)
|
||||
|
||||
int size() nogil:
|
||||
if this._pixels_in_block != NULL:
|
||||
return this._pixels_in_block.size()
|
||||
else:
|
||||
return this._pixels.size()
|
||||
|
||||
cnumpy.int32_t[:] index_array():
|
||||
cdef:
|
||||
int i = 0
|
||||
clist[pixel_t].iterator it_points
|
||||
|
||||
if this._pixels_in_block != NULL:
|
||||
return this._pixels_in_block.index_array()
|
||||
|
||||
data = numpy.empty(this.size(), dtype=numpy.int32)
|
||||
it_points = this._pixels.begin()
|
||||
while it_points != this._pixels.end():
|
||||
data[i] = dereference(it_points).index
|
||||
preincrement(it_points)
|
||||
i += 1
|
||||
return data
|
||||
|
||||
cnumpy.float32_t[:] coef_array():
|
||||
cdef:
|
||||
int i = 0
|
||||
clist[pixel_t].iterator it_points
|
||||
|
||||
if this._pixels_in_block != NULL:
|
||||
return this._pixels_in_block.coef_array()
|
||||
|
||||
data = numpy.empty(this.size(), dtype=numpy.float32)
|
||||
it_points = this._pixels.begin()
|
||||
while it_points != this._pixels.end():
|
||||
data[i] = dereference(it_points).coef
|
||||
preincrement(it_points)
|
||||
i += 1
|
||||
return data
|
||||
|
||||
|
||||
cdef class SparseBuilder(object):
|
||||
|
||||
cdef PixelBin **_bins
|
||||
cdef int _nbin
|
||||
cdef int _block_size
|
||||
|
||||
def __init__(self, nbin, block_size=512):
|
||||
self._block_size = block_size
|
||||
self._nbin = nbin
|
||||
self._bins = <PixelBin **>libc.stdlib.malloc(self._nbin * sizeof(PixelBin *))
|
||||
libc.string.memset(self._bins, 0, self._nbin * sizeof(PixelBin *))
|
||||
|
||||
def __dealloc__(self):
|
||||
cdef:
|
||||
PixelBin *pixel_bin
|
||||
int i
|
||||
for i in range(self._nbin):
|
||||
pixel_bin = self._bins[i]
|
||||
if pixel_bin != NULL:
|
||||
del pixel_bin
|
||||
libc.stdlib.free(self._bins)
|
||||
|
||||
@cython.boundscheck(False)
|
||||
@cython.wraparound(False)
|
||||
@cython.cdivision(True)
|
||||
cdef PixelBin *_create_bin(self) nogil:
|
||||
return new PixelBin(self._block_size)
|
||||
|
||||
@cython.boundscheck(False)
|
||||
@cython.wraparound(False)
|
||||
@cython.cdivision(True)
|
||||
cdef void cinsert(self, int bin_id, int index, cnumpy.float32_t coef) nogil:
|
||||
cdef:
|
||||
pixel_t pixel
|
||||
PixelBin *pixel_bin
|
||||
if bin_id < 0 or bin_id >= self._nbin:
|
||||
return
|
||||
pixel.index = index
|
||||
pixel.coef = coef
|
||||
|
||||
pixel_bin = self._bins[bin_id]
|
||||
if pixel_bin == NULL:
|
||||
pixel_bin = self._create_bin()
|
||||
self._bins[bin_id] = pixel_bin
|
||||
self._bins[bin_id].push(pixel)
|
||||
|
||||
def insert(self, bin_id, index, coef):
|
||||
if bin_id < 0 or bin_id >= self._nbin:
|
||||
raise ValueError("bin_id out of range")
|
||||
self.cinsert(bin_id, index, coef)
|
||||
|
||||
@cython.boundscheck(False)
|
||||
@cython.wraparound(False)
|
||||
@cython.cdivision(True)
|
||||
def get_bin_coefs(self, bin_id):
|
||||
cdef:
|
||||
PixelBin *pixel_bin
|
||||
if bin_id < 0 or bin_id >= self._nbin:
|
||||
raise ValueError("bin_id out of range")
|
||||
pixel_bin = self._bins[bin_id]
|
||||
if pixel_bin == NULL:
|
||||
return numpy.empty(shape=(0, 1), dtype=numpy.float32)
|
||||
return numpy.array(pixel_bin.coef_array())
|
||||
|
||||
@cython.boundscheck(False)
|
||||
@cython.wraparound(False)
|
||||
@cython.cdivision(True)
|
||||
def get_bin_indexes(self, bin_id):
|
||||
cdef:
|
||||
PixelBin *pixel_bin
|
||||
if bin_id < 0 or bin_id >= self._nbin:
|
||||
raise ValueError("bin_id out of range")
|
||||
pixel_bin = self._bins[bin_id]
|
||||
if pixel_bin == NULL:
|
||||
return numpy.empty(shape=(0, 1), dtype=numpy.int32)
|
||||
return numpy.array(pixel_bin.index_array())
|
||||
|
||||
@cython.boundscheck(False)
|
||||
@cython.wraparound(False)
|
||||
@cython.cdivision(True)
|
||||
def get_bin_size(self, bin_id):
|
||||
cdef:
|
||||
PixelBin *pixel_bin
|
||||
if bin_id < 0 or bin_id >= self._nbin:
|
||||
raise ValueError("bin_id out of range")
|
||||
pixel_bin = self._bins[bin_id]
|
||||
if pixel_bin == NULL:
|
||||
return 0
|
||||
return pixel_bin.size()
|
||||
|
||||
@cython.boundscheck(False)
|
||||
@cython.wraparound(False)
|
||||
@cython.cdivision(True)
|
||||
def size(self):
|
||||
cdef:
|
||||
PixelBin *pixel_bin
|
||||
int size
|
||||
int bin_id
|
||||
|
||||
size = 0
|
||||
for bin_id in range(self._nbin):
|
||||
pixel_bin = self._bins[bin_id]
|
||||
if pixel_bin != NULL:
|
||||
size += pixel_bin.size()
|
||||
return size
|
||||
|
||||
@cython.boundscheck(False)
|
||||
@cython.wraparound(False)
|
||||
@cython.cdivision(True)
|
||||
def to_csr(self):
|
||||
cdef:
|
||||
cnumpy.int32_t[:] indexes
|
||||
cnumpy.int32_t[:] indexes2
|
||||
cnumpy.float32_t[:] coefs
|
||||
cnumpy.float32_t[:] coefs2
|
||||
cnumpy.int32_t[:] nbins
|
||||
PixelBin *pixel_bin
|
||||
int size
|
||||
int i
|
||||
int begin, end
|
||||
int bin_id
|
||||
int bin_size
|
||||
|
||||
# indexes of the first and the last+1 elements of each bins
|
||||
size = 0
|
||||
nbins = numpy.empty(self._nbin + 1, dtype=numpy.int32)
|
||||
nbins[0] = size
|
||||
for bin_id in range(self._nbin):
|
||||
pixel_bin = self._bins[bin_id]
|
||||
if pixel_bin != NULL:
|
||||
bin_size = pixel_bin.size()
|
||||
else:
|
||||
bin_size = 0
|
||||
size += bin_size
|
||||
nbins[bin_id + 1] = size
|
||||
|
||||
indexes = numpy.empty(size, dtype=numpy.int32)
|
||||
coefs = numpy.empty(size, dtype=numpy.float32)
|
||||
|
||||
for bin_id in range(self._nbin):
|
||||
pixel_bin = self._bins[bin_id]
|
||||
if pixel_bin == NULL or pixel_bin.size() == 0:
|
||||
continue
|
||||
begin = nbins[bin_id]
|
||||
end = nbins[bin_id + 1]
|
||||
indexes[begin:end] = pixel_bin.index_array()
|
||||
coefs[begin:end] = pixel_bin.coef_array()
|
||||
|
||||
return numpy.array(coefs), numpy.array(indexes), numpy.array(nbins)
|
||||
|
|
@ -0,0 +1,65 @@
|
|||
import numpy
|
||||
cimport numpy as cnumpy
|
||||
|
||||
from libcpp.vector cimport vector
|
||||
from libcpp.list cimport list as clist
|
||||
from libcpp cimport bool
|
||||
from libc.math cimport fabs
|
||||
cimport libc.stdlib
|
||||
cimport libc.string
|
||||
|
||||
from cython.parallel import prange
|
||||
from cython.operator cimport dereference
|
||||
from cython.operator cimport preincrement
|
||||
cimport cython
|
||||
from cython cimport floating
|
||||
|
||||
include "sparse_builder.pxi"
|
||||
|
||||
|
||||
@cython.cdivision(True)
|
||||
cdef inline floating calc_upper_bound(floating maximum_value) nogil:
|
||||
if maximum_value > 0:
|
||||
return maximum_value * EPS32
|
||||
else:
|
||||
return maximum_value / EPS32
|
||||
|
||||
|
||||
@cython.cdivision(True)
|
||||
cdef inline floating get_bin_number(floating x0, floating pos0_min, floating delta) nogil:
|
||||
return (x0 - pos0_min) / delta
|
||||
|
||||
|
||||
@cython.cdivision(True)
|
||||
@cython.boundscheck(False)
|
||||
@cython.wraparound(False)
|
||||
def feed_histogram(SparseBuilder builder not None,
|
||||
cnumpy.ndarray pos not None,
|
||||
cnumpy.ndarray weights not None,
|
||||
int bins=100,
|
||||
double empty=0.0,
|
||||
double normalization_factor=1.0):
|
||||
assert pos.size == weights.size
|
||||
assert bins > 1
|
||||
cdef:
|
||||
int size = pos.size
|
||||
cnumpy.float32_t[::1] cpos = numpy.ascontiguousarray(pos.ravel(), dtype=numpy.float32)
|
||||
cnumpy.float32_t delta, min0, max0, maxin0
|
||||
cnumpy.float32_t a = 0.0
|
||||
cnumpy.float32_t d = 0.0
|
||||
cnumpy.float32_t fbin = 0.0
|
||||
cnumpy.float32_t tmp_count, tmp_data = 0.0
|
||||
cnumpy.float32_t epsilon = 1e-10
|
||||
int bin = 0, i, idx
|
||||
|
||||
min0 = pos.min()
|
||||
maxin0 = pos.max()
|
||||
max0 = calc_upper_bound(maxin0)
|
||||
delta = (max0 - min0) / float(bins)
|
||||
|
||||
with nogil:
|
||||
for i in range(size):
|
||||
a = cpos[i]
|
||||
fbin = get_bin_number(a, min0, delta)
|
||||
bin = < int > fbin
|
||||
builder.cinsert(bin, i, 1.0)
|
||||
Loading…
Reference in New Issue