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

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since ff11b21 was ff11b21, checked in by lewis, 8 years ago

Show extrapolation parameters in GUI

Also change I_obs label to I_obs - Bg and correct Guinier equation in the
docs.

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