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

ESS_GUIESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_opencl
Last change on this file since dc3b34e was a74dd016, checked in by wojciech, 5 years ago

Fixing faling P(r) on files without added errors

  • Property mode set to 100644
File size: 5.7 KB
Line 
1import math
2import logging
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)$"
14PR_PLOT_PTS = 51
15
16logger = logging.getLogger(__name__)
17
18
19class InversionLogic(object):
20    """
21    All the data-related logic. This class deals exclusively with Data1D/2D
22    No QStandardModelIndex here.
23    """
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 = (self._data is not None)
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 = np.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                logger.log(("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 = np.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                    logger.log(("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        new_plot.symbol = 'Line'
114        new_plot.hide_error = True
115
116        return new_plot
117
118    def newPRPlot(self, out, pr, cov=None):
119        """
120        """
121        # Show P(r)
122        x = np.arange(0.0, pr.d_max, pr.d_max / PR_PLOT_PTS)
123
124        y = np.zeros(len(x))
125        dy = np.zeros(len(x))
126
127        total = 0.0
128        pmax = 0.0
129        cov2 = np.ascontiguousarray(cov)
130
131        for i in range(len(x)):
132            if cov2 is None:
133                value = pr.pr(out, x[i])
134            else:
135                (value, dy[i]) = pr.pr_err(out, cov2, x[i])
136            total += value * pr.d_max / len(x)
137
138            # keep track of the maximum P(r) value
139            if value > pmax:
140                pmax = value
141
142            y[i] = value
143
144        if cov2 is None:
145            new_plot = Data1D(x, y)
146        else:
147            new_plot = Data1D(x, y, dy=dy)
148        new_plot.name = PR_FIT_LABEL
149        new_plot.xaxis("\\rm{r}", 'A')
150        new_plot.yaxis("\\rm{P(r)} ", "cm^{-3}")
151        new_plot.title = "P(r) fit"
152        new_plot.id = PR_FIT_LABEL
153        new_plot.scale = "linear"
154        new_plot.group_id = GROUP_ID_PR_FIT
155
156        return new_plot
157
158    def add_errors(self, sigma=0.05):
159        """
160        Adds errors to data set is they are not available.
161        Uses  $\Delta y = \sigma | y |$.
162        """
163        self._data.dy = sigma * np.fabs(self._data.y)
164
165    def computeDataRange(self):
166        """
167        Wrapper for calculating the data range based on local dataset
168        """
169        return self.computeRangeFromData(self.data)
170
171    def computeRangeFromData(self, data):
172        """
173        Compute the minimum and the maximum range of the data
174        return the npts contains in data
175        """
176        qmin, qmax = None, None
177        if isinstance(data, Data1D):
178            try:
179                qmin = min(data.x)
180                qmax = max(data.x)
181            except (ValueError, TypeError):
182                msg = "Unable to find min/max/length of \n data named %s" % \
183                            self.data.filename
184                raise ValueError(msg)
185
186        else:
187            qmin = 0
188            try:
189                x = max(np.fabs(data.xmin), np.fabs(data.xmax))
190                y = max(np.fabs(data.ymin), np.fabs(data.ymax))
191            except (ValueError, TypeError):
192                msg = "Unable to find min/max of \n data named %s" % \
193                            self.data.filename
194                raise ValueError(msg)
195            qmax = np.sqrt(x * x + y * y)
196        return qmin, qmax
Note: See TracBrowser for help on using the repository browser.