source: sasview/src/sas/sascalc/corfunc/corfunc_calculator.py @ c728295

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalcmagnetic_scattrelease-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since c728295 was c728295, checked in by lewis, 7 years ago

Remove unused argument and import

  • Property mode set to 100644
File size: 10.3 KB
Line 
1"""
2This module implements corfunc
3"""
4import warnings
5import numpy as np
6from scipy.optimize import curve_fit
7from scipy.interpolate import interp1d
8from scipy.fftpack import dct
9from scipy.signal import argrelextrema
10from numpy.linalg import lstsq
11from sas.sascalc.dataloader.data_info import Data1D
12from sas.sascalc.corfunc.transform_thread import FourierThread
13from sas.sascalc.corfunc.transform_thread import HilbertThread
14
15class CorfuncCalculator(object):
16
17    class _Interpolator(object):
18        """
19        Interpolates between curve f and curve g over the range start:stop and
20        caches the result of the function when it's called
21
22        :param f: The first curve to interpolate
23        :param g: The second curve to interpolate
24        :param start: The value at which to start the interpolation
25        :param stop: The value at which to stop the interpolation
26        """
27        def __init__(self, f, g, start, stop):
28            self.f = f
29            self.g = g
30            self.start = start
31            self.stop = stop
32            self._lastx = []
33            self._lasty = []
34
35        def __call__(self, x):
36            # If input is a single number, evaluate the function at that number
37            # and return a single number
38            if type(x) == float or type(x) == int:
39                return self._smoothed_function(np.array([x]))[0]
40            # If input is a list, and is different to the last input, evaluate
41            # the function at each point. If the input is the same as last time
42            # the function was called, return the result that was calculated
43            # last time instead of explicity evaluating the function again.
44            elif self._lastx == [] or x.tolist() != self._lastx.tolist():
45                self._lasty = self._smoothed_function(x)
46                self._lastx = x
47            return self._lasty
48
49        def _smoothed_function(self,x):
50            ys = np.zeros(x.shape)
51            ys[x <= self.start] = self.f(x[x <= self.start])
52            ys[x >= self.stop] = self.g(x[x >= self.stop])
53            with warnings.catch_warnings():
54                # Ignore divide by zero error
55                warnings.simplefilter('ignore')
56                h = 1/(1+(x-self.stop)**2/(self.start-x)**2)
57            mask = np.logical_and(x > self.start, x < self.stop)
58            ys[mask] = h[mask]*self.g(x[mask])+(1-h[mask])*self.f(x[mask])
59            return ys
60
61
62    def __init__(self, data=None, lowerq=None, upperq=None, scale=1):
63        """
64        Initialize the class.
65
66        :param data: Data of the type DataLoader.Data1D
67        :param lowerq: The Q value to use as the boundary for
68            Guinier extrapolation
69        :param upperq: A tuple of the form (lower, upper).
70            Values between lower and upper will be used for Porod extrapolation
71        :param scale: Scaling factor for I(q)
72        """
73        self._data = None
74        self.set_data(data, scale)
75        self.lowerq = lowerq
76        self.upperq = upperq
77        self.background = self.compute_background()
78        self._transform_thread = None
79
80    def set_data(self, data, scale=1):
81        """
82        Prepares the data for analysis
83
84        :return: new_data = data * scale - background
85        """
86        if data is None:
87            return
88        # Only process data of the class Data1D
89        if not issubclass(data.__class__, Data1D):
90            raise ValueError, "Data must be of the type DataLoader.Data1D"
91
92        # Prepare the data
93        new_data = Data1D(x=data.x, y=data.y)
94        new_data *= scale
95
96        # Ensure the errors are set correctly
97        if new_data.dy is None or len(new_data.x) != len(new_data.dy) or \
98            (min(new_data.dy) == 0 and max(new_data.dy) == 0):
99            new_data.dy = np.ones(len(new_data.x))
100
101        self._data = new_data
102
103    def compute_background(self, upperq=None):
104        """
105        Compute the background level from the Porod region of the data
106        """
107        if self._data is None: return 0
108        elif upperq is None and self.upperq is not None: upperq = self.upperq
109        elif upperq is None and self.upperq is None: return 0
110        q = self._data.x
111        mask = np.logical_and(q > upperq[0], q < upperq[1])
112        _, _, bg = self._fit_porod(q[mask], self._data.y[mask])
113
114        return bg
115
116    def compute_extrapolation(self):
117        """
118        Extrapolate and interpolate scattering data
119
120        :return: The extrapolated data
121        """
122        q = self._data.x
123        iq = self._data.y
124
125        params, s2 = self._fit_data(q, iq)
126        qs = np.arange(0, q[-1]*100, (q[1]-q[0]))
127        iqs = s2(qs)
128
129        extrapolation = Data1D(qs, iqs)
130
131        return params, extrapolation, s2
132
133    def compute_transform(self, extrapolation, trans_type, background=None,
134        completefn=None, updatefn=None):
135        """
136        Transform an extrapolated scattering curve into a correlation function.
137
138        :param extrapolation: The extrapolated data
139        :param background: The background value (if not provided, previously
140            calculated value will be used)
141        :param extrap_fn: A callable function representing the extraoplated data
142        :param completefn: The function to call when the transform calculation
143            is complete`
144        :param updatefn: The function to call to update the GUI with the status
145            of the transform calculation
146        :return: The transformed data
147        """
148        if self._transform_thread is not None:
149            if self._transform_thread.isrunning(): return
150
151        if background is None: background = self.background
152
153        if trans_type == 'fourier':
154            self._transform_thread = FourierThread(self._data, extrapolation,
155            background, completefn=completefn,
156            updatefn=updatefn)
157        elif trans_type == 'hilbert':
158            self._transform_thread = HilbertThread(self._data, extrapolation,
159            background, completefn=completefn, updatefn=updatefn)
160        else:
161            err = ("Incorrect transform type supplied, must be 'fourier'",
162                " or 'hilbert'")
163            raise ValueError, err
164
165        self._transform_thread.queue()
166
167    def transform_isrunning(self):
168        if self._transform_thread is None: return False
169        return self._transform_thread.isrunning()
170
171    def stop_transform(self):
172        if self._transform_thread.isrunning():
173            self._transform_thread.stop()
174
175    def extract_parameters(self, transformed_data):
176        """
177        Extract the interesting measurements from a correlation function
178        :param transformed_data: Fourier transformation of the
179            extrapolated data
180        """
181        # Calculate indexes of maxima and minima
182        x = transformed_data.x
183        y = transformed_data.y
184        maxs = argrelextrema(y, np.greater)[0]
185        mins = argrelextrema(y, np.less)[0]
186
187        # If there are no maxima, return None
188        if len(maxs) == 0:
189            return None
190
191        GammaMin = y[mins[0]]  # The value at the first minimum
192
193        ddy = (y[:-2]+y[2:]-2*y[1:-1])/(x[2:]-x[:-2])**2  # 2nd derivative of y
194        dy = (y[2:]-y[:-2])/(x[2:]-x[:-2])  # 1st derivative of y
195        # Find where the second derivative goes to zero
196        zeros = argrelextrema(np.abs(ddy), np.less)[0]
197        # locate the first inflection point
198        linear_point = zeros[0]
199
200        # Try to calculate slope around linear_point using 80 data points
201        lower = linear_point - 40
202        upper = linear_point + 40
203
204        # If too few data points to the left, use linear_point*2 data points
205        if lower < 0:
206            lower = 0
207            upper = linear_point * 2
208        # If too few to right, use 2*(dy.size - linear_point) data points
209        elif upper > len(dy):
210            upper = len(dy)
211            width = len(dy) - linear_point
212            lower = 2*linear_point - dy.size
213
214        m = np.mean(dy[lower:upper])  # Linear slope
215        b = y[1:-1][linear_point]-m*x[1:-1][linear_point]  # Linear intercept
216
217        Lc = (GammaMin-b)/# Hard block thickness
218
219        # Find the data points where the graph is linear to within 1%
220        mask = np.where(np.abs((y-(m*x+b))/y) < 0.01)[0]
221        if len(mask) == 0:  # Return garbage for bad fits
222            return { 'max': self._round_sig_figs(x[maxs[0]], 6) }
223        dtr = x[mask[0]]  # Beginning of Linear Section
224        d0 = x[mask[-1]]  # End of Linear Section
225        GammaMax = y[mask[-1]]
226        A = np.abs(GammaMin/GammaMax)  # Normalized depth of minimum
227
228        params = {
229            'max': x[maxs[0]],
230            'dtr': dtr,
231            'Lc': Lc,
232            'd0': d0,
233            'A': A,
234            'fill': Lc/x[maxs[0]]
235        }
236
237        return params
238
239
240    def _porod(self, q, K, sigma, bg):
241        """Equation for the Porod region of the data"""
242        return bg + (K*q**(-4))*np.exp(-q**2*sigma**2)
243
244    def _fit_guinier(self, q, iq):
245        """Fit the Guinier region of the curve"""
246        A = np.vstack([q**2, np.ones(q.shape)]).T
247        return lstsq(A, np.log(iq))
248
249    def _fit_porod(self, q, iq):
250        """Fit the Porod region of the curve"""
251        fitp = curve_fit(lambda q, k, sig, bg: self._porod(q, k, sig, bg)*q**2,
252                         q, iq*q**2, bounds=([-np.inf, 0, -np.inf], [np.inf, np.inf, np.inf]))[0]
253        k, sigma, bg = fitp
254        return k, sigma, bg
255
256    def _fit_data(self, q, iq):
257        """
258        Given a data set, extrapolate out to large q with Porod and
259        to q=0 with Guinier
260        """
261        mask = np.logical_and(q > self.upperq[0], q < self.upperq[1])
262
263        # Returns an array where the 1st and 2nd elements are the values of k
264        # and sigma for the best-fit Porod function
265        k, sigma, _ = self._fit_porod(q[mask], iq[mask])
266        bg = self.background
267
268        # Smooths between the best-fit porod function and the data to produce a
269        # better fitting curve
270        data = interp1d(q, iq)
271        s1 = self._Interpolator(data,
272            lambda x: self._porod(x, k, sigma, bg), self.upperq[0], q[-1])
273
274        mask = np.logical_and(q < self.lowerq, 0 < q)
275
276        # Returns parameters for the best-fit Guinier function
277        g = self._fit_guinier(q[mask], iq[mask])[0]
278
279        # Smooths between the best-fit Guinier function and the Porod curve
280        s2 = self._Interpolator((lambda x: (np.exp(g[1]+g[0]*x**2))), s1, q[0],
281            self.lowerq)
282
283        params = {'A': g[1], 'B': g[0], 'K': k, 'sigma': sigma}
284
285        return params, s2
Note: See TracBrowser for help on using the repository browser.