source: sasview/src/sas/sascalc/pr/c_extensions/Cinvertor.c @ 952ea1f

magnetic_scattrelease-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249
Last change on this file since 952ea1f was 952ea1f, checked in by Paul Kienzle <pkienzle@…>, 11 months ago

move c extensions from core/module.so to _module.so

  • Property mode set to 100644
File size: 31.9 KB
Line 
1/**
2 * C implementation of the P(r) inversion
3 * Cinvertor is the base class for the Invertor class
4 * and provides the underlying computations.
5 *
6 */
7#include <stdio.h>
8#include <stdlib.h>
9#include <math.h>
10#include <time.h>
11
12//#define Py_LIMITED_API 0x03050000
13#include <Python.h>
14#include <structmember.h>
15
16#include "invertor.h"
17
18/// Error object for raised exceptions
19PyObject * CinvertorError;
20
21#define INVECTOR(obj,buf,len)                                                                           \
22    do { \
23        int err = PyObject_AsReadBuffer(obj, (const void **)(&buf), &len); \
24        if (err < 0) return NULL; \
25        len /= sizeof(*buf); \
26    } while (0)
27
28#define OUTVECTOR(obj,buf,len) \
29    do { \
30        int err = PyObject_AsWriteBuffer(obj, (void **)(&buf), &len); \
31        if (err < 0) return NULL; \
32        len /= sizeof(*buf); \
33    } while (0)
34
35
36// Class definition
37/**
38 * C implementation of the P(r) inversion
39 * Cinvertor is the base class for the Invertor class
40 * and provides the underlying computations.
41 *
42 */
43typedef struct {
44    PyObject_HEAD
45    /// Internal data structure
46    Invertor_params params;
47} Cinvertor;
48
49
50static void
51Cinvertor_dealloc(Cinvertor* self)
52{
53    invertor_dealloc(&(self->params));
54
55    Py_TYPE(self)->tp_free((PyObject*)self);
56
57}
58
59static PyObject *
60Cinvertor_new(PyTypeObject *type, PyObject *args, PyObject *kwds)
61{
62    Cinvertor *self;
63
64    self = (Cinvertor *)type->tp_alloc(type, 0);
65
66    return (PyObject *)self;
67}
68
69static int
70Cinvertor_init(Cinvertor *self, PyObject *args, PyObject *kwds)
71{
72    if (self != NULL) {
73        // Create parameters
74        invertor_init(&(self->params));
75    }
76    return 0;
77}
78
79static PyMemberDef Cinvertor_members[] = {
80    //{"params", T_OBJECT, offsetof(Cinvertor, params), 0,
81    // "Parameters"},
82    {NULL}  /* Sentinel */
83};
84
85const char set_x_doc[] =
86        "Function to set the x data\n"
87        "Takes an array of doubles as input.\n"
88        " @return: number of entries found";
89
90/**
91 * Function to set the x data
92 * Takes an array of doubles as input
93 * Returns the number of entries found
94 */
95static PyObject * set_x(Cinvertor *self, PyObject *args) {
96        PyObject *data_obj;
97        Py_ssize_t ndata;
98        double *data;
99        int i;
100
101        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
102        OUTVECTOR(data_obj,data,ndata);
103
104        free(self->params.x);
105        self->params.x = (double*) malloc(ndata*sizeof(double));
106
107        if(self->params.x==NULL) {
108            PyErr_SetString(CinvertorError,
109                "Cinvertor.set_x: problem allocating memory.");
110                return NULL;
111        }
112
113        for (i=0; i<ndata; i++) {
114                self->params.x[i] = data[i];
115        }
116
117        //self->params.x = data;
118        self->params.npoints = (int)ndata;
119        return Py_BuildValue("i", self->params.npoints);
120}
121
122const char get_x_doc[] =
123        "Function to get the x data\n"
124        "Takes an array of doubles as input.\n"
125        " @return: number of entries found";
126
127static PyObject * get_x(Cinvertor *self, PyObject *args) {
128        PyObject *data_obj;
129        Py_ssize_t ndata;
130        double *data;
131    int i;
132
133        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
134        OUTVECTOR(data_obj, data, ndata);
135
136        // Check that the input array is large enough
137        if (ndata < self->params.npoints) {
138            PyErr_SetString(CinvertorError,
139                "Cinvertor.get_x: input array too short for data.");
140                return NULL;
141        }
142
143        for(i=0; i<self->params.npoints; i++){
144                data[i] = self->params.x[i];
145        }
146
147        return Py_BuildValue("i", self->params.npoints);
148}
149
150const char set_y_doc[] =
151        "Function to set the y data\n"
152        "Takes an array of doubles as input.\n"
153        " @return: number of entries found";
154
155/**
156 * Function to set the y data
157 * Takes an array of doubles as input
158 * Returns the number of entries found
159 */
160static PyObject * set_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        free(self->params.y);
170        self->params.y = (double*) malloc(ndata*sizeof(double));
171
172        if(self->params.y==NULL) {
173            PyErr_SetString(CinvertorError,
174                "Cinvertor.set_y: problem allocating memory.");
175                return NULL;
176        }
177
178        for (i=0; i<ndata; i++) {
179                self->params.y[i] = data[i];
180        }
181
182        //self->params.y = data;
183        self->params.ny = (int)ndata;
184        return Py_BuildValue("i", self->params.ny);
185}
186
187const char get_y_doc[] =
188        "Function to get the y data\n"
189        "Takes an array of doubles as input.\n"
190        " @return: number of entries found";
191
192static PyObject * get_y(Cinvertor *self, PyObject *args) {
193        PyObject *data_obj;
194        Py_ssize_t ndata;
195        double *data;
196    int i;
197
198        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
199        OUTVECTOR(data_obj, data, ndata);
200
201        // Check that the input array is large enough
202        if (ndata < self->params.ny) {
203            PyErr_SetString(CinvertorError,
204                "Cinvertor.get_y: input array too short for data.");
205                return NULL;
206        }
207
208        for(i=0; i<self->params.ny; i++){
209                data[i] = self->params.y[i];
210        }
211
212        return Py_BuildValue("i", self->params.npoints);
213}
214
215const char set_err_doc[] =
216        "Function to set the err data\n"
217        "Takes an array of doubles as input.\n"
218        " @return: number of entries found";
219
220/**
221 * Function to set the x data
222 * Takes an array of doubles as input
223 * Returns the number of entries found
224 */
225static PyObject * set_err(Cinvertor *self, PyObject *args) {
226        PyObject *data_obj;
227        Py_ssize_t ndata;
228        double *data;
229        int i;
230
231        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
232        OUTVECTOR(data_obj,data,ndata);
233
234        free(self->params.err);
235        self->params.err = (double*) malloc(ndata*sizeof(double));
236
237        if(self->params.err==NULL) {
238            PyErr_SetString(CinvertorError,
239                "Cinvertor.set_err: problem allocating memory.");
240                return NULL;
241        }
242
243        for (i=0; i<ndata; i++) {
244                self->params.err[i] = data[i];
245        }
246
247        //self->params.err = data;
248        self->params.nerr = (int)ndata;
249        return Py_BuildValue("i", self->params.nerr);
250}
251
252const char get_err_doc[] =
253        "Function to get the err data\n"
254        "Takes an array of doubles as input.\n"
255        " @return: number of entries found";
256
257static PyObject * get_err(Cinvertor *self, PyObject *args) {
258        PyObject *data_obj;
259        Py_ssize_t ndata;
260        double *data;
261    int i;
262
263        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
264        OUTVECTOR(data_obj, data, ndata);
265
266        // Check that the input array is large enough
267        if (ndata < self->params.nerr) {
268            PyErr_SetString(CinvertorError,
269                "Cinvertor.get_err: input array too short for data.");
270                return NULL;
271        }
272
273        for(i=0; i<self->params.nerr; i++){
274                data[i] = self->params.err[i];
275        }
276
277        return Py_BuildValue("i", self->params.npoints);
278}
279
280const char is_valid_doc[] =
281        "Check the validity of the stored data\n"
282        " @return: Returns the number of points if it's all good, -1 otherwise";
283
284/**
285 * Check the validity of the stored data
286 * Returns the number of points if it's all good, -1 otherwise
287 */
288static PyObject * is_valid(Cinvertor *self, PyObject *args) {
289        if(self->params.npoints==self->params.ny &&
290                        self->params.npoints==self->params.nerr) {
291                return Py_BuildValue("i", self->params.npoints);
292        } else {
293                return Py_BuildValue("i", -1);
294        }
295}
296
297const char set_est_bck_doc[] =
298        "Sets background flag\n";
299
300/**
301 * Sets the maximum distance
302 */
303static PyObject * set_est_bck(Cinvertor *self, PyObject *args) {
304        int est_bck;
305
306        if (!PyArg_ParseTuple(args, "i", &est_bck)) return NULL;
307        self->params.est_bck = est_bck;
308        return Py_BuildValue("i", self->params.est_bck);
309}
310
311const char get_est_bck_doc[] =
312        "Gets background flag\n";
313
314/**
315 * Gets the maximum distance
316 */
317static PyObject * get_est_bck(Cinvertor *self, PyObject *args) {
318        return Py_BuildValue("i", self->params.est_bck);
319}
320
321const char set_dmax_doc[] =
322        "Sets the maximum distance\n";
323
324/**
325 * Sets the maximum distance
326 */
327static PyObject * set_dmax(Cinvertor *self, PyObject *args) {
328        double d_max;
329
330        if (!PyArg_ParseTuple(args, "d", &d_max)) return NULL;
331        self->params.d_max = d_max;
332        return Py_BuildValue("d", self->params.d_max);
333}
334
335const char get_dmax_doc[] =
336        "Gets the maximum distance\n";
337
338/**
339 * Gets the maximum distance
340 */
341static PyObject * get_dmax(Cinvertor *self, PyObject *args) {
342        return Py_BuildValue("d", self->params.d_max);
343}
344
345const char set_slit_height_doc[] =
346        "Sets the slit height in units of q [A-1]\n";
347
348/**
349 * Sets the slit height
350 */
351static PyObject * set_slit_height(Cinvertor *self, PyObject *args) {
352        double slit_height;
353
354        if (!PyArg_ParseTuple(args, "d", &slit_height)) return NULL;
355        self->params.slit_height = slit_height;
356        return Py_BuildValue("d", self->params.slit_height);
357}
358
359const char get_slit_height_doc[] =
360        "Gets the slit height\n";
361
362/**
363 * Gets the slit height
364 */
365static PyObject * get_slit_height(Cinvertor *self, PyObject *args) {
366        return Py_BuildValue("d", self->params.slit_height);
367}
368
369const char set_slit_width_doc[] =
370        "Sets the slit width in units of q [A-1]\n";
371
372/**
373 * Sets the slit width
374 */
375static PyObject * set_slit_width(Cinvertor *self, PyObject *args) {
376        double slit_width;
377
378        if (!PyArg_ParseTuple(args, "d", &slit_width)) return NULL;
379        self->params.slit_width = slit_width;
380        return Py_BuildValue("d", self->params.slit_width);
381}
382
383const char get_slit_width_doc[] =
384        "Gets the slit width\n";
385
386/**
387 * Gets the slit width
388 */
389static PyObject * get_slit_width(Cinvertor *self, PyObject *args) {
390        return Py_BuildValue("d", self->params.slit_width);
391}
392
393
394const char set_qmin_doc[] =
395        "Sets the minimum q\n";
396
397/**
398 * Sets the minimum q
399 */
400static PyObject * set_qmin(Cinvertor *self, PyObject *args) {
401        double q_min;
402
403        if (!PyArg_ParseTuple(args, "d", &q_min)) return NULL;
404        self->params.q_min = q_min;
405        return Py_BuildValue("d", self->params.q_min);
406}
407
408const char get_qmin_doc[] =
409        "Gets the minimum q\n";
410
411/**
412 * Gets the minimum q
413 */
414static PyObject * get_qmin(Cinvertor *self, PyObject *args) {
415        return Py_BuildValue("d", self->params.q_min);
416}
417
418const char set_qmax_doc[] =
419        "Sets the maximum q\n";
420
421/**
422 * Sets the maximum q
423 */
424static PyObject * set_qmax(Cinvertor *self, PyObject *args) {
425        double q_max;
426
427        if (!PyArg_ParseTuple(args, "d", &q_max)) return NULL;
428        self->params.q_max = q_max;
429        return Py_BuildValue("d", self->params.q_max);
430}
431
432const char get_qmax_doc[] =
433        "Gets the maximum q\n";
434
435/**
436 * Gets the maximum q
437 */
438static PyObject * get_qmax(Cinvertor *self, PyObject *args) {
439        return Py_BuildValue("d", self->params.q_max);
440}
441
442const char set_alpha_doc[] =
443        "Sets the alpha parameter\n";
444
445static PyObject * set_alpha(Cinvertor *self, PyObject *args) {
446        double alpha;
447
448        if (!PyArg_ParseTuple(args, "d", &alpha)) return NULL;
449        self->params.alpha = alpha;
450        return Py_BuildValue("d", self->params.alpha);
451}
452
453const char get_alpha_doc[] =
454        "Gets the alpha parameter\n";
455
456/**
457 * Gets the maximum distance
458 */
459static PyObject * get_alpha(Cinvertor *self, PyObject *args) {
460        return Py_BuildValue("d", self->params.alpha);
461}
462
463const char get_nx_doc[] =
464        "Gets the number of x points\n";
465
466/**
467 * Gets the number of x points
468 */
469static PyObject * get_nx(Cinvertor *self, PyObject *args) {
470        return Py_BuildValue("i", self->params.npoints);
471}
472
473const char get_ny_doc[] =
474        "Gets the number of y points\n";
475
476/**
477 * Gets the number of y points
478 */
479static PyObject * get_ny(Cinvertor *self, PyObject *args) {
480        return Py_BuildValue("i", self->params.ny);
481}
482
483const char get_nerr_doc[] =
484        "Gets the number of err points\n";
485
486/**
487 * Gets the number of error points
488 */
489static PyObject * get_nerr(Cinvertor *self, PyObject *args) {
490        return Py_BuildValue("i", self->params.nerr);
491}
492
493
494const char residuals_doc[] =
495        "Function to call to evaluate the residuals\n"
496        "for P(r) inversion\n"
497        " @param args: input parameters\n"
498        " @return: list of residuals";
499
500/**
501 * Function to call to evaluate the residuals
502 * @param args: input parameters
503 * @return: list of residuals
504 */
505static PyObject * residuals(Cinvertor *self, PyObject *args) {
506        double *pars;
507        PyObject* residuals;
508        int i;
509        double residual, diff;
510        // Regularization factor
511        double regterm = 0.0;
512        // Number of slices in regularization term estimate
513        int nslice = 25;
514
515        PyObject *data_obj;
516        Py_ssize_t npars;
517
518        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
519
520        OUTVECTOR(data_obj,pars,npars);
521
522    // PyList of residuals
523        // Should create this list only once and refill it
524    residuals = PyList_New(self->params.npoints);
525
526    regterm = reg_term(pars, self->params.d_max, (int)npars, nslice);
527
528    for(i=0; i<self->params.npoints; i++) {
529        diff = self->params.y[i] - iq(pars, self->params.d_max, (int)npars, self->params.x[i]);
530        residual = diff*diff / (self->params.err[i]*self->params.err[i]);
531
532        // regularization term
533        residual += self->params.alpha * regterm;
534
535        if (PyList_SetItem(residuals, i, Py_BuildValue("d",residual) ) < 0){
536            PyErr_SetString(CinvertorError,
537                "Cinvertor.residuals: error setting residual.");
538                return NULL;
539        };
540    }
541        return residuals;
542}
543
544const char pr_residuals_doc[] =
545        "Function to call to evaluate the residuals\n"
546        "for P(r) minimization (for testing purposes)\n"
547        " @param args: input parameters\n"
548        " @return: list of residuals";
549
550/**
551 * Function to call to evaluate the residuals
552 * for P(r) minimization (for testing purposes)
553 * @param args: input parameters
554 * @return: list of residuals
555 */
556static PyObject * pr_residuals(Cinvertor *self, PyObject *args) {
557        double *pars;
558        PyObject* residuals;
559        int i;
560        double residual, diff;
561        // Regularization factor
562        double regterm = 0.0;
563        // Number of slices in regularization term estimate
564        int nslice = 25;
565
566        PyObject *data_obj;
567        Py_ssize_t npars;
568
569        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
570
571        OUTVECTOR(data_obj,pars,npars);
572
573        // Should create this list only once and refill it
574    residuals = PyList_New(self->params.npoints);
575
576    regterm = reg_term(pars, self->params.d_max, (int)npars, nslice);
577
578
579    for(i=0; i<self->params.npoints; i++) {
580        diff = self->params.y[i] - pr(pars, self->params.d_max, (int)npars, self->params.x[i]);
581        residual = diff*diff / (self->params.err[i]*self->params.err[i]);
582
583        // regularization term
584        residual += self->params.alpha * regterm;
585
586        if (PyList_SetItem(residuals, i, Py_BuildValue("d",residual) ) < 0){
587            PyErr_SetString(CinvertorError,
588                "Cinvertor.residuals: error setting residual.");
589                return NULL;
590        };
591    }
592        return residuals;
593}
594
595const char get_iq_doc[] =
596        "Function to call to evaluate the scattering intensity\n"
597        " @param args: c-parameters, and q\n"
598        " @return: I(q)";
599
600/**
601 * Function to call to evaluate the scattering intensity
602 * @param args: c-parameters, and q
603 * @return: I(q)
604 */
605static PyObject * get_iq(Cinvertor *self, PyObject *args) {
606        double *pars;
607        double q, iq_value;
608        PyObject *data_obj;
609        Py_ssize_t npars;
610
611        if (!PyArg_ParseTuple(args, "Od", &data_obj, &q)) return NULL;
612        OUTVECTOR(data_obj,pars,npars);
613
614        iq_value = iq(pars, self->params.d_max, (int)npars, q);
615        return Py_BuildValue("f", iq_value);
616}
617
618const char get_iq_smeared_doc[] =
619        "Function to call to evaluate the scattering intensity.\n"
620        "The scattering intensity is slit-smeared."
621        " @param args: c-parameters, and q\n"
622        " @return: I(q)";
623
624/**
625 * Function to call to evaluate the scattering intensity
626 * The scattering intensity is slit-smeared.
627 * @param args: c-parameters, and q
628 * @return: I(q)
629 */
630static PyObject * get_iq_smeared(Cinvertor *self, PyObject *args) {
631        double *pars;
632        double q, iq_value;
633        PyObject *data_obj;
634        Py_ssize_t npars;
635
636        if (!PyArg_ParseTuple(args, "Od", &data_obj, &q)) return NULL;
637        OUTVECTOR(data_obj,pars,npars);
638
639        iq_value = iq_smeared(pars, self->params.d_max, (int)npars,
640                                                        self->params.slit_height, self->params.slit_width,
641                                                        q, 21);
642        return Py_BuildValue("f", iq_value);
643}
644
645const char get_pr_doc[] =
646        "Function to call to evaluate P(r)\n"
647        " @param args: c-parameters and r\n"
648        " @return: P(r)";
649
650/**
651 * Function to call to evaluate P(r)
652 * @param args: c-parameters and r
653 * @return: P(r)
654 */
655static PyObject * get_pr(Cinvertor *self, PyObject *args) {
656        double *pars;
657        double r, pr_value;
658        PyObject *data_obj;
659        Py_ssize_t npars;
660
661        if (!PyArg_ParseTuple(args, "Od", &data_obj, &r)) return NULL;
662        OUTVECTOR(data_obj,pars,npars);
663
664        pr_value = pr(pars, self->params.d_max, (int)npars, r);
665        return Py_BuildValue("f", pr_value);
666}
667
668const char get_pr_err_doc[] =
669        "Function to call to evaluate P(r) with errors\n"
670        " @param args: c-parameters and r\n"
671        " @return: (P(r),dP(r))";
672
673/**
674 * Function to call to evaluate P(r) with errors
675 * @param args: c-parameters and r
676 * @return: P(r)
677 */
678static PyObject * get_pr_err(Cinvertor *self, PyObject *args) {
679        double *pars;
680        double *pars_err;
681        double pr_err_value;
682        double r, pr_value;
683        PyObject *data_obj;
684        Py_ssize_t npars;
685        PyObject *err_obj;
686        Py_ssize_t npars2;
687
688        if (!PyArg_ParseTuple(args, "OOd", &data_obj, &err_obj, &r)) return NULL;
689        OUTVECTOR(data_obj,pars,npars);
690
691        if (err_obj == Py_None) {
692                pr_value = pr(pars, self->params.d_max, (int)npars, r);
693                pr_err_value = 0.0;
694        } else {
695                OUTVECTOR(err_obj,pars_err,npars2);
696                pr_err(pars, pars_err, self->params.d_max, (int)npars, r, &pr_value, &pr_err_value);
697        }
698        return Py_BuildValue("ff", pr_value, pr_err_value);
699}
700
701const char basefunc_ft_doc[] =
702        "Returns the value of the nth Fourier transofrmed base function\n"
703        " @param args: c-parameters, n and q\n"
704        " @return: nth Fourier transformed base function, evaluated at q";
705
706static PyObject * basefunc_ft(Cinvertor *self, PyObject *args) {
707        double d_max, q;
708        int n;
709
710        if (!PyArg_ParseTuple(args, "did", &d_max, &n, &q)) return NULL;
711        return Py_BuildValue("f", ortho_transformed(d_max, n, q));
712
713}
714
715const char oscillations_doc[] =
716        "Returns the value of the oscillation figure of merit for\n"
717        "the given set of coefficients. For a sphere, the oscillation\n"
718        "figure of merit is 1.1.\n"
719        " @param args: c-parameters\n"
720        " @return: oscillation figure of merit";
721
722static PyObject * oscillations(Cinvertor *self, PyObject *args) {
723        double *pars;
724        PyObject *data_obj;
725        Py_ssize_t npars;
726        double oscill, norm;
727
728        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
729        OUTVECTOR(data_obj,pars,npars);
730
731        oscill = reg_term(pars, self->params.d_max, (int)npars, 100);
732        norm   = int_p2(pars, self->params.d_max, (int)npars, 100);
733        return Py_BuildValue("f", sqrt(oscill/norm)/acos(-1.0)*self->params.d_max );
734
735}
736
737const char get_peaks_doc[] =
738        "Returns the number of peaks in the output P(r) distribution\n"
739        "for the given set of coefficients.\n"
740        " @param args: c-parameters\n"
741        " @return: number of P(r) peaks";
742
743static PyObject * get_peaks(Cinvertor *self, PyObject *args) {
744        double *pars;
745        PyObject *data_obj;
746        Py_ssize_t npars;
747        int count;
748
749        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
750        OUTVECTOR(data_obj,pars,npars);
751
752        count = npeaks(pars, self->params.d_max, (int)npars, 100);
753
754        return Py_BuildValue("i", count );
755
756}
757
758const char get_positive_doc[] =
759        "Returns the fraction of P(r) that is positive over\n"
760        "the full range of r for the given set of coefficients.\n"
761        " @param args: c-parameters\n"
762        " @return: fraction of P(r) that is positive";
763
764static PyObject * get_positive(Cinvertor *self, PyObject *args) {
765        double *pars;
766        PyObject *data_obj;
767        Py_ssize_t npars;
768        double fraction;
769
770        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
771        OUTVECTOR(data_obj,pars,npars);
772
773        fraction = positive_integral(pars, self->params.d_max, (int)npars, 100);
774
775        return Py_BuildValue("f", fraction );
776
777}
778
779const char get_pos_err_doc[] =
780        "Returns the fraction of P(r) that is 1 standard deviation\n"
781        "above zero over the full range of r for the given set of coefficients.\n"
782        " @param args: c-parameters\n"
783        " @return: fraction of P(r) that is positive";
784
785static PyObject * get_pos_err(Cinvertor *self, PyObject *args) {
786        double *pars;
787        double *pars_err;
788        PyObject *data_obj;
789        PyObject *err_obj;
790        Py_ssize_t npars;
791        Py_ssize_t npars2;
792        double fraction;
793
794        if (!PyArg_ParseTuple(args, "OO", &data_obj, &err_obj)) return NULL;
795        OUTVECTOR(data_obj,pars,npars);
796        OUTVECTOR(err_obj,pars_err,npars2);
797
798        fraction = positive_errors(pars, pars_err, self->params.d_max, (int)npars, 51);
799
800        return Py_BuildValue("f", fraction );
801
802}
803
804const char get_rg_doc[] =
805        "Returns the value of the radius of gyration Rg.\n"
806        " @param args: c-parameters\n"
807        " @return: Rg";
808
809static PyObject * get_rg(Cinvertor *self, PyObject *args) {
810        double *pars;
811        PyObject *data_obj;
812        Py_ssize_t npars;
813        double value;
814
815        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
816        OUTVECTOR(data_obj,pars,npars);
817
818        value = rg(pars, self->params.d_max, (int)npars, 101);
819
820        return Py_BuildValue("f", value );
821
822}
823
824const char get_iq0_doc[] =
825        "Returns the value of I(q=0).\n"
826        " @param args: c-parameters\n"
827        " @return: I(q=0)";
828
829static PyObject * get_iq0(Cinvertor *self, PyObject *args) {
830        double *pars;
831        PyObject *data_obj;
832        Py_ssize_t npars;
833        double value;
834
835        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL;
836        OUTVECTOR(data_obj,pars,npars);
837
838        value = 4.0*acos(-1.0)*int_pr(pars, self->params.d_max, (int)npars, 101);
839
840        return Py_BuildValue("f", value );
841
842}
843
844/**
845 * Check whether a q-value is within acceptabel limits
846 * Return 1 if accepted, 0 if rejected.
847 */
848int accept_q(Cinvertor *self, double q) {
849    if (self->params.q_min>0 && q<self->params.q_min) return 0;
850    if (self->params.q_max>0 && q>self->params.q_max) return 0;
851    return 1;
852}
853
854const char get_matrix_doc[] =
855        "Returns A matrix and b vector for least square problem.\n"
856        " @param nfunc: number of base functions\n"
857        " @param nr: number of r-points used when evaluating reg term.\n"
858        " @param a: A array to fill\n"
859        " @param b: b vector to fill\n"
860        " @return: 0";
861
862static PyObject * get_matrix(Cinvertor *self, PyObject *args) {
863        double *a;
864        double *b;
865        PyObject *a_obj;
866        PyObject *b_obj;
867        Py_ssize_t n_a;
868        Py_ssize_t n_b;
869        // Number of bins for regularization term evaluation
870        int nr, nfunc;
871        int i, j, i_r;
872        double r, sqrt_alpha, pi;
873        double tmp;
874        int offset;
875
876        if (!PyArg_ParseTuple(args, "iiOO", &nfunc, &nr, &a_obj, &b_obj)) return NULL;
877        OUTVECTOR(a_obj,a,n_a);
878        OUTVECTOR(b_obj,b,n_b);
879
880        assert(n_b>=nfunc);
881        assert(n_a>=nfunc*(nr+self->params.npoints));
882
883        sqrt_alpha = sqrt(self->params.alpha);
884        pi = acos(-1.0);
885        offset = (self->params.est_bck==1) ? 0 : 1;
886
887    for (j=0; j<nfunc; j++) {
888        for (i=0; i<self->params.npoints; i++) {
889            if (self->params.err[i]==0.0) {
890              PyErr_SetString(CinvertorError,
891                "Cinvertor.get_matrix: Some I(Q) points have no error.");
892              return NULL;
893            }
894            if (accept_q(self, self->params.x[i])){
895                if (self->params.est_bck==1 && j==0) {
896                    a[i*nfunc+j] = 1.0/self->params.err[i];
897                } else {
898                        if (self->params.slit_width>0 || self->params.slit_height>0) {
899                                a[i*nfunc+j] = ortho_transformed_smeared(self->params.d_max,
900                                                j+offset, self->params.slit_height, self->params.slit_width,
901                                                self->params.x[i], 21)/self->params.err[i];
902                        } else {
903                                a[i*nfunc+j] = ortho_transformed(self->params.d_max, j+offset, self->params.x[i])/self->params.err[i];
904                        }
905                }
906            }
907        }
908        for (i_r=0; i_r<nr; i_r++){
909            if (self->params.est_bck==1 && j==0) {
910                a[(i_r+self->params.npoints)*nfunc+j] = 0.0;
911            } else {
912                    r = self->params.d_max/nr*i_r;
913                    tmp = pi*(j+offset)/self->params.d_max;
914                    a[(i_r+self->params.npoints)*nfunc+j] = sqrt_alpha * 1.0/nr*self->params.d_max*2.0*
915                    (2.0*pi*(j+offset)/self->params.d_max*cos(pi*(j+offset)*r/self->params.d_max) +
916                    tmp*tmp*r * sin(pi*(j+offset)*r/self->params.d_max));
917                }
918        }
919    }
920
921    for (i=0; i<self->params.npoints; i++) {
922        if (accept_q(self, self->params.x[i])){
923            b[i] = self->params.y[i]/self->params.err[i];
924         }
925        }
926
927        return Py_BuildValue("i", 0);
928
929}
930
931const char get_invcov_matrix_doc[] =
932        " Compute the inverse covariance matrix, defined as inv_cov = a_transposed x a.\n"
933        " @param nfunc: number of base functions\n"
934        " @param nr: number of r-points used when evaluating reg term.\n"
935        " @param a: A array to fill\n"
936        " @param inv_cov: inverse covariance array to be filled\n"
937        " @return: 0";
938
939static PyObject * get_invcov_matrix(Cinvertor *self, PyObject *args) {
940        double *a;
941        PyObject *a_obj;
942        Py_ssize_t n_a;
943        double *inv_cov;
944        PyObject *cov_obj;
945        Py_ssize_t n_cov;
946        int nr, nfunc;
947        int i, j, k;
948
949        if (!PyArg_ParseTuple(args, "iiOO", &nfunc, &nr, &a_obj, &cov_obj)) return NULL;
950        OUTVECTOR(a_obj,a,n_a);
951        OUTVECTOR(cov_obj,inv_cov,n_cov);
952
953        assert(n_cov>=nfunc*nfunc);
954        assert(n_a>=nfunc*(nr+self->params.npoints));
955
956        for (i=0; i<nfunc; i++) {
957                for (j=0; j<nfunc; j++) {
958                        inv_cov[i*nfunc+j] = 0.0;
959                        for (k=0; k<nr+self->params.npoints; k++) {
960                                inv_cov[i*nfunc+j] += a[k*nfunc+i]*a[k*nfunc+j];
961                        }
962                }
963        }
964        return Py_BuildValue("i", 0);
965}
966
967const char get_reg_size_doc[] =
968        " Compute the covariance matrix, defined as inv_cov = a_transposed x a.\n"
969        " @param nfunc: number of base functions\n"
970        " @param nr: number of r-points used when evaluating reg term.\n"
971        " @param a: A array to fill\n"
972        " @param inv_cov: inverse covariance array to be filled\n"
973        " @return: 0";
974
975static PyObject * get_reg_size(Cinvertor *self, PyObject *args) {
976        double *a;
977        PyObject *a_obj;
978        Py_ssize_t n_a;
979        int nr, nfunc;
980        int i, j;
981        double sum_sig, sum_reg;
982
983        if (!PyArg_ParseTuple(args, "iiO", &nfunc, &nr, &a_obj)) return NULL;
984        OUTVECTOR(a_obj,a,n_a);
985
986        assert(n_a>=nfunc*(nr+self->params.npoints));
987
988    sum_sig = 0.0;
989    sum_reg = 0.0;
990    for (j=0; j<nfunc; j++){
991        for (i=0; i<self->params.npoints; i++){
992            if (accept_q(self, self->params.x[i])==1)
993                sum_sig += (a[i*nfunc+j])*(a[i*nfunc+j]);
994        }
995        for (i=0; i<nr; i++){
996            sum_reg += (a[(i+self->params.npoints)*nfunc+j])*(a[(i+self->params.npoints)*nfunc+j]);
997        }
998    }
999    return Py_BuildValue("ff", sum_sig, sum_reg);
1000}
1001
1002const char eeeget_qmin_doc[] = "\
1003This is a multiline doc string.\n\
1004\n\
1005This is the second line.";
1006const char eeeset_qmin_doc[] =
1007        "This is a multiline doc string.\n"
1008        "\n"
1009        "This is the second line.";
1010
1011static PyMethodDef Cinvertor_methods[] = {
1012                   {"residuals", (PyCFunction)residuals, METH_VARARGS, residuals_doc},
1013                   {"pr_residuals", (PyCFunction)pr_residuals, METH_VARARGS, pr_residuals_doc},
1014                   {"set_x", (PyCFunction)set_x, METH_VARARGS, set_x_doc},
1015                   {"get_x", (PyCFunction)get_x, METH_VARARGS, get_x_doc},
1016                   {"set_y", (PyCFunction)set_y, METH_VARARGS, set_y_doc},
1017                   {"get_y", (PyCFunction)get_y, METH_VARARGS, get_y_doc},
1018                   {"set_err", (PyCFunction)set_err, METH_VARARGS, set_err_doc},
1019                   {"get_err", (PyCFunction)get_err, METH_VARARGS, get_err_doc},
1020                   {"set_dmax", (PyCFunction)set_dmax, METH_VARARGS, set_dmax_doc},
1021                   {"get_dmax", (PyCFunction)get_dmax, METH_VARARGS, get_dmax_doc},
1022                   {"set_qmin", (PyCFunction)set_qmin, METH_VARARGS, set_qmin_doc},
1023                   {"get_qmin", (PyCFunction)get_qmin, METH_VARARGS, get_qmin_doc},
1024                   {"set_qmax", (PyCFunction)set_qmax, METH_VARARGS, set_qmax_doc},
1025                   {"get_qmax", (PyCFunction)get_qmax, METH_VARARGS, get_qmax_doc},
1026                   {"set_alpha", (PyCFunction)set_alpha, METH_VARARGS, set_alpha_doc},
1027                   {"get_alpha", (PyCFunction)get_alpha, METH_VARARGS, get_alpha_doc},
1028                   {"set_slit_width", (PyCFunction)set_slit_width, METH_VARARGS, set_slit_width_doc},
1029                   {"get_slit_width", (PyCFunction)get_slit_width, METH_VARARGS, get_slit_width_doc},
1030                   {"set_slit_height", (PyCFunction)set_slit_height, METH_VARARGS, set_slit_height_doc},
1031                   {"get_slit_height", (PyCFunction)get_slit_height, METH_VARARGS, get_slit_height_doc},
1032                   {"set_est_bck", (PyCFunction)set_est_bck, METH_VARARGS, set_est_bck_doc},
1033                   {"get_est_bck", (PyCFunction)get_est_bck, METH_VARARGS, get_est_bck_doc},
1034                   {"get_nx", (PyCFunction)get_nx, METH_VARARGS, get_nx_doc},
1035                   {"get_ny", (PyCFunction)get_ny, METH_VARARGS, get_ny_doc},
1036                   {"get_nerr", (PyCFunction)get_nerr, METH_VARARGS, get_nerr_doc},
1037                   {"iq", (PyCFunction)get_iq, METH_VARARGS, get_iq_doc},
1038                   {"iq_smeared", (PyCFunction)get_iq_smeared, METH_VARARGS, get_iq_smeared_doc},
1039                   {"pr", (PyCFunction)get_pr, METH_VARARGS, get_pr_doc},
1040                   {"get_pr_err", (PyCFunction)get_pr_err, METH_VARARGS, get_pr_err_doc},
1041                   {"is_valid", (PyCFunction)is_valid, METH_VARARGS, is_valid_doc},
1042                   {"basefunc_ft", (PyCFunction)basefunc_ft, METH_VARARGS, basefunc_ft_doc},
1043                   {"oscillations", (PyCFunction)oscillations, METH_VARARGS, oscillations_doc},
1044                   {"get_peaks", (PyCFunction)get_peaks, METH_VARARGS, get_peaks_doc},
1045                   {"get_positive", (PyCFunction)get_positive, METH_VARARGS, get_positive_doc},
1046                   {"get_pos_err", (PyCFunction)get_pos_err, METH_VARARGS, get_pos_err_doc},
1047                   {"rg", (PyCFunction)get_rg, METH_VARARGS, get_rg_doc},
1048                   {"iq0", (PyCFunction)get_iq0, METH_VARARGS, get_iq0_doc},
1049                   {"_get_matrix", (PyCFunction)get_matrix, METH_VARARGS, get_matrix_doc},
1050                   {"_get_invcov_matrix", (PyCFunction)get_invcov_matrix, METH_VARARGS, get_invcov_matrix_doc},
1051                   {"_get_reg_size", (PyCFunction)get_reg_size, METH_VARARGS, get_reg_size_doc},
1052
1053   {NULL}
1054};
1055
1056static PyTypeObject CinvertorType = {
1057    //PyObject_HEAD_INIT(NULL)
1058    //0,                         /*ob_size*/
1059    PyVarObject_HEAD_INIT(NULL, 0)
1060    "Cinvertor",             /*tp_name*/
1061    sizeof(Cinvertor),             /*tp_basicsize*/
1062    0,                         /*tp_itemsize*/
1063    (destructor)Cinvertor_dealloc, /*tp_dealloc*/
1064    0,                         /*tp_print*/
1065    0,                         /*tp_getattr*/
1066    0,                         /*tp_setattr*/
1067    0,                         /*tp_compare*/
1068    0,                         /*tp_repr*/
1069    0,                         /*tp_as_number*/
1070    0,                         /*tp_as_sequence*/
1071    0,                         /*tp_as_mapping*/
1072    0,                         /*tp_hash */
1073    0,                         /*tp_call*/
1074    0,                         /*tp_str*/
1075    0,                         /*tp_getattro*/
1076    0,                         /*tp_setattro*/
1077    0,                         /*tp_as_buffer*/
1078    Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, /*tp_flags*/
1079    "Cinvertor objects",           /* tp_doc */
1080    0,                         /* tp_traverse */
1081    0,                         /* tp_clear */
1082    0,                         /* tp_richcompare */
1083    0,                         /* tp_weaklistoffset */
1084    0,                         /* tp_iter */
1085    0,                         /* tp_iternext */
1086    Cinvertor_methods,             /* tp_methods */
1087    Cinvertor_members,             /* tp_members */
1088    0,                         /* tp_getset */
1089    0,                         /* tp_base */
1090    0,                         /* tp_dict */
1091    0,                         /* tp_descr_get */
1092    0,                         /* tp_descr_set */
1093    0,                         /* tp_dictoffset */
1094    (initproc)Cinvertor_init,      /* tp_init */
1095    0,                         /* tp_alloc */
1096    Cinvertor_new,                 /* tp_new */
1097};
1098
1099
1100static PyMethodDef module_methods[] = {
1101    {NULL}
1102};
1103
1104/**
1105 * Function used to add the model class to a module
1106 * @param module: module to add the class to
1107 */
1108void addCinvertor(PyObject *module) {
1109        PyObject *d;
1110
1111    if (PyType_Ready(&CinvertorType) < 0)
1112        return;
1113
1114    Py_INCREF(&CinvertorType);
1115    PyModule_AddObject(module, "Cinvertor", (PyObject *)&CinvertorType);
1116
1117    d = PyModule_GetDict(module);
1118    CinvertorError = PyErr_NewException("sas.sascalc.pr.invertor.Cinvertor.InvertorError", PyExc_RuntimeError, NULL);
1119    PyDict_SetItemString(d, "CinvertorError", CinvertorError);
1120}
1121
1122
1123#define MODULE_DOC "C extension module for inversion to P(r)."
1124#define MODULE_NAME "_pr_inversion"
1125#define MODULE_INIT2 init_pr_inversion
1126#define MODULE_INIT3 PyInit__pr_inversion
1127#define MODULE_METHODS module_methods
1128
1129/* ==== boilerplate python 2/3 interface bootstrap ==== */
1130
1131
1132#if defined(WIN32) && !defined(__MINGW32__)
1133    #define DLL_EXPORT __declspec(dllexport)
1134#else
1135    #define DLL_EXPORT
1136#endif
1137
1138#if PY_MAJOR_VERSION >= 3
1139
1140  DLL_EXPORT PyMODINIT_FUNC MODULE_INIT3(void)
1141  {
1142    static struct PyModuleDef moduledef = {
1143      PyModuleDef_HEAD_INIT,
1144      MODULE_NAME,         /* m_name */
1145      MODULE_DOC,          /* m_doc */
1146      -1,                  /* m_size */
1147      MODULE_METHODS,      /* m_methods */
1148      NULL,                /* m_reload */
1149      NULL,                /* m_traverse */
1150      NULL,                /* m_clear */
1151      NULL,                /* m_free */
1152    };
1153        PyObject* m = PyModule_Create(&moduledef);
1154        addCinvertor(m);
1155        return m;
1156  }
1157
1158#else /* !PY_MAJOR_VERSION >= 3 */
1159
1160  DLL_EXPORT PyMODINIT_FUNC MODULE_INIT2(void)
1161  {
1162    PyObject* m = Py_InitModule4(MODULE_NAME,
1163                 MODULE_METHODS,
1164                 MODULE_DOC,
1165                 0,
1166                 PYTHON_API_VERSION
1167                 );
1168        addCinvertor(m);
1169  }
1170
1171#endif /* !PY_MAJOR_VERSION >= 3 */
Note: See TracBrowser for help on using the repository browser.