Changeset 9a23253e in sasview


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
Files:
6 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 
  • pr_inversion/invertor.py

    rb00b487 r9a23253e  
    6363        - get_pos_err(pars): returns the fraction of P(r) that is 1-sigma above zero 
    6464    """ 
     65    #TODO: Allow for slit smearing. Smear each base function once before filling 
     66    # the A matrix. 
     67     
    6568    ## Chisqr of the last computation 
    6669    chi2  = 0 
     
    7578    ## Last errors on output values 
    7679    cov = None 
     80    ## Flag to allow I(q) data with constant background 
     81    #has_bck = False 
     82    ## Background value 
     83    background = 0 
     84     
    7785     
    7886    def __init__(self): 
     
    106114        elif name=='alpha': 
    107115            return self.set_alpha(value) 
     116        elif name=='slit_height': 
     117            return self.set_slit_height(value) 
     118        elif name=='slit_width': 
     119            return self.set_slit_width(value) 
     120        elif name=='has_bck': 
     121            if value==True: 
     122                return self.set_has_bck(1) 
     123            elif value==False: 
     124                return self.set_has_bck(0) 
     125            else: 
     126                raise ValueError, "Invertor: has_bck can only be True or False" 
    108127             
    109128        return Cinvertor.__setattr__(self, name, value) 
     
    142161        elif name=='alpha': 
    143162            return self.get_alpha() 
     163        elif name=='slit_height': 
     164            return self.get_slit_height() 
     165        elif name=='slit_width': 
     166            return self.get_slit_width() 
     167        elif name=='has_bck': 
     168            value = self.get_has_bck() 
     169            if value==1: 
     170                return True 
     171            else: 
     172                return False 
    144173        elif name in self.__dict__: 
    145174            return self.__dict__[name] 
     
    161190        invertor.y = self.y 
    162191        invertor.err = self.err 
     192        invertor.has_bck = self.has_bck 
    163193         
    164194        return invertor 
     
    194224            @return: c_out, c_cov - the coefficients with covariance matrix  
    195225        """ 
    196         #TODO: call the pyhton implementation for now. In the future, translate this to C. 
     226        # Reset the background value before proceeding 
     227        self.background = 0.0 
    197228        return self.lstsq(nfunc, nr=nr) 
     229     
     230    def iq(self, out, q): 
     231        """ 
     232            Function to call to evaluate the scattering intensity 
     233            @param args: c-parameters, and q 
     234            @return: I(q) 
     235        """ 
     236        return Cinvertor.iq(self, out, q)+self.background 
    198237     
    199238    def invert_optimize(self, nfunc=10, nr=20): 
     
    325364 
    326365        """ 
    327         import math 
     366        #TODO: Allow for background by starting at n=0 (since the base function 
     367        # is zero for n=0). 
     368        import math, time 
    328369        from scipy.linalg.basic import lstsq 
    329370         
     
    339380        if sqrt_alpha<0.0: 
    340381            nq = 0 
    341          
     382 
     383        # If we need to fit the background, add a term 
     384        if self.has_bck==True: 
     385            nfunc_0 = nfunc 
     386            nfunc += 1 
     387 
    342388        a = numpy.zeros([npts+nq, nfunc]) 
    343389        b = numpy.zeros(npts+nq) 
    344390        err = numpy.zeros([nfunc, nfunc]) 
    345391         
    346         for j in range(nfunc): 
    347             for i in range(npts): 
    348                 if self._accept_q(self.x[i]): 
    349                     a[i][j] = self.basefunc_ft(self.d_max, j+1, self.x[i])/self.err[i] 
    350                      
    351             #TODO: refactor this: i_q should really be i_r 
    352             for i_q in range(nq): 
    353                 r = self.d_max/nq*i_q 
    354                 #a[i_q+npts][j] = sqrt_alpha * 1.0/nq*self.d_max*2.0*math.fabs(math.sin(math.pi*(j+1)*r/self.d_max) + math.pi*(j+1)*r/self.d_max * math.cos(math.pi*(j+1)*r/self.d_max))      
    355                 a[i_q+npts][j] = sqrt_alpha * 1.0/nq*self.d_max*2.0*(2.0*math.pi*(j+1)/self.d_max*math.cos(math.pi*(j+1)*r/self.d_max) + math.pi**2*(j+1)**2*r/self.d_max**2 * math.sin(math.pi*(j+1)*r/self.d_max))      
    356          
    357         for i in range(npts): 
    358             if self._accept_q(self.x[i]): 
    359                 b[i] = self.y[i]/self.err[i] 
    360              
     392        # Construct the a matrix and b vector that represent the problem 
     393        self._get_matrix(nfunc, nq, a, b) 
     394              
     395        # Perform the inversion (least square fit) 
    361396        c, chi2, rank, n = lstsq(a, b) 
    362397        # Sanity check 
     
    367402        self.chi2 = chi2 
    368403                 
    369         at = numpy.transpose(a) 
    370404        inv_cov = numpy.zeros([nfunc,nfunc]) 
    371         for i in range(nfunc): 
    372             for j in range(nfunc): 
    373                 inv_cov[i][j] = 0.0 
    374                 for k in range(npts+nr): 
    375                     #if self._accept_q(self.x[i]): 
    376                     inv_cov[i][j] += at[i][k]*a[k][j] 
     405        # Get the covariance matrix, defined as inv_cov = a_transposed * a 
     406        self._get_invcov_matrix(nfunc, nr, a, inv_cov) 
    377407                     
    378408        # Compute the reg term size for the output 
    379         sum_sig = 0.0 
    380         sum_reg = 0.0 
    381         for j in range(nfunc): 
    382             for i in range(npts): 
    383                 if self._accept_q(self.x[i]): 
    384                     sum_sig += (a[i][j])**2 
    385             for i in range(nq): 
    386                 sum_reg += (a[i+npts][j])**2 
     409        sum_sig, sum_reg = self._get_reg_size(nfunc, nr, a) 
    387410                     
    388411        if math.fabs(self.alpha)>0: 
     
    403426             
    404427        # Keep a copy of the last output 
    405         self.out = c 
    406         self.cov = err 
    407          
    408         return c, err 
     428        if self.has_bck==False: 
     429            self.background = 0 
     430            self.out = c 
     431            self.cov = err 
     432        else: 
     433            self.background = c[0] 
     434             
     435            err_0 = numpy.zeros([nfunc, nfunc]) 
     436            c_0 = numpy.zeros(nfunc) 
     437             
     438            for i in range(nfunc_0): 
     439                c_0[i] = c[i+1] 
     440                for j in range(nfunc_0): 
     441                    err_0[i][j] = err[i+1][j+1] 
     442                     
     443            self.out = c_0 
     444            self.cov = err_0 
     445             
     446        return self.out, self.cov 
     447         
     448    def lstsq_bck(self, nfunc=5, nr=20): 
     449        #TODO: Allow for background by starting at n=0 (since the base function 
     450        # is zero for n=0). 
     451        import math 
     452        from scipy.linalg.basic import lstsq 
     453         
     454        if self.is_valid()<0: 
     455            raise RuntimeError, "Invertor: invalid data; incompatible data lengths." 
     456         
     457        self.nfunc = nfunc 
     458        # a -- An M x N matrix. 
     459        # b -- An M x nrhs matrix or M vector. 
     460        npts = len(self.x) 
     461        nq   = nr 
     462        sqrt_alpha = math.sqrt(math.fabs(self.alpha)) 
     463        if sqrt_alpha<0.0: 
     464            nq = 0 
     465         
     466        err_0 = numpy.zeros([nfunc, nfunc]) 
     467        c_0 = numpy.zeros(nfunc) 
     468        nfunc_0 = nfunc 
     469        nfunc += 1 
     470         
     471        a = numpy.zeros([npts+nq, nfunc]) 
     472        b = numpy.zeros(npts+nq) 
     473        err = numpy.zeros([nfunc, nfunc]) 
     474         
     475        # Construct the a matrix and b vector that represent the problem 
     476        self._get_matrix(nfunc, nq, a, b) 
     477             
     478        c, chi2, rank, n = lstsq(a, b) 
     479        # Sanity check 
     480        try: 
     481            float(chi2) 
     482        except: 
     483            chi2 = -1.0 
     484        self.chi2 = chi2 
     485 
     486        inv_cov = numpy.zeros([nfunc,nfunc]) 
     487        # Get the covariance matrix, defined as inv_cov = a_transposed * a 
     488        self._get_invcov_matrix(nfunc, nr, a, inv_cov) 
     489                     
     490        # Compute the reg term size for the output 
     491        sum_sig, sum_reg = self._get_reg_size(nfunc, nr, a) 
     492         
     493        if math.fabs(self.alpha)>0: 
     494            new_alpha = sum_sig/(sum_reg/self.alpha) 
     495        else: 
     496            new_alpha = 0.0 
     497        self.suggested_alpha = new_alpha 
     498         
     499        try: 
     500            cov = numpy.linalg.pinv(inv_cov) 
     501            err = math.fabs(chi2/float(npts-nfunc)) * cov 
     502        except: 
     503            # We were not able to estimate the errors, 
     504            # returns an empty covariance matrix 
     505            print "lstsq:", sys.exc_value 
     506            print chi2 
     507            pass 
     508             
     509        # Keep a copy of the last output 
     510         
     511        print "BACKGROUND =", c[0] 
     512        self.background = c[0] 
     513         
     514        for i in range(nfunc_0): 
     515            c_0[i] = c[i+1] 
     516            for j in range(nfunc_0): 
     517                err_0[i][j] = err[i+1][j+1] 
     518                 
     519        self.out = c_0 
     520        self.cov = err_0 
     521         
     522        return c_0, err_0 
    409523         
    410524    def estimate_alpha(self, nfunc): 
     
    433547                  
    434548            # Perform inversion to find the largest alpha 
    435             out, cov = pr.lstsq(nfunc) 
     549            out, cov = pr.invert(nfunc) 
    436550            elapsed = time.time()-starttime 
    437551            initial_alpha = pr.alpha 
     
    440554            # Try the inversion with the estimated alpha 
    441555            pr.alpha = pr.suggested_alpha 
    442             out, cov = pr.lstsq(nfunc) 
     556            out, cov = pr.invert(nfunc) 
    443557     
    444558            npeaks = pr.get_peaks(out) 
     
    458572                for i in range(10): 
    459573                    pr.alpha = (0.33)**(i+1)*alpha 
    460                     out, cov = pr.lstsq(nfunc) 
     574                    out, cov = pr.invert(nfunc) 
    461575                     
    462576                    peaks = pr.get_peaks(out) 
  • pr_inversion/test/test_output.txt

    rb00b487 r9a23253e  
    55#elapsed=0 
    66#alpha_estimate=4.58991e-006 
    7 #C_0=217.876478953+-0.0394946442786 
    8 #C_1=0.30164617344+-0.407412922728 
    9 #C_2=7.15694574584+-2.73726481828 
    10 #C_3=-0.79690812244+-5.7183647083 
    11 #C_4=0.948333847403+-16.2086432642 
    12 #C_5=-0.751612279429+-28.2108893047 
    13 #C_6=-0.0796650783501+-52.0568306543 
    14 #C_7=-0.518109559543+-88.6068891309 
    15 #C_8=-0.129910831984+-141.7178308 
    16 #C_9=-0.258134624641+-215.767430098 
     7#C_0=217.876478953+-2522.42314999 
     8#C_1=0.30164617344+-226.644556603 
     9#C_2=7.15694574584+-43.4919208221 
     10#C_3=-0.79690812244+-25.8031634616 
     11#C_4=0.948333847403+-11.3794625443 
     12#C_5=-0.751612279429+-8.64083336005 
     13#C_6=-0.0796650783501+-5.73820483268 
     14#C_7=-0.518109559543+-3.56660274557 
     15#C_8=-0.129910831984+-1.99841316045 
     16#C_9=-0.258134624641+-0.756509301101 
    1717<r>  <Pr>  <dPr> 
    18180  0  0 
    19 1.6  22.981  20.5137 
    20 3.2  91.986  78.6612 
    21 4.8  207.187  164.741 
    22 6.4  368.819  264.129 
    23 8  577.117  359.427 
    24 9.6  832.224  433.162 
    25 11.2  1134.11  470.89 
    26 12.8  1482.46  464.814 
    27 14.4  1876.63  419.086 
    28 16  2315.57  360.693 
    29 17.6  2797.81  353.245 
    30 19.2  3321.46  448.357 
    31 20.8  3884.25  611.481 
    32 22.4  4483.6  786.464 
    33 24  5116.72  934.559 
    34 25.6  5780.7  1030.78 
    35 27.2  6472.66  1062.21 
    36 28.8  7189.82  1030.28 
    37 30.4  7929.59  955.681 
    38 32  8689.7  883.97 
    39 33.6  9468.17  878.871 
    40 35.2  10263.4  977.553 
    41 36.8  11073.9  1152 
    42 38.4  11898.7  1342.59 
    43 40  12737  1498.3 
    44 41.6  13587.7  1587.08 
    45 43.2  14450.2  1596.95 
    46 44.8  15323.5  1537.99 
    47 46.4  16206.5  1445.58 
    48 48  17098  1379.39 
    49 49.6  17996.3  1401.86 
    50 51.2  18899.8  1532.06 
    51 52.8  19806.5  1728.86 
    52 54.4  20714.3  1927.37 
    53 56  21620.8  2073.39 
    54 57.6  22524  2135.21 
    55 59.2  23421.6  2106.97 
    56 60.8  24311.7  2011.28 
    57 62.4  25192.3  1900.75 
    58 64  26061.8  1848.67 
    59 65.6  26918.7  1912.83 
    60 67.2  27761.7  2089.96 
    61 68.8  28589.6  2319.48 
    62 70.4  29401.3  2527.71 
    63 72  30195.6  2658.76 
    64 73.6  30971.4  2685.17 
    65 75.2  31727.2  2611.68 
    66 76.8  32461.5  2478.05 
    67 78.4  33172.7  2357.19 
    68 80  33858.8  2333.81 
    69 81.6  34517.8  2453.43 
    70 83.2  35147.6  2684.61 
    71 84.8  35745.9  2945.13 
    72 86.4  36310.5  3152.18 
    73 88  36839.5  3249.59 
    74 89.6  37330.9  3217.04 
    75 91.2  37783.1  3074.61 
    76 92.8  38194.6  2886.09 
    77 94.4  38564.3  2752.58 
    78 96  38891.1  2773.03 
    79 97.6  39174.1  2973.02 
    80 99.2  39412.6  3281.62 
    81 100.8  39605.8  3587.59 
    82 102.4  39752.6  3795.95 
    83 104  39851.9  3850.56 
    84 105.6  39902.4  3741.32 
    85 107.2  39902.5  3510.03 
    86 108.8  39850.1  3255.58 
    87 110.4  39743.2  3120.25 
    88 112  39579.6  3218.61 
    89 113.6  39356.7  3539.26 
    90 115.2  39072.3  3954.92 
    91 116.8  38724.1  4317.73 
    92 118.4  38310.3  4517.08 
    93 120  37829  4494.86 
    94 121.6  37278.9  4251.44 
    95 123.2  36659  3856.58 
    96 124.8  35968.4  3465.78 
    97 126.4  35206.6  3306.22 
    98 128  34373  3535.97 
    99 129.6  33466.9  4074.35 
    100 131.2  32487.5  4693.64 
    101 132.8  31433.4  5186.14 
    102 134.4  30303.1  5411.6 
    103 136  29094.1  5295.75 
    104 137.6  27803.8  4830.98 
    105 139.2  26428.9  4092.02 
    106 140.8  24965.9  3287.78 
    107 142.4  23411.1  2858.06 
    108 144  21761.3  3246.41 
    109 145.6  20013.4  4237.67 
    110 147.2  18165.4  5345.59 
    111 148.8  16216.4  6256.27 
    112 150.4  14166.9  6786.26 
    113 152  12019  6829.04 
    114 153.6  9776.61  6338.91 
    115 155.2  7445.34  5326.31 
    116 156.8  5032.57  3854 
    117 158.4  2547.22  2030.4 
     191.6  22.981  7.48836 
     203.2  91.986  29.737 
     214.8  207.187  66.1302 
     226.4  368.819  115.749 
     238  577.117  177.508 
     249.6  832.224  250.315 
     2511.2  1134.11  333.205 
     2612.8  1482.46  425.442 
     2714.4  1876.63  526.549 
     2816  2315.57  636.272 
     2917.6  2797.81  754.499 
     3019.2  3321.46  881.145 
     3120.8  3884.25  1016.06 
     3222.4  4483.6  1158.97 
     3324  5116.72  1309.47 
     3425.6  5780.7  1467.05 
     3527.2  6472.66  1631.18 
     3628.8  7189.82  1801.33 
     3730.4  7929.59  1977.02 
     3832  8689.7  2157.88 
     3933.6  9468.17  2343.55 
     4035.2  10263.4  2533.72 
     4136.8  11073.9  2728.06 
     4238.4  11898.7  2926.2 
     4340  12737  3127.75 
     4441.6  13587.7  3332.27 
     4543.2  14450.2  3539.34 
     4644.8  15323.5  3748.54 
     4746.4  16206.5  3959.49 
     4848  17098  4171.83 
     4949.6  17996.3  4385.18 
     5051.2  18899.8  4599.19 
     5152.8  19806.5  4813.46 
     5254.4  20714.3  5027.56 
     5356  21620.8  5241.07 
     5457.6  22524  5453.53 
     5559.2  23421.6  5664.53 
     5660.8  24311.7  5873.68 
     5762.4  25192.3  6080.65 
     5864  26061.8  6285.13 
     5965.6  26918.7  6486.85 
     6067.2  27761.7  6685.53 
     6168.8  28589.6  6880.91 
     6270.4  29401.3  7072.73 
     6372  30195.6  7260.72 
     6473.6  30971.4  7444.63 
     6575.2  31727.2  7624.23 
     6676.8  32461.5  7799.3 
     6778.4  33172.7  7969.61 
     6880  33858.8  8134.91 
     6981.6  34517.8  8294.9 
     7083.2  35147.6  8449.24 
     7184.8  35745.9  8597.54 
     7286.4  36310.5  8739.34 
     7388  36839.5  8874.2 
     7489.6  37330.9  9001.64 
     7591.2  37783.1  9121.2 
     7692.8  38194.6  9232.39 
     7794.4  38564.3  9334.72 
     7896  38891.1  9427.68 
     7997.6  39174.1  9510.74 
     8099.2  39412.6  9583.35 
     81100.8  39605.8  9644.98 
     82102.4  39752.6  9695.14 
     83104  39851.9  9733.38 
     84105.6  39902.4  9759.36 
     85107.2  39902.5  9772.74 
     86108.8  39850.1  9773.24 
     87110.4  39743.2  9760.52 
     88112  39579.6  9734.22 
     89113.6  39356.7  9693.89 
     90115.2  39072.3  9639.07 
     91116.8  38724.1  9569.28 
     92118.4  38310.3  9484.1 
     93120  37829  9383.18 
     94121.6  37278.9  9266.24 
     95123.2  36659  9133 
     96124.8  35968.4  8983.13 
     97126.4  35206.6  8816.15 
     98128  34373  8631.45 
     99129.6  33466.9  8428.29 
     100131.2  32487.5  8205.96 
     101132.8  31433.4  7963.92 
     102134.4  30303.1  7701.95 
     103136  29094.1  7420.23 
     104137.6  27803.8  7119.29 
     105139.2  26428.9  6799.69 
     106140.8  24965.9  6461.63 
     107142.4  23411.1  6104.48 
     108144  21761.3  5726.35 
     109145.6  20013.4  5323.88 
     110147.2  18165.4  4892.46 
     111148.8  16216.4  4426.74 
     112150.4  14166.9  3921.46 
     113152  12019  3372.52 
     114153.6  9776.61  2777.82 
     115155.2  7445.34  2138.05 
     116156.8  5032.57  1456.94 
     117158.4  2547.22  741.175 
  • pr_inversion/test/utest_invertor.py

    rb00b487 r9a23253e  
    6363            self.x_in[i] = 1.0*(i+1) 
    6464 
     65    def test_has_bck_flag(self): 
     66        """ 
     67            Tests the has_bck flag operations 
     68        """ 
     69        self.assertEqual(self.invertor.has_bck, False) 
     70        self.invertor.has_bck=True 
     71        self.assertEqual(self.invertor.has_bck, True) 
     72        def doit_float(): 
     73            self.invertor.has_bck  = 2.0 
     74        def doit_str(): 
     75            self.invertor.has_bck  = 'a' 
     76 
     77        self.assertRaises(ValueError, doit_float) 
     78        self.assertRaises(ValueError, doit_str) 
     79         
    6580 
    6681    def testset_dmax(self): 
     
    177192        # Perform inversion 
    178193        out, cov = self.invertor.invert_optimize(10) 
     194        #out, cov = self.invertor.invert(10) 
    179195        # This is a very specific case 
    180196        # We should make sure it always passes 
Note: See TracChangeset for help on using the changeset viewer.