Changeset 9a23253e in sasview for pr_inversion/c_extensions


Ignore:
Timestamp:
Jun 30, 2008 2:18:00 PM (16 years ago)
Author:
Mathieu Doucet <doucetm@…>
Branches:
master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
Children:
a4bd2ac
Parents:
7bb61da
Message:

Started to add slit smearing. Added Rg and I(q=0) as outputs.

Location:
pr_inversion/c_extensions
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • pr_inversion/c_extensions/Cinvertor.c

    r896abb3 r9a23253e  
    33 * Cinvertor is the base class for the Invertor class 
    44 * and provides the underlying computations. 
    5  *  
     5 * 
    66 */ 
    77#include <Python.h> 
     
    3838 * Cinvertor is the base class for the Invertor class 
    3939 * and provides the underlying computations. 
    40  *  
     40 * 
    4141 */ 
    4242typedef struct { 
    43     PyObject_HEAD     
     43    PyObject_HEAD 
    4444    /// Internal data structure 
    45     Invertor_params params;  
     45    Invertor_params params; 
    4646} Cinvertor; 
    4747 
     
    5151{ 
    5252    invertor_dealloc(&(self->params)); 
    53       
     53 
    5454    self->ob_type->tp_free((PyObject*)self); 
    5555 
     
    6060{ 
    6161    Cinvertor *self; 
    62      
     62 
    6363    self = (Cinvertor *)type->tp_alloc(type, 0); 
    64     
     64 
    6565    return (PyObject *)self; 
    6666} 
     
    6969Cinvertor_init(Cinvertor *self, PyObject *args, PyObject *kwds) 
    7070{ 
    71     if (self != NULL) {          
     71    if (self != NULL) { 
    7272        // Create parameters 
    7373        invertor_init(&(self->params)); 
     
    8282}; 
    8383 
    84 const char set_x_doc[] =  
     84const char set_x_doc[] = 
    8585        "Function to set the x data\n" 
    8686        "Takes an array of doubles as input.\n" 
     
    9797        double *data; 
    9898        int i; 
    99    
     99 
    100100        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 
    101101        OUTVECTOR(data_obj,data,ndata); 
    102          
     102 
    103103        free(self->params.x); 
    104104        self->params.x = (double*) malloc(ndata*sizeof(double)); 
    105          
     105 
    106106        if(self->params.x==NULL) { 
    107             PyErr_SetString(CinvertorError,  
     107            PyErr_SetString(CinvertorError, 
    108108                "Cinvertor.set_x: problem allocating memory."); 
    109                 return NULL;             
    110         } 
    111          
     109                return NULL; 
     110        } 
     111 
    112112        for (i=0; i<ndata; i++) { 
    113113                self->params.x[i] = data[i]; 
    114114        } 
    115          
     115 
    116116        //self->params.x = data; 
    117117        self->params.npoints = ndata; 
    118         return Py_BuildValue("i", self->params.npoints);         
    119 } 
    120  
    121 const char get_x_doc[] =  
     118        return Py_BuildValue("i", self->params.npoints); 
     119} 
     120 
     121const char get_x_doc[] = 
    122122        "Function to get the x data\n" 
    123123        "Takes an array of doubles as input.\n" 
     
    129129        double *data; 
    130130    int i; 
    131      
     131 
    132132        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 
    133133        OUTVECTOR(data_obj, data, ndata); 
    134          
     134 
    135135        // Check that the input array is large enough 
    136136        if (ndata < self->params.npoints) { 
    137             PyErr_SetString(CinvertorError,  
     137            PyErr_SetString(CinvertorError, 
    138138                "Cinvertor.get_x: input array too short for data."); 
    139                 return NULL;             
    140         } 
    141          
     139                return NULL; 
     140        } 
     141 
    142142        for(i=0; i<self->params.npoints; i++){ 
    143143                data[i] = self->params.x[i]; 
    144144        } 
    145          
    146         return Py_BuildValue("i", self->params.npoints);         
    147 } 
    148  
    149 const char set_y_doc[] =  
     145 
     146        return Py_BuildValue("i", self->params.npoints); 
     147} 
     148 
     149const char set_y_doc[] = 
    150150        "Function to set the y data\n" 
    151151        "Takes an array of doubles as input.\n" 
     
    162162        double *data; 
    163163        int i; 
    164    
     164 
    165165        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 
    166166        OUTVECTOR(data_obj,data,ndata); 
    167          
     167 
    168168        free(self->params.y); 
    169169        self->params.y = (double*) malloc(ndata*sizeof(double)); 
    170          
     170 
    171171        if(self->params.y==NULL) { 
    172             PyErr_SetString(CinvertorError,  
     172            PyErr_SetString(CinvertorError, 
    173173                "Cinvertor.set_y: problem allocating memory."); 
    174                 return NULL;             
    175         } 
    176          
     174                return NULL; 
     175        } 
     176 
    177177        for (i=0; i<ndata; i++) { 
    178178                self->params.y[i] = data[i]; 
    179         }        
    180          
     179        } 
     180 
    181181        //self->params.y = data; 
    182182        self->params.ny = ndata; 
    183         return Py_BuildValue("i", self->params.ny);      
    184 } 
    185  
    186 const char get_y_doc[] =  
     183        return Py_BuildValue("i", self->params.ny); 
     184} 
     185 
     186const char get_y_doc[] = 
    187187        "Function to get the y data\n" 
    188188        "Takes an array of doubles as input.\n" 
     
    194194        double *data; 
    195195    int i; 
    196      
     196 
    197197        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 
    198198        OUTVECTOR(data_obj, data, ndata); 
    199          
     199 
    200200        // Check that the input array is large enough 
    201201        if (ndata < self->params.ny) { 
    202             PyErr_SetString(CinvertorError,  
     202            PyErr_SetString(CinvertorError, 
    203203                "Cinvertor.get_y: input array too short for data."); 
    204                 return NULL;             
    205         } 
    206          
     204                return NULL; 
     205        } 
     206 
    207207        for(i=0; i<self->params.ny; i++){ 
    208208                data[i] = self->params.y[i]; 
    209209        } 
    210          
    211         return Py_BuildValue("i", self->params.npoints);         
    212 } 
    213  
    214 const char set_err_doc[] =  
     210 
     211        return Py_BuildValue("i", self->params.npoints); 
     212} 
     213 
     214const char set_err_doc[] = 
    215215        "Function to set the err data\n" 
    216216        "Takes an array of doubles as input.\n" 
     
    227227        double *data; 
    228228        int i; 
    229    
     229 
    230230        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 
    231231        OUTVECTOR(data_obj,data,ndata); 
    232          
     232 
    233233        free(self->params.err); 
    234234        self->params.err = (double*) malloc(ndata*sizeof(double)); 
    235          
     235 
    236236        if(self->params.err==NULL) { 
    237             PyErr_SetString(CinvertorError,  
     237            PyErr_SetString(CinvertorError, 
    238238                "Cinvertor.set_err: problem allocating memory."); 
    239                 return NULL;             
    240         } 
    241          
     239                return NULL; 
     240        } 
     241 
    242242        for (i=0; i<ndata; i++) { 
    243243                self->params.err[i] = data[i]; 
    244244        } 
    245          
     245 
    246246        //self->params.err = data; 
    247247        self->params.nerr = ndata; 
    248         return Py_BuildValue("i", self->params.nerr);    
    249 } 
    250  
    251 const char get_err_doc[] =  
     248        return Py_BuildValue("i", self->params.nerr); 
     249} 
     250 
     251const char get_err_doc[] = 
    252252        "Function to get the err data\n" 
    253253        "Takes an array of doubles as input.\n" 
     
    259259        double *data; 
    260260    int i; 
    261      
     261 
    262262        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 
    263263        OUTVECTOR(data_obj, data, ndata); 
    264          
     264 
    265265        // Check that the input array is large enough 
    266266        if (ndata < self->params.nerr) { 
    267             PyErr_SetString(CinvertorError,  
     267            PyErr_SetString(CinvertorError, 
    268268                "Cinvertor.get_err: input array too short for data."); 
    269                 return NULL;             
    270         } 
    271          
     269                return NULL; 
     270        } 
     271 
    272272        for(i=0; i<self->params.nerr; i++){ 
    273273                data[i] = self->params.err[i]; 
    274274        } 
    275          
    276         return Py_BuildValue("i", self->params.npoints);         
    277 } 
    278  
    279 const char is_valid_doc[] =  
     275 
     276        return Py_BuildValue("i", self->params.npoints); 
     277} 
     278 
     279const char is_valid_doc[] = 
    280280        "Check the validity of the stored data\n" 
    281281        " @return: Returns the number of points if it's all good, -1 otherwise"; 
     
    291291        } else { 
    292292                return Py_BuildValue("i", -1); 
    293         }        
    294 } 
    295  
    296 const char set_dmax_doc[] =  
     293        } 
     294} 
     295 
     296const char set_has_bck_doc[] = 
     297        "Sets background flag\n"; 
     298 
     299/** 
     300 * Sets the maximum distance 
     301 */ 
     302static PyObject * set_has_bck(Cinvertor *self, PyObject *args) { 
     303        int has_bck; 
     304 
     305        if (!PyArg_ParseTuple(args, "i", &has_bck)) return NULL; 
     306        self->params.has_bck = has_bck; 
     307        return Py_BuildValue("i", self->params.has_bck); 
     308} 
     309 
     310const char get_has_bck_doc[] = 
     311        "Gets background flag\n"; 
     312 
     313/** 
     314 * Gets the maximum distance 
     315 */ 
     316static PyObject * get_has_bck(Cinvertor *self, PyObject *args) { 
     317        return Py_BuildValue("i", self->params.has_bck); 
     318} 
     319 
     320const char set_dmax_doc[] = 
    297321        "Sets the maximum distance\n"; 
    298322 
     
    302326static PyObject * set_dmax(Cinvertor *self, PyObject *args) { 
    303327        double d_max; 
    304    
     328 
    305329        if (!PyArg_ParseTuple(args, "d", &d_max)) return NULL; 
    306330        self->params.d_max = d_max; 
    307         return Py_BuildValue("d", self->params.d_max);   
    308 } 
    309  
    310 const char get_dmax_doc[] =  
     331        return Py_BuildValue("d", self->params.d_max); 
     332} 
     333 
     334const char get_dmax_doc[] = 
    311335        "Gets the maximum distance\n"; 
    312336 
     
    315339 */ 
    316340static PyObject * get_dmax(Cinvertor *self, PyObject *args) { 
    317         return Py_BuildValue("d", self->params.d_max);   
    318 } 
    319  
    320 const char set_qmin_doc[] =  
     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[] = 
    321394        "Sets the minimum q\n"; 
    322395 
     
    326399static PyObject * set_qmin(Cinvertor *self, PyObject *args) { 
    327400        double q_min; 
    328    
     401 
    329402        if (!PyArg_ParseTuple(args, "d", &q_min)) return NULL; 
    330403        self->params.q_min = q_min; 
    331         return Py_BuildValue("d", self->params.q_min);   
    332 } 
    333  
    334 const char get_qmin_doc[] =  
     404        return Py_BuildValue("d", self->params.q_min); 
     405} 
     406 
     407const char get_qmin_doc[] = 
    335408        "Gets the minimum q\n"; 
    336409 
     
    339412 */ 
    340413static PyObject * get_qmin(Cinvertor *self, PyObject *args) { 
    341         return Py_BuildValue("d", self->params.q_min);   
    342 } 
    343  
    344 const char set_qmax_doc[] =  
     414        return Py_BuildValue("d", self->params.q_min); 
     415} 
     416 
     417const char set_qmax_doc[] = 
    345418        "Sets the maximum q\n"; 
    346419 
     
    350423static PyObject * set_qmax(Cinvertor *self, PyObject *args) { 
    351424        double q_max; 
    352    
     425 
    353426        if (!PyArg_ParseTuple(args, "d", &q_max)) return NULL; 
    354427        self->params.q_max = q_max; 
    355         return Py_BuildValue("d", self->params.q_max);   
    356 } 
    357  
    358 const char get_qmax_doc[] =  
     428        return Py_BuildValue("d", self->params.q_max); 
     429} 
     430 
     431const char get_qmax_doc[] = 
    359432        "Gets the maximum q\n"; 
    360433 
     
    363436 */ 
    364437static PyObject * get_qmax(Cinvertor *self, PyObject *args) { 
    365         return Py_BuildValue("d", self->params.q_max);   
    366 } 
    367  
    368 const char set_alpha_doc[] =  
     438        return Py_BuildValue("d", self->params.q_max); 
     439} 
     440 
     441const char set_alpha_doc[] = 
    369442        "Sets the alpha parameter\n"; 
    370443 
    371444static PyObject * set_alpha(Cinvertor *self, PyObject *args) { 
    372445        double alpha; 
    373    
     446 
    374447        if (!PyArg_ParseTuple(args, "d", &alpha)) return NULL; 
    375448        self->params.alpha = alpha; 
    376         return Py_BuildValue("d", self->params.alpha);   
    377 } 
    378  
    379 const char get_alpha_doc[] =  
     449        return Py_BuildValue("d", self->params.alpha); 
     450} 
     451 
     452const char get_alpha_doc[] = 
    380453        "Gets the alpha parameter\n"; 
    381454 
     
    384457 */ 
    385458static PyObject * get_alpha(Cinvertor *self, PyObject *args) { 
    386         return Py_BuildValue("d", self->params.alpha);   
    387 } 
    388  
    389 const char get_nx_doc[] =  
     459        return Py_BuildValue("d", self->params.alpha); 
     460} 
     461 
     462const char get_nx_doc[] = 
    390463        "Gets the number of x points\n"; 
    391464 
     
    394467 */ 
    395468static PyObject * get_nx(Cinvertor *self, PyObject *args) { 
    396         return Py_BuildValue("i", self->params.npoints);         
    397 } 
    398  
    399 const char get_ny_doc[] =  
     469        return Py_BuildValue("i", self->params.npoints); 
     470} 
     471 
     472const char get_ny_doc[] = 
    400473        "Gets the number of y points\n"; 
    401474 
     
    404477 */ 
    405478static PyObject * get_ny(Cinvertor *self, PyObject *args) { 
    406         return Py_BuildValue("i", self->params.ny);      
    407 } 
    408  
    409 const char get_nerr_doc[] =  
     479        return Py_BuildValue("i", self->params.ny); 
     480} 
     481 
     482const char get_nerr_doc[] = 
    410483        "Gets the number of err points\n"; 
    411484 
     
    414487 */ 
    415488static PyObject * get_nerr(Cinvertor *self, PyObject *args) { 
    416         return Py_BuildValue("i", self->params.nerr);    
    417 } 
    418  
    419  
    420 const char residuals_doc[] =  
     489        return Py_BuildValue("i", self->params.nerr); 
     490} 
     491 
     492 
     493const char residuals_doc[] = 
    421494        "Function to call to evaluate the residuals\n" 
    422495        "for P(r) inversion\n" 
     
    441514        // Number of slices in regularization term estimate 
    442515        int nslice = 25; 
    443          
     516 
    444517        PyObject *data_obj; 
    445518        Py_ssize_t npars; 
    446            
    447         if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 
    448          
     519 
     520        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 
     521 
    449522        OUTVECTOR(data_obj,pars,npars); 
    450                  
     523 
    451524    // PyList of residuals 
    452525        // Should create this list only once and refill it 
     
    454527 
    455528    regterm = reg_term(pars, self->params.d_max, npars, nslice); 
    456      
     529 
    457530    for(i=0; i<self->params.npoints; i++) { 
    458531        diff = self->params.y[i] - iq(pars, self->params.d_max, npars, self->params.x[i]); 
    459532        residual = diff*diff / (self->params.err[i]*self->params.err[i]); 
    460533        tmp = residual; 
    461          
     534 
    462535        // regularization term 
    463536        residual += self->params.alpha * regterm; 
    464          
     537 
    465538        if (PyList_SetItem(residuals, i, Py_BuildValue("d",residual) ) < 0){ 
    466             PyErr_SetString(CinvertorError,  
     539            PyErr_SetString(CinvertorError, 
    467540                "Cinvertor.residuals: error setting residual."); 
    468541                return NULL; 
    469542        }; 
    470                  
     543 
    471544    } 
    472      
     545 
    473546        return residuals; 
    474547} 
    475548 
    476 const char pr_residuals_doc[] =  
     549const char pr_residuals_doc[] = 
    477550        "Function to call to evaluate the residuals\n" 
    478551        "for P(r) minimization (for testing purposes)\n" 
     
    498571        // Number of slices in regularization term estimate 
    499572        int nslice = 25; 
    500          
     573 
    501574        PyObject *data_obj; 
    502575        Py_ssize_t npars; 
    503            
    504         if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 
    505          
     576 
     577        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 
     578 
    506579        OUTVECTOR(data_obj,pars,npars); 
    507                  
     580 
    508581        // Should create this list only once and refill it 
    509582    residuals = PyList_New(self->params.npoints); 
     
    511584    regterm = reg_term(pars, self->params.d_max, npars, nslice); 
    512585 
    513      
     586 
    514587    for(i=0; i<self->params.npoints; i++) { 
    515588        diff = self->params.y[i] - pr(pars, self->params.d_max, npars, self->params.x[i]); 
    516589        residual = diff*diff / (self->params.err[i]*self->params.err[i]); 
    517590        tmp = residual; 
    518          
     591 
    519592        // regularization term 
    520593        residual += self->params.alpha * regterm; 
    521          
     594 
    522595        if (PyList_SetItem(residuals, i, Py_BuildValue("d",residual) ) < 0){ 
    523             PyErr_SetString(CinvertorError,  
     596            PyErr_SetString(CinvertorError, 
    524597                "Cinvertor.residuals: error setting residual."); 
    525598                return NULL; 
    526599        }; 
    527                  
     600 
    528601    } 
    529      
     602 
    530603        return residuals; 
    531604} 
    532605 
    533 const char get_iq_doc[] =  
     606const char get_iq_doc[] = 
    534607        "Function to call to evaluate the scattering intensity\n" 
    535608        " @param args: c-parameters, and q\n" 
     
    546619        PyObject *data_obj; 
    547620        Py_ssize_t npars; 
    548            
     621 
    549622        if (!PyArg_ParseTuple(args, "Od", &data_obj, &q)) return NULL; 
    550623        OUTVECTOR(data_obj,pars,npars); 
    551                  
     624 
    552625        iq_value = iq(pars, self->params.d_max, npars, q); 
    553         return Py_BuildValue("f", iq_value);     
    554 } 
    555  
    556 const char get_pr_doc[] =  
     626        return Py_BuildValue("f", iq_value); 
     627} 
     628 
     629const char get_pr_doc[] = 
    557630        "Function to call to evaluate P(r)\n" 
    558631        " @param args: c-parameters and r\n" 
     
    569642        PyObject *data_obj; 
    570643        Py_ssize_t npars; 
    571            
     644 
    572645        if (!PyArg_ParseTuple(args, "Od", &data_obj, &r)) return NULL; 
    573646        OUTVECTOR(data_obj,pars,npars); 
    574                  
     647 
    575648        pr_value = pr(pars, self->params.d_max, npars, r); 
    576         return Py_BuildValue("f", pr_value);     
    577 } 
    578  
    579 const char get_pr_err_doc[] =  
     649        return Py_BuildValue("f", pr_value); 
     650} 
     651 
     652const char get_pr_err_doc[] = 
    580653        "Function to call to evaluate P(r) with errors\n" 
    581654        " @param args: c-parameters and r\n" 
     
    596669        PyObject *err_obj; 
    597670        Py_ssize_t npars2; 
    598         int i;  
    599          
     671        int i; 
     672 
    600673        if (!PyArg_ParseTuple(args, "OOd", &data_obj, &err_obj, &r)) return NULL; 
    601         OUTVECTOR(data_obj,pars,npars);  
     674        OUTVECTOR(data_obj,pars,npars); 
    602675        OUTVECTOR(err_obj,pars_err,npars2); 
    603676 
    604677        pr_err(pars, pars_err, self->params.d_max, npars, r, &pr_value, &pr_err_value); 
    605         return Py_BuildValue("ff", pr_value, pr_err_value);      
    606 } 
    607  
    608 const char basefunc_ft_doc[] =  
     678        return Py_BuildValue("ff", pr_value, pr_err_value); 
     679} 
     680 
     681const char basefunc_ft_doc[] = 
    609682        "Returns the value of the nth Fourier transofrmed base function\n" 
    610683        " @param args: c-parameters, n and q\n" 
     
    614687        double d_max, q; 
    615688        int n; 
    616          
     689 
    617690        if (!PyArg_ParseTuple(args, "did", &d_max, &n, &q)) return NULL; 
    618         return Py_BuildValue("f", ortho_transformed(d_max, n, q));       
    619          
    620 } 
    621  
    622 const char oscillations_doc[] =  
     691        return Py_BuildValue("f", ortho_transformed(d_max, n, q)); 
     692 
     693} 
     694 
     695const char oscillations_doc[] = 
    623696        "Returns the value of the oscillation figure of merit for\n" 
    624697        "the given set of coefficients. For a sphere, the oscillation\n" 
     
    632705        Py_ssize_t npars; 
    633706        double oscill, norm; 
    634          
     707 
    635708        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 
    636709        OUTVECTOR(data_obj,pars,npars); 
    637          
     710 
    638711        oscill = reg_term(pars, self->params.d_max, npars, 100); 
    639712        norm   = int_p2(pars, self->params.d_max, npars, 100); 
    640         return Py_BuildValue("f", sqrt(oscill/norm)/acos(-1.0)*self->params.d_max );     
    641          
    642 } 
    643  
    644 const char get_peaks_doc[] =  
     713        return Py_BuildValue("f", sqrt(oscill/norm)/acos(-1.0)*self->params.d_max ); 
     714 
     715} 
     716 
     717const char get_peaks_doc[] = 
    645718        "Returns the number of peaks in the output P(r) distrubution\n" 
    646719        "for the given set of coefficients.\n" 
     
    653726        Py_ssize_t npars; 
    654727        int count; 
    655          
     728 
    656729        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 
    657730        OUTVECTOR(data_obj,pars,npars); 
    658          
     731 
    659732        count = npeaks(pars, self->params.d_max, npars, 100); 
    660733 
    661         return Py_BuildValue("i", count );       
    662          
    663 } 
    664  
    665 const char get_positive_doc[] =  
     734        return Py_BuildValue("i", count ); 
     735 
     736} 
     737 
     738const char get_positive_doc[] = 
    666739        "Returns the fraction of P(r) that is positive over\n" 
    667740        "the full range of r for the given set of coefficients.\n" 
     
    674747        Py_ssize_t npars; 
    675748        double fraction; 
    676           
     749 
    677750        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 
    678751        OUTVECTOR(data_obj,pars,npars); 
    679          
     752 
    680753        fraction = positive_integral(pars, self->params.d_max, npars, 100); 
    681754 
    682         return Py_BuildValue("f", fraction );    
    683          
    684 } 
    685  
    686 const char get_pos_err_doc[] =  
     755        return Py_BuildValue("f", fraction ); 
     756 
     757} 
     758 
     759const char get_pos_err_doc[] = 
    687760        "Returns the fraction of P(r) that is 1 standard deviation\n" 
    688761        "above zero over the full range of r for the given set of coefficients.\n" 
     
    698771        Py_ssize_t npars2; 
    699772        double fraction; 
    700          
     773 
    701774        if (!PyArg_ParseTuple(args, "OO", &data_obj, &err_obj)) return NULL; 
    702         OUTVECTOR(data_obj,pars,npars);  
     775        OUTVECTOR(data_obj,pars,npars); 
    703776        OUTVECTOR(err_obj,pars_err,npars2); 
    704          
     777 
    705778        fraction = positive_errors(pars, pars_err, self->params.d_max, npars, 51); 
    706779 
    707         return Py_BuildValue("f", fraction );    
    708          
    709 } 
    710  
     780        return Py_BuildValue("f", fraction ); 
     781 
     782} 
     783 
     784const char get_rg_doc[] = 
     785        "Returns the value of the radius of gyration Rg.\n" 
     786        " @param args: c-parameters\n" 
     787        " @return: Rg"; 
     788 
     789static PyObject * get_rg(Cinvertor *self, PyObject *args) { 
     790        double *pars; 
     791        double *pars_err; 
     792        PyObject *data_obj; 
     793        Py_ssize_t npars; 
     794        Py_ssize_t npars2; 
     795        double value; 
     796 
     797        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 
     798        OUTVECTOR(data_obj,pars,npars); 
     799 
     800        value = rg(pars, self->params.d_max, npars, 101); 
     801 
     802        return Py_BuildValue("f", value ); 
     803 
     804} 
     805 
     806const char get_iq0_doc[] = 
     807        "Returns the value of I(q=0).\n" 
     808        " @param args: c-parameters\n" 
     809        " @return: I(q=0)"; 
     810 
     811static PyObject * get_iq0(Cinvertor *self, PyObject *args) { 
     812        double *pars; 
     813        double *pars_err; 
     814        PyObject *data_obj; 
     815        Py_ssize_t npars; 
     816        Py_ssize_t npars2; 
     817        double value; 
     818 
     819        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 
     820        OUTVECTOR(data_obj,pars,npars); 
     821 
     822        value = 4.0*acos(-1.0)*int_pr(pars, self->params.d_max, npars, 101); 
     823 
     824        return Py_BuildValue("f", value ); 
     825 
     826} 
     827 
     828/** 
     829 * Check whether a q-value is within acceptabel limits 
     830 * Return 1 if accepted, 0 if rejected. 
     831 */ 
     832int accept_q(Cinvertor *self, double q) { 
     833    if (self->params.q_min>0 && q<self->params.q_min) return 0; 
     834    if (self->params.q_max>0 && q>self->params.q_max) return 0; 
     835    return 1; 
     836} 
     837 
     838const char get_matrix_doc[] = 
     839        "Returns A matrix and b vector for least square problem.\n" 
     840        " @param nfunc: number of base functions\n" 
     841        " @param nr: number of r-points used when evaluating reg term.\n" 
     842        " @param a: A array to fill\n" 
     843        " @param b: b vector to fill\n" 
     844        " @return: 0"; 
     845 
     846static PyObject * get_matrix(Cinvertor *self, PyObject *args) { 
     847        double *a; 
     848        double *b; 
     849        PyObject *a_obj; 
     850        PyObject *b_obj; 
     851        Py_ssize_t n_a; 
     852        Py_ssize_t n_b; 
     853        // Number of bins for regularization term evaluation 
     854        int nr, nfunc; 
     855        int i, j, i_r; 
     856        double r, sqrt_alpha, pi; 
     857        double tmp; 
     858        int offset; 
     859 
     860        if (!PyArg_ParseTuple(args, "iiOO", &nfunc, &nr, &a_obj, &b_obj)) return NULL; 
     861        OUTVECTOR(a_obj,a,n_a); 
     862        OUTVECTOR(b_obj,b,n_b); 
     863 
     864        assert(n_b>=nfunc); 
     865        assert(n_a>=nfunc*(nr*self->params.npoints)); 
     866 
     867        sqrt_alpha = sqrt(self->params.alpha); 
     868        pi = acos(-1.0); 
     869        offset = (self->params.has_bck==1) ? 0 : 1; 
     870 
     871    for (j=0; j<nfunc; j++) { 
     872        for (i=0; i<self->params.npoints; i++) { 
     873            if (accept_q(self, self->params.x[i])){ 
     874                if (self->params.has_bck==1 && j==0) { 
     875                    a[i*nfunc+j] = 1.0/self->params.err[i]; 
     876                } else { 
     877                        a[i*nfunc+j] = ortho_transformed(self->params.d_max, j+offset, self->params.x[i])/self->params.err[i]; 
     878                } 
     879            } 
     880        } 
     881        for (i_r=0; i_r<nr; i_r++){ 
     882            if (self->params.has_bck==1 && j==0) { 
     883                a[(i_r+self->params.npoints)*nfunc+j] = 0.0; 
     884            } else { 
     885                    r = self->params.d_max/nr*i_r; 
     886                    tmp = pi*(j+offset)/self->params.d_max; 
     887                    a[(i_r+self->params.npoints)*nfunc+j] = sqrt_alpha * 1.0/nr*self->params.d_max*2.0* 
     888                    (2.0*pi*(j+offset)/self->params.d_max*cos(pi*(j+offset)*r/self->params.d_max) + 
     889                    tmp*tmp*r * sin(pi*(j+offset)*r/self->params.d_max)); 
     890                } 
     891        } 
     892    } 
     893 
     894    for (i=0; i<self->params.npoints; i++) { 
     895        if (accept_q(self, self->params.x[i])){ 
     896            b[i] = self->params.y[i]/self->params.err[i]; 
     897         } 
     898        } 
     899 
     900        return Py_BuildValue("i", 0); 
     901 
     902} 
     903 
     904const char get_invcov_matrix_doc[] = 
     905        " Compute the inverse covariance matrix, defined as inv_cov = a_transposed x a.\n" 
     906        " @param nfunc: number of base functions\n" 
     907        " @param nr: number of r-points used when evaluating reg term.\n" 
     908        " @param a: A array to fill\n" 
     909        " @param inv_cov: inverse covariance array to be filled\n" 
     910        " @return: 0"; 
     911 
     912static PyObject * get_invcov_matrix(Cinvertor *self, PyObject *args) { 
     913        double *a; 
     914        PyObject *a_obj; 
     915        Py_ssize_t n_a; 
     916        double *inv_cov; 
     917        PyObject *cov_obj; 
     918        Py_ssize_t n_cov; 
     919        int nr, nfunc; 
     920        int i, j, k; 
     921 
     922        if (!PyArg_ParseTuple(args, "iiOO", &nfunc, &nr, &a_obj, &cov_obj)) return NULL; 
     923        OUTVECTOR(a_obj,a,n_a); 
     924        OUTVECTOR(cov_obj,inv_cov,n_cov); 
     925 
     926        assert(n_cov>=nfunc*nfunc); 
     927        assert(n_a>=nfunc*(nr+self->params.npoints)); 
     928 
     929        for (i=0; i<nfunc; i++) { 
     930                for (j=0; j<nfunc; j++) { 
     931                        inv_cov[i*nfunc+j] = 0.0; 
     932                        for (k=0; k<nr+self->params.npoints; k++) { 
     933                                inv_cov[i*nfunc+j] += a[k*nfunc+i]*a[k*nfunc+j]; 
     934                        } 
     935                } 
     936        } 
     937        return Py_BuildValue("i", 0); 
     938} 
     939 
     940const char get_reg_size_doc[] = 
     941        " Compute the covariance matrix, defined as inv_cov = a_transposed x a.\n" 
     942        " @param nfunc: number of base functions\n" 
     943        " @param nr: number of r-points used when evaluating reg term.\n" 
     944        " @param a: A array to fill\n" 
     945        " @param inv_cov: inverse covariance array to be filled\n" 
     946        " @return: 0"; 
     947 
     948static PyObject * get_reg_size(Cinvertor *self, PyObject *args) { 
     949        double *a; 
     950        PyObject *a_obj; 
     951        Py_ssize_t n_a; 
     952        int nr, nfunc; 
     953        int i, j; 
     954        double sum_sig, sum_reg; 
     955 
     956        if (!PyArg_ParseTuple(args, "iiO", &nfunc, &nr, &a_obj)) return NULL; 
     957        OUTVECTOR(a_obj,a,n_a); 
     958 
     959        assert(n_a>=nfunc*(nr+self->params.npoints)); 
     960 
     961    sum_sig = 0.0; 
     962    sum_reg = 0.0; 
     963    for (j=0; j<nfunc; j++){ 
     964        for (i=0; i<self->params.npoints; i++){ 
     965            if (accept_q(self, self->params.x[i])==1) 
     966                sum_sig += (a[i*nfunc+j])*(a[i*nfunc+j]); 
     967        } 
     968        for (i=0; i<nr; i++){ 
     969            sum_reg += (a[(i+self->params.npoints)*nfunc+j])*(a[(i+self->params.npoints)*nfunc+j]); 
     970        } 
     971    } 
     972    return Py_BuildValue("ff", sum_sig, sum_reg); 
     973} 
    711974 
    712975const char eeeget_qmin_doc[] = "\ 
     
    714977\n\ 
    715978This is the second line."; 
    716 const char eeeset_qmin_doc[] =  
     979const char eeeset_qmin_doc[] = 
    717980        "This is a multiline doc string.\n" 
    718981        "\n" 
     
    736999                   {"set_alpha", (PyCFunction)set_alpha, METH_VARARGS, set_alpha_doc}, 
    7371000                   {"get_alpha", (PyCFunction)get_alpha, METH_VARARGS, get_alpha_doc}, 
     1001                   {"set_slit_width", (PyCFunction)set_slit_width, METH_VARARGS, set_slit_width_doc}, 
     1002                   {"get_slit_width", (PyCFunction)get_slit_width, METH_VARARGS, get_slit_width_doc}, 
     1003                   {"set_slit_height", (PyCFunction)set_slit_height, METH_VARARGS, set_slit_height_doc}, 
     1004                   {"get_slit_height", (PyCFunction)get_slit_height, METH_VARARGS, get_slit_height_doc}, 
     1005                   {"set_has_bck", (PyCFunction)set_has_bck, METH_VARARGS, set_has_bck_doc}, 
     1006                   {"get_has_bck", (PyCFunction)get_has_bck, METH_VARARGS, get_has_bck_doc}, 
    7381007                   {"get_nx", (PyCFunction)get_nx, METH_VARARGS, get_nx_doc}, 
    7391008                   {"get_ny", (PyCFunction)get_ny, METH_VARARGS, get_ny_doc}, 
     
    7481017                   {"get_positive", (PyCFunction)get_positive, METH_VARARGS, get_positive_doc}, 
    7491018                   {"get_pos_err", (PyCFunction)get_pos_err, METH_VARARGS, get_pos_err_doc}, 
    750     
     1019                   {"rg", (PyCFunction)get_rg, METH_VARARGS, get_rg_doc}, 
     1020                   {"iq0", (PyCFunction)get_iq0, METH_VARARGS, get_iq0_doc}, 
     1021                   {"_get_matrix", (PyCFunction)get_matrix, METH_VARARGS, get_matrix_doc}, 
     1022                   {"_get_invcov_matrix", (PyCFunction)get_invcov_matrix, METH_VARARGS, get_invcov_matrix_doc}, 
     1023                   {"_get_reg_size", (PyCFunction)get_reg_size, METH_VARARGS, get_reg_size_doc}, 
     1024 
    7511025   {NULL} 
    7521026}; 
     
    7961070 
    7971071static PyMethodDef module_methods[] = { 
    798     {NULL}  
     1072    {NULL} 
    7991073}; 
    8001074 
     
    8021076 * Function used to add the model class to a module 
    8031077 * @param module: module to add the class to 
    804  */  
     1078 */ 
    8051079void addCinvertor(PyObject *module) { 
    8061080        PyObject *d; 
    807          
     1081 
    8081082    if (PyType_Ready(&CinvertorType) < 0) 
    8091083        return; 
     
    8111085    Py_INCREF(&CinvertorType); 
    8121086    PyModule_AddObject(module, "Cinvertor", (PyObject *)&CinvertorType); 
    813      
     1087 
    8141088    d = PyModule_GetDict(module); 
    8151089    CinvertorError = PyErr_NewException("Cinvertor.error", NULL, NULL); 
     
    8221096#endif 
    8231097PyMODINIT_FUNC 
    824 initpr_inversion(void)  
     1098initpr_inversion(void) 
    8251099{ 
    8261100    PyObject* m; 
     
    8281102    m = Py_InitModule3("pr_inversion", module_methods, 
    8291103                       "C extension module for inversion to P(r)."); 
    830                         
     1104 
    8311105    addCinvertor(m); 
    8321106} 
  • pr_inversion/c_extensions/invertor.c

    rc61228f r9a23253e  
    2020        pars->q_min = -1.0; 
    2121        pars->q_max = -1.0; 
     22        pars->has_bck = 0; 
    2223} 
    2324 
     
    2526/** 
    2627 * P(r) of a sphere, for test purposes 
    27  *  
     28 * 
    2829 * @param R: radius of the sphere 
    2930 * @param r: distance, in the same units as the radius 
     
    4142 * Orthogonal functions: 
    4243 * B(r) = 2r sin(pi*nr/d) 
    43  *  
     44 * 
    4445 */ 
    4546double ortho(double d_max, int n, double r) { 
     
    4950/** 
    5051 * Fourier transform of the nth orthogonal function 
    51  *  
     52 * 
    5253 */ 
    5354double ortho_transformed(double d_max, int n, double q) { 
     
    5758 
    5859/** 
     60 * Slit-smeared Fourier transform of the nth orthogonal function. 
     61 * Smearing follows Lake, Acta Cryst. (1967) 23, 191. 
     62 */ 
     63double ortho_transformed_smeared(double d_max, int n, double height, double width, double q, int npts) { 
     64        double sum, value, y, z; 
     65        int i, j; 
     66        double fnpts; 
     67        sum = 0.0; 
     68        fnpts = (float)npts-1.0; 
     69 
     70        for(i=0; i<npts; i++) { 
     71                y = -width/2.0+width/fnpts*(float)i; 
     72                for(j=0; j<npts; j++) { 
     73                        z = height/fnpts*(float)j; 
     74                        sum += ortho_transformed(d_max, n, sqrt((q-y)*(q-y)+z*z)); 
     75                } 
     76        } 
     77 
     78        return sum/npts/npts/height/width; 
     79} 
     80 
     81/** 
    5982 * First derivative in of the orthogonal function dB(r)/dr 
    60  *  
     83 * 
    6184 */ 
    6285double ortho_derived(double d_max, int n, double r) { 
     
    80103 */ 
    81104double pr(double *pars, double d_max, int n_c, double r) { 
    82     double sum = 0.0;   
     105    double sum = 0.0; 
    83106        int i; 
    84107    for (i=0; i<n_c; i++) { 
     
    91114 * P(r) calculated from the expansion, with errors 
    92115 */ 
    93 void pr_err(double *pars, double *err, double d_max, int n_c,  
     116void pr_err(double *pars, double *err, double d_max, int n_c, 
    94117                double r, double *pr_value, double *pr_value_err) { 
    95     double sum = 0.0;  
     118    double sum = 0.0; 
    96119    double sum_err = 0.0; 
    97120    double func_value; 
     
    109132        *pr_value_err = sum; 
    110133    } 
    111 }  
     134} 
    112135 
    113136/** 
     
    115138 */ 
    116139double dprdr(double *pars, double d_max, int n_c, double r) { 
    117     double sum = 0.0;  
    118         int i;  
     140    double sum = 0.0; 
     141        int i; 
    119142    for (i=0; i<n_c; i++) { 
    120143        sum += pars[i] * 2.0*(sin(pi*(i+1)*r/d_max) + pi*(i+1)*r/d_max * cos(pi*(i+1)*r/d_max)); 
     
    127150 */ 
    128151double reg_term(double *pars, double d_max, int n_c, int nslice) { 
    129     double sum = 0.0;  
     152    double sum = 0.0; 
    130153    double r; 
    131154    double deriv; 
     
    143166 */ 
    144167double int_p2(double *pars, double d_max, int n_c, int nslice) { 
    145     double sum = 0.0;  
    146     double r;  
     168    double sum = 0.0; 
     169    double r; 
    147170    double value; 
    148171        int i; 
     
    156179 
    157180/** 
     181 * Integral of P(r) 
     182 */ 
     183double int_pr(double *pars, double d_max, int n_c, int nslice) { 
     184    double sum = 0.0; 
     185    double r; 
     186    double value; 
     187        int i; 
     188    for (i=0; i<nslice; i++) { 
     189        r = d_max/(1.0*nslice)*i; 
     190        value = pr(pars, d_max, n_c, r); 
     191        sum += value; 
     192    } 
     193    return sum/(1.0*nslice)*d_max; 
     194} 
     195 
     196/** 
    158197 * Get the number of P(r) peaks. 
    159198 */ 
    160199int npeaks(double *pars, double d_max, int n_c, int nslice) { 
    161     double r;  
     200    double r; 
    162201    double value; 
    163202        int i; 
     
    187226 */ 
    188227double positive_integral(double *pars, double d_max, int n_c, int nslice) { 
    189     double r;  
     228    double r; 
    190229    double value; 
    191230        int i; 
    192231        double sum_pos = 0.0; 
    193232        double sum = 0.0; 
    194          
     233 
    195234    for (i=0; i<nslice; i++) { 
    196235        r = d_max/(1.0*nslice)*i; 
     
    207246 */ 
    208247double positive_errors(double *pars, double *err, double d_max, int n_c, int nslice) { 
    209     double r;  
    210     double value; 
    211         int i;  
     248    double r; 
     249    double value; 
     250        int i; 
    212251        double sum_pos = 0.0; 
    213252        double sum = 0.0; 
    214253        double pr_val; 
    215254        double pr_val_err; 
    216          
     255 
    217256    for (i=0; i<nslice; i++) { 
    218257        r = d_max/(1.0*nslice)*i; 
     
    220259        if (pr_val>pr_val_err) sum_pos += pr_val; 
    221260        sum += fabs(pr_val); 
    222          
     261 
    223262 
    224263    } 
    225264    return sum_pos/sum; 
    226265} 
     266 
     267/** 
     268 * R_g radius of gyration calculation 
     269 * 
     270 * R_g**2 = integral[r**2 * p(r) dr] /  (2.0 * integral[p(r) dr]) 
     271 */ 
     272double rg(double *pars, double d_max, int n_c, int nslice) { 
     273    double sum_r2 = 0.0; 
     274    double sum    = 0.0; 
     275    double r; 
     276    double value; 
     277        int i; 
     278    for (i=0; i<nslice; i++) { 
     279        r = d_max/(1.0*nslice)*i; 
     280        value = pr(pars, d_max, n_c, r); 
     281        sum += value; 
     282        sum_r2 += r*r*value; 
     283    } 
     284    return sqrt(sum_r2/(2.0*sum)); 
     285} 
     286 
  • pr_inversion/c_extensions/invertor.h

    r896abb3 r9a23253e  
    1515    double *err; 
    1616    /// Number of q points 
    17     int npoints;     
     17    int npoints; 
    1818    /// Number of I(q) points 
    19     int ny;     
     19    int ny; 
    2020    /// Number of dI(q) points 
    21     int nerr;   
     21    int nerr; 
    2222    /// Alpha value 
    2323    double alpha; 
     
    2626    /// Maximum q to include in inversion 
    2727    double q_max; 
    28 } Invertor_params;  
     28    /// Flag for whether or not to evalute a constant background while inverting 
     29    int has_bck; 
     30    /// Slit height in units of q [A-1] 
     31    double slit_height; 
     32    /// Slit width in units of q [A-1] 
     33    double slit_width; 
     34} Invertor_params; 
    2935 
    3036void invertor_dealloc(Invertor_params *pars); 
    3137 
    3238void invertor_init(Invertor_params *pars); 
    33  
    3439 
    3540double pr_sphere(double R, double r); 
     
    4247double reg_term(double *pars, double d_max, int n_c, int nslice); 
    4348double int_p2(double *pars, double d_max, int n_c, int nslice); 
    44 void pr_err(double *pars, double *err, double d_max, int n_c,  
     49void pr_err(double *pars, double *err, double d_max, int n_c, 
    4550                double r, double *pr_value, double *pr_value_err); 
    4651int npeaks(double *pars, double d_max, int n_c, int nslice); 
    4752double positive_integral(double *pars, double d_max, int n_c, int nslice); 
    4853double positive_errors(double *pars, double *err, double d_max, int n_c, int nslice); 
     54double rg(double *pars, double d_max, int n_c, int nslice); 
     55double int_pr(double *pars, double d_max, int n_c, int nslice); 
     56double ortho_transformed_smeared(double d_max, int n, double heigth, double width, double q, int npts); 
    4957 
    5058#endif 
Note: See TracChangeset for help on using the changeset viewer.