Changeset eba9885 in sasview for sansmodels/src/sans/models


Ignore:
Timestamp:
Aug 5, 2009 5:39:03 PM (15 years ago)
Author:
Gervaise Alina <gervyh@…>
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:
f88624d
Parents:
8344c50
Message:

code for evalDistribution

Location:
sansmodels/src/sans/models
Files:
6 added
10 edited

Legend:

Unmodified
Added
Removed

  • sansmodels/src/sans/models/c_extensions/c_models.c

    r812b901 reba9885  
    5252        addDisperser(m); 
    5353        addCGaussian(m); 
     54        addCSchulz(m); 
     55        addCLogNormal(m); 
    5456        addCLorentzian(m); 
    5557        addCHollowCylinderModel(m); 

  • sansmodels/src/sans/models/c_models/c_models.cpp

    r812b901 reba9885  
    4848        void addCGaussian(PyObject *module); 
    4949        void addCLorentzian(PyObject *module); 
     50        void addCLogNormal(PyObject *module); 
     51        void addCSchulz(PyObject *module); 
    5052} 
    5153 
    … …  
    6971 
    7072/** 
     73 * Delete a lognormal dispersion model object 
     74 */ 
     75void del_lognormal_dispersion(void *ptr){ 
     76        LogNormalDispersion * disp = static_cast<LogNormalDispersion *>(ptr); 
     77        delete disp; 
     78        return; 
     79} 
     80 
     81/** 
     82 * Create a lognormal dispersion model as a python object 
     83 */ 
     84PyObject * new_lognormal_dispersion(PyObject *, PyObject *args) { 
     85        LogNormalDispersion *disp = new LogNormalDispersion(); 
     86        return PyCObject_FromVoidPtr(disp, del_lognormal_dispersion); 
     87} 
     88 
     89/** 
    7190 * Delete a gaussian dispersion model object 
    7291 */ 
    … …  
    84103        return PyCObject_FromVoidPtr(disp, del_gaussian_dispersion); 
    85104} 
     105 
     106/** 
     107 * Delete a schulz dispersion model object 
     108 */ 
     109void del_schulz_dispersion(void *ptr){ 
     110        SchulzDispersion * disp = static_cast<SchulzDispersion *>(ptr); 
     111        delete disp; 
     112        return; 
     113} 
     114/** 
     115 * Create a schulz dispersion model as a python object 
     116 */ 
     117PyObject * new_schulz_dispersion(PyObject *, PyObject *args) { 
     118        SchulzDispersion *disp = new SchulzDispersion(); 
     119        return PyCObject_FromVoidPtr(disp, del_schulz_dispersion); 
     120} 
     121 
    86122 
    87123/** 
    … …  
    147183        {"new_gaussian_model",   (PyCFunction)new_gaussian_dispersion, METH_VARARGS, 
    148184                  "Create a new GaussianDispersion object"}, 
     185    {"new_lognormal_model",   (PyCFunction)new_lognormal_dispersion, METH_VARARGS, 
     186                  "Create a new LogNormalDispersion object"}, 
     187    {"new_schulz_model",   (PyCFunction)new_schulz_dispersion, METH_VARARGS, 
     188                  "Create a new SchulzDispersion object"}, 
    149189        {"new_array_model",      (PyCFunction)new_array_dispersion  , METH_VARARGS, 
    150190                  "Create a new ArrayDispersion object"}, 
    … …  
    194234        addDisperser(m); 
    195235        addCGaussian(m); 
     236        addCSchulz(m); 
     237        addCLogNormal(m); 
    196238        addCLorentzian(m); 
    197239        addCVesicleModel(m); 
    198240 
    199  
    200 } 
     241} 

  • sansmodels/src/sans/models/c_models/dispersion_visitor.cpp

    rfca6936 reba9885  
    2525 
    2626        PyDict_SetItemString(dict, "type",  Py_BuildValue("s", "gaussian")); 
     27    PyDict_SetItemString(dict, "npts",  Py_BuildValue("i", disp->npts)); 
     28    PyDict_SetItemString(dict, "width", Py_BuildValue("d", disp->width)); 
     29    PyDict_SetItemString(dict, "nsigmas", Py_BuildValue("i", disp->nsigmas)); 
     30#endif 
     31} 
     32 
     33void DispersionVisitor:: lognormal_to_dict(void* dispersion, void* dictionary) { 
     34#ifndef __MODELS_STANDALONE__ 
     35        LogNormalDispersion * disp = (LogNormalDispersion*)dispersion; 
     36        PyObject * dict = (PyObject*)dictionary; 
     37 
     38        PyDict_SetItemString(dict, "type",  Py_BuildValue("s", "lognormal")); 
     39    PyDict_SetItemString(dict, "npts",  Py_BuildValue("i", disp->npts)); 
     40    PyDict_SetItemString(dict, "width", Py_BuildValue("d", disp->width)); 
     41    PyDict_SetItemString(dict, "nsigmas", Py_BuildValue("i", disp->nsigmas)); 
     42#endif 
     43} 
     44 
     45 
     46void DispersionVisitor:: schulz_to_dict(void* dispersion, void* dictionary) { 
     47#ifndef __MODELS_STANDALONE__ 
     48        SchulzDispersion * disp = (SchulzDispersion*)dispersion; 
     49        PyObject * dict = (PyObject*)dictionary; 
     50 
     51        PyDict_SetItemString(dict, "type",  Py_BuildValue("s", "schulz")); 
    2752    PyDict_SetItemString(dict, "npts",  Py_BuildValue("i", disp->npts)); 
    2853    PyDict_SetItemString(dict, "width", Py_BuildValue("d", disp->width)); 
    … …  
    6186} 
    6287 
     88void DispersionVisitor:: lognormal_from_dict(void* dispersion, void* dictionary) { 
     89#ifndef __MODELS_STANDALONE__ 
     90        LogNormalDispersion * disp = (LogNormalDispersion*)dispersion; 
     91        PyObject * dict = (PyObject*)dictionary; 
     92 
     93        disp->npts    = PyInt_AsLong( PyDict_GetItemString(dict, "npts") ); 
     94        disp->width   = PyFloat_AsDouble( PyDict_GetItemString(dict, "width") ); 
     95        disp->nsigmas = PyFloat_AsDouble( PyDict_GetItemString(dict, "nsigmas") ); 
     96#endif 
     97} 
     98void DispersionVisitor:: schulz_from_dict(void* dispersion, void* dictionary) { 
     99#ifndef __MODELS_STANDALONE__ 
     100        SchulzDispersion * disp = (SchulzDispersion*)dispersion; 
     101        PyObject * dict = (PyObject*)dictionary; 
     102 
     103        disp->npts    = PyInt_AsLong( PyDict_GetItemString(dict, "npts") ); 
     104        disp->width   = PyFloat_AsDouble( PyDict_GetItemString(dict, "width") ); 
     105        disp->nsigmas = PyFloat_AsDouble( PyDict_GetItemString(dict, "nsigmas") ); 
     106#endif 
     107} 
     108 
    63109void DispersionVisitor:: array_from_dict(void* dispersion, void* dictionary) {} 

  • sansmodels/src/sans/models/c_models/dispersion_visitor.hh

    rfca6936 reba9885  
    1919        void dispersion_to_dict(void *, void *); 
    2020        void gaussian_to_dict(void *, void *); 
     21        void lognormal_to_dict(void *, void *); 
     22        void schulz_to_dict(void*, void *); 
    2123        void array_to_dict(void *, void *); 
    2224 
    2325        void dispersion_from_dict(void*, void *); 
    2426        void gaussian_from_dict(void*, void *); 
     27        void lognormal_from_dict(void*, void *); 
     28        void schulz_from_dict(void*, void *); 
    2529        void array_from_dict(void*, void *); 
    2630 

  • sansmodels/src/sans/models/c_models/parameters.cpp

    r07da749 reba9885  
    136136        } 
    137137} 
     138 
     139 
     140/** 
     141 * LogNormal dispersion 
     142 */ 
     143 
     144LogNormalDispersion :: LogNormalDispersion() { 
     145        npts  = 1; 
     146        width = 0.0; 
     147        nsigmas = 2; 
     148}; 
     149 
     150void LogNormalDispersion :: accept_as_source(DispersionVisitor* visitor, void* from, void* to) { 
     151        visitor->lognormal_to_dict(from, to); 
     152} 
     153void LogNormalDispersion :: accept_as_destination(DispersionVisitor* visitor, void* from, void* to) { 
     154        visitor->lognormal_from_dict(from, to); 
     155} 
     156 
     157double lognormal_weight(double mean, double sigma, double x) { 
     158         
     159        double sigma2 = pow(sigma, 2);   
     160        return 1/(x*sigma2) * exp( -pow((log(x) -mean), 2) / (2*sigma2)); 
     161  
     162} 
     163 
     164/** 
     165 * Lognormal dispersion 
     166 * @param mean: mean value of the LogNormal 
     167 * @param sigma: standard deviation of the LogNormal 
     168 * @param x: value at which the LogNormal is evaluated 
     169 * @return: value of the LogNormal 
     170 */ 
     171void LogNormalDispersion :: operator() (void *param, vector<WeightPoint> &weights){ 
     172        // Check against zero width 
     173        if (width<=0) { 
     174                width = 0.0; 
     175                npts  = 1; 
     176                nsigmas = 3; 
     177        } 
     178 
     179        Parameter* par = (Parameter*)param; 
     180        double value = (*par)(); 
     181 
     182        if (npts<2) { 
     183                weights.insert(weights.end(), WeightPoint(value, 1.0)); 
     184        } else { 
     185                for(int i=0; i<npts; i++) { 
     186                        // We cover n(nsigmas) times sigmas on each side of the mean 
     187                        double val = value + width * (2.0*nsigmas*i/float(npts-1) - nsigmas); 
     188 
     189                        if ( ((*par).has_min==false || val>(*par).min) 
     190                          && ((*par).has_max==false || val<(*par).max)  ) { 
     191                                double _w = lognormal_weight(value, width, val); 
     192                                weights.insert(weights.end(), WeightPoint(val, _w)); 
     193                        } 
     194                } 
     195        } 
     196} 
     197 
     198 
     199 
     200/** 
     201 * Schulz dispersion 
     202 */ 
     203 
     204SchulzDispersion :: SchulzDispersion() { 
     205        npts  = 1; 
     206        width = 0.0; 
     207        nsigmas = 2; 
     208}; 
     209 
     210void SchulzDispersion :: accept_as_source(DispersionVisitor* visitor, void* from, void* to) { 
     211        visitor->schulz_to_dict(from, to); 
     212} 
     213void SchulzDispersion :: accept_as_destination(DispersionVisitor* visitor, void* from, void* to) { 
     214        visitor->schulz_from_dict(from, to); 
     215} 
     216 
     217double schulz_weight(double mean, double sigma, double x) { 
     218        double vary, expo_value; 
     219    double z = pow(mean/ sigma, 2)-1;    
     220        double R= x/mean; 
     221        double zz= z+1; 
     222        return  pow(zz,zz) * pow(R,z) * exp(-1*R*zz)/((mean) * tgamma(zz)) ; 
     223} 
     224 
     225/** 
     226 * Schulz dispersion 
     227 * @param mean: mean value of the Schulz 
     228 * @param sigma: standard deviation of the Schulz 
     229 * @param x: value at which the Schulz is evaluated 
     230 * @return: value of the Schulz 
     231 */ 
     232void SchulzDispersion :: operator() (void *param, vector<WeightPoint> &weights){ 
     233        // Check against zero width 
     234        if (width<=0) { 
     235                width = 0.0; 
     236                npts  = 1; 
     237                nsigmas = 3; 
     238        } 
     239 
     240        Parameter* par = (Parameter*)param; 
     241        double value = (*par)(); 
     242 
     243        if (npts<2) { 
     244                weights.insert(weights.end(), WeightPoint(value, 1.0)); 
     245        } else { 
     246                for(int i=0; i<npts; i++) { 
     247                        // We cover n(nsigmas) times sigmas on each side of the mean 
     248                        double val = value + width * (2.0*nsigmas*i/float(npts-1) - nsigmas); 
     249 
     250                        if ( ((*par).has_min==false || val>(*par).min) 
     251                          && ((*par).has_max==false || val<(*par).max)  ) { 
     252                                double _w = schulz_weight(value, width, val); 
     253                                weights.insert(weights.end(), WeightPoint(val, _w)); 
     254                        } 
     255                } 
     256        } 
     257} 
     258 
     259 
     260 
    138261 
    139262/** 

  • sansmodels/src/sans/models/c_models/parameters.hh

    rfca6936 reba9885  
    7777 
    7878/** 
     79 * Schulz dispersion model 
     80 */ 
     81class SchulzDispersion: public DispersionModel { 
     82public: 
     83        /// Number of sigmas on each side of the mean 
     84        int nsigmas; 
     85 
     86        SchulzDispersion(); 
     87        void operator()(void *, vector<WeightPoint>&); 
     88        void accept_as_source(DispersionVisitor*, void*, void*); 
     89        void accept_as_destination(DispersionVisitor*, void*, void*); 
     90}; 
     91 
     92/** 
     93 * LogNormal dispersion model 
     94 */ 
     95class LogNormalDispersion: public DispersionModel { 
     96public: 
     97        /// Number of sigmas on each side of the mean 
     98        int nsigmas; 
     99 
     100        LogNormalDispersion(); 
     101        void operator()(void *, vector<WeightPoint>&); 
     102        void accept_as_source(DispersionVisitor*, void*, void*); 
     103        void accept_as_destination(DispersionVisitor*, void*, void*); 
     104}; 
     105 
     106 
     107/** 
    79108 * Dispersion model based on arrays provided by the user 
    80109 */ 

  • sansmodels/src/sans/models/dispersion_models.py

    r988130c6 reba9885  
    7474            Set the weights of an array dispersion 
    7575        """ 
    76         message = "set_weights is not available for GaussiantDispersion.\n" 
     76        message = "set_weights is not available for GaussianDispersion.\n" 
     77        message += "  Solution: Use an ArrayDispersion object" 
     78        raise "RuntimeError", message 
     79     
     80class SchulzDispersion(DispersionModel): 
     81    """ 
     82        Python bridge class for a dispersion model based  
     83        on a Schulz distribution. 
     84    """ 
     85    def __init__(self): 
     86        self.cdisp = c_models.new_schulz_model() 
     87         
     88    def set_weights(self, values, weights): 
     89        """ 
     90            Set the weights of an array dispersion 
     91        """ 
     92        message = "set_weights is not available for SchulzDispersion.\n" 
     93        message += "  Solution: Use an ArrayDispersion object" 
     94        raise "RuntimeError", message 
     95     
     96class LogNormalDispersion(DispersionModel): 
     97    """ 
     98        Python bridge class for a dispersion model based  
     99        on a Log Normal distribution. 
     100    """ 
     101    def __init__(self): 
     102        self.cdisp = c_models.new_lognormal_model() 
     103         
     104    def set_weights(self, values, weights): 
     105        """ 
     106            Set the weights of an array dispersion 
     107        """ 
     108        message = "set_weights is not available for LogNormalDispersion.\n" 
    77109        message += "  Solution: Use an ArrayDispersion object" 
    78110        raise "RuntimeError", message 
    … …  
    100132        c_models.set_dispersion_weights(self.cdisp, values, weights) 
    101133  
    102 models = {GaussianDispersion:"GaussianModel", ArrayDispersion:"MyModel"}        
     134models = {GaussianDispersion:"GaussianModel", ArrayDispersion:"MyModel", 
     135          SchulzDispersion: "Schulz", LogNormalDispersion: "LogNormal"}        
    103136         

  • sansmodels/src/sans/models/test/utest_dispersity.py

    r8809e48 reba9885  
    5757         
    5858        new_model = self.model.clone() 
     59        print "gaussian",self.model.run(0.001) 
    5960        self.assertAlmostEqual(new_model.run(0.001), 4723.32213339, 3) 
    6061        self.assertAlmostEqual(new_model.runXY([0.001,0.001]), 4743.56, 2) 
     62         
     63    def test_schulz_zero(self): 
     64        from sans.models.dispersion_models import SchulzDispersion 
     65        disp = SchulzDispersion() 
     66        self.model.set_dispersion('radius', disp) 
     67        self.model.dispersion['radius']['width'] = 5.0 
     68        #self.model.dispersion['radius']['width'] = 0.0 
     69        self.model.dispersion['radius']['npts'] = 100 
     70        #self.model.setParam('scale', 1.0) 
     71        self.model.setParam('scale', 10.0) 
     72        print "schulz",self.model.run(0.001), self.model.dispersion 
     73        self.assertAlmostEqual(self.model.run(0.001), 450.355, 3) 
     74        self.assertAlmostEqual(self.model.runXY([0.001,0.001]), 452.299, 3) 
     75         
     76    def test_lognormal_zero(self): 
     77        from sans.models.dispersion_models import LogNormalDispersion 
     78        disp = LogNormalDispersion() 
     79        self.model.set_dispersion('radius', disp) 
     80        self.model.dispersion['radius']['width'] = 5.0 
     81        #self.model.dispersion['radius']['width'] = 0.0 
     82        self.model.dispersion['radius']['npts'] = 100 
     83        #self.model.setParam('scale', 1.0) 
     84        self.model.setParam('scale', 10.0) 
     85        print "model dispersion",self.model.dispersion 
     86        print "lognormal",self.model.run(0.001) 
     87        self.assertAlmostEqual(self.model.run(0.001), 450.355, 3) 
     88        self.assertAlmostEqual(self.model.runXY([0.001,0.001]), 452.299, 3) 
    6189         
    6290    def test_gaussian_zero(self): 

  • sansmodels/src/sans/models/test/utest_models.py

    rae60f86 reba9885  
    2424    def test1D(self): 
    2525        """ Test 1D model for a sphere """ 
    26         self.assertAlmostEqual(self.comp.run(1.0), 56.3878, 4) 
     26        self.assertAlmostEqual(self.comp.run(1.0), 5.6387e-5, 4) 
    2727         
    2828    def test1D_2(self): 
    2929        """ Test 2D model for a sphere """ 
    30         self.assertAlmostEqual(self.comp.run([1.0, 1.3]), 56.3878, 4) 
     30        self.assertAlmostEqual(self.comp.run([1.0, 1.3]), 5.63878e-5, 4) 
    3131 
    3232class TestCyl(unittest.TestCase): 

  • sansmodels/src/sans/models/test/utest_models_array.py

    r812b901 reba9885  
    88         
    99class TestSphere(unittest.TestCase): 
    10     """ Unit tests for sphere model """ 
     10    """ Unit tests for sphere model using evalDistribution function """ 
    1111     
    1212    def setUp(self): 
    … …  
    1818    def test1D(self): 
    1919        """ Test 1D model for a sphere  with vector as input""" 
    20         answer=numpy.array([5.63877831e-05,2.57231782e-06,2.73704050e-07,2.54229069e-08]) 
     20        answer = numpy.array([5.63877831e-05,2.57231782e-06,2.73704050e-07,2.54229069e-08]) 
    2121        
    22         testvector= self.comp.run(self.x) 
     22        
     23        testvector= self.comp.evalDistribution(self.x) 
    2324        
    2425        self.assertAlmostEqual(len(testvector),4) 
    … …  
    3334        """ Test 2D model for a sphere for 2 scalar """ 
    3435        self.assertAlmostEqual(self.comp.run([1.0, 1.3]), 56.3878e-06, 4) 
    35     
     36         
     37    def test1D_3(self): 
     38        """ Test 2D model for a Shpere for 2 vectors as input """ 
     39        x= numpy.reshape(self.x, [len(self.x),1]) 
     40        y= numpy.reshape(self.y, [1,len(self.y)]) 
     41        vect = self.comp.evalDistribution([x,y]) 
     42        self.assertAlmostEqual(vect[0][0],9.2985e-07, 4) 
     43        self.assertAlmostEqual(vect[len(self.x)-1][len(self.y)-1],1.3871e-08, 4) 
     44         
    3645         
    3746class TestCylinder(unittest.TestCase): 
    38     """ Unit tests for sphere model """ 
     47    """ Unit tests for Cylinder model using evalDistribution function """ 
    3948     
    4049    def setUp(self): 
    … …  
    4251        self.comp = CylinderModel() 
    4352        self.x = numpy.array([1.0,2.0,3.0, 4.0]) 
    44         self.y = numpy.array([1.0,2.0,3.0, 4.0]) 
     53        self.y = self.x +1 
    4554         
    4655    def test1D(self): 
    47         """ Test 1D model for a Cylinder  with vector as input""" 
     56        """ Test 1D model for a cylinder with vector as input""" 
     57         
    4858        answer = numpy.array([1.98860592e-04,7.03686335e-05,2.89144683e-06,2.04282827e-06]) 
    49         testvector= self.comp.run(self.x) 
    50         
     59 
     60        testvector= self.comp.evalDistribution(self.x) 
    5161        self.assertAlmostEqual(len(testvector),4) 
    5262        for i in xrange(len(answer)): 
    … …  
    5464        
    5565    def test1D_1(self): 
    56         """ Test 2D model for a Cylinder with scalar as input""" 
    57         self.assertAlmostEqual(self.comp.run(1.0),1.9886e-04, 4) 
     66        """ Test 2D model for a cylinder with scalar as input""" 
     67        self.assertAlmostEqual(self.comp.run(0.2), 0.041761386790780453, 4) 
    5868          
    5969    def test1D_2(self): 
    60         """ Test 2D model for a Cylinder for 2 scalar """ 
    61         self.assertAlmostEqual(self.comp.run([1.0, 1.3]), 56.3878e-06, 4) 
     70        """ Test 2D model of a cylinder """  
     71        self.comp.setParam('cyl_theta', 1.0) 
     72        self.comp.setParam('cyl_phi', 1.0) 
     73        self.assertAlmostEqual(self.comp.run([0.2, 2.5]),  
     74                               0.038176446608393366, 4) 
    6275         
    6376    def test1D_3(self): 
    64         """ Test 2D model for a Cylinder for 2 vectors as input """ 
    65         ans_input = numpy.zeros(len(self.x)) 
    66         temp_x = numpy.zeros(len(self.x)) 
    67         temp_y = numpy.zeros(len(self.y)) 
     77        """ Test 2D model for a cylinder for 2 vectors as input """ 
     78        x= numpy.reshape(self.x, [len(self.x),1]) 
     79        y= numpy.reshape(self.y, [1,len(self.y)]) 
     80        vect = self.comp.evalDistribution([x,y]) 
     81       
     82        self.assertAlmostEqual(vect[0][0],5.06121018e-08,4) 
     83        self.assertAlmostEqual(vect[len(self.x)-1][len(self.y)-1],2.5978e-11, 4) 
    6884         
    69         for i in xrange(len(self.x)): 
    70             qx = self.x[i] 
    71             qy = self.y[i] 
    72            
    73             temp_x[i]= qx*math.cos(qy) 
    74             temp_y[i]= qx*math.sin(qy) 
    75              
    76             value = math.sqrt(temp_x[i]*temp_x[i]+ temp_y[i]*temp_y[i] )#qx*qx +qy*qy) 
    77             ans_input[i]= value 
    78           
    79         vect_runXY_qx_qy = self.comp.runXY([temp_x, temp_y]) 
    80         vect_run_x_y = self.comp.run([self.x, self.y]) 
    8185         
    82         for i in xrange(len(vect_runXY_qx_qy)): 
    83             self.assertAlmostEqual(vect_runXY_qx_qy[i], vect_run_x_y[i]) 
    84              
    85         vect_run_x = self.comp.run(self.x) 
    86         vect_run_answer = self.comp.run(ans_input) 
    87          
    88         for i in xrange(len(vect_run_x )): 
    89             self.assertAlmostEqual(vect_run_x [i], vect_run_answer[i]) 
    90              
    91  
    9286 
    9387if __name__ == '__main__': 
Note: See TracChangeset for help on using the changeset viewer.