source: sasview/pr_inversion/c_extensions/Cinvertor.c @ 43c0a8e

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 43c0a8e was 43c0a8e, checked in by Mathieu Doucet <doucetm@…>, 16 years ago

Add a couple of figures of merit
  • Property mode set to 100644
File size: 18.4 KB
Line 
1
2#include <Python.h>
3#include "structmember.h"
4#include <stdio.h>
5#include <stdlib.h>
6#include <math.h>
7#include <time.h>
8
9#include "invertor.h"
10
11
12/// Error object for raised exceptions
13static PyObject * CinvertorError = NULL;
14
15#define INVECTOR(obj,buf,len)                                                                           \
16    do { \
17        int err = PyObject_AsReadBuffer(obj, (const void **)(&buf), &len); \
18        if (err < 0) return NULL; \
19        len /= sizeof(*buf); \
20    } while (0)
21   
22#define OUTVECTOR(obj,buf,len) \
23    do { \
24        int err = PyObject_AsWriteBuffer(obj, (void **)(&buf), &len); \
25        if (err < 0) return NULL; \
26        len /= sizeof(*buf); \
27    } while (0)
28
29
30// Class definition
31typedef struct {
32    PyObject_HEAD   
33    Invertor_params params; 
34} Cinvertor;
35
36
37static void
38Cinvertor_dealloc(Cinvertor* self)
39{
40    invertor_dealloc(&(self->params));
41     
42    self->ob_type->tp_free((PyObject*)self);
43
44}
45
46static PyObject *
47Cinvertor_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
48{
49    Cinvertor *self;
50   
51    self = (Cinvertor *)type->tp_alloc(type, 0);
52   
53    return (PyObject *)self;
54}
55
56static int
57Cinvertor_init(Cinvertor *self, PyObject *args, PyObject *kwds)
58{
59    if (self != NULL) {         
60        // Create parameters
61        invertor_init(&(self->params));
62    }
63    return 0;
64}
65
66static PyMemberDef Cinvertor_members[] = {
67    //{"params", T_OBJECT, offsetof(Cinvertor, params), 0,
68    // "Parameters"},
69    {NULL}  /* Sentinel */
70};
71
72
73/**
74 * Function to set the x data
75 * Takes an array of doubles as input
76 * Returns the number of entries found
77 */
78static PyObject * set_x(Cinvertor *self, PyObject *args) {
79        PyObject *data_obj;
80        Py_ssize_t ndata;
81        double *data;
82        int i;
83 
84        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
85        OUTVECTOR(data_obj,data,ndata);
86       
87        free(self->params.x);
88        self->params.x = (double*) malloc(ndata*sizeof(double));
89       
90        if(self->params.x==NULL) {
91            PyErr_SetString(CinvertorError, 
92                "Cinvertor.set_x: problem allocating memory.");
93                return NULL;           
94        }
95       
96        for (i=0; i<ndata; i++) {
97                self->params.x[i] = data[i];
98        }
99       
100        //self->params.x = data;
101        self->params.npoints = ndata;
102        return Py_BuildValue("i", self->params.npoints);       
103}
104
105static PyObject * get_x(Cinvertor *self, PyObject *args) {
106        PyObject *data_obj;
107        Py_ssize_t ndata;
108        double *data;
109    int i;
110   
111        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
112        OUTVECTOR(data_obj, data, ndata);
113       
114        // Check that the input array is large enough
115        if (ndata < self->params.npoints) {
116            PyErr_SetString(CinvertorError, 
117                "Cinvertor.get_x: input array too short for data.");
118                return NULL;           
119        }
120       
121        for(i=0; i<self->params.npoints; i++){
122                data[i] = self->params.x[i];
123        }
124       
125        return Py_BuildValue("i", self->params.npoints);       
126}
127
128/**
129 * Function to set the y data
130 * Takes an array of doubles as input
131 * Returns the number of entries found
132 */
133static PyObject * set_y(Cinvertor *self, PyObject *args) {
134        PyObject *data_obj;
135        Py_ssize_t ndata;
136        double *data;
137        int i;
138 
139        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
140        OUTVECTOR(data_obj,data,ndata);
141       
142        free(self->params.y);
143        self->params.y = (double*) malloc(ndata*sizeof(double));
144       
145        if(self->params.y==NULL) {
146            PyErr_SetString(CinvertorError, 
147                "Cinvertor.set_y: problem allocating memory.");
148                return NULL;           
149        }
150       
151        for (i=0; i<ndata; i++) {
152                self->params.y[i] = data[i];
153        }       
154       
155        //self->params.y = data;
156        self->params.ny = ndata;
157        return Py_BuildValue("i", self->params.ny);     
158}
159
160static PyObject * get_y(Cinvertor *self, PyObject *args) {
161        PyObject *data_obj;
162        Py_ssize_t ndata;
163        double *data;
164    int i;
165   
166        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
167        OUTVECTOR(data_obj, data, ndata);
168       
169        // Check that the input array is large enough
170        if (ndata < self->params.ny) {
171            PyErr_SetString(CinvertorError, 
172                "Cinvertor.get_y: input array too short for data.");
173                return NULL;           
174        }
175       
176        for(i=0; i<self->params.ny; i++){
177                data[i] = self->params.y[i];
178        }
179       
180        return Py_BuildValue("i", self->params.npoints);       
181}
182
183/**
184 * Function to set the x data
185 * Takes an array of doubles as input
186 * Returns the number of entries found
187 */
188static PyObject * set_err(Cinvertor *self, PyObject *args) {
189        PyObject *data_obj;
190        Py_ssize_t ndata;
191        double *data;
192        int i;
193 
194        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
195        OUTVECTOR(data_obj,data,ndata);
196       
197        free(self->params.err);
198        self->params.err = (double*) malloc(ndata*sizeof(double));
199       
200        if(self->params.err==NULL) {
201            PyErr_SetString(CinvertorError, 
202                "Cinvertor.set_err: problem allocating memory.");
203                return NULL;           
204        }
205       
206        for (i=0; i<ndata; i++) {
207                self->params.err[i] = data[i];
208        }
209       
210        //self->params.err = data;
211        self->params.nerr = ndata;
212        return Py_BuildValue("i", self->params.nerr);   
213}
214
215static PyObject * get_err(Cinvertor *self, PyObject *args) {
216        PyObject *data_obj;
217        Py_ssize_t ndata;
218        double *data;
219    int i;
220   
221        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
222        OUTVECTOR(data_obj, data, ndata);
223       
224        // Check that the input array is large enough
225        if (ndata < self->params.nerr) {
226            PyErr_SetString(CinvertorError, 
227                "Cinvertor.get_err: input array too short for data.");
228                return NULL;           
229        }
230       
231        for(i=0; i<self->params.nerr; i++){
232                data[i] = self->params.err[i];
233        }
234       
235        return Py_BuildValue("i", self->params.npoints);       
236}
237
238/**
239 * Check the validity of the stored data
240 * Returns the number of points if it's all good, -1 otherwise
241 */
242static PyObject * is_valid(Cinvertor *self, PyObject *args) {
243        if(self->params.npoints==self->params.ny &&
244                        self->params.npoints==self->params.nerr) {
245                return Py_BuildValue("i", self->params.npoints);
246        } else {
247                return Py_BuildValue("i", -1);
248        }       
249}
250
251/**
252 * Sets the maximum distance
253 */
254static PyObject * set_dmax(Cinvertor *self, PyObject *args) {
255        double d_max;
256 
257        if (!PyArg_ParseTuple(args, "d", &d_max)) return NULL;
258        self->params.d_max = d_max;
259        return Py_BuildValue("d", self->params.d_max); 
260}
261
262/**
263 * Gets the maximum distance
264 */
265static PyObject * get_dmax(Cinvertor *self, PyObject *args) {
266        return Py_BuildValue("d", self->params.d_max); 
267}
268
269/**
270 * Sets the minimum q
271 */
272static PyObject * set_qmin(Cinvertor *self, PyObject *args) {
273        double q_min;
274 
275        if (!PyArg_ParseTuple(args, "d", &q_min)) return NULL;
276        self->params.q_min = q_min;
277        return Py_BuildValue("d", self->params.q_min); 
278}
279
280/**
281 * Gets the minimum q
282 */
283static PyObject * get_qmin(Cinvertor *self, PyObject *args) {
284        return Py_BuildValue("d", self->params.q_min); 
285}
286
287
288/**
289 * Sets the maximum q
290 */
291static PyObject * set_qmax(Cinvertor *self, PyObject *args) {
292        double q_max;
293 
294        if (!PyArg_ParseTuple(args, "d", &q_max)) return NULL;
295        self->params.q_max = q_max;
296        return Py_BuildValue("d", self->params.q_max); 
297}
298
299/**
300 * Gets the maximum q
301 */
302static PyObject * get_qmax(Cinvertor *self, PyObject *args) {
303        return Py_BuildValue("d", self->params.q_max); 
304}
305
306
307static PyObject * set_alpha(Cinvertor *self, PyObject *args) {
308        double alpha;
309 
310        if (!PyArg_ParseTuple(args, "d", &alpha)) return NULL;
311        self->params.alpha = alpha;
312        return Py_BuildValue("d", self->params.alpha); 
313}
314
315/**
316 * Gets the maximum distance
317 */
318static PyObject * get_alpha(Cinvertor *self, PyObject *args) {
319        return Py_BuildValue("d", self->params.alpha); 
320}
321
322/**
323 * Gets the number of x points
324 */
325static PyObject * get_nx(Cinvertor *self, PyObject *args) {
326        return Py_BuildValue("i", self->params.npoints);       
327}
328
329/**
330 * Gets the number of y points
331 */
332static PyObject * get_ny(Cinvertor *self, PyObject *args) {
333        return Py_BuildValue("i", self->params.ny);     
334}
335
336/**
337 * Gets the number of error points
338 */
339static PyObject * get_nerr(Cinvertor *self, PyObject *args) {
340        return Py_BuildValue("i", self->params.nerr);   
341}
342
343
344/**
345 * Function to call to evaluate the residuals
346 * @param args: input parameters
347 * @return: list of residuals
348 */
349static PyObject * residuals(Cinvertor *self, PyObject *args) {
350        double *pars;
351        PyObject* residuals;
352        PyObject* temp;
353        double *res;
354        int i;
355        double residual, diff;
356        // Regularization factor
357        double regterm = 0.0;
358        double tmp = 0.0;
359        // Number of slices in regularization term estimate
360        int nslice = 25;
361       
362        PyObject *data_obj;
363        Py_ssize_t npars;
364         
365        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
366       
367        OUTVECTOR(data_obj,pars,npars);
368               
369    // PyList of residuals
370        // Should create this list only once and refill it
371    residuals = PyList_New(self->params.npoints);
372
373    regterm = reg_term(pars, self->params.d_max, npars, nslice);
374   
375    for(i=0; i<self->params.npoints; i++) {
376        diff = self->params.y[i] - iq(pars, self->params.d_max, npars, self->params.x[i]);
377        residual = diff*diff / (self->params.err[i]*self->params.err[i]);
378        tmp = residual;
379       
380        // regularization term
381        residual += self->params.alpha * regterm;
382       
383        if (PyList_SetItem(residuals, i, Py_BuildValue("d",residual) ) < 0){
384            PyErr_SetString(CinvertorError, 
385                "Cinvertor.residuals: error setting residual.");
386                return NULL;
387        };
388               
389    }
390   
391        return residuals;
392}
393/**
394 * Function to call to evaluate the residuals
395 * for P(r) minimization (for testing purposes)
396 * @param args: input parameters
397 * @return: list of residuals
398 */
399static PyObject * pr_residuals(Cinvertor *self, PyObject *args) {
400        double *pars;
401        PyObject* residuals;
402        PyObject* temp;
403        double *res;
404        int i;
405        double residual, diff;
406        // Regularization factor
407        double regterm = 0.0;
408        double tmp = 0.0;
409        // Number of slices in regularization term estimate
410        int nslice = 25;
411       
412        PyObject *data_obj;
413        Py_ssize_t npars;
414         
415        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
416       
417        OUTVECTOR(data_obj,pars,npars);
418               
419        // Should create this list only once and refill it
420    residuals = PyList_New(self->params.npoints);
421
422    regterm = reg_term(pars, self->params.d_max, npars, nslice);
423
424   
425    for(i=0; i<self->params.npoints; i++) {
426        diff = self->params.y[i] - pr(pars, self->params.d_max, npars, self->params.x[i]);
427        residual = diff*diff / (self->params.err[i]*self->params.err[i]);
428        tmp = residual;
429       
430        // regularization term
431        residual += self->params.alpha * regterm;
432       
433        if (PyList_SetItem(residuals, i, Py_BuildValue("d",residual) ) < 0){
434            PyErr_SetString(CinvertorError, 
435                "Cinvertor.residuals: error setting residual.");
436                return NULL;
437        };
438               
439    }
440   
441        return residuals;
442}
443
444/**
445 * Function to call to evaluate the scattering intensity
446 * @param args: c-parameters, and q
447 * @return: I(q)
448 */
449static PyObject * get_iq(Cinvertor *self, PyObject *args) {
450        double *pars;
451        double q, iq_value;
452        PyObject *data_obj;
453        Py_ssize_t npars;
454         
455        if (!PyArg_ParseTuple(args, "Od", &data_obj, &q)) return NULL;
456        OUTVECTOR(data_obj,pars,npars);
457               
458        iq_value = iq(pars, self->params.d_max, npars, q);
459        return Py_BuildValue("f", iq_value);   
460}
461
462/**
463 * Function to call to evaluate P(r)
464 * @param args: c-parameters and r
465 * @return: P(r)
466 */
467static PyObject * get_pr(Cinvertor *self, PyObject *args) {
468        double *pars;
469        double r, pr_value;
470        PyObject *data_obj;
471        Py_ssize_t npars;
472         
473        if (!PyArg_ParseTuple(args, "Od", &data_obj, &r)) return NULL;
474        OUTVECTOR(data_obj,pars,npars);
475               
476        pr_value = pr(pars, self->params.d_max, npars, r);
477        return Py_BuildValue("f", pr_value);   
478}
479
480/**
481 * Function to call to evaluate P(r) with errors
482 * @param args: c-parameters and r
483 * @return: P(r)
484 */
485static PyObject * get_pr_err(Cinvertor *self, PyObject *args) {
486        double *pars;
487        double *pars_err;
488        double pr_err_value;
489        double r, pr_value;
490        PyObject *data_obj;
491        Py_ssize_t npars;
492        PyObject *err_obj;
493        Py_ssize_t npars2;
494        int i; 
495       
496        if (!PyArg_ParseTuple(args, "OOd", &data_obj, &err_obj, &r)) return NULL;
497        OUTVECTOR(data_obj,pars,npars); 
498        OUTVECTOR(err_obj,pars_err,npars2);
499
500        pr_err(pars, pars_err, self->params.d_max, npars, r, &pr_value, &pr_err_value);
501        return Py_BuildValue("ff", pr_value, pr_err_value);     
502}
503
504static PyObject * basefunc_ft(Cinvertor *self, PyObject *args) {
505        double d_max, q;
506        int n;
507       
508        if (!PyArg_ParseTuple(args, "did", &d_max, &n, &q)) return NULL;
509        return Py_BuildValue("f", ortho_transformed(d_max, n, q));     
510       
511}
512
513static PyObject * oscillations(Cinvertor *self, PyObject *args) {
514        double *pars;
515        PyObject *data_obj;
516        Py_ssize_t npars;
517        double oscill, norm;
518       
519        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
520        OUTVECTOR(data_obj,pars,npars);
521       
522        oscill = reg_term(pars, self->params.d_max, npars, 100);
523        norm   = int_p2(pars, self->params.d_max, npars, 100);
524        return Py_BuildValue("f", sqrt(oscill/norm)/acos(-1.0)*self->params.d_max );   
525       
526}
527
528static PyObject * get_peaks(Cinvertor *self, PyObject *args) {
529        double *pars;
530        PyObject *data_obj;
531        Py_ssize_t npars;
532        int count;
533       
534        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
535        OUTVECTOR(data_obj,pars,npars);
536       
537        count = npeaks(pars, self->params.d_max, npars, 100);
538
539        return Py_BuildValue("i", count );     
540       
541}
542
543static PyObject * get_positive(Cinvertor *self, PyObject *args) {
544        double *pars;
545        PyObject *data_obj;
546        Py_ssize_t npars;
547        double fraction;
548         
549        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
550        OUTVECTOR(data_obj,pars,npars);
551       
552        fraction = positive_integral(pars, self->params.d_max, npars, 100);
553
554        return Py_BuildValue("f", fraction );   
555       
556}
557
558static PyObject * get_pos_err(Cinvertor *self, PyObject *args) {
559        double *pars;
560        double *pars_err;
561        PyObject *data_obj;
562        PyObject *err_obj;
563        Py_ssize_t npars;
564        Py_ssize_t npars2;
565        double fraction;
566       
567        if (!PyArg_ParseTuple(args, "OO", &data_obj, &err_obj)) return NULL;
568        OUTVECTOR(data_obj,pars,npars); 
569        OUTVECTOR(err_obj,pars_err,npars2);
570       
571        fraction = positive_errors(pars, pars_err, self->params.d_max, npars, 51);
572
573        return Py_BuildValue("f", fraction );   
574       
575}
576
577static PyMethodDef Cinvertor_methods[] = {
578                   {"residuals", (PyCFunction)residuals, METH_VARARGS, "Get the list of residuals"},
579                   {"pr_residuals", (PyCFunction)pr_residuals, METH_VARARGS, "Get the list of residuals"},
580                   {"set_x", (PyCFunction)set_x, METH_VARARGS, ""},
581                   {"get_x", (PyCFunction)get_x, METH_VARARGS, ""},
582                   {"set_y", (PyCFunction)set_y, METH_VARARGS, ""},
583                   {"get_y", (PyCFunction)get_y, METH_VARARGS, ""},
584                   {"set_err", (PyCFunction)set_err, METH_VARARGS, ""},
585                   {"get_err", (PyCFunction)get_err, METH_VARARGS, ""},
586                   {"set_dmax", (PyCFunction)set_dmax, METH_VARARGS, ""},
587                   {"get_dmax", (PyCFunction)get_dmax, METH_VARARGS, ""},
588                   {"set_qmin", (PyCFunction)set_qmin, METH_VARARGS, ""},
589                   {"get_qmin", (PyCFunction)get_qmin, METH_VARARGS, ""},
590                   {"set_qmax", (PyCFunction)set_qmax, METH_VARARGS, ""},
591                   {"get_qmax", (PyCFunction)get_qmax, METH_VARARGS, ""},
592                   {"set_alpha", (PyCFunction)set_alpha, METH_VARARGS, ""},
593                   {"get_alpha", (PyCFunction)get_alpha, METH_VARARGS, ""},
594                   {"get_nx", (PyCFunction)get_nx, METH_VARARGS, ""},
595                   {"get_ny", (PyCFunction)get_ny, METH_VARARGS, ""},
596                   {"get_nerr", (PyCFunction)get_nerr, METH_VARARGS, ""},
597                   {"iq", (PyCFunction)get_iq, METH_VARARGS, ""},
598                   {"pr", (PyCFunction)get_pr, METH_VARARGS, ""},
599                   {"get_pr_err", (PyCFunction)get_pr_err, METH_VARARGS, ""},
600                   {"is_valid", (PyCFunction)is_valid, METH_VARARGS, ""},
601                   {"basefunc_ft", (PyCFunction)basefunc_ft, METH_VARARGS, ""},
602                   {"oscillations", (PyCFunction)oscillations, METH_VARARGS, ""},
603                   {"get_peaks", (PyCFunction)get_peaks, METH_VARARGS, ""},
604                   {"get_positive", (PyCFunction)get_positive, METH_VARARGS, ""},
605                   {"get_pos_err", (PyCFunction)get_pos_err, METH_VARARGS, ""},
606   
607   {NULL}
608};
609
610static PyTypeObject CinvertorType = {
611    PyObject_HEAD_INIT(NULL)
612    0,                         /*ob_size*/
613    "Cinvertor",             /*tp_name*/
614    sizeof(Cinvertor),             /*tp_basicsize*/
615    0,                         /*tp_itemsize*/
616    (destructor)Cinvertor_dealloc, /*tp_dealloc*/
617    0,                         /*tp_print*/
618    0,                         /*tp_getattr*/
619    0,                         /*tp_setattr*/
620    0,                         /*tp_compare*/
621    0,                         /*tp_repr*/
622    0,                         /*tp_as_number*/
623    0,                         /*tp_as_sequence*/
624    0,                         /*tp_as_mapping*/
625    0,                         /*tp_hash */
626    0,                         /*tp_call*/
627    0,                         /*tp_str*/
628    0,                         /*tp_getattro*/
629    0,                         /*tp_setattro*/
630    0,                         /*tp_as_buffer*/
631    Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /*tp_flags*/
632    "Cinvertor objects",           /* tp_doc */
633    0,                         /* tp_traverse */
634    0,                         /* tp_clear */
635    0,                         /* tp_richcompare */
636    0,                         /* tp_weaklistoffset */
637    0,                         /* tp_iter */
638    0,                         /* tp_iternext */
639    Cinvertor_methods,             /* tp_methods */
640    Cinvertor_members,             /* tp_members */
641    0,                         /* tp_getset */
642    0,                         /* tp_base */
643    0,                         /* tp_dict */
644    0,                         /* tp_descr_get */
645    0,                         /* tp_descr_set */
646    0,                         /* tp_dictoffset */
647    (initproc)Cinvertor_init,      /* tp_init */
648    0,                         /* tp_alloc */
649    Cinvertor_new,                 /* tp_new */
650};
651
652
653static PyMethodDef module_methods[] = {
654    {NULL} 
655};
656
657/**
658 * Function used to add the model class to a module
659 * @param module: module to add the class to
660 */ 
661void addCinvertor(PyObject *module) {
662        PyObject *d;
663       
664    if (PyType_Ready(&CinvertorType) < 0)
665        return;
666
667    Py_INCREF(&CinvertorType);
668    PyModule_AddObject(module, "Cinvertor", (PyObject *)&CinvertorType);
669   
670    d = PyModule_GetDict(module);
671    CinvertorError = PyErr_NewException("Cinvertor.error", NULL, NULL);
672    PyDict_SetItemString(d, "CinvertorError", CinvertorError);
673}
674
675
676#ifndef PyMODINIT_FUNC  /* declarations for DLL import/export */
677#define PyMODINIT_FUNC void
678#endif
679PyMODINIT_FUNC
680initpr_inversion(void) 
681{
682    PyObject* m;
683
684    m = Py_InitModule3("pr_inversion", module_methods,
685                       "C extension module for inversion to P(r).");
686                       
687    addCinvertor(m);
688}
Note: See TracBrowser for help on using the repository browser.