##// END OF EJS Templates
schain-jc
RSS
Update NPY_ARRAY definition in C extension
Juan C. Espinoza -
r1169:668c33de071f
parent child
Show More
Side by Side Unified
Context file:
r1169:668c33de071f Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,139 +1,139
1 #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
1 #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
2 #define NUM_CPY_THREADS 8
2 #define NUM_CPY_THREADS 8
3 #include <Python.h>
3 #include <Python.h>
4 #include <numpy/arrayobject.h>
4 #include <numpy/arrayobject.h>
5 #include <math.h>
5 #include <math.h>
6 #include <complex.h>
6 #include <complex.h>
7 #include <time.h>
7 #include <time.h>
8
8
9 // void printArr(int *array);
9 // void printArr(int *array);
10 static PyObject *hildebrand_sekhon(PyObject *self, PyObject *args);
10 static PyObject *hildebrand_sekhon(PyObject *self, PyObject *args);
11 static PyObject *correlateByBlock(PyObject *self, PyObject *args);
11 static PyObject *correlateByBlock(PyObject *self, PyObject *args);
12 #ifndef PyMODINIT_FUNC /* declarations for DLL import/export */
12 #ifndef PyMODINIT_FUNC /* declarations for DLL import/export */
13 #define PyMODINIT_FUNC void
13 #define PyMODINIT_FUNC void
14 #endif
14 #endif
15
15
16 static PyMethodDef extensionsMethods[] = {
16 static PyMethodDef extensionsMethods[] = {
17 { "correlateByBlock", (PyCFunction)correlateByBlock, METH_VARARGS, "get correlation by block" },
17 { "correlateByBlock", (PyCFunction)correlateByBlock, METH_VARARGS, "get correlation by block" },
18 { "hildebrand_sekhon", (PyCFunction)hildebrand_sekhon, METH_VARARGS, "get noise with hildebrand_sekhon" },
18 { "hildebrand_sekhon", (PyCFunction)hildebrand_sekhon, METH_VARARGS, "get noise with hildebrand_sekhon" },
19 { NULL, NULL, 0, NULL }
19 { NULL, NULL, 0, NULL }
20 };
20 };
21
21
22 PyMODINIT_FUNC initcSchain() {
22 PyMODINIT_FUNC initcSchain() {
23 Py_InitModule("cSchain", extensionsMethods);
23 Py_InitModule("cSchain", extensionsMethods);
24 import_array();
24 import_array();
25 }
25 }
26
26
27 static PyObject *correlateByBlock(PyObject *self, PyObject *args) {
27 static PyObject *correlateByBlock(PyObject *self, PyObject *args) {
28
28
29 // int *x = (int*) malloc(4000000 * 216 * sizeof(int));;
29 // int *x = (int*) malloc(4000000 * 216 * sizeof(int));;
30 // int a = 5;
30 // int a = 5;
31 // x = &a;
31 // x = &a;
32 // int b = 6;
32 // int b = 6;
33 // x = &b;
33 // x = &b;
34 // printf("Antes de imprimir x \n");
34 // printf("Antes de imprimir x \n");
35 // printf("%d \n", x[0]);
35 // printf("%d \n", x[0]);
36
36
37 PyObject *data_obj1, *data_obj2;
37 PyObject *data_obj1, *data_obj2;
38 PyArrayObject *data_array1, *data_array2, *correlateRow, *out, *dataRow, *codeRow; //, ,
38 PyArrayObject *data_array1, *data_array2, *correlateRow, *out, *dataRow, *codeRow; //, ,
39 int mode;
39 int mode;
40
40
41 if (!PyArg_ParseTuple(args, "OOi", &data_obj1, &data_obj2, &mode)) return NULL;
41 if (!PyArg_ParseTuple(args, "OOi", &data_obj1, &data_obj2, &mode)) return NULL;
42
42
43 data_array1 = (PyArrayObject *) PyArray_FROM_OTF(data_obj1, NPY_COMPLEX128, NPY_ARRAY_DEFAULT);
43 data_array1 = (PyArrayObject *) PyArray_FROM_OTF(data_obj1, NPY_COMPLEX128, NPY_ARRAY_IN_ARRAY);
44 data_array2 = (PyArrayObject *) PyArray_FROM_OTF(data_obj2, NPY_FLOAT64, NPY_ARRAY_DEFAULT);
44 data_array2 = (PyArrayObject *) PyArray_FROM_OTF(data_obj2, NPY_FLOAT64, NPY_ARRAY_IN_ARRAY);
45
45
46 npy_intp dims[1];
46 npy_intp dims[1];
47 dims[0] = 200;
47 dims[0] = 200;
48 npy_intp dims_code[1];
48 npy_intp dims_code[1];
49 dims_code[0] = 16;
49 dims_code[0] = 16;
50
50
51 double complex * dataRaw;
51 double complex * dataRaw;
52 double * codeRaw;
52 double * codeRaw;
53 dataRaw = (double complex*) PyArray_DATA(data_array1);
53 dataRaw = (double complex*) PyArray_DATA(data_array1);
54 codeRaw = (double *) PyArray_DATA(data_array2);
54 codeRaw = (double *) PyArray_DATA(data_array2);
55 double complex ** outC = malloc(40000*200*sizeof(double complex));
55 double complex ** outC = malloc(40000*200*sizeof(double complex));
56 int i;
56 int i;
57
57
58 clock_t start = clock();
58 clock_t start = clock();
59 for(i=0; i<40000; i++){
59 for(i=0; i<40000; i++){
60 // codeRow = PyArray_SimpleNewFromData(1, dims_code, NPY_FLOAT64, codeRaw + 16 * i);
60 // codeRow = PyArray_SimpleNewFromData(1, dims_code, NPY_FLOAT64, codeRaw + 16 * i);
61 // dataRow = PyArray_SimpleNewFromData(1, dims, NPY_COMPLEX128, dataRaw + 200 * i);
61 // dataRow = PyArray_SimpleNewFromData(1, dims, NPY_COMPLEX128, dataRaw + 200 * i);
62 // Py_INCREF(codeRow);
62 // Py_INCREF(codeRow);
63 // Py_INCREF(dataRow);
63 // Py_INCREF(dataRow);
64 // PyArray_ENABLEFLAGS(codeRow, NPY_ARRAY_OWNDATA);
64 // PyArray_ENABLEFLAGS(codeRow, NPY_ARRAY_OWNDATA);
65 // PyArray_ENABLEFLAGS(dataRow, NPY_ARRAY_OWNDATA);
65 // PyArray_ENABLEFLAGS(dataRow, NPY_ARRAY_OWNDATA);
66 correlateRow = (PyArrayObject *) PyArray_Correlate2(PyArray_SimpleNewFromData(1, dims_code, NPY_FLOAT64, codeRaw + 16 * i), PyArray_SimpleNewFromData(1, dims, NPY_COMPLEX128, dataRaw + 200 * i), (npy_intp) 2);
66 correlateRow = (PyArrayObject *) PyArray_Correlate2(PyArray_SimpleNewFromData(1, dims_code, NPY_FLOAT64, codeRaw + 16 * i), PyArray_SimpleNewFromData(1, dims, NPY_COMPLEX128, dataRaw + 200 * i), (npy_intp) 2);
67 //Py_INCREF(correlateRow);
67 //Py_INCREF(correlateRow);
68 // PyArray_ENABLEFLAGS(correlateRow, NPY_ARRAY_OWNDATA);
68 // PyArray_ENABLEFLAGS(correlateRow, NPY_ARRAY_OWNDATA);
69 memcpy(outC + 200*i, (double complex*) PyArray_DATA(correlateRow), 200 * sizeof(double complex));
69 memcpy(outC + 200*i, (double complex*) PyArray_DATA(correlateRow), 200 * sizeof(double complex));
70
70
71 Py_DECREF(correlateRow);
71 Py_DECREF(correlateRow);
72 // Py_DECREF(codeRow);
72 // Py_DECREF(codeRow);
73 // Py_DECREF(dataRow);
73 // Py_DECREF(dataRow);
74 }
74 }
75 clock_t end = clock();
75 clock_t end = clock();
76 float seconds = (float)(end - start) / CLOCKS_PER_SEC;
76 float seconds = (float)(end - start) / CLOCKS_PER_SEC;
77 printf("%f", seconds);
77 printf("%f", seconds);
78 //
78 //
79 npy_intp dimsret[2];
79 npy_intp dimsret[2];
80 dimsret[0] = 40000;
80 dimsret[0] = 40000;
81 dimsret[1] = 200;
81 dimsret[1] = 200;
82 out = PyArray_SimpleNewFromData(2, dimsret, NPY_COMPLEX128, outC);
82 out = PyArray_SimpleNewFromData(2, dimsret, NPY_COMPLEX128, outC);
83 PyArray_ENABLEFLAGS(out, NPY_ARRAY_OWNDATA);
83 PyArray_ENABLEFLAGS(out, NPY_ARRAY_OWNDATA);
84 //Py_INCREF(out);
84 //Py_INCREF(out);
85 Py_DECREF(data_array1);
85 Py_DECREF(data_array1);
86 Py_DECREF(data_array2);
86 Py_DECREF(data_array2);
87 // PyArray_DebugPrint(out);
87 // PyArray_DebugPrint(out);
88 // Py_DECREF(data_obj2);
88 // Py_DECREF(data_obj2);
89 // Py_DECREF(data_obj1);
89 // Py_DECREF(data_obj1);
90 // Py_DECREF(codeRow);
90 // Py_DECREF(codeRow);
91 // Py_DECREF(dataRow);
91 // Py_DECREF(dataRow);
92 // free(dataRaw);
92 // free(dataRaw);
93 // free(codeRaw);
93 // free(codeRaw);
94
94
95 return PyArray_Return(out);
95 return PyArray_Return(out);
96 }
96 }
97
97
98 static PyObject *hildebrand_sekhon(PyObject *self, PyObject *args) {
98 static PyObject *hildebrand_sekhon(PyObject *self, PyObject *args) {
99 double navg;
99 double navg;
100 PyObject *data_obj, *data_array;
100 PyObject *data_obj, *data_array;
101
101
102 if (!PyArg_ParseTuple(args, "Od", &data_obj, &navg)) return NULL;
102 if (!PyArg_ParseTuple(args, "Od", &data_obj, &navg)) return NULL;
103 data_array = PyArray_FROM_OTF(data_obj, NPY_FLOAT64, NPY_ARRAY_DEFAULT);
103 data_array = PyArray_FROM_OTF(data_obj, NPY_FLOAT64, NPY_ARRAY_IN_ARRAY);
104 if (data_array == NULL) {
104 if (data_array == NULL) {
105 Py_XDECREF(data_array);
105 Py_XDECREF(data_array);
106 Py_XDECREF(data_obj);
106 Py_XDECREF(data_obj);
107 return NULL;
107 return NULL;
108 }
108 }
109 double *sortdata = (double*)PyArray_DATA(data_array);
109 double *sortdata = (double*)PyArray_DATA(data_array);
110 int lenOfData = (int)PyArray_SIZE(data_array) ;
110 int lenOfData = (int)PyArray_SIZE(data_array) ;
111 double nums_min = lenOfData*0.2;
111 double nums_min = lenOfData*0.2;
112 if (nums_min <= 5) nums_min = 5;
112 if (nums_min <= 5) nums_min = 5;
113 double sump = 0;
113 double sump = 0;
114 double sumq = 0;
114 double sumq = 0;
115 int j = 0;
115 int j = 0;
116 int cont = 1;
116 int cont = 1;
117 double rtest = 0;
117 double rtest = 0;
118 while ((cont == 1) && (j < lenOfData)) {
118 while ((cont == 1) && (j < lenOfData)) {
119 sump = sump + sortdata[j];
119 sump = sump + sortdata[j];
120 sumq = sumq + pow(sortdata[j], 2);
120 sumq = sumq + pow(sortdata[j], 2);
121 if (j > nums_min) {
121 if (j > nums_min) {
122 rtest = (double)j/(j-1) + 1/navg;
122 rtest = (double)j/(j-1) + 1/navg;
123 if ((sumq*j) > (rtest*pow(sump, 2))) {
123 if ((sumq*j) > (rtest*pow(sump, 2))) {
124 j = j - 1;
124 j = j - 1;
125 sump = sump - sortdata[j];
125 sump = sump - sortdata[j];
126 sumq = sumq - pow(sortdata[j],2);
126 sumq = sumq - pow(sortdata[j],2);
127 cont = 0;
127 cont = 0;
128 }
128 }
129 }
129 }
130 j = j + 1;
130 j = j + 1;
131 }
131 }
132
132
133 double lnoise = sump / j;
133 double lnoise = sump / j;
134
134
135 Py_DECREF(data_array);
135 Py_DECREF(data_array);
136
136
137 return Py_BuildValue("d", lnoise);
137 return Py_BuildValue("d", lnoise);
138 }
138 }
139
139
General Comments 0
You need to be logged in to leave comments. Login now