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

magnetic_scattrelease-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since a52f32f was a52f32f, checked in by Stuart Prescott <stuart@…>, 20 months ago

Fix typo in P(r) peaks description

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