source: sasview/pr_inversion/c_extensions/Cinvertor.c @ b43a009

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 b43a009 was 2d06beb, checked in by Mathieu Doucet <doucetm@…>, 17 years ago

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