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

Last change on this file since 46ca1f4 was fa81e94, checked in by Piotr Rozyczko <rozyczko@…>, 7 years ago

Initial commit of the P(r) inversion perspective.
Code merged from Jeff Krzywon's ESS_GUI_Pr branch.
Also, minor 2to3 mods to sascalc/sasgui to enble error free setup.

  • Property mode set to 100755
File size: 11.6 KB
Line 
1"""
2    Calculation thread for modeling
3"""
4
5import time
6import numpy as np
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 is None:
56            self.qmin = 0
57        if self.qmax is None:
58            if self.data is not 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 = np.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 & np.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            # Calculate smeared Intensity
85            #(by Gaussian averaging): DataLoader/smearing2d/Smearer2D()
86            value = fn.get_value()
87        else:
88            # calculation w/o smearing
89            value = self.model.evalDistribution([
90                self.data.qx_data[index_model],
91                self.data.qy_data[index_model]
92            ])
93        output = np.zeros(len(self.data.qx_data))
94        # output default is None
95        # This method is to distinguish between masked
96        #point(nan) and data point = 0.
97        output = output / output
98        # set value for self.mask==True, else still None to Plottools
99        output[index_model] = value
100        elapsed = time.time() - self.starttime
101        #self.complete(image=output,
102        #               data=self.data,
103        #               page_id=self.page_id,
104        #               model=self.model,
105        #               state=self.state,
106        #               toggle_mode_on=self.toggle_mode_on,
107        #               elapsed=elapsed,
108        #               index=index_model,
109        #               fid=self.fid,
110        #               qmin=self.qmin,
111        #               qmax=self.qmax,
112        #               weight=self.weight,
113        #               #qstep=self.qstep,
114        #               update_chisqr=self.update_chisqr,
115        #               source=self.source)
116        return (output,
117                self.data,
118                self.page_id,
119                self.model,
120                self.state,
121                self.toggle_mode_on,
122                elapsed,
123                index_model,
124                self.fid,
125                self.qmin,
126                self.qmax,
127                self.weight,
128                self.update_chisqr,
129                self.source)
130
131
132class Calc1D(CalcThread):
133    """
134    Compute 1D data
135    """
136    def __init__(self, model,
137                 page_id,
138                 data,
139                 fid=None,
140                 qmin=None,
141                 qmax=None,
142                 weight=None,
143                 smearer=None,
144                 toggle_mode_on=False,
145                 state=None,
146                 completefn=None,
147                 update_chisqr=True,
148                 source='model',
149                 updatefn=None,
150                 yieldtime=0.01,
151                 worktime=0.01,
152                 exception_handler=None,
153                 ):
154        """
155        """
156        CalcThread.__init__(self, completefn, updatefn, yieldtime, worktime,
157                            exception_handler=exception_handler)
158        self.fid = fid
159        self.data = data
160        self.qmin = qmin
161        self.qmax = qmax
162        self.model = model
163        self.weight = weight
164        self.toggle_mode_on = toggle_mode_on
165        self.state = state
166        self.page_id = page_id
167        self.smearer = smearer
168        self.starttime = 0
169        self.update_chisqr = update_chisqr
170        self.source = source
171        self.out = None
172        self.index = None
173
174    def compute(self):
175        """
176        Compute model 1d value given qmin , qmax , x value
177        """
178        self.starttime = time.time()
179        output = np.zeros((len(self.data.x)))
180        index = (self.qmin <= self.data.x) & (self.data.x <= self.qmax)
181
182        # If we use a smearer, also return the unsmeared model
183        unsmeared_output = None
184        unsmeared_data = None
185        unsmeared_error = None
186        ##smearer the ouput of the plot
187        if self.smearer is not None:
188            first_bin, last_bin = self.smearer.get_bin_range(self.qmin,
189                                                             self.qmax)
190            mask = self.data.x[first_bin:last_bin+1]
191            unsmeared_output = np.zeros((len(self.data.x)))
192            unsmeared_output[first_bin:last_bin+1] = self.model.evalDistribution(mask)
193            self.smearer.model = self.model
194            output = self.smearer(unsmeared_output, first_bin, last_bin)
195
196            # Rescale data to unsmeared model
197            # Check that the arrays are compatible. If we only have a model but no data,
198            # the length of data.y will be zero.
199            if isinstance(self.data.y, np.ndarray) and output.shape == self.data.y.shape:
200                unsmeared_data = np.zeros((len(self.data.x)))
201                unsmeared_error = np.zeros((len(self.data.x)))
202                unsmeared_data[first_bin:last_bin+1] = self.data.y[first_bin:last_bin+1]\
203                                                        * unsmeared_output[first_bin:last_bin+1]\
204                                                        / output[first_bin:last_bin+1]
205                unsmeared_error[first_bin:last_bin+1] = self.data.dy[first_bin:last_bin+1]\
206                                                        * unsmeared_output[first_bin:last_bin+1]\
207                                                        / output[first_bin:last_bin+1]
208                unsmeared_output=unsmeared_output[index]
209                unsmeared_data=unsmeared_data[index]
210                unsmeared_error=unsmeared_error
211        else:
212            output[index] = self.model.evalDistribution(self.data.x[index])
213
214        sq_values = None
215        pq_values = None
216        s_model = None
217        p_model = None
218        if isinstance(self.model, MultiplicationModel):
219            s_model = self.model.s_model
220            p_model = self.model.p_model
221        elif hasattr(self.model, "get_composition_models"):
222            p_model, s_model = self.model.get_composition_models()
223
224        if p_model is not None and s_model is not None:
225            sq_values = np.zeros((len(self.data.x)))
226            pq_values = np.zeros((len(self.data.x)))
227            sq_values[index] = s_model.evalDistribution(self.data.x[index])
228            pq_values[index] = p_model.evalDistribution(self.data.x[index])
229
230        elapsed = time.time() - self.starttime
231
232        #self.complete(x=self.data.x[index], y=output[index],
233        #              page_id=self.page_id,
234        #              state=self.state,
235        #              weight=self.weight,
236        #              fid=self.fid,
237        #              toggle_mode_on=self.toggle_mode_on,
238        #              elapsed=elapsed, index=index, model=self.model,
239        #              data=self.data,
240        #              update_chisqr=self.update_chisqr,
241        #              source=self.source)
242        return (self.data.x[index], output[index],
243                self.page_id,
244                self.state,
245                self.weight,
246                self.fid,
247                self.toggle_mode_on,
248                elapsed, index, self.model,
249                self.data,
250                self.update_chisqr,
251                self.source)
252
253        # TODO: as of 4.1, the output contains more items:
254        # unsmeared_* and pq_model/sq_model
255        # Need to add these too
256
257        #self.complete(x=self.data.x[index], y=output[index],
258        #              page_id=self.page_id,
259        #              state=self.state,
260        #              weight=self.weight,
261        #              fid=self.fid,
262        #              toggle_mode_on=self.toggle_mode_on,
263        #              elapsed=elapsed, index=index, model=self.model,
264        #              data=self.data,
265        #              update_chisqr=self.update_chisqr,
266        #              source=self.source,
267        #              unsmeared_model=unsmeared_output,
268        #              unsmeared_data=unsmeared_data,
269        #              unsmeared_error=unsmeared_error,
270        #              pq_model=pq_values,
271        #              sq_model=sq_values)
272
273    def results(self):
274        """
275        Send resuts of the computation
276        """
277        return [self.out, self.index]
278
279"""
280Example: ::
281
282    class CalcCommandline:
283        def __init__(self, n=20000):
284            #print thread.get_ident()
285            from sas.models.CylinderModel import CylinderModel
286
287            model = CylinderModel()
288
289
290            print model.runXY([0.01, 0.02])
291
292            qmax = 0.01
293            qstep = 0.0001
294            self.done = False
295
296            x = numpy.arange(-qmax, qmax+qstep*0.01, qstep)
297            y = numpy.arange(-qmax, qmax+qstep*0.01, qstep)
298
299
300            calc_thread_2D = Calc2D(x, y, None, model.clone(),None,
301                                    -qmax, qmax,qstep,
302                                            completefn=self.complete,
303                                            updatefn=self.update ,
304                                            yieldtime=0.0)
305
306            calc_thread_2D.queue()
307            calc_thread_2D.ready(2.5)
308
309            while not self.done:
310                time.sleep(1)
311
312        def update(self,output):
313            print "update"
314
315        def complete(self, image, data, model, elapsed, qmin, qmax,index, qstep ):
316            print "complete"
317            self.done = True
318
319    if __name__ == "__main__":
320        CalcCommandline()
321"""
Note: See TracBrowser for help on using the repository browser.