source: sasview/sansmodels/src/sans/models/test/SmearList.py @ f99a6b1

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since f99a6b1 was ae3ce4e, checked in by Mathieu Doucet <doucetm@…>, 17 years ago

Moving sansmodels to trunk

  • Property mode set to 100644
File size: 6.6 KB
Line 
1
2from sans.models.ModelFactory import ModelFactory
3from sans.models.BaseComponent import BaseComponent
4import math
5
6class Smear:
7   
8    def __init__(self, model, param, sigma):
9        """
10            @param model: model to smear
11            @param param: parameter to smear
12            @param sigma: std deviations for parameter
13        """
14
15       
16        ## Model to smear
17        #self.model = factory.getModel(model_name)
18        self.model = model
19        ## Standard deviation of the smearing
20        self.sigmas = sigma
21        ## Parameter to smear
22        self.params = param
23        ## Nominal value of the smeared parameter
24        self.centers = []
25        ## Error on last evaluation
26        self.error = 0.0
27        ## Error flag
28        self.doErrors = False
29        for par in self.params:
30            self.centers.append(self.model.getParam(par))
31       
32    def smearParam(self, id, x):
33        """
34            @param x: input value
35        """
36        # from random import random
37        # If we exhausted the parameter array, simply evaluate
38        # the model
39        if id < len(self.params):
40            #print "smearing", self.params[id]
41           
42            # Average over Gaussian distribution (2 sigmas)
43            value_sum = 0.0
44            gauss_sum = 0.0
45           
46            min_value = self.centers[id] - 2*self.sigmas[id]
47            max_value = self.centers[id] + 2*self.sigmas[id]
48            n_pts = 25
49            step = (max_value - min_value)/(n_pts-1)
50            #print min_value, max_value, step, max_value-min_value
51            if step == 0.0:
52                return self.smearParam(id+1,x)
53           
54            # Gaussian function used to weigh points
55            #gaussian = ModelFactory().getModel("Gaussian")
56            gaussian = Gaussian()
57            gaussian.setParam('sigma', self.sigmas[id])
58            gaussian.setParam('mean', self.centers[id])
59                   
60            # Compute average
61            prev_value = None
62            error_sys = 0.0
63            for i in range(n_pts):
64                # Set the parameter value           
65                value = min_value + i*step
66                # value = random()*4.0*self.sigmas[id] + min_value
67                # print value
68                self.model.setParam(self.params[id], value)
69                gauss_value = gaussian.run(value)
70                #gauss_value = 1.0
71                #if id==0: print value, gauss_value
72                func_value, error_1 = self.smearParam(id+1, x)
73                if self.doErrors:
74                    if not prev_value == None:
75                        error_sys += (func_value-prev_value)*(func_value-prev_value)/4.0
76                    prev_value = func_value
77
78                value_sum += gauss_value * func_value
79                gauss_sum += gauss_value
80               
81            #print "Error", math.sqrt(error)
82            return value_sum/gauss_sum, math.sqrt(error_sys)
83       
84        else:
85            return self.model.run(x), 0.0
86       
87    def run(self, x):
88        """
89            @param x: input
90        """
91       
92        # Go through the list of parameters
93        n_par = len(self.params)
94       
95        # Check array lengths
96        if not len(self.centers) == n_par or\
97            not len(self.sigmas) == n_par:
98            raise ValueError, "Incompatible array lengths"
99       
100        # Smear first parameters
101        if n_par > 0:
102            value, error = self.smearParam(0, x)
103            self.error = error
104           
105            # Put back original values
106            for i in range(len(self.centers)):
107                self.model.setParam(self.params[i], self.centers[i])
108           
109           
110            return value
111       
112class Gaussian(BaseComponent):
113    """ Gaussian function """
114   
115    def __init__(self):
116        """ Initialization """
117        BaseComponent.__init__(self)
118       
119        ## Name of the model
120        self.name = "Gaussian"
121        ## Scale factor
122        self.params['scale']  = 1.0
123        ## Mean value
124        self.params['mean']   = 0.0
125        ## Standard deviation
126        self.params['sigma']  = 1.0
127        ## Internal log
128        self.log = {}
129        return
130   
131    def run(self, x=0):
132        """ Evaluate the function
133            @param x: input q or [q,phi]
134            @return: scattering function
135        """
136        if(type(x)==type([])):
137            # vector input
138            if(len(x)==2):
139                return self.analytical_2D(x)
140            else:
141                raise ValueError, "Gaussian takes a scalar or a 2D point"
142        else:
143            return self.analytical_1D(x)
144
145    def analytical_2D(self, x):
146        """ Evaluate 2D model
147            @param x: input [q,phi]
148            @return: scattering function
149        """
150       
151        # 2D sphere is the same as 1D sphere
152        return self.analytical_1D(x[0])
153
154    def analytical_1D(self, x):
155        """ Evaluate 1D model
156            @param x: input q
157            @return: scattering function
158        """
159        vary = x-self.params['mean']
160        expo_value = -vary*vary/(2*self.params['sigma']*self.params['sigma'])
161        value = self.params['scale'] *  math.exp(expo_value)
162       
163        return value
164       
165# main
166if __name__ == '__main__':
167    import math
168    if True:
169        cyl = ModelFactory().getModel("CylinderModel")
170        cyl.setParam('length', 2000.0)
171        cyl.setParam('cyl_phi', .45)
172        cyl.setParam('cyl_theta', 2.)
173        #app = Smear(cyl,['cyl_theta', 'cyl_phi'], [.01, .01])
174        #app = Smear(cyl,['cyl_theta', 'cyl_phi'], [2.0, 2.0])
175        #app = Smear(cyl,['cyl_theta'], [2.5])
176        app = Smear(cyl, ['cyl_phi', 'cyl_theta'], [math.pi, math.pi])
177        #app = Smear(cyl,['scale', 'scale'],[.5,1.0])
178        val_no = cyl.run([.001, 1.0])
179        print "Cylinder (no smear) f(.1):", val_no
180        val_sm =  app.run([.001, 1.0])
181        print "Cylinder (smear)    f(.1):", val_sm, app.error, app.error/val_sm
182        print "Cylinder (1D)    f(.1):", cyl.run(.001), cyl.run(.001)/val_sm
183     
184     
185    print "--------------"
186    sinsin = ModelFactory().getModel("SinSin")
187    sinsin.setParam('a', -1.7)
188    sinsin.setParam('b', 5.45)
189    #app = Smear(sinsin, ['A'], [3.1])
190    #app = Smear(sinsin, ['A', 'b'], [math.pi/2.0, .0])
191    app = Smear(sinsin, ['A'], [math.pi/2.0])
192    #app = Smear(cyl,['scale', 'scale'],[.5,1.0])
193    val_no = sinsin.run(1.0)
194    print "SinSin (no smear) :", val_no
195    val_sm =  app.run(1.0)
196    print "SinSin (smear)    : %g +- %g" % (val_sm, app.error)
197     
Note: See TracBrowser for help on using the repository browser.