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

Last change on this file since 5a2bb75 was 6da860a, checked in by krzywon, 7 years ago

Responding to P(r) code review.

  • Property mode set to 100644
File size: 5.4 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        return new_plot
114
115    def newPRPlot(self, out, pr, cov=None):
116        """
117        """
118        # Show P(r)
119        x = np.arange(0.0, pr.d_max, pr.d_max / PR_PLOT_PTS)
120
121        y = np.zeros(len(x))
122        dy = np.zeros(len(x))
123
124        total = 0.0
125        pmax = 0.0
126        cov2 = np.ascontiguousarray(cov)
127
128        for i in range(len(x)):
129            if cov2 is None:
130                value = pr.pr(out, x[i])
131            else:
132                (value, dy[i]) = pr.pr_err(out, cov2, x[i])
133            total += value * pr.d_max / len(x)
134
135            # keep track of the maximum P(r) value
136            if value > pmax:
137                pmax = value
138
139            y[i] = value
140
141        if cov2 is None:
142            new_plot = Data1D(x, y)
143        else:
144            new_plot = Data1D(x, y, dy=dy)
145        new_plot.name = PR_FIT_LABEL
146        new_plot.xaxis("\\rm{r}", 'A')
147        new_plot.yaxis("\\rm{P(r)} ", "cm^{-3}")
148        new_plot.title = "P(r) fit"
149        new_plot.id = PR_FIT_LABEL
150        new_plot.scale = "linear"
151        new_plot.group_id = GROUP_ID_PR_FIT
152
153        return new_plot
154
155    def computeDataRange(self):
156        """
157        Wrapper for calculating the data range based on local dataset
158        """
159        return self.computeRangeFromData(self.data)
160
161    def computeRangeFromData(self, data):
162        """
163        Compute the minimum and the maximum range of the data
164        return the npts contains in data
165        """
166        qmin, qmax = None, None
167        if isinstance(data, Data1D):
168            try:
169                qmin = min(data.x)
170                qmax = max(data.x)
171            except (ValueError, TypeError):
172                msg = "Unable to find min/max/length of \n data named %s" % \
173                            self.data.filename
174                raise ValueError(msg)
175
176        else:
177            qmin = 0
178            try:
179                x = max(np.fabs(data.xmin), np.fabs(data.xmax))
180                y = max(np.fabs(data.ymin), np.fabs(data.ymax))
181            except (ValueError, TypeError):
182                msg = "Unable to find min/max of \n data named %s" % \
183                            self.data.filename
184                raise ValueError(msg)
185            qmax = np.sqrt(x * x + y * y)
186        return qmin, qmax
Note: See TracBrowser for help on using the repository browser.