source: sasview/src/sas/sasgui/perspectives/fitting/model_thread.py @ d7d0182

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.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since d7d0182 was 286c757, checked in by mathieu, 8 years ago

Make sure we can compute models w/o data. Fixes #714
  • Property mode set to 100644
File size: 10.0 KB
Line 
1"""
2    Calculation thread for modeling
3"""
4
5import time
6import numpy
7import math
8from sas.sascalc.data_util.calcthread import CalcThread
9from sas.sascalc.fit.MultiplicationModel import MultiplicationModel
10
11class Calc2D(CalcThread):
12    """
13    Compute 2D model
14    This calculation assumes a 2-fold symmetry of the model
15    where points are computed for one half of the detector
16    and I(qx, qy) = I(-qx, -qy) is assumed.
17    """
18    def __init__(self, data, model, smearer, qmin, qmax, page_id,
19                 state=None,
20                 weight=None,
21                 fid=None,
22                 toggle_mode_on=False,
23                 completefn=None,
24                 updatefn=None,
25                 update_chisqr=True,
26                 source='model',
27                 yieldtime=0.04,
28                 worktime=0.04,
29                 exception_handler=None,
30                 ):
31        CalcThread.__init__(self, completefn, updatefn, yieldtime, worktime,
32                            exception_handler=exception_handler)
33        self.qmin = qmin
34        self.qmax = qmax
35        self.weight = weight
36        self.fid = fid
37        #self.qstep = qstep
38        self.toggle_mode_on = toggle_mode_on
39        self.data = data
40        self.page_id = page_id
41        self.state = None
42        # the model on to calculate
43        self.model = model
44        self.smearer = smearer
45        self.starttime = 0
46        self.update_chisqr = update_chisqr
47        self.source = source
48
49    def compute(self):
50        """
51        Compute the data given a model function
52        """
53        self.starttime = time.time()
54        # Determine appropriate q range
55        if self.qmin == None:
56            self.qmin = 0
57        if self.qmax == None:
58            if self.data != None:
59                newx = math.pow(max(math.fabs(self.data.xmax),
60                                   math.fabs(self.data.xmin)), 2)
61                newy = math.pow(max(math.fabs(self.data.ymax),
62                                   math.fabs(self.data.ymin)), 2)
63                self.qmax = math.sqrt(newx + newy)
64
65        if self.data is None:
66            msg = "Compute Calc2D receive data = %s.\n" % str(self.data)
67            raise ValueError, msg
68
69        # Define matrix where data will be plotted
70        radius = numpy.sqrt((self.data.qx_data * self.data.qx_data) + \
71                    (self.data.qy_data * self.data.qy_data))
72
73        # For theory, qmax is based on 1d qmax
74        # so that must be mulitified by sqrt(2) to get actual max for 2d
75        index_model = (self.qmin <= radius) & (radius <= self.qmax)
76        index_model = index_model & self.data.mask
77        index_model = index_model & numpy.isfinite(self.data.data)
78
79        if self.smearer is not None:
80            # Set smearer w/ data, model and index.
81            fn = self.smearer
82            fn.set_model(self.model)
83            fn.set_index(index_model)
84            # Get necessary data from self.data and set the data for smearing
85            fn.get_data()
86            # Calculate smeared Intensity
87            #(by Gaussian averaging): DataLoader/smearing2d/Smearer2D()
88            value = fn.get_value()
89        else:
90            # calculation w/o smearing
91            value = self.model.evalDistribution(\
92                [self.data.qx_data[index_model],
93                 self.data.qy_data[index_model]])
94        output = numpy.zeros(len(self.data.qx_data))
95        # output default is None
96        # This method is to distinguish between masked
97        #point(nan) and data point = 0.
98        output = output / output
99        # set value for self.mask==True, else still None to Plottools
100        output[index_model] = value
101        elapsed = time.time() - self.starttime
102        self.complete(image=output,
103                       data=self.data,
104                       page_id=self.page_id,
105                       model=self.model,
106                       state=self.state,
107                       toggle_mode_on=self.toggle_mode_on,
108                       elapsed=elapsed,
109                       index=index_model,
110                       fid=self.fid,
111                       qmin=self.qmin,
112                       qmax=self.qmax,
113                       weight=self.weight,
114                       #qstep=self.qstep,
115                       update_chisqr=self.update_chisqr,
116                       source=self.source)
117
118
119class Calc1D(CalcThread):
120    """
121    Compute 1D data
122    """
123    def __init__(self, model,
124                 page_id,
125                 data,
126                 fid=None,
127                 qmin=None,
128                 qmax=None,
129                 weight=None,
130                 smearer=None,
131                 toggle_mode_on=False,
132                 state=None,
133                 completefn=None,
134                 update_chisqr=True,
135                 source='model',
136                 updatefn=None,
137                 yieldtime=0.01,
138                 worktime=0.01,
139                 exception_handler=None,
140                 ):
141        """
142        """
143        CalcThread.__init__(self, completefn, updatefn, yieldtime, worktime,
144                            exception_handler=exception_handler)
145        self.fid = fid
146        self.data = data
147        self.qmin = qmin
148        self.qmax = qmax
149        self.model = model
150        self.weight = weight
151        self.toggle_mode_on = toggle_mode_on
152        self.state = state
153        self.page_id = page_id
154        self.smearer = smearer
155        self.starttime = 0
156        self.update_chisqr = update_chisqr
157        self.source = source
158        self.out = None
159        self.index = None
160
161    def compute(self):
162        """
163        Compute model 1d value given qmin , qmax , x value
164        """
165        self.starttime = time.time()
166        output = numpy.zeros((len(self.data.x)))
167        index = (self.qmin <= self.data.x) & (self.data.x <= self.qmax)
168
169        # If we use a smearer, also return the unsmeared model
170        unsmeared_output = None
171        unsmeared_data = None
172        unsmeared_error = None
173        ##smearer the ouput of the plot
174        if self.smearer is not None:
175            first_bin, last_bin = self.smearer.get_bin_range(self.qmin,
176                                                             self.qmax)
177            mask = self.data.x[first_bin:last_bin+1]
178            unsmeared_output = numpy.zeros((len(self.data.x)))
179            unsmeared_output[first_bin:last_bin+1] = self.model.evalDistribution(mask)
180            output = self.smearer(unsmeared_output, first_bin, last_bin)
181
182            # Rescale data to unsmeared model
183            # Check that the arrays are compatible. If we only have a model but no data,
184            # the length of data.y will be zero.
185            if isinstance(self.data.y, numpy.ndarray) and output.shape == self.data.y.shape:
186                unsmeared_data = numpy.zeros((len(self.data.x)))
187                unsmeared_error = numpy.zeros((len(self.data.x)))
188                unsmeared_data[first_bin:last_bin+1] = self.data.y[first_bin:last_bin+1]\
189                                                        * unsmeared_output[first_bin:last_bin+1]\
190                                                        / output[first_bin:last_bin+1]
191                unsmeared_error[first_bin:last_bin+1] = self.data.dy[first_bin:last_bin+1]\
192                                                        * unsmeared_output[first_bin:last_bin+1]\
193                                                        / output[first_bin:last_bin+1]
194                unsmeared_output=unsmeared_output[index]
195                unsmeared_data=unsmeared_data[index]
196                unsmeared_error=unsmeared_error
197        else:
198            output[index] = self.model.evalDistribution(self.data.x[index])
199
200        sq_model = None
201        pq_model = None
202        if isinstance(self.model, MultiplicationModel):
203            sq_model = numpy.zeros((len(self.data.x)))
204            pq_model = numpy.zeros((len(self.data.x)))
205            sq_model[index] = self.model.s_model.evalDistribution(self.data.x[index])
206            pq_model[index] = self.model.p_model.evalDistribution(self.data.x[index])
207
208        elapsed = time.time() - self.starttime
209
210        self.complete(x=self.data.x[index], y=output[index],
211                      page_id=self.page_id,
212                      state=self.state,
213                      weight=self.weight,
214                      fid=self.fid,
215                      toggle_mode_on=self.toggle_mode_on,
216                      elapsed=elapsed, index=index, model=self.model,
217                      data=self.data,
218                      update_chisqr=self.update_chisqr,
219                      source=self.source,
220                      unsmeared_model=unsmeared_output,
221                      unsmeared_data=unsmeared_data,
222                      unsmeared_error=unsmeared_error,
223                      pq_model=pq_model,
224                      sq_model=sq_model)
225
226    def results(self):
227        """
228        Send resuts of the computation
229        """
230        return [self.out, self.index]
231
232"""
233Example: ::
234
235    class CalcCommandline:
236        def __init__(self, n=20000):
237            #print thread.get_ident()
238            from sas.models.CylinderModel import CylinderModel
239
240            model = CylinderModel()
241
242
243            print model.runXY([0.01, 0.02])
244
245            qmax = 0.01
246            qstep = 0.0001
247            self.done = False
248
249            x = numpy.arange(-qmax, qmax+qstep*0.01, qstep)
250            y = numpy.arange(-qmax, qmax+qstep*0.01, qstep)
251
252
253            calc_thread_2D = Calc2D(x, y, None, model.clone(),None,
254                                    -qmax, qmax,qstep,
255                                            completefn=self.complete,
256                                            updatefn=self.update ,
257                                            yieldtime=0.0)
258
259            calc_thread_2D.queue()
260            calc_thread_2D.ready(2.5)
261
262            while not self.done:
263                time.sleep(1)
264
265        def update(self,output):
266            print "update"
267
268        def complete(self, image, data, model, elapsed, qmin, qmax,index, qstep ):
269            print "complete"
270            self.done = True
271
272    if __name__ == "__main__":
273        CalcCommandline()
274"""
Note: See TracBrowser for help on using the repository browser.