source: sasview/pr_inversion/c_extensions/Cinvertor.c @ 34ae302

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 34ae302 was eca05c8, checked in by Mathieu Doucet <doucetm@…>, 17 years ago

Added P(r) fit test
  • Property mode set to 100644
File size: 14.1 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 
83        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
84        OUTVECTOR(data_obj,data,ndata);
85        self->params.x = data;
86        self->params.npoints = ndata;
87        return Py_BuildValue("i", self->params.npoints);       
88}
89
90static PyObject * get_x(Cinvertor *self, PyObject *args) {
91        PyObject *data_obj;
92        Py_ssize_t ndata;
93        double *data;
94    int i;
95   
96        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
97        INVECTOR(data_obj, data, ndata);
98       
99        // Check that the input array is large enough
100        if (ndata < self->params.npoints) {
101            PyErr_SetString(CinvertorError, 
102                "Cinvertor.get_x: input array too short for data.");
103                return NULL;           
104        }
105       
106        for(i=0; i<self->params.npoints; i++){
107                data[i] = self->params.x[i];
108        }
109       
110        return Py_BuildValue("i", self->params.npoints);       
111}
112
113/**
114 * Function to set the y data
115 * Takes an array of doubles as input
116 * Returns the number of entries found
117 */
118static PyObject * set_y(Cinvertor *self, PyObject *args) {
119        PyObject *data_obj;
120        Py_ssize_t ndata;
121        double *data;
122 
123        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
124        OUTVECTOR(data_obj,data,ndata);
125        self->params.y = data;
126        self->params.ny = ndata;
127        return Py_BuildValue("i", self->params.ny);     
128}
129
130static PyObject * get_y(Cinvertor *self, PyObject *args) {
131        PyObject *data_obj;
132        Py_ssize_t ndata;
133        double *data;
134    int i;
135   
136        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
137        INVECTOR(data_obj, data, ndata);
138       
139        // Check that the input array is large enough
140        if (ndata < self->params.ny) {
141            PyErr_SetString(CinvertorError, 
142                "Cinvertor.get_y: input array too short for data.");
143                return NULL;           
144        }
145       
146        for(i=0; i<self->params.ny; i++){
147                data[i] = self->params.y[i];
148        }
149       
150        return Py_BuildValue("i", self->params.npoints);       
151}
152
153/**
154 * Function to set the x data
155 * Takes an array of doubles as input
156 * Returns the number of entries found
157 */
158static PyObject * set_err(Cinvertor *self, PyObject *args) {
159        PyObject *data_obj;
160        Py_ssize_t ndata;
161        double *data;
162 
163        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
164        OUTVECTOR(data_obj,data,ndata);
165        self->params.err = data;
166        self->params.nerr = ndata;
167        return Py_BuildValue("i", self->params.nerr);   
168}
169
170static PyObject * get_err(Cinvertor *self, PyObject *args) {
171        PyObject *data_obj;
172        Py_ssize_t ndata;
173        double *data;
174    int i;
175   
176        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
177        INVECTOR(data_obj, data, ndata);
178       
179        // Check that the input array is large enough
180        if (ndata < self->params.nerr) {
181            PyErr_SetString(CinvertorError, 
182                "Cinvertor.get_err: input array too short for data.");
183                return NULL;           
184        }
185       
186        for(i=0; i<self->params.nerr; i++){
187                data[i] = self->params.err[i];
188        }
189       
190        return Py_BuildValue("i", self->params.npoints);       
191}
192
193/**
194 * Check the validity of the stored data
195 * Returns the number of points if it's all good, -1 otherwise
196 */
197static PyObject * is_valid(Cinvertor *self, PyObject *args) {
198        if(self->params.npoints==self->params.ny &&
199                        self->params.npoints==self->params.nerr) {
200                return Py_BuildValue("i", self->params.npoints);
201        } else {
202                return Py_BuildValue("i", -1);
203        }       
204}
205
206/**
207 * Sets the maximum distance
208 */
209static PyObject * set_dmax(Cinvertor *self, PyObject *args) {
210        double d_max;
211 
212        if (!PyArg_ParseTuple(args, "d", &d_max)) return NULL;
213        self->params.d_max = d_max;
214        return Py_BuildValue("d", self->params.d_max); 
215}
216
217/**
218 * Gets the maximum distance
219 */
220static PyObject * get_dmax(Cinvertor *self, PyObject *args) {
221        return Py_BuildValue("d", self->params.d_max); 
222}
223
224static PyObject * set_alpha(Cinvertor *self, PyObject *args) {
225        double alpha;
226 
227        if (!PyArg_ParseTuple(args, "d", &alpha)) return NULL;
228        self->params.alpha = alpha;
229        return Py_BuildValue("d", self->params.alpha); 
230}
231
232/**
233 * Gets the maximum distance
234 */
235static PyObject * get_alpha(Cinvertor *self, PyObject *args) {
236        return Py_BuildValue("d", self->params.alpha); 
237}
238
239/**
240 * Gets the number of x points
241 */
242static PyObject * get_nx(Cinvertor *self, PyObject *args) {
243        return Py_BuildValue("i", self->params.npoints);       
244}
245
246/**
247 * Gets the number of y points
248 */
249static PyObject * get_ny(Cinvertor *self, PyObject *args) {
250        return Py_BuildValue("i", self->params.ny);     
251}
252
253/**
254 * Gets the number of error points
255 */
256static PyObject * get_nerr(Cinvertor *self, PyObject *args) {
257        return Py_BuildValue("i", self->params.nerr);   
258}
259
260
261/**
262 * Function to call to evaluate the residuals
263 * @param args: input parameters
264 * @return: list of residuals
265 */
266static PyObject * residuals(Cinvertor *self, PyObject *args) {
267        double *pars;
268        PyObject* residuals;
269        PyObject* temp;
270        double *res;
271        int i;
272        double residual, diff;
273        // Regularization factor
274        double regterm = 0.0;
275        double tmp = 0.0;
276       
277        PyObject *data_obj;
278        Py_ssize_t npars;
279         
280        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
281       
282        OUTVECTOR(data_obj,pars,npars);
283               
284    // PyList of residuals
285        // Should create this list only once and refill it
286    residuals = PyList_New(self->params.npoints);
287
288    regterm = reg_term(pars, self->params.d_max, npars);
289   
290    for(i=0; i<self->params.npoints; i++) {
291        diff = self->params.y[i] - iq(pars, self->params.d_max, npars, self->params.x[i]);
292        residual = diff*diff / (self->params.err[i]*self->params.err[i]);
293        tmp = residual;
294       
295        // regularization term
296        residual += self->params.alpha * regterm;
297       
298        if (PyList_SetItem(residuals, i, Py_BuildValue("d",residual) ) < 0){
299            PyErr_SetString(CinvertorError, 
300                "Cinvertor.residuals: error setting residual.");
301                return NULL;
302        };
303               
304    }
305   
306        return residuals;
307}
308/**
309 * Function to call to evaluate the residuals
310 * for P(r) minimization (for testing purposes)
311 * @param args: input parameters
312 * @return: list of residuals
313 */
314static PyObject * pr_residuals(Cinvertor *self, PyObject *args) {
315        double *pars;
316        PyObject* residuals;
317        PyObject* temp;
318        double *res;
319        int i;
320        double residual, diff;
321        // Regularization factor
322        double regterm = 0.0;
323        double tmp = 0.0;
324       
325        PyObject *data_obj;
326        Py_ssize_t npars;
327         
328        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
329       
330        OUTVECTOR(data_obj,pars,npars);
331               
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);
336
337   
338    for(i=0; i<self->params.npoints; i++) {
339        diff = self->params.y[i] - pr(pars, self->params.d_max, npars, self->params.x[i]);
340        residual = diff*diff / (self->params.err[i]*self->params.err[i]);
341        tmp = residual;
342       
343        // regularization term
344        residual += self->params.alpha * regterm;
345       
346        if (PyList_SetItem(residuals, i, Py_BuildValue("d",residual) ) < 0){
347            PyErr_SetString(CinvertorError, 
348                "Cinvertor.residuals: error setting residual.");
349                return NULL;
350        };
351               
352    }
353   
354        return residuals;
355}
356
357/**
358 * Function to call to evaluate the scattering intensity
359 * @param args: c-parameters, and q
360 * @return: I(q)
361 */
362static PyObject * get_iq(Cinvertor *self, PyObject *args) {
363        double *pars;
364        double q, iq_value;
365        PyObject *data_obj;
366        Py_ssize_t npars;
367         
368        if (!PyArg_ParseTuple(args, "Od", &data_obj, &q)) return NULL;
369        OUTVECTOR(data_obj,pars,npars);
370               
371        iq_value = iq(pars, self->params.d_max, npars, q);
372        return Py_BuildValue("f", iq_value);   
373}
374
375/**
376 * Function to call to evaluate P(r)
377 * @param args: c-parameters and r
378 * @return: P(r)
379 */
380static PyObject * get_pr(Cinvertor *self, PyObject *args) {
381        double *pars;
382        double r, pr_value;
383        PyObject *data_obj;
384        Py_ssize_t npars;
385         
386        if (!PyArg_ParseTuple(args, "Od", &data_obj, &r)) return NULL;
387        OUTVECTOR(data_obj,pars,npars);
388               
389        pr_value = pr(pars, self->params.d_max, npars, r);
390        return Py_BuildValue("f", pr_value);   
391}
392
393/**
394 * Function to call to evaluate P(r) with errors
395 * @param args: c-parameters and r
396 * @return: P(r)
397 */
398static PyObject * get_pr_err(Cinvertor *self, PyObject *args) {
399        double *pars;
400        double *pars_err;
401        double pr_err_value;
402        double r, pr_value;
403        PyObject *data_obj;
404        Py_ssize_t npars;
405        PyObject *err_obj;
406        Py_ssize_t npars2;
407         
408        if (!PyArg_ParseTuple(args, "OOd", &data_obj, &err_obj, &r)) return NULL;
409        OUTVECTOR(data_obj,pars,npars);
410        OUTVECTOR(err_obj,pars_err,npars2);
411               
412        pr_err(pars, pars_err, self->params.d_max, npars, r, &pr_value, &pr_err_value);
413        return Py_BuildValue("ff", pr_value, pr_err_value);     
414}
415
416
417static PyMethodDef Cinvertor_methods[] = {
418                   {"residuals", (PyCFunction)residuals, METH_VARARGS, "Get the list of residuals"},
419                   {"pr_residuals", (PyCFunction)pr_residuals, METH_VARARGS, "Get the list of residuals"},
420                   {"set_x", (PyCFunction)set_x, METH_VARARGS, ""},
421                   {"get_x", (PyCFunction)get_x, METH_VARARGS, ""},
422                   {"set_y", (PyCFunction)set_y, METH_VARARGS, ""},
423                   {"get_y", (PyCFunction)get_y, METH_VARARGS, ""},
424                   {"set_err", (PyCFunction)set_err, METH_VARARGS, ""},
425                   {"get_err", (PyCFunction)get_err, METH_VARARGS, ""},
426                   {"set_dmax", (PyCFunction)set_dmax, METH_VARARGS, ""},
427                   {"get_dmax", (PyCFunction)get_dmax, METH_VARARGS, ""},
428                   {"set_alpha", (PyCFunction)set_alpha, METH_VARARGS, ""},
429                   {"get_alpha", (PyCFunction)get_alpha, METH_VARARGS, ""},
430                   {"get_nx", (PyCFunction)get_nx, METH_VARARGS, ""},
431                   {"get_ny", (PyCFunction)get_ny, METH_VARARGS, ""},
432                   {"get_nerr", (PyCFunction)get_nerr, METH_VARARGS, ""},
433                   {"iq", (PyCFunction)get_iq, METH_VARARGS, ""},
434                   {"pr", (PyCFunction)get_pr, METH_VARARGS, ""},
435                   {"get_pr_err", (PyCFunction)get_pr_err, METH_VARARGS, ""},
436                   {"is_valid", (PyCFunction)is_valid, METH_VARARGS, ""},
437   
438   {NULL}
439};
440
441static PyTypeObject CinvertorType = {
442    PyObject_HEAD_INIT(NULL)
443    0,                         /*ob_size*/
444    "Cinvertor",             /*tp_name*/
445    sizeof(Cinvertor),             /*tp_basicsize*/
446    0,                         /*tp_itemsize*/
447    (destructor)Cinvertor_dealloc, /*tp_dealloc*/
448    0,                         /*tp_print*/
449    0,                         /*tp_getattr*/
450    0,                         /*tp_setattr*/
451    0,                         /*tp_compare*/
452    0,                         /*tp_repr*/
453    0,                         /*tp_as_number*/
454    0,                         /*tp_as_sequence*/
455    0,                         /*tp_as_mapping*/
456    0,                         /*tp_hash */
457    0,                         /*tp_call*/
458    0,                         /*tp_str*/
459    0,                         /*tp_getattro*/
460    0,                         /*tp_setattro*/
461    0,                         /*tp_as_buffer*/
462    Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /*tp_flags*/
463    "Cinvertor objects",           /* tp_doc */
464    0,                         /* tp_traverse */
465    0,                         /* tp_clear */
466    0,                         /* tp_richcompare */
467    0,                         /* tp_weaklistoffset */
468    0,                         /* tp_iter */
469    0,                         /* tp_iternext */
470    Cinvertor_methods,             /* tp_methods */
471    Cinvertor_members,             /* tp_members */
472    0,                         /* tp_getset */
473    0,                         /* tp_base */
474    0,                         /* tp_dict */
475    0,                         /* tp_descr_get */
476    0,                         /* tp_descr_set */
477    0,                         /* tp_dictoffset */
478    (initproc)Cinvertor_init,      /* tp_init */
479    0,                         /* tp_alloc */
480    Cinvertor_new,                 /* tp_new */
481};
482
483
484static PyMethodDef module_methods[] = {
485    {NULL} 
486};
487
488/**
489 * Function used to add the model class to a module
490 * @param module: module to add the class to
491 */ 
492void addCinvertor(PyObject *module) {
493        PyObject *d;
494       
495    if (PyType_Ready(&CinvertorType) < 0)
496        return;
497
498    Py_INCREF(&CinvertorType);
499    PyModule_AddObject(module, "Cinvertor", (PyObject *)&CinvertorType);
500   
501    d = PyModule_GetDict(module);
502    CinvertorError = PyErr_NewException("Cinvertor.error", NULL, NULL);
503    PyDict_SetItemString(d, "CinvertorError", CinvertorError);
504}
505
506
507#ifndef PyMODINIT_FUNC  /* declarations for DLL import/export */
508#define PyMODINIT_FUNC void
509#endif
510PyMODINIT_FUNC
511initpr_inversion(void) 
512{
513    PyObject* m;
514
515    m = Py_InitModule3("pr_inversion", module_methods,
516                       "C extension module for inversion to P(r).");
517                       
518    addCinvertor(m);
519}
Note: See TracBrowser for help on using the repository browser.