source: sasview/src/sas/qtgui/Perspectives/Inversion/InversionLogic.py @ 9c0ce68

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalc
Last change on this file since 9c0ce68 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: 8.1 KB
Line 
1import math
2import pylab
3import numpy as np
4
5from sas.qtgui.Plotting.PlotterData import Data1D
6
7PR_FIT_LABEL = r"$P_{fit}(r)$"
8PR_LOADED_LABEL = r"$P_{loaded}(r)$"
9IQ_DATA_LABEL = r"$I_{obs}(q)$"
10IQ_FIT_LABEL = r"$I_{fit}(q)$"
11IQ_SMEARED_LABEL = r"$I_{smeared}(q)$"
12GROUP_ID_IQ_DATA = r"$I_{obs}(q)$"
13GROUP_ID_PR_FIT = r"$P_{fit}(r)$"
14
15
16class InversionLogic(object):
17    """
18    All the data-related logic. This class deals exclusively with Data1D/2D
19    No QStandardModelIndex here.
20    """
21
22    # TODO: Add way to change this value
23    _pr_n_pts = 51
24
25    def __init__(self, data=None):
26        self._data = data
27        self.data_is_loaded = False
28        if data is not None:
29            self.data_is_loaded = True
30
31    @property
32    def data(self):
33        return self._data
34
35    @data.setter
36    def data(self, value):
37        """ data setter """
38        self._data = value
39        self.data_is_loaded = True
40
41    def isLoadedData(self):
42        """ accessor """
43        return self.data_is_loaded
44
45    def new1DPlot(self, out, pr, q=None):
46        """
47        Create a new 1D data instance based on fitting results
48        """
49
50        qtemp = pr.x
51        if q is not None:
52            qtemp = q
53
54        # Make a plot
55        maxq = max(qtemp)
56
57        minq = min(qtemp)
58
59        # Check for user min/max
60        if pr.q_min is not None and maxq >= pr.q_min >= minq:
61            minq = pr.q_min
62        if pr.q_max is not None and maxq >= pr.q_max >= minq:
63            maxq = pr.q_max
64
65        x = pylab.arange(minq, maxq, maxq / 301.0)
66        y = np.zeros(len(x))
67        err = np.zeros(len(x))
68        for i in range(len(x)):
69            value = pr.iq(out, x[i])
70            y[i] = value
71            try:
72                err[i] = math.sqrt(math.fabs(value))
73            except:
74                err[i] = 1.0
75                print(("Error getting error", value, x[i]))
76
77        new_plot = Data1D(x, y)
78        new_plot.name = IQ_FIT_LABEL
79        new_plot.xaxis("\\rm{Q}", 'A^{-1}')
80        new_plot.yaxis("\\rm{Intensity} ", "cm^{-1}")
81        title = "I(q)"
82        new_plot.title = title
83
84        # If we have a group ID, use it
85        if 'plot_group_id' in pr.info:
86            new_plot.group_id = pr.info["plot_group_id"]
87        new_plot.id = IQ_FIT_LABEL
88
89        # If we have used slit smearing, plot the smeared I(q) too
90        if pr.slit_width > 0 or pr.slit_height > 0:
91            x = pylab.arange(minq, maxq, maxq / 301.0)
92            y = np.zeros(len(x))
93            err = np.zeros(len(x))
94            for i in range(len(x)):
95                value = pr.iq_smeared(pr.out, x[i])
96                y[i] = value
97                try:
98                    err[i] = math.sqrt(math.fabs(value))
99                except:
100                    err[i] = 1.0
101                    print(("Error getting error", value, x[i]))
102
103            new_plot = Data1D(x, y)
104            new_plot.name = IQ_SMEARED_LABEL
105            new_plot.xaxis("\\rm{Q}", 'A^{-1}')
106            new_plot.yaxis("\\rm{Intensity} ", "cm^{-1}")
107            # If we have a group ID, use it
108            if 'plot_group_id' in pr.info:
109                new_plot.group_id = pr.info["plot_group_id"]
110            new_plot.id = IQ_SMEARED_LABEL
111            new_plot.title = title
112
113        return new_plot
114
115    def update1DPlot(self, plot, out, pr, q=None):
116        """
117        Create a new 1D data instance based on fitting results
118        """
119
120        qtemp = pr.x
121        if q is not None:
122            qtemp = q
123
124        # Make a plot
125        maxq = max(qtemp)
126
127        minq = min(qtemp)
128
129        # Check for user min/max
130        if pr.q_min is not None and maxq >= pr.q_min >= minq:
131            minq = pr.q_min
132        if pr.q_max is not None and maxq >= pr.q_max >= minq:
133            maxq = pr.q_max
134
135        x = pylab.arange(minq, maxq, maxq / 301.0)
136        y = np.zeros(len(x))
137        err = np.zeros(len(x))
138        for i in range(len(x)):
139            value = pr.iq(out, x[i])
140            y[i] = value
141            try:
142                err[i] = math.sqrt(math.fabs(value))
143            except:
144                err[i] = 1.0
145                print(("Error getting error", value, x[i]))
146
147        plot.x = x
148        plot.y = y
149
150        # If we have used slit smearing, plot the smeared I(q) too
151        if pr.slit_width > 0 or pr.slit_height > 0:
152            x = pylab.arange(minq, maxq, maxq / 301.0)
153            y = np.zeros(len(x))
154            err = np.zeros(len(x))
155            for i in range(len(x)):
156                value = pr.iq_smeared(pr.out, x[i])
157                y[i] = value
158                try:
159                    err[i] = math.sqrt(math.fabs(value))
160                except:
161                    err[i] = 1.0
162                    print(("Error getting error", value, x[i]))
163
164            plot.x = x
165            plot.y = y
166
167        return plot
168
169    def newPRPlot(self, out, pr, cov=None):
170        """
171        """
172        # Show P(r)
173        x = pylab.arange(0.0, pr.d_max, pr.d_max / self._pr_n_pts)
174
175        y = np.zeros(len(x))
176        dy = np.zeros(len(x))
177
178        total = 0.0
179        pmax = 0.0
180        cov2 = np.ascontiguousarray(cov)
181
182        for i in range(len(x)):
183            if cov2 is None:
184                value = pr.pr(out, x[i])
185            else:
186                (value, dy[i]) = pr.pr_err(out, cov2, x[i])
187            total += value * pr.d_max / len(x)
188
189            # keep track of the maximum P(r) value
190            if value > pmax:
191                pmax = value
192
193            y[i] = value
194
195        # if self._normalize_output == True:
196        #     y = y / total
197        #     dy = dy / total
198        # elif self._scale_output_unity == True:
199        #     y = y / pmax
200        #     dy = dy / pmax
201
202        if cov2 is None:
203            new_plot = Data1D(x, y)
204        else:
205            new_plot = Data1D(x, y, dy=dy)
206        new_plot.name = PR_FIT_LABEL
207        new_plot.xaxis("\\rm{r}", 'A')
208        new_plot.yaxis("\\rm{P(r)} ", "cm^{-3}")
209        new_plot.title = "P(r) fit"
210        new_plot.id = PR_FIT_LABEL
211        # Make sure that the plot is linear
212        new_plot.xtransform = "x"
213        new_plot.ytransform = "y"
214        new_plot.group_id = GROUP_ID_PR_FIT
215
216        return new_plot
217
218    def updatePRPlot(self, plot, out, pr, cov=None):
219        x = pylab.arange(0.0, pr.d_max, pr.d_max / self._pr_n_pts)
220
221        y = np.zeros(len(x))
222        dy = np.zeros(len(x))
223
224        total = 0.0
225        pmax = 0.0
226        cov2 = np.ascontiguousarray(cov)
227
228        for i in range(len(x)):
229            if cov2 is None:
230                value = pr.pr(out, x[i])
231            else:
232                (value, dy[i]) = pr.pr_err(out, cov2, x[i])
233            total += value * pr.d_max / len(x)
234
235            # keep track of the maximum P(r) value
236            if value > pmax:
237                pmax = value
238
239            y[i] = value
240
241        # if self._normalize_output == True:
242        #     y = y / total
243        #     dy = dy / total
244        # elif self._scale_output_unity == True:
245        #     y = y / pmax
246        #     dy = dy / pmax
247        plot.x = x
248        plot.y = y
249
250        if cov2 is not None:
251            plot.dy = dy
252
253        return plot
254
255    def computeDataRange(self):
256        """
257        Wrapper for calculating the data range based on local dataset
258        """
259        return self.computeRangeFromData(self.data)
260
261    def computeRangeFromData(self, data):
262        """
263        Compute the minimum and the maximum range of the data
264        return the npts contains in data
265        """
266        qmin, qmax = None, None
267        if isinstance(data, Data1D):
268            try:
269                qmin = min(data.x)
270                qmax = max(data.x)
271            except (ValueError, TypeError):
272                msg = "Unable to find min/max/length of \n data named %s" % \
273                            self.data.filename
274                raise ValueError(msg)
275
276        else:
277            qmin = 0
278            try:
279                x = max(np.fabs(data.xmin), np.fabs(data.xmax))
280                y = max(np.fabs(data.ymin), np.fabs(data.ymax))
281            except (ValueError, TypeError):
282                msg = "Unable to find min/max of \n data named %s" % \
283                            self.data.filename
284                raise ValueError(msg)
285            qmax = np.sqrt(x * x + y * y)
286        return qmin, qmax
Note: See TracBrowser for help on using the repository browser.