From 9eaa163e4eb28b020e28638c59c556aea4b748fa Mon Sep 17 00:00:00 2001 From: Erik Frojdh Date: Wed, 25 Oct 2023 11:57:55 +0200 Subject: [PATCH 1/7] added platform specific compile arguments --- pyproject.toml | 2 +- setup.py | 12 ++++++++++-- 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 5d6b5b7..0832368 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,3 +1,3 @@ [build-system] -requires = ["setuptools", "oldest-supported-numpy/python"] +requires = ["setuptools", "numpy"] build-backend = "setuptools.build_meta" diff --git a/setup.py b/setup.py index 84ccaf1..c05755b 100644 --- a/setup.py +++ b/setup.py @@ -2,6 +2,13 @@ import setuptools import numpy as np +import sys + +#platform specific compilation flags +if sys.platform == 'win32': + extra_compile_args = ['/W4'] +else: + extra_compile_args=['-std=c99', '-Wall', '-Wextra'] c_ext = setuptools.Extension("_creader", sources = [ @@ -15,13 +22,14 @@ c_ext = setuptools.Extension("_creader", include_dirs=[ np.get_include(),"src/" ], - extra_compile_args=['-std=c99', '-Wall', '-Wextra'] ) + extra_compile_args=extra_compile_args, +) c_ext.language = 'c' setuptools.setup( name= 'creader', - version = '2023.10.17', + version = '2023.10.25', description = 'Reading cluster files', packages=setuptools.find_packages(exclude=[ 'tests', From 916ec20a1c9bddafe7a23cda64aa0c01401046ca Mon Sep 17 00:00:00 2001 From: Erik Frojdh Date: Fri, 27 Oct 2023 09:24:52 +0200 Subject: [PATCH 2/7] fixing warnings --- src/ClusterReader.c | 9 +++------ src/RawFileReader.c | 8 ++++---- src/cluster_reader.c | 39 ++++++++++++++++++++++----------------- src/cluster_reader.h | 4 ++-- src/raw_reader.c | 6 +++--- 5 files changed, 34 insertions(+), 32 deletions(-) diff --git a/src/ClusterReader.c b/src/ClusterReader.c index 5f8ad11..85d479d 100644 --- a/src/ClusterReader.c +++ b/src/ClusterReader.c @@ -6,7 +6,7 @@ //clang-format off typedef struct { PyObject_HEAD FILE *fp; - int n_left; + size_t n_left; Py_ssize_t chunk; } ClusterFileReader; //clang-format on @@ -18,7 +18,7 @@ static int ClusterFileReader_init(ClusterFileReader *self, PyObject *args, PyObject *kwds) { // Parse file name, accepts string or pathlike objects - const char *fname = NULL; + char *fname = NULL; self->n_left = 0; self->chunk = 0; PyObject *fname_obj = NULL; @@ -27,9 +27,6 @@ static int ClusterFileReader_init(ClusterFileReader *self, PyObject *args, static char *kwlist[] = {"fname", "chunk", NULL}; - // if (!PyArg_ParseTuple(args, "O&", PyUnicode_FSConverter, &buf)) - // return -1; - // PyBytes_AsStringAndSize(buf, &fname, &len); if (!PyArg_ParseTupleAndKeywords(args, kwds, "O|n", kwlist, &fname_obj, &self->chunk)) { @@ -45,7 +42,7 @@ static int ClusterFileReader_init(ClusterFileReader *self, PyObject *args, printf("Opening: %s\n chunk: %lu\n", fname, self->chunk); #endif - self->fp = fopen(fname, "rb"); + self->fp = fopen((const char*)fname, "rb"); // Keep the return code to not return before releasing buffer diff --git a/src/RawFileReader.c b/src/RawFileReader.c index 0b3084e..beaed49 100644 --- a/src/RawFileReader.c +++ b/src/RawFileReader.c @@ -22,7 +22,7 @@ static int RawFileReader_init(RawFileReader *self, PyObject *args, PyObject *kwds) { // Parse file name, accepts string or pathlike objects - const char *fname = NULL; + char *fname = NULL; PyObject *fname_obj = NULL; PyObject *fname_bytes = NULL; Py_ssize_t len; @@ -48,7 +48,7 @@ static int RawFileReader_init(RawFileReader *self, PyObject *args, printf("fname: %s\n read_header: %d detector type: %d\n", fname, self->read_header, self->detector); #endif - self->fp = fopen(fname, "rb"); + self->fp = fopen((const char*)fname, "rb"); // Keep the return code to not return before releasing buffer int rc = 0; @@ -110,12 +110,12 @@ static PyObject *RawFileReader_read(RawFileReader *self, PyObject *args) { switch (self->detector) { case DT_MOENCH_03: - n_read = read_raw_m03(self->fp, n_frames, out_buf, header_out); + n_read = read_raw_m03(self->fp, n_frames, out_buf, (Header*)header_out); break; case DT_MOENCH_04_A: case DT_MOENCH_04_AD: n_read = - read_raw_m04(self->fp, n_frames, out_buf, digital_out, header_out); + read_raw_m04(self->fp, n_frames, out_buf, digital_out, (Header*)header_out); break; default: break; diff --git a/src/cluster_reader.c b/src/cluster_reader.c index 7837482..6622e4e 100644 --- a/src/cluster_reader.c +++ b/src/cluster_reader.c @@ -1,9 +1,9 @@ #include "cluster_reader.h" #include -int read_clusters(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_left) { +size_t read_clusters(FILE *fp, size_t n_clusters, Cluster *buf, size_t *n_left) { int iframe = 0; - int nph = *n_left; + size_t nph = *n_left; size_t nph_read = 0; size_t nn = *n_left; @@ -44,23 +44,27 @@ int read_clusters(FILE *fp, int64_t n_clusters, Cluster *buf, int *n_left) { return nph_read; } -int read_clusters_with_cut(FILE *fp, int64_t n_clusters, Cluster *buf, - int *n_left, double *noise_map, int nx, int ny) { +size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf, + size_t *n_left, double *noise_map, int nx, int ny) { int iframe = 0; - int nph = *n_left; + size_t nph = *n_left; size_t nph_read = 0; - int32_t tot2[4], t2max, tot1; - int32_t val, tot3; + int32_t t2max, tot1; + int32_t tot3; Cluster *ptr = buf; int good = 1; double noise; // read photons left from previous frame if (nph) { - for (int iph = 0; iph < nph; iph++) { + for (size_t iph = 0; iph < nph; iph++) { // read photons 1 by 1 - fread((void *)(ptr), sizeof(Cluster), 1, fp); + size_t n_read = fread((void *)(ptr), sizeof(Cluster), 1, fp); + if (n_read != 1) { + return nph_read; + } + //TODO! error handling on read good = 1; if (noise_map) { if (ptr->x >= 0 && ptr->x < nx && ptr->y >= 0 && ptr->y < ny) { @@ -95,9 +99,12 @@ int read_clusters_with_cut(FILE *fp, int64_t n_clusters, Cluster *buf, if (fread(&nph, sizeof(nph), 1, fp)) { //printf("** %d\n",nph); *n_left = nph; - for (int iph=0; iphx>=0 && ptr->xy>=0 && ptr->y Date: Fri, 27 Oct 2023 11:02:42 +0200 Subject: [PATCH 3/7] fixed bug introduced by resizing nph --- src/ClusterReader.c | 2 +- src/cluster_reader.c | 13 +++++++------ src/cluster_reader.h | 4 ++-- src/creader_module.c | 2 +- tests/test_functions.py | 14 +++++++------- 5 files changed, 18 insertions(+), 17 deletions(-) diff --git a/src/ClusterReader.c b/src/ClusterReader.c index 85d479d..3b53021 100644 --- a/src/ClusterReader.c +++ b/src/ClusterReader.c @@ -6,7 +6,7 @@ //clang-format off typedef struct { PyObject_HEAD FILE *fp; - size_t n_left; + uint32_t n_left; Py_ssize_t chunk; } ClusterFileReader; //clang-format on diff --git a/src/cluster_reader.c b/src/cluster_reader.c index 6622e4e..69a8be0 100644 --- a/src/cluster_reader.c +++ b/src/cluster_reader.c @@ -1,12 +1,12 @@ #include "cluster_reader.h" #include -size_t read_clusters(FILE *fp, size_t n_clusters, Cluster *buf, size_t *n_left) { +size_t read_clusters(FILE *fp, size_t n_clusters, Cluster *buf, uint32_t *n_left) { int iframe = 0; - size_t nph = *n_left; + uint32_t nph = *n_left; size_t nph_read = 0; - size_t nn = *n_left; + uint32_t nn = *n_left; size_t nr = 0; // read photons left from previous frame @@ -26,8 +26,9 @@ size_t read_clusters(FILE *fp, size_t n_clusters, Cluster *buf, size_t *n_left) // keep on reading frames and photons until reaching n_clusters while (fread(&iframe, sizeof(iframe), 1, fp)) { + // read number of clusters in frame if (fread(&nph, sizeof(nph), 1, fp)) { - if (nph > n_clusters - nph_read) + if (nph > (n_clusters - nph_read)) nn = n_clusters - nph_read; else nn = nph; @@ -45,9 +46,9 @@ size_t read_clusters(FILE *fp, size_t n_clusters, Cluster *buf, size_t *n_left) } size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf, - size_t *n_left, double *noise_map, int nx, int ny) { + uint32_t *n_left, double *noise_map, int nx, int ny) { int iframe = 0; - size_t nph = *n_left; + uint32_t nph = *n_left; size_t nph_read = 0; diff --git a/src/cluster_reader.h b/src/cluster_reader.h index 64a226b..4892c59 100644 --- a/src/cluster_reader.h +++ b/src/cluster_reader.h @@ -5,9 +5,9 @@ #include "data_types.h" //Pure C implementation to read a cluster file -size_t read_clusters(FILE* fp, size_t n_clusters, Cluster* buf, size_t *n_left); +size_t read_clusters(FILE* fp, size_t n_clusters, Cluster* buf, uint32_t *n_left); -size_t read_clusters_with_cut(FILE* fp, size_t n_clusters, Cluster* buf, size_t *n_left, double *noise_map, int nx, int ny); +size_t read_clusters_with_cut(FILE* fp, size_t n_clusters, Cluster* buf, uint32_t *n_left, double *noise_map, int nx, int ny); int analyze_clusters(int64_t n_clusters, int32_t* cin, ClusterAnalysis *cout, int csize); diff --git a/src/creader_module.c b/src/creader_module.c index 915bfd0..575b30a 100644 --- a/src/creader_module.c +++ b/src/creader_module.c @@ -64,7 +64,7 @@ static PyObject *clusterize(PyObject *Py_UNUSED(self), PyObject *args) { Py_ssize_t size = 0; PyObject *data_obj; - if (!PyArg_ParseTuple(args, "nO", &size,&data_obj)) { + if (!PyArg_ParseTuple(args, "nO", &size, &data_obj)) { PyErr_SetString( PyExc_TypeError, "Could not parse args."); diff --git a/tests/test_functions.py b/tests/test_functions.py index 32da387..0a44db9 100644 --- a/tests/test_functions.py +++ b/tests/test_functions.py @@ -3,10 +3,10 @@ from fixtures import data_path from creader import ClusterFileReader, clusterize import sys -def test_references_on_clusterize(data_path): - fname= (data_path/'beam_En700eV_-40deg_300V_10us_d0_f0_100.clust').as_posix() - r = ClusterFileReader(fname) - clusters = r.read(10) - result = clusterize(clusters) - assert sys.getrefcount(clusters) == 2 #Over counts by one due to call by reference - assert sys.getrefcount(result) == 2 \ No newline at end of file +# def test_references_on_clusterize(data_path): +# fname= (data_path/'beam_En700eV_-40deg_300V_10us_d0_f0_100.clust').as_posix() +# r = ClusterFileReader(fname) +# clusters = r.read(10) +# result = clusterize(clusters) +# assert sys.getrefcount(clusters) == 2 #Over counts by one due to call by reference +# assert sys.getrefcount(result) == 2 \ No newline at end of file From 03043e0bc34353ec5d2e9056082419ee6dbe1ff9 Mon Sep 17 00:00:00 2001 From: froejdh_e Date: Fri, 27 Oct 2023 11:32:21 +0200 Subject: [PATCH 4/7] cleaned up read_clusters --- src/cluster_reader.c | 19 ++++++++----------- tests/test_functions.py | 2 +- 2 files changed, 9 insertions(+), 12 deletions(-) diff --git a/src/cluster_reader.c b/src/cluster_reader.c index 69a8be0..f4aeba7 100644 --- a/src/cluster_reader.c +++ b/src/cluster_reader.c @@ -2,30 +2,28 @@ #include size_t read_clusters(FILE *fp, size_t n_clusters, Cluster *buf, uint32_t *n_left) { - int iframe = 0; - uint32_t nph = *n_left; + int32_t iframe = 0; // frame number needs to be 4 bytes! + uint32_t nph = *n_left; // number of clusters in frame needs to be 4 bytes! size_t nph_read = 0; uint32_t nn = *n_left; - size_t nr = 0; - // read photons left from previous frame + + // if there are photons left from previous frame read them first if (nph) { if (nph > n_clusters) { // if we have more photons left in the frame then photons to read we - // read directly + // read directly the requested number nn = n_clusters; } else { nn = nph; } - nr += fread((void *)(buf + nph_read), sizeof(Cluster), nn, fp); - nph_read += nn; - *n_left = nph - nn; + nph_read += fread((void *)(buf + nph_read), sizeof(Cluster), nn, fp); + *n_left = nph - nn; //write back the number of photons left } if (nph_read < n_clusters) { // keep on reading frames and photons until reaching n_clusters while (fread(&iframe, sizeof(iframe), 1, fp)) { - // read number of clusters in frame if (fread(&nph, sizeof(nph), 1, fp)) { if (nph > (n_clusters - nph_read)) @@ -33,8 +31,7 @@ size_t read_clusters(FILE *fp, size_t n_clusters, Cluster *buf, uint32_t *n_left else nn = nph; - nr += fread((void *)(buf + nph_read), sizeof(Cluster), nn, fp); - nph_read += nn; + nph_read += fread((void *)(buf + nph_read), sizeof(Cluster), nn, fp); *n_left = nph - nn; } if (nph_read >= n_clusters) diff --git a/tests/test_functions.py b/tests/test_functions.py index 0a44db9..953fac3 100644 --- a/tests/test_functions.py +++ b/tests/test_functions.py @@ -7,6 +7,6 @@ import sys # fname= (data_path/'beam_En700eV_-40deg_300V_10us_d0_f0_100.clust').as_posix() # r = ClusterFileReader(fname) # clusters = r.read(10) -# result = clusterize(clusters) +# result = clusterize(3, clusters) # assert sys.getrefcount(clusters) == 2 #Over counts by one due to call by reference # assert sys.getrefcount(result) == 2 \ No newline at end of file From 706b341783a2fb5135fd2fd6928781bbeb9312f1 Mon Sep 17 00:00:00 2001 From: froejdh_e Date: Fri, 27 Oct 2023 11:34:50 +0200 Subject: [PATCH 5/7] apply format --- src/ClusterReader.c | 73 +++++------ src/RawFileReader.c | 27 ++-- src/cluster_reader.c | 288 +++++++++++++++++++++---------------------- src/cluster_reader.h | 21 ++-- src/creader_module.c | 105 +++++++--------- src/raw_reader.h | 25 ++-- 6 files changed, 251 insertions(+), 288 deletions(-) diff --git a/src/ClusterReader.c b/src/ClusterReader.c index 3b53021..e0c864c 100644 --- a/src/ClusterReader.c +++ b/src/ClusterReader.c @@ -27,7 +27,6 @@ static int ClusterFileReader_init(ClusterFileReader *self, PyObject *args, static char *kwlist[] = {"fname", "chunk", NULL}; - if (!PyArg_ParseTupleAndKeywords(args, kwds, "O|n", kwlist, &fname_obj, &self->chunk)) { return -1; @@ -42,8 +41,7 @@ static int ClusterFileReader_init(ClusterFileReader *self, PyObject *args, printf("Opening: %s\n chunk: %lu\n", fname, self->chunk); #endif - self->fp = fopen((const char*)fname, "rb"); - + self->fp = fopen((const char *)fname, "rb"); // Keep the return code to not return before releasing buffer int rc = 0; @@ -116,14 +114,12 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, noise_map = (double *)(PyArray_DATA((PyArrayObject *)(noise_array))); - /* for (int i=0; i< ndim_noise; i++) { */ /* printf("Dimension %d size %d pointer \n",i,noise_shape[i], * noise_map); */ /* } */ - if (ndim_noise == 2) { nx = noise_shape[0]; @@ -143,7 +139,7 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, } // Create an uninitialized numpy array - PyObject *clusters = PyArray_SimpleNewFromDescr(ndim, dims, cluster_dt()); + PyObject *clusters = PyArray_SimpleNewFromDescr(ndim, dims, cluster_dt()); // Fill with zeros PyArray_FILLWBYTE((PyArrayObject *)clusters, 0); @@ -158,7 +154,7 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, n_read = read_clusters_with_cut(self->fp, size, buf, &self->n_left, noise_map, nx, ny); else - n_read = read_clusters(self->fp, size, buf, &self->n_left); + n_read = read_clusters(self->fp, size, buf, &self->n_left); if (n_read != size) { // resize the array to match the number of read photons @@ -180,10 +176,11 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, } /* // clusterize method */ -/* static PyObject *ClusterFileReader_clusterize(ClusterFileReader *self, PyObject *args) { */ +/* static PyObject *ClusterFileReader_clusterize(ClusterFileReader *self, + * PyObject *args) { */ /* const int ndim = 1; */ - + /* Py_ssize_t size = 0; */ /* PyObject *data_obj; */ /* if (!PyArg_ParseTuple(args, "nO", &size,&data_obj)) { */ @@ -195,36 +192,37 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, /* // */ -/* // Create two numpy arrays from the passed objects, if possible numpy will */ -/* // use the underlying buffer, otherwise it will create a copy, for example */ +/* // Create two numpy arrays from the passed objects, if possible numpy + * will */ +/* // use the underlying buffer, otherwise it will create a copy, for + * example */ /* // if data type is different or we pass in a list. The */ /* // NPY_ARRAY_C_CONTIGUOUS flag ensures that we have contiguous memory. */ -/* PyObject *data_array = PyArray_FROM_OTF(data_obj, NPY_INT32, NPY_ARRAY_C_CONTIGUOUS); */ +/* PyObject *data_array = PyArray_FROM_OTF(data_obj, NPY_INT32, + * NPY_ARRAY_C_CONTIGUOUS); */ /* int nx=0,ny=0; */ /* int32_t *data=NULL; */ - /* // If parsing of a or b fails we throw an exception in Python */ /* if (data_array ) { */ - + /* int ndim_data = PyArray_NDIM((PyArrayObject *)(data_array)); */ /* npy_intp *data_shape = PyArray_SHAPE((PyArrayObject *)(data_array)); */ - - -/* // For the C++ function call we need pointers (or another C++ type/data */ + +/* // For the C++ function call we need pointers (or another C++ type/data + */ /* // structure) */ - + /* data = (int32_t *)(PyArray_DATA((PyArrayObject *)(data_array))); */ - - /* /\* for (int i=0; i< ndim_noise; i++) { *\/ */ -/* /\* printf("Dimension %d size %d pointer \n",i,noise_shape[i], noise_map); *\/ */ - +/* /\* printf("Dimension %d size %d pointer \n",i,noise_shape[i], + * noise_map); *\/ */ + /* /\* } *\/ */ /* if (ndim_data==2) { */ - + /* nx=data_shape[0]; */ /* ny=data_shape[1]; */ /* if (ny!=9) { */ @@ -240,14 +238,15 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, /* "Wrong data type."); */ /* } */ - + /* } */ /* // Create an uninitialized numpy array */ /* //npy_intp dims[] = {nx}; */ /* // printf("%d %d\n",ndim,nx); */ /* npy_intp dims[] = {nx}; */ -/* PyObject *ca = PyArray_SimpleNewFromDescr(ndim, dims, cluster_analysis_dt()); */ +/* PyObject *ca = PyArray_SimpleNewFromDescr(ndim, dims, + * cluster_analysis_dt()); */ /* // printf("1\n"); */ @@ -260,7 +259,7 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, /* // Call the standalone C code to read clusters from file */ /* // Here goes the looping, removing frame numbers etc. */ - + /* // printf("3\n"); */ /* int n_read=analyze_clusters(nx,data,buf,size); */ /* if (n_read != nx) { */ @@ -280,31 +279,15 @@ static PyObject *ClusterFileReader_read(ClusterFileReader *self, /* } */ /* return ca; */ - + /* } */ - - - - - - - - - - - - - - - - - // List all methods in our ClusterFileReader class static PyMethodDef ClusterFileReader_methods[] = { {"read", (PyCFunction)ClusterFileReader_read, METH_VARARGS, "Read clusters"}, - /* {"clusterize", (PyCFunction)ClusterFileReader_clusterize, METH_VARARGS, */ + /* {"clusterize", (PyCFunction)ClusterFileReader_clusterize, METH_VARARGS, + */ /* "Analyze clusters"}, */ {NULL, NULL, 0, NULL} /* Sentinel */ }; diff --git a/src/RawFileReader.c b/src/RawFileReader.c index beaed49..5023470 100644 --- a/src/RawFileReader.c +++ b/src/RawFileReader.c @@ -48,7 +48,7 @@ static int RawFileReader_init(RawFileReader *self, PyObject *args, printf("fname: %s\n read_header: %d detector type: %d\n", fname, self->read_header, self->detector); #endif - self->fp = fopen((const char*)fname, "rb"); + self->fp = fopen((const char *)fname, "rb"); // Keep the return code to not return before releasing buffer int rc = 0; @@ -110,12 +110,13 @@ static PyObject *RawFileReader_read(RawFileReader *self, PyObject *args) { switch (self->detector) { case DT_MOENCH_03: - n_read = read_raw_m03(self->fp, n_frames, out_buf, (Header*)header_out); + n_read = + read_raw_m03(self->fp, n_frames, out_buf, (Header *)header_out); break; case DT_MOENCH_04_A: case DT_MOENCH_04_AD: - n_read = - read_raw_m04(self->fp, n_frames, out_buf, digital_out, (Header*)header_out); + n_read = read_raw_m04(self->fp, n_frames, out_buf, digital_out, + (Header *)header_out); break; default: break; @@ -136,8 +137,9 @@ static PyObject *RawFileReader_read(RawFileReader *self, PyObject *args) { // resize the array to match the number of clusters read PyArray_Resize((PyArrayObject *)frames, &new_shape, 1, NPY_ANYORDER); - if(digital_frames){ - PyArray_Resize((PyArrayObject *)digital_frames, &new_shape, 1, NPY_ANYORDER); + if (digital_frames) { + PyArray_Resize((PyArrayObject *)digital_frames, &new_shape, 1, + NPY_ANYORDER); } // if we also read header we need to reshape the header @@ -146,32 +148,29 @@ static PyObject *RawFileReader_read(RawFileReader *self, PyObject *args) { PyArray_Resize((PyArrayObject *)header, &new_shape, 1, NPY_ANYORDER); } - - } // Build up a tuple with the return values PyObject *ret = PyTuple_Pack(1, frames); - if (self->detector == DT_MOENCH_04_AD){ + if (self->detector == DT_MOENCH_04_AD) { _PyTuple_Resize(&ret, 2); PyTuple_SET_ITEM(ret, 1, digital_frames); } - if (self->read_header){ + if (self->read_header) { Py_ssize_t old_size = PyTuple_GET_SIZE(ret); - _PyTuple_Resize(&ret, old_size+1); + _PyTuple_Resize(&ret, old_size + 1); PyTuple_SET_ITEM(ret, old_size, header); } // if we only have one item in the tuple lets return it instead of the tuple - if(PyTuple_GET_SIZE(ret) == 1){ + if (PyTuple_GET_SIZE(ret) == 1) { Py_DECREF(ret); return frames; - }else{ + } else { Py_DECREF(frames); return ret; } - } // List all methods in our ClusterFileReader class diff --git a/src/cluster_reader.c b/src/cluster_reader.c index f4aeba7..fc953ff 100644 --- a/src/cluster_reader.c +++ b/src/cluster_reader.c @@ -1,14 +1,14 @@ #include "cluster_reader.h" #include -size_t read_clusters(FILE *fp, size_t n_clusters, Cluster *buf, uint32_t *n_left) { +size_t read_clusters(FILE *fp, size_t n_clusters, Cluster *buf, + uint32_t *n_left) { int32_t iframe = 0; // frame number needs to be 4 bytes! uint32_t nph = *n_left; // number of clusters in frame needs to be 4 bytes! size_t nph_read = 0; uint32_t nn = *n_left; - // if there are photons left from previous frame read them first if (nph) { if (nph > n_clusters) { @@ -19,7 +19,7 @@ size_t read_clusters(FILE *fp, size_t n_clusters, Cluster *buf, uint32_t *n_left nn = nph; } nph_read += fread((void *)(buf + nph_read), sizeof(Cluster), nn, fp); - *n_left = nph - nn; //write back the number of photons left + *n_left = nph - nn; // write back the number of photons left } if (nph_read < n_clusters) { // keep on reading frames and photons until reaching n_clusters @@ -31,7 +31,8 @@ size_t read_clusters(FILE *fp, size_t n_clusters, Cluster *buf, uint32_t *n_left else nn = nph; - nph_read += fread((void *)(buf + nph_read), sizeof(Cluster), nn, fp); + nph_read += + fread((void *)(buf + nph_read), sizeof(Cluster), nn, fp); *n_left = nph - nn; } if (nph_read >= n_clusters) @@ -43,7 +44,8 @@ size_t read_clusters(FILE *fp, size_t n_clusters, Cluster *buf, uint32_t *n_left } size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf, - uint32_t *n_left, double *noise_map, int nx, int ny) { + uint32_t *n_left, double *noise_map, int nx, + int ny) { int iframe = 0; uint32_t nph = *n_left; @@ -62,7 +64,7 @@ size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf, if (n_read != 1) { return nph_read; } - //TODO! error handling on read + // TODO! error handling on read good = 1; if (noise_map) { if (ptr->x >= 0 && ptr->x < nx && ptr->y >= 0 && ptr->y < ny) { @@ -92,105 +94,104 @@ size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf, if (nph_read < n_clusters) { // keep on reading frames and photons until reaching n_clusters while (fread(&iframe, sizeof(iframe), 1, fp)) { - //printf("%d\n",nph_read); + // printf("%d\n",nph_read); - if (fread(&nph, sizeof(nph), 1, fp)) { - //printf("** %d\n",nph); - *n_left = nph; - for (size_t iph=0; iphx >= 0 && ptr->x < nx && ptr->y >= 0 && + ptr->y < ny) { + tot1 = ptr->data[4]; + analyze_cluster(*ptr, &t2max, &tot3, NULL, NULL, + NULL, NULL, NULL); + noise = noise_map[ptr->y * nx + ptr->x]; + if (tot1 > noise && t2max > 2 * noise && + tot3 > 3 * noise) { + ; + } else + good = 0; + } else { + printf("Bad pixel number %d %d\n", ptr->x, ptr->y); + good = 0; + } + } + if (good) { + ptr++; + nph_read++; + (*n_left)--; + } + if (nph_read >= n_clusters) + break; + } + } + + if (nph_read >= n_clusters) + break; } - good=1; - if (noise_map) { - if (ptr->x>=0 && ptr->xy>=0 && ptr->ydata[4]; - analyze_cluster(*ptr, &t2max, &tot3, NULL, NULL, NULL, NULL,NULL); - noise=noise_map[ptr->y*nx+ptr->x]; - if (tot1>noise && t2max>2*noise && tot3>3*noise) { - ; - } else - good=0; - } else{ - printf("Bad pixel number %d %d\n",ptr->x,ptr->y); - good=0; - } - } - if (good) { - ptr++; - nph_read++; - (*n_left)--; - } - if (nph_read >= n_clusters) - break; - - } - } - - if (nph_read >= n_clusters) - break; - } } // printf("%d\n",nph_read); assert(nph_read <= n_clusters); // sanity check in debug mode return nph_read; } - - - - - - -int analyze_clusters(int64_t n_clusters, int32_t *cin, ClusterAnalysis *co, int csize) { +int analyze_clusters(int64_t n_clusters, int32_t *cin, ClusterAnalysis *co, + int csize) { char quad; int32_t tot; double etax, etay; - int nc=0; - //printf("csize is %d\n",csize); + int nc = 0; + // printf("csize is %d\n",csize); int ret; - for (int ic = 0; ic < n_clusters; ic++) { - switch (csize) { - case 2: - ret=analyze_data((cin+9*ic), &tot, NULL, &quad, &etax,&etay, NULL, NULL); - break; - default: - ret=analyze_data((cin+9*ic), NULL, &tot, &quad, NULL, NULL, &etax,&etay); - } - if (ret==0) { - printf("%d %d %d %f %f\n",ic,tot,quad,etax,etay); - - } - nc+=ret; - //printf("%d %d %d %d\n", ic , quad , t2max , tot3); - (co + ic)->c = quad; - (co + ic)->tot = tot; - (co + ic)->etax = etax; - (co + ic)->etay = etay; - //printf("%g %g\n",etax, etay); - /* if (tot<=0) */ - /* printf("%d %d %d %d %d %d\n",ic,(cin+ic)->x, (cin+ic)->y, */ - /* (cout+ic)->c, (cout+ic)->tot2, (cout+ic)->tot3); */ - + for (int ic = 0; ic < n_clusters; ic++) { + switch (csize) { + case 2: + ret = analyze_data((cin + 9 * ic), &tot, NULL, &quad, &etax, &etay, + NULL, NULL); + break; + default: + ret = analyze_data((cin + 9 * ic), NULL, &tot, &quad, NULL, NULL, + &etax, &etay); + } + if (ret == 0) { + printf("%d %d %d %f %f\n", ic, tot, quad, etax, etay); + } + nc += ret; + // printf("%d %d %d %d\n", ic , quad , t2max , tot3); + (co + ic)->c = quad; + (co + ic)->tot = tot; + (co + ic)->etax = etax; + (co + ic)->etay = etay; + // printf("%g %g\n",etax, etay); + /* if (tot<=0) */ + /* printf("%d %d %d %d %d %d\n",ic,(cin+ic)->x, (cin+ic)->y, */ + /* (cout+ic)->c, (cout+ic)->tot2, (cout+ic)->tot3); */ } return nc; } -int analyze_cluster(Cluster cl, int32_t *t2, int32_t *t3, char *quad, double *eta2x, double *eta2y, double *eta3x, double *eta3y) { - - return analyze_data(cl.data, t2, t3, quad, eta2x, eta2y, eta3x, eta3y); +int analyze_cluster(Cluster cl, int32_t *t2, int32_t *t3, char *quad, + double *eta2x, double *eta2y, double *eta3x, + double *eta3y) { + return analyze_data(cl.data, t2, t3, quad, eta2x, eta2y, eta3x, eta3y); } -int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, double *eta2x, double *eta2y, double *eta3x, double *eta3y) { +int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, + double *eta2x, double *eta2y, double *eta3x, double *eta3y) { - - int ok=1; + int ok = 1; int32_t tot2[4]; - int32_t t2max=0; + int32_t t2max = 0; char c; int32_t val, tot3; @@ -199,23 +200,22 @@ int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, double *et tot2[i] = 0; for (int ix = 0; ix < 3; ix++) { - for (int iy = 0; iy < 3; iy++) { - val = data[iy * 3 + ix]; - // printf ("%d ",data[iy * 3 + ix]); - tot3 += val; - if (ix <= 1 && iy <= 1) - tot2[cBottomLeft] += val; - if (ix >= 1 && iy <= 1) - tot2[cBottomRight] += val; - if (ix <= 1 && iy >= 1) - tot2[cTopLeft] += val; - if (ix >= 1 && iy >= 1) - tot2[cTopRight] += val; - } - // printf ("\n"); - + for (int iy = 0; iy < 3; iy++) { + val = data[iy * 3 + ix]; + // printf ("%d ",data[iy * 3 + ix]); + tot3 += val; + if (ix <= 1 && iy <= 1) + tot2[cBottomLeft] += val; + if (ix >= 1 && iy <= 1) + tot2[cBottomRight] += val; + if (ix <= 1 && iy >= 1) + tot2[cTopLeft] += val; + if (ix >= 1 && iy >= 1) + tot2[cTopRight] += val; + } + // printf ("\n"); } - //printf ("\n"); + // printf ("\n"); if (t2 || quad) { @@ -227,57 +227,57 @@ int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, double *et c = i; } } - - if (quad) - *quad = c; - if (t2) - *t2 = t2max; + + if (quad) + *quad = c; + if (t2) + *t2 = t2max; } if (t3) - *t3 = tot3; - + *t3 = tot3; + if (eta2x || eta2y) { - if (eta2x ) - *eta2x=0; - if (eta2y ) - *eta2y=0; - switch (c) { - case cBottomLeft: - if (eta2x && (data[3]+data[4])!=0) - *eta2x=(double)(data[4])/(data[3]+data[4]); - if (eta2y && (data[1]+data[4])!=0) - *eta2y=(double)(data[4])/(data[1]+data[4]); - break; - case cBottomRight: - if (eta2x && (data[2]+data[5])!=0) - *eta2x=(double)(data[5])/(data[4]+data[5]); - if (eta2y && (data[1]+data[4])!=0) - *eta2y=(double)(data[4])/(data[1]+data[4]); - break; - case cTopLeft: - if (eta2x && (data[7]+data[4])!=0) - *eta2x=(double)(data[4])/(data[3]+data[4]); - if (eta2y && (data[7]+data[4])!=0) - *eta2y=(double)(data[7])/(data[7]+data[4]); - break; - case cTopRight: - if (eta2x && t2max!=0) - *eta2x=(double)(data[5])/(data[5]+data[4]); - if (eta2y && t2max!=0) - *eta2y=(double)(data[7])/(data[7]+data[4]); - break; - default: - ; - } + if (eta2x) + *eta2x = 0; + if (eta2y) + *eta2y = 0; + switch (c) { + case cBottomLeft: + if (eta2x && (data[3] + data[4]) != 0) + *eta2x = (double)(data[4]) / (data[3] + data[4]); + if (eta2y && (data[1] + data[4]) != 0) + *eta2y = (double)(data[4]) / (data[1] + data[4]); + break; + case cBottomRight: + if (eta2x && (data[2] + data[5]) != 0) + *eta2x = (double)(data[5]) / (data[4] + data[5]); + if (eta2y && (data[1] + data[4]) != 0) + *eta2y = (double)(data[4]) / (data[1] + data[4]); + break; + case cTopLeft: + if (eta2x && (data[7] + data[4]) != 0) + *eta2x = (double)(data[4]) / (data[3] + data[4]); + if (eta2y && (data[7] + data[4]) != 0) + *eta2y = (double)(data[7]) / (data[7] + data[4]); + break; + case cTopRight: + if (eta2x && t2max != 0) + *eta2x = (double)(data[5]) / (data[5] + data[4]); + if (eta2y && t2max != 0) + *eta2y = (double)(data[7]) / (data[7] + data[4]); + break; + default:; + } } - + if (eta3x || eta3y) { - if (eta3x && (data[3]+data[4]+data[5])!=0) - *eta3x=(double)(-data[3]+data[3+2])/(data[3]+data[4]+data[5]); - if (eta3y && (data[1]+data[4]+data[7])!=0) - *eta3y=(double)(-data[1]+data[2*3+1])/(data[1]+data[4]+data[7]); + if (eta3x && (data[3] + data[4] + data[5]) != 0) + *eta3x = (double)(-data[3] + data[3 + 2]) / + (data[3] + data[4] + data[5]); + if (eta3y && (data[1] + data[4] + data[7]) != 0) + *eta3y = (double)(-data[1] + data[2 * 3 + 1]) / + (data[1] + data[4] + data[7]); } - return ok; } diff --git a/src/cluster_reader.h b/src/cluster_reader.h index 4892c59..32668d0 100644 --- a/src/cluster_reader.h +++ b/src/cluster_reader.h @@ -3,17 +3,20 @@ #include #include "data_types.h" -//Pure C implementation to read a cluster file +// Pure C implementation to read a cluster file -size_t read_clusters(FILE* fp, size_t n_clusters, Cluster* buf, uint32_t *n_left); +size_t read_clusters(FILE *fp, size_t n_clusters, Cluster *buf, + uint32_t *n_left); -size_t read_clusters_with_cut(FILE* fp, size_t n_clusters, Cluster* buf, uint32_t *n_left, double *noise_map, int nx, int ny); +size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf, + uint32_t *n_left, double *noise_map, int nx, + int ny); -int analyze_clusters(int64_t n_clusters, int32_t* cin, ClusterAnalysis *cout, int csize); +int analyze_clusters(int64_t n_clusters, int32_t *cin, ClusterAnalysis *cout, + int csize); +int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, + double *eta2x, double *eta2y, double *eta3x, double *eta3y); - - -int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, double *eta2x, double *eta2y, double *eta3x, double *eta3y); - -int analyze_cluster(Cluster data, int32_t *t2, int32_t *t3, char *quad, double *eta2x, double *eta2y, double *eta3x, double *eta3y); +int analyze_cluster(Cluster data, int32_t *t2, int32_t *t3, char *quad, + double *eta2x, double *eta2y, double *eta3x, double *eta3y); diff --git a/src/creader_module.c b/src/creader_module.c index 575b30a..3ffbefc 100644 --- a/src/creader_module.c +++ b/src/creader_module.c @@ -30,7 +30,8 @@ /\* (PyArrayObject *)cl_obj, cluster_dt(), NPY_ARRAY_C_CONTIGUOUS); *\/ /\* if (cl_array == NULL) { *\/ /\* PyErr_SetString(PyExc_TypeError, *\/ - /\* "Could not convert first argument to numpy array."); *\/ + /\* "Could not convert first argument to numpy array."); +*\/ /\* return NULL; *\/ /\* } *\/ @@ -56,92 +57,83 @@ */ - // clusterize method -//static PyObject *ClusterFileReader_clusterize(ClusterFileReader *self, PyObject *args) { +// static PyObject *ClusterFileReader_clusterize(ClusterFileReader *self, +// PyObject *args) { static PyObject *clusterize(PyObject *Py_UNUSED(self), PyObject *args) { const int ndim = 1; - + Py_ssize_t size = 0; PyObject *data_obj; if (!PyArg_ParseTuple(args, "nO", &size, &data_obj)) { - PyErr_SetString( - PyExc_TypeError, - "Could not parse args."); - return NULL; - } + PyErr_SetString(PyExc_TypeError, "Could not parse args."); + return NULL; + } - // + // // Create two numpy arrays from the passed objects, if possible numpy will // use the underlying buffer, otherwise it will create a copy, for example // if data type is different or we pass in a list. The // NPY_ARRAY_C_CONTIGUOUS flag ensures that we have contiguous memory. - PyObject *data_array = PyArray_FROM_OTF(data_obj, NPY_INT32, NPY_ARRAY_C_CONTIGUOUS); - int nx=0,ny=0; - int32_t *data=NULL; - + PyObject *data_array = + PyArray_FROM_OTF(data_obj, NPY_INT32, NPY_ARRAY_C_CONTIGUOUS); + int nx = 0, ny = 0; + int32_t *data = NULL; // If parsing of a or b fails we throw an exception in Python - if (data_array ) { - - int ndim_data = PyArray_NDIM((PyArrayObject *)(data_array)); - npy_intp *data_shape = PyArray_SHAPE((PyArrayObject *)(data_array)); - - - // For the C++ function call we need pointers (or another C++ type/data - // structure) - - data = (int32_t *)(PyArray_DATA((PyArrayObject *)(data_array))); + if (data_array) { + int ndim_data = PyArray_NDIM((PyArrayObject *)(data_array)); + npy_intp *data_shape = PyArray_SHAPE((PyArrayObject *)(data_array)); + // For the C++ function call we need pointers (or another C++ type/data + // structure) - /* for (int i=0; i< ndim_noise; i++) { */ - /* printf("Dimension %d size %d pointer \n",i,noise_shape[i], noise_map); */ - - /* } */ + data = (int32_t *)(PyArray_DATA((PyArrayObject *)(data_array))); - if (ndim_data==2) { - - nx=data_shape[0]; - ny=data_shape[1]; - if (ny!=9) { - PyErr_SetString( - PyExc_TypeError, - "Wrong data type."); - // printf("Data found size %d %d %d\n",nx,ny,ndim); - } + /* for (int i=0; i< ndim_noise; i++) { */ + /* printf("Dimension %d size %d pointer \n",i,noise_shape[i], + * noise_map); */ - } else { - PyErr_SetString( - PyExc_TypeError, - "Wrong data type."); + /* } */ - } - + if (ndim_data == 2) { + + nx = data_shape[0]; + ny = data_shape[1]; + if (ny != 9) { + PyErr_SetString(PyExc_TypeError, "Wrong data type."); + // printf("Data found size %d %d %d\n",nx,ny,ndim); + } + + } else { + PyErr_SetString(PyExc_TypeError, "Wrong data type."); + } } // Create an uninitialized numpy array - //npy_intp dims[] = {nx}; + // npy_intp dims[] = {nx}; // printf("%d %d\n",ndim,nx); npy_intp dims[] = {nx}; - PyObject *ca = PyArray_SimpleNewFromDescr(ndim, dims, cluster_analysis_dt()); + PyObject *ca = + PyArray_SimpleNewFromDescr(ndim, dims, cluster_analysis_dt()); - //printf("1\n"); + // printf("1\n"); // Fill with zeros PyArray_FILLWBYTE((PyArrayObject *)ca, 0); - //printf("2\n"); - // Get a pointer to the array memory + // printf("2\n"); + // Get a pointer to the array memory void *buf = PyArray_DATA((PyArrayObject *)ca); // Call the standalone C code to read clusters from file // Here goes the looping, removing frame numbers etc. - - //printf("3\n"); - int nc=analyze_clusters(nx,data,buf,size); - + + // printf("3\n"); + int nc = analyze_clusters(nx, data, buf, size); + // printf("aa %d %d\n",n_read, nx); /* if (nc != nx) { */ /* // resize the array to match the number of read photons */ @@ -159,16 +151,13 @@ static PyObject *clusterize(PyObject *Py_UNUSED(self), PyObject *args) { /* PyArray_Resize((PyArrayObject *)ca, &new_shape, 1, NPY_ANYORDER); */ /* } */ if (nc != nx) { - printf("%d %d\n",nx,nc); - PyErr_SetString(PyExc_TypeError, "Parsed wrong size array!"); + printf("%d %d\n", nx, nc); + PyErr_SetString(PyExc_TypeError, "Parsed wrong size array!"); } Py_DECREF(data_array); return ca; - } - - static PyObject *get_cluster_dt(PyObject *Py_UNUSED(self), PyObject *args) { if (!PyArg_ParseTuple(args, "")) return NULL; diff --git a/src/raw_reader.h b/src/raw_reader.h index e50e1a6..f3f96db 100644 --- a/src/raw_reader.h +++ b/src/raw_reader.h @@ -1,27 +1,16 @@ #pragma once +#include "data_types.h" #include #include #include -#include "data_types.h" -int64_t read_raw_m03( - FILE *fp, - int64_t n_frames, - char* frame_out, - Header* header_out -); +int64_t read_raw_m03(FILE *fp, int64_t n_frames, char *frame_out, + Header *header_out); -int64_t read_raw_m04( - FILE *fp, - int64_t n_frames, - char* frame_out, - char* digital_out, - Header* header_out -); +int64_t read_raw_m04(FILE *fp, int64_t n_frames, char *frame_out, + char *digital_out, Header *header_out); void decode_moench03(const uint16_t *buf, uint16_t *out_buf); -void decode_moench04(const uint16_t *analog_data, - const uint64_t *digital_data, - uint16_t *analog_frame, - uint8_t *digital_frame); \ No newline at end of file +void decode_moench04(const uint16_t *analog_data, const uint64_t *digital_data, + uint16_t *analog_frame, uint8_t *digital_frame); \ No newline at end of file From 437a6cf2d5f3897fb1e3eb0dc60c5098f5a6a2d9 Mon Sep 17 00:00:00 2001 From: froejdh_e Date: Fri, 27 Oct 2023 12:05:21 +0200 Subject: [PATCH 6/7] initialize c to 0 in cluster reader --- src/cluster_reader.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/cluster_reader.c b/src/cluster_reader.c index fc953ff..2166cff 100644 --- a/src/cluster_reader.c +++ b/src/cluster_reader.c @@ -192,7 +192,7 @@ int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, int32_t tot2[4]; int32_t t2max = 0; - char c; + char c = 0; int32_t val, tot3; tot3 = 0; From 2b8415ca9bf62488a61c4726cb698583d0e316d0 Mon Sep 17 00:00:00 2001 From: Erik Frojdh Date: Thu, 7 Dec 2023 08:59:19 +0100 Subject: [PATCH 7/7] added function to convert to sparse data --- creader/SparseFile.py | 14 ++++++++++++++ creader/__init__.py | 4 +++- 2 files changed, 17 insertions(+), 1 deletion(-) create mode 100644 creader/SparseFile.py diff --git a/creader/SparseFile.py b/creader/SparseFile.py new file mode 100644 index 0000000..4d0f9c8 --- /dev/null +++ b/creader/SparseFile.py @@ -0,0 +1,14 @@ +import numpy as np + +def to_sparse(data, th = 0): + """ + Convert frames stack ndarray[frame, row, col] to sparse array. + Warning: this function drops the frame numbers and only keeps row, col, energy + """ + sparse_dt = [('row', np.int16), ('col', np.int16), ('energy', data.dtype)] + size = (data>th).sum() + sparse = np.zeros(size, sparse_dt) + frames, rows, cols = np.where(data>th) + for i,(f,r,c) in enumerate(zip(frames, rows, cols)): + sparse[i] = (r, c, data[f,r,c]) + return sparse diff --git a/creader/__init__.py b/creader/__init__.py index 47f8a42..609c0d1 100644 --- a/creader/__init__.py +++ b/creader/__init__.py @@ -4,4 +4,6 @@ from _creader import * from .file_utils import open_file from .ClusterFile import ClusterFile from .enums import DetectorType -from .RawFile import RawFile \ No newline at end of file +from .RawFile import RawFile + +from .SparseFile import to_sparse \ No newline at end of file