cluster reader seems to work

setup.py

@@ -1,14 +1,20 @@
-import setuptools
-import numpy
 
-c_ext = setuptools.Extension("cluster_reader_cpp",
+
-                             sources=["src/cluster_reader_module.cpp", "src/cluster_reader.cpp"],
+import setuptools
-                             extra_compile_args=['-std=c++11', '-Wall'],
+import numpy as np
-                             include_dirs=[numpy.get_include()])
+
+c_ext = setuptools.Extension("creader",
+                    sources = ["src/creader_module.c", "src/cluster_reader.c"],
+                    include_dirs=[
+                            np.get_include(),"src/"
+                            ],
+                    extra_compile_args=['-std=c++11', '-Wall', '-Wextra'] )
+                    
+
 c_ext.language = 'c++'
 setuptools.setup(
-    name='cluster_reader_cpp',
+    name= 'creader',
-    version='0.1',
+    version = '0.1',
-    description='Tools to read cluster files.',
+    description = 'Reading cluster files',
     ext_modules=[c_ext],
 )

src/Makefile

@@ -0,0 +1,2 @@
+test: main.cpp cluster_reader.c cluster_reader.h
+	g++ -o test main.cpp cluster_reader.c -std=c++11

src/Makefile.test

@@ -1,5 +0,0 @@
-
-
-
-test: main.cpp cluster_reader.cpp cluster_reader.h
-      g++ -o test main.cpp cluster_reader.cpp -std=c++11

src/cluster_reader.c

@@ -0,0 +1,45 @@
+#include "cluster_reader.h"
+
+int read_clusters(FILE* fp, int64_t n_clusters, Cluster* buf, int *n_left){
+  printf("Item size: %d n_clusters: %d, n_left: %d\n", sizeof(Cluster), n_clusters,*n_left);
+    int iframe=0, nph=*n_left;
+    size_t n_read=0, nph_read=0, nn=*n_left, nr=0;
+    //n_left=n_clusters;
+    
+    //read photons left from previous frame
+    if (nph) {
+	if (nph>n_clusters-nph_read)
+	  nn=n_clusters-nph_read;
+	else 
+	  nn=nph;
+	//printf("* %d %d %d %d\n",iframe,nph,nn,n_left);
+      	nr += fread((void*)(buf+nph_read), sizeof(Cluster), nn, fp);
+	n_read+=nr/sizeof(Cluster);
+	nph_read+=nn;
+	*n_left=nph-nn;
+    } 
+    if (nph_read<n_clusters) {
+    //keep on reading frames and photons until reaching n_clusters
+    while (fread(&iframe, sizeof(iframe), 1, fp)) {
+      if (fread(&nph, sizeof(nph), 1, fp) ) {
+	if (nph>n_clusters-nph_read)
+	  nn=n_clusters-nph_read;
+	else 
+	  nn=nph;
+	
+	//printf("%d %d %d %d\n",iframe,nph,nr,n_left);
+	nr += fread((void*)(buf+nph_read), sizeof(Cluster), nn, fp);
+	//printf("%d %d %d %d\n",iframe,nph,nr,n_left);
+	n_read+=nr;
+	nph_read+=nn;
+	*n_left=nph-nn;
+      }
+      if (nph_read>=n_clusters)
+	break;
+    }
+    }
+    // size_t n_read = fread(buf, sizeof(Cluster), n_clusters, fp);
+    printf("Read: %d items %d left %d\n", nph_read, n_read,*n_left);
+    return nph_read;
+}
+  

src/cluster_reader.cpp

@@ -1,67 +0,0 @@
-
-#include "cluster_reader.h"
-#include <iostream>
-#include <vector>
-
-FILE *f;
-int ph_left;
-
-void cpp_open_cluster_file(char *fname, int &ok){
-  f=fopen(fname,"r");
-  ph_left=0;
-  if (f)
-    ok=1;
-  else
-    ok=0;
-}
-
-void cpp_close_cluster_file(int &ok){
-  if (f){
-    fclose(f);
-    ph_left=0;
-    f=NULL;
-    ok=1;
-  } else
-    ok=0;
-  
-}
-void  cpp_read_clusters(struct cluster *clust, int &iiph, int &nframes, int &ok, int maxph){
-  iiph=0;
-  int iframe, iff, nph;
-
-  ph_left=0;
-  iff=0;
-  ok=0;
-  iiph=0;
-  int nn;
-  if (f) {
-    ok=1;
-    while (fread((void*)&iframe, 1,  sizeof(int), f)) {
-      if (fread((void*)&nph, 1,  sizeof(int), f)) {
-	ph_left=nph;
-	if (nph>0) {
-	  //  for (int iph=0; iph<nph; iph++){ 
-	  std::cout << iframe << " " << nph << std::endl;
-	    if ((maxph-iiph)>nph) nn=nph;
-	    else nn=(maxph-iiph);
-	      
-	    if ( fread( (void*)(clust+iiph), sizeof(struct cluster), nn, f) ) { 
-	      ph_left=nph-nn;
-	    // if ( fread( (void*)&x, sizeof(int16_t), 1, myFile) ) {     //reads x
-	    //   if ( fread( (void*)&y, sizeof(int16_t), 1, myFile ) )    //reads y
-	    // 	return fread( (void*)data, sizeof(int), dx*dy, myFile );
-
-	   
-	    
-	      iiph+=nn;
-	      nframes=iff;
-	      if (iiph>=maxph)
-		break;
-	    }
-	}
-      iff++;
-      }
-    }
-    //return v;
-  }
-}

src/cluster_reader.h

@@ -1,22 +1,8 @@
-#ifndef CLUSTER_READER_H
+#pragma once
-#define CLUSTER_READER_H
 #include <stdint.h>
-#include <cstdio>
+#include <stdio.h>
 
-#define MAXPH 100000000
+#include "data_types.h"
+//Pure C implementation to read a cluster file
 
-struct cluster {
+int read_clusters(FILE* fp, int64_t n_clusters, Cluster* buf, int *n_left);
-  int16_t x;
-  int16_t y;
-  //data[iy * 3 + ix ]
-  int32_t data[9]; 
-} ;
-
-void  cpp_read_clusters(struct cluster *clust, int &iiph, int &nframes, int &ok, int maxph=MAXPH);
-
-
-void cpp_open_cluster_file(char *fname, int &ok);
-
-void cpp_close_cluster_file(int &ok);
-
-#endif

src/cluster_reader_module.cpp

@@ -1,147 +0,0 @@
-#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
-#include "cluster_reader.h" //this is the function that we want to call
-#include <Python.h>
-#include <numpy/arrayobject.h>
-#include<iostream>
-using namespace std;
-
-
-// Docstring
-static char module_docstring[] = "Cluster reader";
-
-
-PyDoc_STRVAR(read_clusters_doc, "Read clusters from open file.\n\n"
-                      "Parameters\n"
-                      "----------\n"
-                      "f: file name\n"
-
-                      "Returns\n"
-                      "----------\n"
-                      "nothing\n\n");
-
-
-
-
-PyDoc_STRVAR(open_cluster_file_doc, "Read cluster file.\n\n"
-                      "Parameters\n"
-                      "----------\n"
-                      "maxframes: maximum number of frames to be read\n"
-
-                      "Returns\n"
-                      "----------\n"
-                      "v: numpy_array\n"
-                      "   clusters\n\n");
-
-
-PyDoc_STRVAR(close_cluster_file_doc, "Read cluster file.\n\n"
-                      "Parameters\n"
-                      "----------\n"
-                      "none\n"
-                      "Returns\n"
-                      "----------\n"
-                      "nothing\n\n");
-
-
-// Declare functions that should go into the module
-static PyObject *read_clusters(PyObject *self, PyObject *args);
-static PyObject *open_cluster_file(PyObject *self, PyObject *args);
-static PyObject *close_cluster_file(PyObject *self, PyObject *args);
-
-// This is the module itself
-static PyMethodDef module_methods[] = {{"open_cluster_file", (PyCFunction)open_cluster_file, METH_VARARGS, open_cluster_file_doc},{"close_cluster_file", (PyCFunction)close_cluster_file, METH_VARARGS, close_cluster_file_doc},{"read_clusters", (PyCFunction)read_clusters, METH_VARARGS, read_clusters_doc},
-    {NULL, NULL, 0, NULL}};
-
-// Don't touch except the name
-static struct PyModuleDef cluster_reader_module_def = {PyModuleDef_HEAD_INIT, "cluster_reader_module",
-                                       module_docstring, -1, module_methods};
-
-// Don't touch
-PyMODINIT_FUNC PyInit_cluster_reader_cpp(void) {
-    PyObject *m = PyModule_Create(&cluster_reader_module_def);
-    if (m == NULL)
-        return NULL;
-
-    import_array();
-    return m;
-}
-
-static PyObject *open_cluster_file(PyObject *self, PyObject *args) {
-
-//    cout << "Hello there" << endl; // General Kenobi
-    // Python offers a lot of flexibility so the first thing we need to do
-    // is to parse the arguments and make sure we are called with correct
-    // parameters
-  char *fname="";
-    if (!PyArg_ParseTuple(args, "s", fname))
-        return NULL;
-
-    // Now call add wih pointers and size
-    int ok;
-    cpp_open_cluster_file(fname,ok);
-
-    return Py_BuildValue("i",ok);
-}
-
-static PyObject *close_cluster_file(PyObject *self, PyObject *args) {
-
-//    cout << "Hello there" << endl; // General Kenobi
-    // Python offers a lot of flexibility so the first thing we need to do
-    // is to parse the arguments and make sure we are called with correct
-    // parameters
- 
-    // Now call add wih pointers and size
-    int ok;
-    cpp_close_cluster_file(ok);
-
-    return Py_BuildValue("i",ok);
-}
-
-
-static PyObject *read_clusters(PyObject *self, PyObject *args) {
-
-//    cout << "Hello there" << endl; // General Kenobi
-    // Python offers a lot of flexibility so the first thing we need to do
-    // is to parse the arguments and make sure we are called with correct
-    // parameters
-    PyObject *clust_obj;
-    //PyObject *maxframes_obj;
-    if (!PyArg_ParseTuple(args, "O", &clust_obj))
-        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 *clust_array =
-        PyArray_FROM_OTF(clust_obj, NPY_OBJECT, NPY_ARRAY_C_CONTIGUOUS);
-    // npy_intp ind=0;
-    // void* clust_ptr = PyArray_GetPtr((PyArrayObject*)clust_obj, &ind);
-    
-      // If parsing of a or b fails we throw an exception in Python
-    if (clust_array == NULL) {
-        PyErr_SetString(
-            PyExc_TypeError,
-            "Could not convert one of the arguments to a numpy array.");
-        return NULL;
-    }
-    int N=MAXPH;
-    //  int N = (int)PyArray_DIM(clust_array, 0);
-    //npy_intp dims[2] = { 400, 400 };
-
-    // Create array for return values
-  
-    // For the C++ function call we need pointers (or another C++ type/data
-    // structure)
-    struct cluster *clust = reinterpret_cast<struct cluster *>(
-        PyArray_DATA(reinterpret_cast<PyArrayObject *>(clust_array)));
-
-    int nph, nframes,ok;
-    // Now call add wih pointers and size
-    cpp_read_clusters(clust,nph,nframes,ok,N);
-
-    // Clean up
-    Py_DECREF(clust_array);
-
-//    return result;
-    return Py_BuildValue("iii",nph,nframes,ok);
-}

src/data_types.h

@@ -0,0 +1,12 @@
+#ifndef CLUSTER_H
+#define CLUSTER_H
+#include <stdint.h>
+// Data type for the cluster
+// has to be extended
+typedef struct {
+  int16_t x;
+  int16_t y;
+  int32_t data[9];
+} Cluster ;
+
+#endif

src/main.cpp

@@ -1,22 +1,34 @@
+#include "data_types.h"
 #include "cluster_reader.h"
 #include <iostream>
 
+#define MAXPH 100000000
 using namespace std;
 
 int main(int argc, char* argv[]) {
 
   int nph, nframes;
   int ok;
-  struct cluster *clust=new struct cluster[MAXPH];
+  Cluster *clust=new Cluster[MAXPH];
-  cout << "Clusters allocated " << endl;
+  //Cluster c;
-  cpp_open_cluster_file("/mnt/moench_data/Moench_LGAD_SIM_Nov22/moenchLGAD202211/clustW17new/beam_En800eV_-40deg_300V_10us_d0_f5_0.clust",ok);
+  cout << "Clusters allocated " << sizeof(clust) << endl;
+  FILE *fp=fopen("/mnt/moench_data/Moench_LGAD_SIM_Nov22/moenchLGAD202211/clustW17new/beam_En800eV_-40deg_300V_10us_d0_f5_0.clust","rb");
+  //cpp_open_cluster_file("/mnt/moench_data/Moench_LGAD_SIM_Nov22/moenchLGAD202211/clustW17new/beam_En800eV_-40deg_300V_10us_d0_f5_0.clust",ok);
   cout << "File opened " << endl;
-  cpp_read_clusters(clust, nph,nframes,ok);
+  int nr,ntot;
+  int nleft=0;
+  while (nr =read_clusters(fp,MAXPH, clust+ntot,&nleft)){
+    ntot+=nr;
+    printf("** %d ",nr);
+  }
+  // cpp_read_clusters(clust, nph,nframes,ok);
   cout << "Clusters read " << endl;
-  cpp_close_cluster_file(ok);
+  //cpp_close_cluster_file(ok);
+  fclose(fp);
   cout << "File closed " << endl;
   int iph=0; 
-  // for (int iph=0; iph<nph; iph++){
+  printf("[");
+  for (int iph=0; iph<ntot; iph+=100000){
   //   int tot=0;
   //   for (int iy=0; iy<3; iy++) {
   //     for (int ix=0; ix<3; ix++) 
@@ -27,9 +39,23 @@ int main(int argc, char* argv[]) {
     
   //   // if (tot>500)
   //   //   std::cout << iph << "\t" << (clust+iph)->x << "\t" << (clust+iph)->y << "\t" << tot << endl;
+    printf("(");
+    printf("%3d, ",(clust+iph)->x);
+    printf("%3d, ",(clust+iph)->y);
  
-  // }
+    printf("[");
-  cout << nph << "\t"<< nframes<<endl;
+    // printf("\n \t");
+    for (int iy=0; iy<9; iy++) {
+	printf("%3d",(clust+iph)->data[iy]);
+	if (iy<8)printf(", "); 
+      //printf("\n \t");
+    }
+    printf("]");
+    printf(")");
+    printf("\n");
+  }
+  printf("]");
+  cout << ntot <<endl;
 
   return 0;
 }

test.py

@@ -0,0 +1,8 @@
+from creader import ClusterFileReader
+import numpy as np
+
+maxph=100000000
+
+r = ClusterFileReader("/mnt/moench_data/Moench_LGAD_SIM_Nov22/moenchLGAD202211/clustW17new/beam_En800eV_-40deg_300V_10us_d0_f5_0.clust")
+a=r.read(maxph)
+print(a[::100000])

tester.py

@@ -1,52 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-import pytest
-import numpy as np
-import mymod
-
-
-def test_add_two_int():
-    a = 5
-    b = 7
-    #Check for types to make sure that the test is ok
-    assert type(a) is int
-    assert type(b) is int
-
-    #This is actually what we want to test
-    assert mymod.add(a,b) == 12
-
-def test_add_int_float():
-    a = 723
-    b = 4.98
-    assert type(a) is int
-    assert type(b) is float
-    assert mymod.add(a,b) == pytest.approx(727.98)
-
-def test_add_two_float():
-    a = 309.123
-    b = 4.98
-    assert type(a) is float
-    assert type(b) is float
-    assert mymod.add(a,b) == pytest.approx(314.103)
-
-def test_add_two_arrays():
-    a = np.array((3,4,5))
-    b = np.array((9,1,3))
-    c = mymod.add(a,b) 
-    for aa,bb,cc in zip(a,b,c):
-        assert cc == aa+bb
-
-def test_throws_on_string():
-    with pytest.raises(TypeError):
-        mymod.add('a', 7)
-
-def test_throws_when_dims_disagree():
-    with pytest.raises(ValueError):
-        mymod.add(3, (4,5))
-
-def test_throws_when_dims_disagree():
-    a = np.arange(10).reshape(2,5)
-    b = np.arange(10).reshape(5,2)
-    with pytest.raises(ValueError):
-        mymod.add(a, b)