source: sasview/pr_inversion/c_extensions/Cinvertor.c @ 4241e48

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 4241e48 was abad620, checked in by Mathieu Doucet <doucetm@…>, 16 years ago

Added alpha estimate
  • Property mode set to 100644
File size: 16.0 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
269static PyObject * set_alpha(Cinvertor *self, PyObject *args) {
270        double alpha;
271 
272        if (!PyArg_ParseTuple(args, "d", &alpha)) return NULL;
273        self->params.alpha = alpha;
274        return Py_BuildValue("d", self->params.alpha); 
275}
276
277/**
278 * Gets the maximum distance
279 */
280static PyObject * get_alpha(Cinvertor *self, PyObject *args) {
281        return Py_BuildValue("d", self->params.alpha); 
282}
283
284/**
285 * Gets the number of x points
286 */
287static PyObject * get_nx(Cinvertor *self, PyObject *args) {
288        return Py_BuildValue("i", self->params.npoints);       
289}
290
291/**
292 * Gets the number of y points
293 */
294static PyObject * get_ny(Cinvertor *self, PyObject *args) {
295        return Py_BuildValue("i", self->params.ny);     
296}
297
298/**
299 * Gets the number of error points
300 */
301static PyObject * get_nerr(Cinvertor *self, PyObject *args) {
302        return Py_BuildValue("i", self->params.nerr);   
303}
304
305
306/**
307 * Function to call to evaluate the residuals
308 * @param args: input parameters
309 * @return: list of residuals
310 */
311static PyObject * residuals(Cinvertor *self, PyObject *args) {
312        double *pars;
313        PyObject* residuals;
314        PyObject* temp;
315        double *res;
316        int i;
317        double residual, diff;
318        // Regularization factor
319        double regterm = 0.0;
320        double tmp = 0.0;
321        // Number of slices in regularization term estimate
322        int nslice = 25;
323       
324        PyObject *data_obj;
325        Py_ssize_t npars;
326         
327        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
328       
329        OUTVECTOR(data_obj,pars,npars);
330               
331    // PyList of residuals
332        // Should create this list only once and refill it
333    residuals = PyList_New(self->params.npoints);
334
335    regterm = reg_term(pars, self->params.d_max, npars, nslice);
336   
337    for(i=0; i<self->params.npoints; i++) {
338        diff = self->params.y[i] - iq(pars, self->params.d_max, npars, self->params.x[i]);
339        residual = diff*diff / (self->params.err[i]*self->params.err[i]);
340        tmp = residual;
341       
342        // regularization term
343        residual += self->params.alpha * regterm;
344       
345        if (PyList_SetItem(residuals, i, Py_BuildValue("d",residual) ) < 0){
346            PyErr_SetString(CinvertorError, 
347                "Cinvertor.residuals: error setting residual.");
348                return NULL;
349        };
350               
351    }
352   
353        return residuals;
354}
355/**
356 * Function to call to evaluate the residuals
357 * for P(r) minimization (for testing purposes)
358 * @param args: input parameters
359 * @return: list of residuals
360 */
361static PyObject * pr_residuals(Cinvertor *self, PyObject *args) {
362        double *pars;
363        PyObject* residuals;
364        PyObject* temp;
365        double *res;
366        int i;
367        double residual, diff;
368        // Regularization factor
369        double regterm = 0.0;
370        double tmp = 0.0;
371        // Number of slices in regularization term estimate
372        int nslice = 25;
373       
374        PyObject *data_obj;
375        Py_ssize_t npars;
376         
377        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
378       
379        OUTVECTOR(data_obj,pars,npars);
380               
381        // Should create this list only once and refill it
382    residuals = PyList_New(self->params.npoints);
383
384    regterm = reg_term(pars, self->params.d_max, npars, nslice);
385
386   
387    for(i=0; i<self->params.npoints; i++) {
388        diff = self->params.y[i] - pr(pars, self->params.d_max, npars, self->params.x[i]);
389        residual = diff*diff / (self->params.err[i]*self->params.err[i]);
390        tmp = residual;
391       
392        // regularization term
393        residual += self->params.alpha * regterm;
394       
395        if (PyList_SetItem(residuals, i, Py_BuildValue("d",residual) ) < 0){
396            PyErr_SetString(CinvertorError, 
397                "Cinvertor.residuals: error setting residual.");
398                return NULL;
399        };
400               
401    }
402   
403        return residuals;
404}
405
406/**
407 * Function to call to evaluate the scattering intensity
408 * @param args: c-parameters, and q
409 * @return: I(q)
410 */
411static PyObject * get_iq(Cinvertor *self, PyObject *args) {
412        double *pars;
413        double q, iq_value;
414        PyObject *data_obj;
415        Py_ssize_t npars;
416         
417        if (!PyArg_ParseTuple(args, "Od", &data_obj, &q)) return NULL;
418        OUTVECTOR(data_obj,pars,npars);
419               
420        iq_value = iq(pars, self->params.d_max, npars, q);
421        return Py_BuildValue("f", iq_value);   
422}
423
424/**
425 * Function to call to evaluate P(r)
426 * @param args: c-parameters and r
427 * @return: P(r)
428 */
429static PyObject * get_pr(Cinvertor *self, PyObject *args) {
430        double *pars;
431        double r, pr_value;
432        PyObject *data_obj;
433        Py_ssize_t npars;
434         
435        if (!PyArg_ParseTuple(args, "Od", &data_obj, &r)) return NULL;
436        OUTVECTOR(data_obj,pars,npars);
437               
438        pr_value = pr(pars, self->params.d_max, npars, r);
439        return Py_BuildValue("f", pr_value);   
440}
441
442/**
443 * Function to call to evaluate P(r) with errors
444 * @param args: c-parameters and r
445 * @return: P(r)
446 */
447static PyObject * get_pr_err(Cinvertor *self, PyObject *args) {
448        double *pars;
449        double *pars_err;
450        double pr_err_value;
451        double r, pr_value;
452        PyObject *data_obj;
453        Py_ssize_t npars;
454        PyObject *err_obj;
455        Py_ssize_t npars2;
456         
457        if (!PyArg_ParseTuple(args, "OOd", &data_obj, &err_obj, &r)) return NULL;
458        OUTVECTOR(data_obj,pars,npars);
459        OUTVECTOR(err_obj,pars_err,npars2);
460               
461        pr_err(pars, pars_err, self->params.d_max, npars, r, &pr_value, &pr_err_value);
462        return Py_BuildValue("ff", pr_value, pr_err_value);     
463}
464
465static PyObject * basefunc_ft(Cinvertor *self, PyObject *args) {
466        double d_max, q;
467        int n;
468       
469        if (!PyArg_ParseTuple(args, "did", &d_max, &n, &q)) return NULL;
470        return Py_BuildValue("f", ortho_transformed(d_max, n, q));     
471       
472}
473
474static PyObject * oscillations(Cinvertor *self, PyObject *args) {
475        double *pars;
476        PyObject *data_obj;
477        Py_ssize_t npars;
478        double oscill, norm;
479       
480        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
481        OUTVECTOR(data_obj,pars,npars);
482       
483        oscill = reg_term(pars, self->params.d_max, npars, 100);
484        norm   = int_p2(pars, self->params.d_max, npars, 100);
485        return Py_BuildValue("f", sqrt(oscill/norm)/acos(-1.0)*self->params.d_max );   
486       
487}
488
489static PyMethodDef Cinvertor_methods[] = {
490                   {"residuals", (PyCFunction)residuals, METH_VARARGS, "Get the list of residuals"},
491                   {"pr_residuals", (PyCFunction)pr_residuals, METH_VARARGS, "Get the list of residuals"},
492                   {"set_x", (PyCFunction)set_x, METH_VARARGS, ""},
493                   {"get_x", (PyCFunction)get_x, METH_VARARGS, ""},
494                   {"set_y", (PyCFunction)set_y, METH_VARARGS, ""},
495                   {"get_y", (PyCFunction)get_y, METH_VARARGS, ""},
496                   {"set_err", (PyCFunction)set_err, METH_VARARGS, ""},
497                   {"get_err", (PyCFunction)get_err, METH_VARARGS, ""},
498                   {"set_dmax", (PyCFunction)set_dmax, METH_VARARGS, ""},
499                   {"get_dmax", (PyCFunction)get_dmax, METH_VARARGS, ""},
500                   {"set_alpha", (PyCFunction)set_alpha, METH_VARARGS, ""},
501                   {"get_alpha", (PyCFunction)get_alpha, METH_VARARGS, ""},
502                   {"get_nx", (PyCFunction)get_nx, METH_VARARGS, ""},
503                   {"get_ny", (PyCFunction)get_ny, METH_VARARGS, ""},
504                   {"get_nerr", (PyCFunction)get_nerr, METH_VARARGS, ""},
505                   {"iq", (PyCFunction)get_iq, METH_VARARGS, ""},
506                   {"pr", (PyCFunction)get_pr, METH_VARARGS, ""},
507                   {"get_pr_err", (PyCFunction)get_pr_err, METH_VARARGS, ""},
508                   {"is_valid", (PyCFunction)is_valid, METH_VARARGS, ""},
509                   {"basefunc_ft", (PyCFunction)basefunc_ft, METH_VARARGS, ""},
510                   {"oscillations", (PyCFunction)oscillations, METH_VARARGS, ""},
511   
512   {NULL}
513};
514
515static PyTypeObject CinvertorType = {
516    PyObject_HEAD_INIT(NULL)
517    0,                         /*ob_size*/
518    "Cinvertor",             /*tp_name*/
519    sizeof(Cinvertor),             /*tp_basicsize*/
520    0,                         /*tp_itemsize*/
521    (destructor)Cinvertor_dealloc, /*tp_dealloc*/
522    0,                         /*tp_print*/
523    0,                         /*tp_getattr*/
524    0,                         /*tp_setattr*/
525    0,                         /*tp_compare*/
526    0,                         /*tp_repr*/
527    0,                         /*tp_as_number*/
528    0,                         /*tp_as_sequence*/
529    0,                         /*tp_as_mapping*/
530    0,                         /*tp_hash */
531    0,                         /*tp_call*/
532    0,                         /*tp_str*/
533    0,                         /*tp_getattro*/
534    0,                         /*tp_setattro*/
535    0,                         /*tp_as_buffer*/
536    Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /*tp_flags*/
537    "Cinvertor objects",           /* tp_doc */
538    0,                         /* tp_traverse */
539    0,                         /* tp_clear */
540    0,                         /* tp_richcompare */
541    0,                         /* tp_weaklistoffset */
542    0,                         /* tp_iter */
543    0,                         /* tp_iternext */
544    Cinvertor_methods,             /* tp_methods */
545    Cinvertor_members,             /* tp_members */
546    0,                         /* tp_getset */
547    0,                         /* tp_base */
548    0,                         /* tp_dict */
549    0,                         /* tp_descr_get */
550    0,                         /* tp_descr_set */
551    0,                         /* tp_dictoffset */
552    (initproc)Cinvertor_init,      /* tp_init */
553    0,                         /* tp_alloc */
554    Cinvertor_new,                 /* tp_new */
555};
556
557
558static PyMethodDef module_methods[] = {
559    {NULL} 
560};
561
562/**
563 * Function used to add the model class to a module
564 * @param module: module to add the class to
565 */ 
566void addCinvertor(PyObject *module) {
567        PyObject *d;
568       
569    if (PyType_Ready(&CinvertorType) < 0)
570        return;
571
572    Py_INCREF(&CinvertorType);
573    PyModule_AddObject(module, "Cinvertor", (PyObject *)&CinvertorType);
574   
575    d = PyModule_GetDict(module);
576    CinvertorError = PyErr_NewException("Cinvertor.error", NULL, NULL);
577    PyDict_SetItemString(d, "CinvertorError", CinvertorError);
578}
579
580
581#ifndef PyMODINIT_FUNC  /* declarations for DLL import/export */
582#define PyMODINIT_FUNC void
583#endif
584PyMODINIT_FUNC
585initpr_inversion(void) 
586{
587    PyObject* m;
588
589    m = Py_InitModule3("pr_inversion", module_methods,
590                       "C extension module for inversion to P(r).");
591                       
592    addCinvertor(m);
593}
Note: See TracBrowser for help on using the repository browser.