source: sasview/src/sas/sascalc/pr/c_extensions/Cinvertor.c @ 7ba6470

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

remove deprecation warnings for old buffer protocol from py37 build

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