Changeset eba9885 in sasview for sansmodels/src

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:

• Schulz.py (added)
• c_extensions/CSchulz.c (added)
• c_extensions/c_models.c (modified) (1 diff)
• c_extensions/logNormal.c (added)
• c_extensions/logNormal.h (added)
• c_extensions/schulz.c (added)
• c_extensions/schulz.h (added)
• c_models/c_models.cpp (modified) (5 diffs)
• c_models/dispersion_visitor.cpp (modified) (2 diffs)
• c_models/dispersion_visitor.hh (modified) (1 diff)
• c_models/parameters.cpp (modified) (1 diff)
• c_models/parameters.hh (modified) (1 diff)
• dispersion_models.py (modified) (2 diffs)
• test/utest_dispersity.py (modified) (1 diff)
• test/utest_models.py (modified) (1 diff)
• test/utest_models_array.py (modified) (5 diffs)

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

@@ -52 +52 @@
         addDisperser(m);
         addCGaussian(m);
+        addCSchulz(m);
+        addCLogNormal(m);
         addCLorentzian(m);
         addCHollowCylinderModel(m);

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

@@ -48 +48 @@
         void addCGaussian(PyObject *module);
         void addCLorentzian(PyObject *module);
+        void addCLogNormal(PyObject *module);
+        void addCSchulz(PyObject *module);
 }
 
@@ -69 +71 @@
 
 /**
+ * Delete a lognormal dispersion model object
+ */
+void del_lognormal_dispersion(void *ptr){
+        LogNormalDispersion * disp = static_cast<LogNormalDispersion *>(ptr);
+        delete disp;
+        return;
+}
+
+/**
+ * Create a lognormal dispersion model as a python object
+ */
+PyObject * new_lognormal_dispersion(PyObject *, PyObject *args) {
+        LogNormalDispersion *disp = new LogNormalDispersion();
+        return PyCObject_FromVoidPtr(disp, del_lognormal_dispersion);
+}
+
+/**
  * Delete a gaussian dispersion model object
  */
@@ -84 +103 @@
         return PyCObject_FromVoidPtr(disp, del_gaussian_dispersion);
 }
+
+/**
+ * Delete a schulz dispersion model object
+ */
+void del_schulz_dispersion(void *ptr){
+        SchulzDispersion * disp = static_cast<SchulzDispersion *>(ptr);
+        delete disp;
+        return;
+}
+/**
+ * Create a schulz dispersion model as a python object
+ */
+PyObject * new_schulz_dispersion(PyObject *, PyObject *args) {
+        SchulzDispersion *disp = new SchulzDispersion();
+        return PyCObject_FromVoidPtr(disp, del_schulz_dispersion);
+}
+
 
 /**
@@ -147 +183 @@
         {"new_gaussian_model",   (PyCFunction)new_gaussian_dispersion, METH_VARARGS,
                   "Create a new GaussianDispersion object"},
+    {"new_lognormal_model",   (PyCFunction)new_lognormal_dispersion, METH_VARARGS,
+                  "Create a new LogNormalDispersion object"},
+    {"new_schulz_model",   (PyCFunction)new_schulz_dispersion, METH_VARARGS,
+                  "Create a new SchulzDispersion object"},
         {"new_array_model",      (PyCFunction)new_array_dispersion  , METH_VARARGS,
                   "Create a new ArrayDispersion object"},
@@ -194 +234 @@
         addDisperser(m);
         addCGaussian(m);
+        addCSchulz(m);
+        addCLogNormal(m);
         addCLorentzian(m);
         addCVesicleModel(m);
 
-
-}
+}

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

@@ -25 +25 @@
 
         PyDict_SetItemString(dict, "type",  Py_BuildValue("s", "gaussian"));
+    PyDict_SetItemString(dict, "npts",  Py_BuildValue("i", disp->npts));
+    PyDict_SetItemString(dict, "width", Py_BuildValue("d", disp->width));
+    PyDict_SetItemString(dict, "nsigmas", Py_BuildValue("i", disp->nsigmas));
+#endif
+}
+
+void DispersionVisitor:: lognormal_to_dict(void* dispersion, void* dictionary) {
+#ifndef __MODELS_STANDALONE__
+        LogNormalDispersion * disp = (LogNormalDispersion*)dispersion;
+        PyObject * dict = (PyObject*)dictionary;
+
+        PyDict_SetItemString(dict, "type",  Py_BuildValue("s", "lognormal"));
+    PyDict_SetItemString(dict, "npts",  Py_BuildValue("i", disp->npts));
+    PyDict_SetItemString(dict, "width", Py_BuildValue("d", disp->width));
+    PyDict_SetItemString(dict, "nsigmas", Py_BuildValue("i", disp->nsigmas));
+#endif
+}
+
+
+void DispersionVisitor:: schulz_to_dict(void* dispersion, void* dictionary) {
+#ifndef __MODELS_STANDALONE__
+        SchulzDispersion * disp = (SchulzDispersion*)dispersion;
+        PyObject * dict = (PyObject*)dictionary;
+
+        PyDict_SetItemString(dict, "type",  Py_BuildValue("s", "schulz"));
     PyDict_SetItemString(dict, "npts",  Py_BuildValue("i", disp->npts));
     PyDict_SetItemString(dict, "width", Py_BuildValue("d", disp->width));
@@ -61 +86 @@
 }
 
+void DispersionVisitor:: lognormal_from_dict(void* dispersion, void* dictionary) {
+#ifndef __MODELS_STANDALONE__
+        LogNormalDispersion * disp = (LogNormalDispersion*)dispersion;
+        PyObject * dict = (PyObject*)dictionary;
+
+        disp->npts    = PyInt_AsLong( PyDict_GetItemString(dict, "npts") );
+        disp->width   = PyFloat_AsDouble( PyDict_GetItemString(dict, "width") );
+        disp->nsigmas = PyFloat_AsDouble( PyDict_GetItemString(dict, "nsigmas") );
+#endif
+}
+void DispersionVisitor:: schulz_from_dict(void* dispersion, void* dictionary) {
+#ifndef __MODELS_STANDALONE__
+        SchulzDispersion * disp = (SchulzDispersion*)dispersion;
+        PyObject * dict = (PyObject*)dictionary;
+
+        disp->npts    = PyInt_AsLong( PyDict_GetItemString(dict, "npts") );
+        disp->width   = PyFloat_AsDouble( PyDict_GetItemString(dict, "width") );
+        disp->nsigmas = PyFloat_AsDouble( PyDict_GetItemString(dict, "nsigmas") );
+#endif
+}
+
 void DispersionVisitor:: array_from_dict(void* dispersion, void* dictionary) {}

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

@@ -19 +19 @@
         void dispersion_to_dict(void *, void *);
         void gaussian_to_dict(void *, void *);
+        void lognormal_to_dict(void *, void *);
+        void schulz_to_dict(void*, void *);
         void array_to_dict(void *, void *);
 
         void dispersion_from_dict(void*, void *);
         void gaussian_from_dict(void*, void *);
+        void lognormal_from_dict(void*, void *);
+        void schulz_from_dict(void*, void *);
         void array_from_dict(void*, void *);
 

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

@@ -136 +136 @@
         }
 }
+
+
+/**
+ * LogNormal dispersion
+ */
+
+LogNormalDispersion :: LogNormalDispersion() {
+        npts  = 1;
+        width = 0.0;
+        nsigmas = 2;
+};
+
+void LogNormalDispersion :: accept_as_source(DispersionVisitor* visitor, void* from, void* to) {
+        visitor->lognormal_to_dict(from, to);
+}
+void LogNormalDispersion :: accept_as_destination(DispersionVisitor* visitor, void* from, void* to) {
+        visitor->lognormal_from_dict(from, to);
+}
+
+double lognormal_weight(double mean, double sigma, double x) {
+        
+        double sigma2 = pow(sigma, 2);  
+        return 1/(x*sigma2) * exp( -pow((log(x) -mean), 2) / (2*sigma2));
+ 
+}
+
+/**
+ * Lognormal dispersion
+ * @param mean: mean value of the LogNormal
+ * @param sigma: standard deviation of the LogNormal
+ * @param x: value at which the LogNormal is evaluated
+ * @return: value of the LogNormal
+ */
+void LogNormalDispersion :: operator() (void *param, vector<WeightPoint> &weights){
+        // Check against zero width
+        if (width<=0) {
+                width = 0.0;
+                npts  = 1;
+                nsigmas = 3;
+        }
+
+        Parameter* par = (Parameter*)param;
+        double value = (*par)();
+
+        if (npts<2) {
+                weights.insert(weights.end(), WeightPoint(value, 1.0));
+        } else {
+                for(int i=0; i<npts; i++) {
+                        // We cover n(nsigmas) times sigmas on each side of the mean
+                        double val = value + width * (2.0*nsigmas*i/float(npts-1) - nsigmas);
+
+                        if ( ((*par).has_min==false || val>(*par).min)
+                          && ((*par).has_max==false || val<(*par).max)  ) {
+                                double _w = lognormal_weight(value, width, val);
+                                weights.insert(weights.end(), WeightPoint(val, _w));
+                        }
+                }
+        }
+}
+
+
+
+/**
+ * Schulz dispersion
+ */
+
+SchulzDispersion :: SchulzDispersion() {
+        npts  = 1;
+        width = 0.0;
+        nsigmas = 2;
+};
+
+void SchulzDispersion :: accept_as_source(DispersionVisitor* visitor, void* from, void* to) {
+        visitor->schulz_to_dict(from, to);
+}
+void SchulzDispersion :: accept_as_destination(DispersionVisitor* visitor, void* from, void* to) {
+        visitor->schulz_from_dict(from, to);
+}
+
+double schulz_weight(double mean, double sigma, double x) {
+        double vary, expo_value;
+    double z = pow(mean/ sigma, 2)-1;   
+        double R= x/mean;
+        double zz= z+1;
+        return  pow(zz,zz) * pow(R,z) * exp(-1*R*zz)/((mean) * tgamma(zz)) ;
+}
+
+/**
+ * Schulz dispersion
+ * @param mean: mean value of the Schulz
+ * @param sigma: standard deviation of the Schulz
+ * @param x: value at which the Schulz is evaluated
+ * @return: value of the Schulz
+ */
+void SchulzDispersion :: operator() (void *param, vector<WeightPoint> &weights){
+        // Check against zero width
+        if (width<=0) {
+                width = 0.0;
+                npts  = 1;
+                nsigmas = 3;
+        }
+
+        Parameter* par = (Parameter*)param;
+        double value = (*par)();
+
+        if (npts<2) {
+                weights.insert(weights.end(), WeightPoint(value, 1.0));
+        } else {
+                for(int i=0; i<npts; i++) {
+                        // We cover n(nsigmas) times sigmas on each side of the mean
+                        double val = value + width * (2.0*nsigmas*i/float(npts-1) - nsigmas);
+
+                        if ( ((*par).has_min==false || val>(*par).min)
+                          && ((*par).has_max==false || val<(*par).max)  ) {
+                                double _w = schulz_weight(value, width, val);
+                                weights.insert(weights.end(), WeightPoint(val, _w));
+                        }
+                }
+        }
+}
+
+
+
 
 /**

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

@@ -77 +77 @@
 
 /**
+ * Schulz dispersion model
+ */
+class SchulzDispersion: public DispersionModel {
+public:
+        /// Number of sigmas on each side of the mean
+        int nsigmas;
+
+        SchulzDispersion();
+        void operator()(void *, vector<WeightPoint>&);
+        void accept_as_source(DispersionVisitor*, void*, void*);
+        void accept_as_destination(DispersionVisitor*, void*, void*);
+};
+
+/**
+ * LogNormal dispersion model
+ */
+class LogNormalDispersion: public DispersionModel {
+public:
+        /// Number of sigmas on each side of the mean
+        int nsigmas;
+
+        LogNormalDispersion();
+        void operator()(void *, vector<WeightPoint>&);
+        void accept_as_source(DispersionVisitor*, void*, void*);
+        void accept_as_destination(DispersionVisitor*, void*, void*);
+};
+
+
+/**
  * Dispersion model based on arrays provided by the user
  */

sansmodels/src/sans/models/dispersion_models.py

@@ -74 +74 @@
             Set the weights of an array dispersion
         """
-        message = "set_weights is not available for GaussiantDispersion.\n"
+        message = "set_weights is not available for GaussianDispersion.\n"
+        message += "  Solution: Use an ArrayDispersion object"
+        raise "RuntimeError", message
+    
+class SchulzDispersion(DispersionModel):
+    """
+        Python bridge class for a dispersion model based 
+        on a Schulz distribution.
+    """
+    def __init__(self):
+        self.cdisp = c_models.new_schulz_model()
+        
+    def set_weights(self, values, weights):
+        """
+            Set the weights of an array dispersion
+        """
+        message = "set_weights is not available for SchulzDispersion.\n"
+        message += "  Solution: Use an ArrayDispersion object"
+        raise "RuntimeError", message
+    
+class LogNormalDispersion(DispersionModel):
+    """
+        Python bridge class for a dispersion model based 
+        on a Log Normal distribution.
+    """
+    def __init__(self):
+        self.cdisp = c_models.new_lognormal_model()
+        
+    def set_weights(self, values, weights):
+        """
+            Set the weights of an array dispersion
+        """
+        message = "set_weights is not available for LogNormalDispersion.\n"
         message += "  Solution: Use an ArrayDispersion object"
         raise "RuntimeError", message
@@ -100 +132 @@
         c_models.set_dispersion_weights(self.cdisp, values, weights)
  
-models = {GaussianDispersion:"GaussianModel", ArrayDispersion:"MyModel"}       
+models = {GaussianDispersion:"GaussianModel", ArrayDispersion:"MyModel",
+          SchulzDispersion: "Schulz", LogNormalDispersion: "LogNormal"}       
         

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

@@ -57 +57 @@
         
         new_model = self.model.clone()
+        print "gaussian",self.model.run(0.001)
         self.assertAlmostEqual(new_model.run(0.001), 4723.32213339, 3)
         self.assertAlmostEqual(new_model.runXY([0.001,0.001]), 4743.56, 2)
+        
+    def test_schulz_zero(self):
+        from sans.models.dispersion_models import SchulzDispersion
+        disp = SchulzDispersion()
+        self.model.set_dispersion('radius', disp)
+        self.model.dispersion['radius']['width'] = 5.0
+        #self.model.dispersion['radius']['width'] = 0.0
+        self.model.dispersion['radius']['npts'] = 100
+        #self.model.setParam('scale', 1.0)
+        self.model.setParam('scale', 10.0)
+        print "schulz",self.model.run(0.001), self.model.dispersion
+        self.assertAlmostEqual(self.model.run(0.001), 450.355, 3)
+        self.assertAlmostEqual(self.model.runXY([0.001,0.001]), 452.299, 3)
+        
+    def test_lognormal_zero(self):
+        from sans.models.dispersion_models import LogNormalDispersion
+        disp = LogNormalDispersion()
+        self.model.set_dispersion('radius', disp)
+        self.model.dispersion['radius']['width'] = 5.0
+        #self.model.dispersion['radius']['width'] = 0.0
+        self.model.dispersion['radius']['npts'] = 100
+        #self.model.setParam('scale', 1.0)
+        self.model.setParam('scale', 10.0)
+        print "model dispersion",self.model.dispersion
+        print "lognormal",self.model.run(0.001)
+        self.assertAlmostEqual(self.model.run(0.001), 450.355, 3)
+        self.assertAlmostEqual(self.model.runXY([0.001,0.001]), 452.299, 3)
         
     def test_gaussian_zero(self):

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

@@ -24 +24 @@
     def test1D(self):
         """ Test 1D model for a sphere """
-        self.assertAlmostEqual(self.comp.run(1.0), 56.3878, 4)
+        self.assertAlmostEqual(self.comp.run(1.0), 5.6387e-5, 4)
         
     def test1D_2(self):
         """ Test 2D model for a sphere """
-        self.assertAlmostEqual(self.comp.run([1.0, 1.3]), 56.3878, 4)
+        self.assertAlmostEqual(self.comp.run([1.0, 1.3]), 5.63878e-5, 4)
 
 class TestCyl(unittest.TestCase):

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

@@ -8 +8 @@
         
 class TestSphere(unittest.TestCase):
-    """ Unit tests for sphere model """
+    """ Unit tests for sphere model using evalDistribution function """
     
     def setUp(self):
@@ -18 +18 @@
     def test1D(self):
         """ Test 1D model for a sphere  with vector as input"""
-        answer=numpy.array([5.63877831e-05,2.57231782e-06,2.73704050e-07,2.54229069e-08])
+        answer = numpy.array([5.63877831e-05,2.57231782e-06,2.73704050e-07,2.54229069e-08])
        
-        testvector= self.comp.run(self.x)
+       
+        testvector= self.comp.evalDistribution(self.x)
        
         self.assertAlmostEqual(len(testvector),4)
@@ -33 +34 @@
         """ Test 2D model for a sphere for 2 scalar """
         self.assertAlmostEqual(self.comp.run([1.0, 1.3]), 56.3878e-06, 4)
-   
+        
+    def test1D_3(self):
+        """ Test 2D model for a Shpere for 2 vectors as input """
+        x= numpy.reshape(self.x, [len(self.x),1])
+        y= numpy.reshape(self.y, [1,len(self.y)])
+        vect = self.comp.evalDistribution([x,y])
+        self.assertAlmostEqual(vect[0][0],9.2985e-07, 4)
+        self.assertAlmostEqual(vect[len(self.x)-1][len(self.y)-1],1.3871e-08, 4)
+        
         
 class TestCylinder(unittest.TestCase):
-    """ Unit tests for sphere model """
+    """ Unit tests for Cylinder model using evalDistribution function """
     
     def setUp(self):
@@ -42 +51 @@
         self.comp = CylinderModel()
         self.x = numpy.array([1.0,2.0,3.0, 4.0])
-        self.y = numpy.array([1.0,2.0,3.0, 4.0])
+        self.y = self.x +1
         
     def test1D(self):
-        """ Test 1D model for a Cylinder  with vector as input"""
+        """ Test 1D model for a cylinder with vector as input"""
+        
         answer = numpy.array([1.98860592e-04,7.03686335e-05,2.89144683e-06,2.04282827e-06])
-        testvector= self.comp.run(self.x)
-       
+
+        testvector= self.comp.evalDistribution(self.x)
         self.assertAlmostEqual(len(testvector),4)
         for i in xrange(len(answer)):
@@ -54 +64 @@
        
     def test1D_1(self):
-        """ Test 2D model for a Cylinder with scalar as input"""
-        self.assertAlmostEqual(self.comp.run(1.0),1.9886e-04, 4)
+        """ Test 2D model for a cylinder  with scalar as input"""
+        self.assertAlmostEqual(self.comp.run(0.2), 0.041761386790780453, 4)
          
     def test1D_2(self):
-        """ Test 2D model for a Cylinder for 2 scalar """
-        self.assertAlmostEqual(self.comp.run([1.0, 1.3]), 56.3878e-06, 4)
+        """ Test 2D model of a cylinder """ 
+        self.comp.setParam('cyl_theta', 1.0)
+        self.comp.setParam('cyl_phi', 1.0)
+        self.assertAlmostEqual(self.comp.run([0.2, 2.5]), 
+                               0.038176446608393366, 4)
         
     def test1D_3(self):
-        """ Test 2D model for a Cylinder for 2 vectors as input """
-        ans_input = numpy.zeros(len(self.x))
-        temp_x = numpy.zeros(len(self.x))
-        temp_y = numpy.zeros(len(self.y))
+        """ Test 2D model for a cylinder for 2 vectors as input """
+        x= numpy.reshape(self.x, [len(self.x),1])
+        y= numpy.reshape(self.y, [1,len(self.y)])
+        vect = self.comp.evalDistribution([x,y])
+      
+        self.assertAlmostEqual(vect[0][0],5.06121018e-08,4)
+        self.assertAlmostEqual(vect[len(self.x)-1][len(self.y)-1],2.5978e-11, 4)
         
-        for i in xrange(len(self.x)):
-            qx = self.x[i]
-            qy = self.y[i]
-          
-            temp_x[i]= qx*math.cos(qy)
-            temp_y[i]= qx*math.sin(qy)
-            
-            value = math.sqrt(temp_x[i]*temp_x[i]+ temp_y[i]*temp_y[i] )#qx*qx +qy*qy)
-            ans_input[i]= value
-         
-        vect_runXY_qx_qy = self.comp.runXY([temp_x, temp_y])
-        vect_run_x_y = self.comp.run([self.x, self.y])
         
-        for i in xrange(len(vect_runXY_qx_qy)):
-            self.assertAlmostEqual(vect_runXY_qx_qy[i], vect_run_x_y[i])
-            
-        vect_run_x = self.comp.run(self.x)
-        vect_run_answer = self.comp.run(ans_input)
-        
-        for i in xrange(len(vect_run_x )):
-            self.assertAlmostEqual(vect_run_x [i], vect_run_answer[i])
-            
-
 
 if __name__ == '__main__':