##// END OF EJS Templates
fix line 461 y 665
fix line 461 y 665

File last commit:

r1543:a3a1a4029bf5
r1791:e4984c469c9f isr
Show More
_HS_algorithm.c
86 lines | 1.9 KiB | text/x-c | CLexer
#include <Python.h>
#include <numpy/arrayobject.h>
#include <math.h>
static PyObject *HS_algorithm(PyObject *self, PyObject *args) {
double navg;
PyObject *data_obj, *data_array;
if (!PyArg_ParseTuple(args, "Od", &data_obj, &navg)) {
return NULL;
}
data_array = PyArray_FROM_OTF(data_obj, NPY_FLOAT64, NPY_IN_ARRAY);
if (data_array == NULL) {
Py_XDECREF(data_array);
Py_XDECREF(data_obj);
return NULL;
}
double *sortdata = (double*)PyArray_DATA(data_array);
int lenOfData = (int)PyArray_SIZE(data_array) ;
double nums_min = lenOfData*0.75;
if (nums_min <= 5) nums_min = 5;
double sump = 0;
double sumq = 0;
long j = 0;
int cont = 1;
double rtest = 0;
while ((cont == 1) && (j < lenOfData)) {
sump = sump + sortdata[j];
sumq = sumq + pow(sortdata[j], 2);
if (j > nums_min) {
rtest = (double)j/(j-1) + 1/navg;
//printf("%ld\n", j);
//printf("%f \n", rtest);
//printf("%f \n", sump);
if ((sumq*j) > (rtest*pow(sump, 2))) {
j = j - 1;
sump = sump - sortdata[j];
sumq = sumq - pow(sortdata[j],2);
cont = 0;
}
}
//printf("%ld\n", j);
j = j + 1;
}
double lnoise = sump / j;
Py_DECREF(data_array);
// return PyLong_FromLong(lnoise);
return PyFloat_FromDouble(j);
}
static PyMethodDef noiseMethods[] = {
{ "HS_algorithm", HS_algorithm, METH_VARARGS, "Applies hildebrand_sekhon algorithm" },
{ NULL, NULL, 0, NULL }
};
#if PY_MAJOR_VERSION >= 3
static struct PyModuleDef noisemodule = {
PyModuleDef_HEAD_INIT,
"_HS_algorithm",
"Applies hildebrand_sekhon algorithm",
-1,
noiseMethods
};
#endif
#if PY_MAJOR_VERSION >= 3
PyMODINIT_FUNC PyInit__HS_algorithm(void) {
Py_Initialize();
import_array();
return PyModule_Create(&noisemodule);
}
#else
PyMODINIT_FUNC init_HS_algorithm() {
Py_InitModule("_HS_algorithm", noiseMethods);
import_array();
}
#endif