| @@ -1,55 +1,58 | |||
|
|
1 | 1 | #include <Python.h> |
|
|
2 | 2 | #include <numpy/arrayobject.h> |
|
|
3 | 3 | #include <math.h> |
|
|
4 | 4 | |
|
|
5 | 5 | static PyObject *hildebrand_sekhon(PyObject *self, PyObject *args); |
|
|
6 | 6 | |
|
|
7 | 7 | static PyMethodDef extensionsMethods[] = { |
|
|
8 | 8 | { "hildebrand_sekhon", (PyCFunction)hildebrand_sekhon, METH_VARARGS, "get noise with" }, |
|
|
9 | 9 | { NULL, NULL, 0, NULL } |
|
|
10 | 10 | }; |
|
|
11 | 11 | |
|
|
12 | 12 | PyMODINIT_FUNC initcSchain() { |
|
|
13 | 13 | Py_InitModule("cSchain", extensionsMethods); |
|
|
14 | 14 | import_array(); |
|
|
15 | 15 | } |
|
|
16 | 16 | |
|
|
17 | 17 | static PyObject *hildebrand_sekhon(PyObject *self, PyObject *args) { |
|
|
18 | 18 | /* Do your stuff here. */ |
|
|
19 | 19 | double navg; |
|
|
20 | 20 | PyObject *data_obj, *data_array; |
|
|
21 | 21 | |
|
|
22 | 22 | if (!PyArg_ParseTuple(args, "Od", &data_obj, &navg)) return NULL; |
|
|
23 | 23 | data_array = PyArray_FROM_OTF(data_obj, NPY_FLOAT64, NPY_IN_ARRAY); |
|
|
24 | 24 | if (data_array == NULL) { |
|
|
25 | 25 | Py_XDECREF(data_array); |
|
|
26 | 26 | Py_XDECREF(data_obj); |
|
|
27 | 27 | return NULL; |
|
|
28 | 28 | } |
|
|
29 | 29 | double *sortdata = (double*)PyArray_DATA(data_array); |
|
|
30 | 30 | int lenOfData = (int)PyArray_SIZE(data_array) ; |
|
|
31 | 31 | double nums_min = lenOfData*0.2; |
|
|
32 | 32 | if (nums_min <= 5) nums_min = 5; |
|
|
33 | 33 | double sump = 0; |
|
|
34 | 34 | double sumq = 0; |
|
|
35 | 35 | int j = 0; |
|
|
36 | 36 | int cont = 1; |
|
|
37 | 37 | double rtest = 0; |
|
|
38 | 38 | while ((cont == 1) && (j < lenOfData)) { |
|
|
39 | 39 | sump = sump + sortdata[j]; |
|
|
40 | 40 | sumq = sumq + pow(sortdata[j], 2); |
|
|
41 | 41 | if (j > nums_min) { |
|
|
42 | 42 | rtest = (double)j/(j-1) + 1/navg; |
|
|
43 | 43 | if ((sumq*j) > (rtest*pow(sump, 2))) { |
|
|
44 | 44 | j = j - 1; |
|
|
45 | 45 | sump = sump - sortdata[j]; |
|
|
46 | 46 | sumq = sumq - pow(sortdata[j],2); |
|
|
47 | 47 | cont = 0; |
|
|
48 | 48 | } |
|
|
49 | 49 | } |
|
|
50 | 50 | j = j + 1; |
|
|
51 | 51 | } |
|
|
52 | 52 | |
|
|
53 | 53 | double lnoise = sump / j; |
|
|
54 | ||
|
|
55 | Py_DECREF(data_array); | |
|
|
56 | ||
|
|
54 | 57 | return Py_BuildValue("d", lnoise); |
|
|
55 | 58 | } |
General Comments 0
You need to be logged in to leave comments.
Login now