Changeset eca05c8 in sasview for pr_inversion


Ignore:
Timestamp:
Apr 30, 2008 4:29:36 PM (17 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:
34ae302
Parents:
47f695c9
Message:

Added P(r) fit test

Location:
pr_inversion
Files:
1 added
5 edited

Legend:

Unmodified
Added
Removed

  • pr_inversion/c_extensions/Cinvertor.c

    r9e8dc22 reca05c8  
    222222} 
    223223 
     224static PyObject * set_alpha(Cinvertor *self, PyObject *args) { 
     225        double alpha; 
     226   
     227        if (!PyArg_ParseTuple(args, "d", &alpha)) return NULL; 
     228        self->params.alpha = alpha; 
     229        return Py_BuildValue("d", self->params.alpha);   
     230} 
     231 
     232/** 
     233 * Gets the maximum distance 
     234 */ 
     235static PyObject * get_alpha(Cinvertor *self, PyObject *args) { 
     236        return Py_BuildValue("d", self->params.alpha);   
     237} 
     238 
    224239/** 
    225240 * Gets the number of x points 
    … …  
    246261/** 
    247262 * Function to call to evaluate the residuals 
    248  * @param args: input q or [q,phi] 
    249  * @return: function value 
     263 * @param args: input parameters 
     264 * @return: list of residuals 
    250265 */ 
    251266static PyObject * residuals(Cinvertor *self, PyObject *args) { 
    … …  
    257272        double residual, diff; 
    258273        // Regularization factor 
    259         double lmult = 0.0; 
     274        double regterm = 0.0; 
     275        double tmp = 0.0; 
    260276         
    261277        PyObject *data_obj; 
    … …  
    266282        OUTVECTOR(data_obj,pars,npars); 
    267283                 
    268         //printf("npars: %i", npars); 
    269284    // PyList of residuals 
    270285        // Should create this list only once and refill it 
    271286    residuals = PyList_New(self->params.npoints); 
     287 
     288    regterm = reg_term(pars, self->params.d_max, npars); 
    272289     
    273290    for(i=0; i<self->params.npoints; i++) { 
    274         //printf("\n%i: ", i); 
    275291        diff = self->params.y[i] - iq(pars, self->params.d_max, npars, self->params.x[i]); 
    276292        residual = diff*diff / (self->params.err[i]*self->params.err[i]); 
    277         //printf("%i %g\n", i, residual); 
     293        tmp = residual; 
     294         
     295        // regularization term 
     296        residual += self->params.alpha * regterm; 
     297         
    278298        if (PyList_SetItem(residuals, i, Py_BuildValue("d",residual) ) < 0){ 
    279299            PyErr_SetString(CinvertorError,  
    … …  
    286306        return residuals; 
    287307} 
     308/** 
     309 * Function to call to evaluate the residuals 
     310 * for P(r) minimization (for testing purposes) 
     311 * @param args: input parameters 
     312 * @return: list of residuals 
     313 */ 
     314static PyObject * pr_residuals(Cinvertor *self, PyObject *args) { 
     315        double *pars; 
     316        PyObject* residuals; 
     317        PyObject* temp; 
     318        double *res; 
     319        int i; 
     320        double residual, diff; 
     321        // Regularization factor 
     322        double regterm = 0.0; 
     323        double tmp = 0.0; 
     324         
     325        PyObject *data_obj; 
     326        Py_ssize_t npars; 
     327           
     328        if (!PyArg_ParseTuple(args, "O", &data_obj)) return NULL; 
     329         
     330        OUTVECTOR(data_obj,pars,npars); 
     331                 
     332        // Should create this list only once and refill it 
     333    residuals = PyList_New(self->params.npoints); 
     334 
     335    regterm = reg_term(pars, self->params.d_max, npars); 
     336 
     337     
     338    for(i=0; i<self->params.npoints; i++) { 
     339        diff = self->params.y[i] - pr(pars, self->params.d_max, npars, self->params.x[i]); 
     340        residual = diff*diff / (self->params.err[i]*self->params.err[i]); 
     341        tmp = residual; 
     342         
     343        // regularization term 
     344        residual += self->params.alpha * regterm; 
     345         
     346        if (PyList_SetItem(residuals, i, Py_BuildValue("d",residual) ) < 0){ 
     347            PyErr_SetString(CinvertorError,  
     348                "Cinvertor.residuals: error setting residual."); 
     349                return NULL; 
     350        }; 
     351                 
     352    } 
     353     
     354        return residuals; 
     355} 
    288356 
    289357/** 
    … …  
    323391} 
    324392 
     393/** 
     394 * Function to call to evaluate P(r) with errors 
     395 * @param args: c-parameters and r 
     396 * @return: P(r) 
     397 */ 
     398static PyObject * get_pr_err(Cinvertor *self, PyObject *args) { 
     399        double *pars; 
     400        double *pars_err; 
     401        double pr_err_value; 
     402        double r, pr_value; 
     403        PyObject *data_obj; 
     404        Py_ssize_t npars; 
     405        PyObject *err_obj; 
     406        Py_ssize_t npars2; 
     407           
     408        if (!PyArg_ParseTuple(args, "OOd", &data_obj, &err_obj, &r)) return NULL; 
     409        OUTVECTOR(data_obj,pars,npars); 
     410        OUTVECTOR(err_obj,pars_err,npars2); 
     411                 
     412        pr_err(pars, pars_err, self->params.d_max, npars, r, &pr_value, &pr_err_value); 
     413        return Py_BuildValue("ff", pr_value, pr_err_value);      
     414} 
     415 
    325416 
    326417static PyMethodDef Cinvertor_methods[] = { 
    327418                   {"residuals", (PyCFunction)residuals, METH_VARARGS, "Get the list of residuals"}, 
     419                   {"pr_residuals", (PyCFunction)pr_residuals, METH_VARARGS, "Get the list of residuals"}, 
    328420                   {"set_x", (PyCFunction)set_x, METH_VARARGS, ""}, 
    329421                   {"get_x", (PyCFunction)get_x, METH_VARARGS, ""}, 
    … …  
    334426                   {"set_dmax", (PyCFunction)set_dmax, METH_VARARGS, ""}, 
    335427                   {"get_dmax", (PyCFunction)get_dmax, METH_VARARGS, ""}, 
     428                   {"set_alpha", (PyCFunction)set_alpha, METH_VARARGS, ""}, 
     429                   {"get_alpha", (PyCFunction)get_alpha, METH_VARARGS, ""}, 
    336430                   {"get_nx", (PyCFunction)get_nx, METH_VARARGS, ""}, 
    337431                   {"get_ny", (PyCFunction)get_ny, METH_VARARGS, ""}, 
    … …  
    339433                   {"iq", (PyCFunction)get_iq, METH_VARARGS, ""}, 
    340434                   {"pr", (PyCFunction)get_pr, METH_VARARGS, ""}, 
     435                   {"get_pr_err", (PyCFunction)get_pr_err, METH_VARARGS, ""}, 
    341436                   {"is_valid", (PyCFunction)is_valid, METH_VARARGS, ""}, 
    342437    

  • pr_inversion/c_extensions/invertor.c

    r9e8dc22 reca05c8  
    4040double ortho(double d_max, int n, double r) { 
    4141        return 2.0*r*sin(pi*n*r/d_max); 
    42          
    4342} 
    4443 
    … …  
    7675 */ 
    7776double pr(double *pars, double d_max, int n_c, double r) { 
    78     double sum = 0.0;  
     77    double sum = 0.0;   
    7978        int i; 
    8079    for (i=0; i<n_c; i++) { 
    … …  
    8483} 
    8584 
     85/** 
     86 * P(r) calculated from the expansion, with errors 
     87 */ 
     88void pr_err(double *pars, double *err, double d_max, int n_c,  
     89                double r, double *pr_value, double *pr_value_err) { 
     90    double sum = 0.0;  
     91    double sum_err = 0.0; 
     92    double func_value; 
     93        int i; 
     94    for (i=0; i<n_c; i++) { 
     95        func_value = ortho(d_max, i+1, r); 
     96        sum += pars[i] * func_value; 
     97        sum_err += err[i]*err[i]*func_value*func_value; 
     98    } 
     99    *pr_value = sum; 
     100    if (sum_err>0) { 
     101        *pr_value_err = sqrt(sum_err); 
     102    } else { 
     103        *pr_value_err = sum; 
     104    } 
     105}  
    86106 
     107/** 
     108 * dP(r)/dr calculated from the expansion. 
     109 */ 
     110double dprdr(double *pars, double d_max, int n_c, double r) { 
     111    double sum = 0.0;  
     112        int i; 
     113    for (i=0; i<n_c; i++) { 
     114        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)); 
     115    } 
     116    return sum; 
     117} 
     118 
     119/** 
     120 * regularization term calculated from the expansion. 
     121 */ 
     122double reg_term(double *pars, double d_max, int n_c) { 
     123    double sum = 0.0;  
     124    double r; 
     125    double deriv; 
     126        int i; 
     127    for (i=0; i<25; i++) { 
     128        r = d_max/25.0*i; 
     129        deriv = dprdr(pars, d_max, n_c, r); 
     130        sum += deriv*deriv; 
     131    } 
     132    return sum/25.0*d_max; 
     133} 
     134 
     135 

  • pr_inversion/c_extensions/invertor.h

    r9e8dc22 reca05c8  
    1010    int npoints;     
    1111    int ny;     
    12     int nerr;     
    13 } Invertor_params; 
     12    int nerr;   
     13    double alpha; 
     14} Invertor_params;  
    1415 
    1516void invertor_dealloc(Invertor_params *pars); 
    … …  
    2425double iq(double *c, double d_max, int n_c, double q); 
    2526double pr(double *c, double d_max, int n_c, double r); 
     27double dprdr(double *pars, double d_max, int n_c, double r); 
     28double reg_term(double *pars, double d_max, int n_c); 
     29void pr_err(double *pars, double *err, double d_max, int n_c,  
     30                double r, double *pr_value, double *pr_value_err); 
    2631 
    2732#endif 

  • pr_inversion/invertor.py

    r9e8dc22 reca05c8  
    33 
    44class Invertor(Cinvertor): 
     5     
     6    ## Chisqr of the last computation 
     7    chisqr = 0 
    58     
    69    def __init__(self): 
    … …  
    1417        """ 
    1518        if   name=='x': 
     19            if 0.0 in value: 
     20                raise ValueError, "Invertor: one of your q-values is zero. Delete that entry before proceeding" 
    1621            return self.set_x(value) 
    1722        elif name=='y': 
    … …  
    2126        elif name=='d_max': 
    2227            return self.set_dmax(value) 
     28        elif name=='alpha': 
     29            return self.set_alpha(value) 
    2330             
    2431        return Cinvertor.__setattr__(self, name, value) 
    … …  
    4552        elif name=='d_max': 
    4653            return self.get_dmax() 
     54        elif name=='alpha': 
     55            return self.get_alpha() 
    4756        elif name in self.__dict__: 
    4857            return self.__dict__[name] 
    … …  
    6271        out, cov_x, info, mesg, success = optimize.leastsq(self.residuals, p, full_output=1, warning=True) 
    6372         
     73        # Compute chi^2 
     74        res = self.residuals(out) 
     75        chisqr = 0 
     76        for i in range(len(res)): 
     77            chisqr += res[i] 
     78         
     79        self.chi2 = chisqr 
     80         
    6481        return out, cov_x 
     82     
     83    def pr_fit(self, nfunc=5): 
     84        """ 
     85            Perform inversion to P(r) 
     86        """ 
     87        from scipy import optimize 
     88         
     89        # First, check that the current data is valid 
     90        if self.is_valid()<=0: 
     91            raise RuntimeError, "Invertor.invert: Data arrays are of different length" 
     92         
     93        p = numpy.ones(nfunc) 
     94        out, cov_x, info, mesg, success = optimize.leastsq(self.pr_residuals, p, full_output=1, warning=True) 
     95         
     96        # Compute chi^2 
     97        res = self.pr_residuals(out) 
     98        chisqr = 0 
     99        for i in range(len(res)): 
     100            chisqr += res[i] 
     101         
     102        self.chisqr = chisqr 
     103         
     104        return out, cov_x 
     105     
     106    def pr_err(self, c, c_cov, r): 
     107        import math 
     108        c_err = numpy.zeros(len(c)) 
     109        for i in range(len(c)): 
     110            try: 
     111                c_err[i] = math.sqrt(math.fabs(c_cov[i][i])) 
     112            except: 
     113                import sys 
     114                print sys.exc_value 
     115                print "oups", c_cov[i][i] 
     116                c_err[i] = c[i] 
     117 
     118        return self.get_pr_err(c, c_err, r) 
     119         
    65120     
    66121if __name__ == "__main__": 

  • pr_inversion/test/utest_invertor.py

    r9e8dc22 reca05c8  
    88 
    99 
    10 import unittest, math, numpy 
     10import unittest, math, numpy, pylab 
    1111from sans.pr.invertor import Invertor 
    1212         
    … …  
    2121        self.x_in = numpy.ones(self.ntest) 
    2222        for i in range(self.ntest): 
    23             self.x_in[i] = 1.0*i 
     23            self.x_in[i] = 1.0*(i+1) 
    2424 
    2525 
    … …  
    3131        self.invertor.d_max = value 
    3232        self.assertEqual(self.invertor.d_max, value) 
     33         
     34    def testset_alpha(self): 
     35        """ 
     36            Set and read alpha 
     37        """ 
     38        value = 15.0 
     39        self.invertor.alpha = value 
     40        self.assertEqual(self.invertor.alpha, value) 
    3341         
    3442    def testset_x_1(self): 
    … …  
    113121        self.assertEqual(self.invertor.test_no_data, None) 
    114122         
     123    def test_inversion(self): 
     124        """ 
     125            Test an inversion for which we know the answer 
     126        """ 
     127        x, y, err = load("sphere_80.txt") 
     128 
     129        # Choose the right d_max... 
     130        self.invertor.d_max = 160.0 
     131        # Set a small alpha 
     132        self.invertor.alpha = 1e-7 
     133        # Set data 
     134        self.invertor.x   = x 
     135        self.invertor.y   = y 
     136        self.invertor.err = err 
     137        # Perform inversion 
     138        out, cov = self.invertor.invert(10) 
     139        # This is a very specific case 
     140        # We should make sure it always passes 
     141        self.assertTrue(self.invertor.chi2/len(x)<200.00) 
     142         
     143        # Check the computed P(r) with the theory 
     144        # for shpere of radius 80 
     145        x = pylab.arange(0.01, self.invertor.d_max, self.invertor.d_max/51.0) 
     146        y = numpy.zeros(len(x)) 
     147        dy = numpy.zeros(len(x)) 
     148        y_true = numpy.zeros(len(x)) 
     149 
     150        sum = 0.0 
     151        sum_true = 0.0 
     152        for i in range(len(x)): 
     153            #y[i] = self.invertor.pr(out, x[i]) 
     154            (y[i], dy[i]) = self.invertor.pr_err(out, cov, x[i]) 
     155            sum += y[i] 
     156            if x[i]<80.0: 
     157                y_true[i] = pr_theory(x[i], 80.0) 
     158            else: 
     159                y_true[i] = 0 
     160            sum_true += y_true[i] 
     161             
     162        y = y/sum*self.invertor.d_max/len(x) 
     163        dy = dy/sum*self.invertor.d_max/len(x) 
     164        y_true = y_true/sum_true*self.invertor.d_max/len(x) 
     165         
     166        chi2 = 0.0 
     167        for i in range(len(x)): 
     168            res = (y[i]-y_true[i])/dy[i] 
     169            chi2 += res*res 
     170             
     171        try: 
     172            self.assertTrue(chi2/51.0<10.0) 
     173        except: 
     174            print "chi2 =", chi2/51.0 
     175            raise 
     176             
     177    def test_q_zero(self): 
     178        """ 
     179            Test error condition where a point has q=0 
     180        """ 
     181        x, y, err = load("sphere_80.txt") 
     182        x[0] = 0.0 
     183         
     184        # Choose the right d_max... 
     185        self.invertor.d_max = 160.0 
     186        # Set a small alpha 
     187        self.invertor.alpha = 1e-7 
     188        # Set data 
     189        def doit(): 
     190            self.invertor.x   = x 
     191        self.assertRaises(ValueError, doit) 
     192         
     193                             
     194    def test_q_neg(self): 
     195        """ 
     196            Test error condition where a point has q<0 
     197        """ 
     198        x, y, err = load("sphere_80.txt") 
     199        x[0] = -0.2 
     200         
     201        # Choose the right d_max... 
     202        self.invertor.d_max = 160.0 
     203        # Set a small alpha 
     204        self.invertor.alpha = 1e-7 
     205        # Set data 
     206        self.invertor.x   = x 
     207        self.invertor.y   = y 
     208        self.invertor.err = err 
     209        # Perform inversion 
     210        out, cov = self.invertor.invert(4) 
     211         
     212        try: 
     213            self.assertTrue(self.invertor.chi2>0) 
     214        except: 
     215            print "Chi2 =", self.invertor.chi2 
     216            raise 
     217                             
     218    def test_Iq_zero(self): 
     219        """ 
     220            Test error condition where a point has q<0 
     221        """ 
     222        x, y, err = load("sphere_80.txt") 
     223        y[0] = 0.0 
     224         
     225        # Choose the right d_max... 
     226        self.invertor.d_max = 160.0 
     227        # Set a small alpha 
     228        self.invertor.alpha = 1e-7 
     229        # Set data 
     230        self.invertor.x   = x 
     231        self.invertor.y   = y 
     232        self.invertor.err = err 
     233        # Perform inversion 
     234        out, cov = self.invertor.invert(4) 
     235         
     236        try: 
     237            self.assertTrue(self.invertor.chi2>0) 
     238        except: 
     239            print "Chi2 =", self.invertor.chi2 
     240            raise 
     241                             
     242         
     243 
     244def pr_theory(r, R): 
     245    """ 
     246       P(r) for a sphere 
     247    """ 
     248    if r<=2*R: 
     249        return 12.0* ((0.5*r/R)**2) * ((1.0-0.5*r/R)**2) * ( 2.0 + 0.5*r/R ) 
     250    else: 
     251        return 0.0 
     252     
     253def load(path = "sphere_60_q0_2.txt"): 
     254    import numpy, math, sys 
     255    # Read the data from the data file 
     256    data_x   = numpy.zeros(0) 
     257    data_y   = numpy.zeros(0) 
     258    data_err = numpy.zeros(0) 
     259    if not path == None: 
     260        input_f = open(path,'r') 
     261        buff    = input_f.read() 
     262        lines   = buff.split('\n') 
     263        for line in lines: 
     264            try: 
     265                toks = line.split() 
     266                x = float(toks[0]) 
     267                y = float(toks[1]) 
     268                data_x = numpy.append(data_x, x) 
     269                data_y = numpy.append(data_y, y) 
     270                # Set the error of the first point to 5% 
     271                # to make theory look like data 
     272                scale = 0.1/math.sqrt(data_x[0]) 
     273                data_err = numpy.append(data_err, scale*math.sqrt(y)) 
     274            except: 
     275                pass 
     276                
     277    return data_x, data_y, data_err  
     278        
    115279if __name__ == '__main__': 
    116280    unittest.main() 
Note: See TracChangeset for help on using the changeset viewer.