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