[3b8efec] | 1 | """ |
---|
| 2 | This module implements corfunc |
---|
| 3 | """ |
---|
| 4 | import warnings |
---|
| 5 | import numpy as np |
---|
| 6 | from scipy.optimize import curve_fit |
---|
| 7 | from scipy.interpolate import interp1d |
---|
| 8 | from scipy.fftpack import dct |
---|
| 9 | from scipy.signal import argrelextrema |
---|
| 10 | from numpy.linalg import lstsq |
---|
| 11 | from sas.sascalc.dataloader.data_info import Data1D |
---|
[d03228e] | 12 | from sas.sascalc.corfunc.transform_thread import FourierThread |
---|
| 13 | from sas.sascalc.corfunc.transform_thread import HilbertThread |
---|
[3b8efec] | 14 | |
---|
| 15 | class 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): |
---|
[412c509] | 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(): |
---|
[3b8efec] | 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 | |
---|
[911dbe4] | 62 | def __init__(self, data=None, lowerq=None, upperq=None, scale=1): |
---|
[3b8efec] | 63 | """ |
---|
| 64 | Initialize the class. |
---|
| 65 | |
---|
| 66 | :param data: Data of the type DataLoader.Data1D |
---|
[911dbe4] | 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 |
---|
[3b8efec] | 71 | :param scale: Scaling factor for I(q) |
---|
| 72 | """ |
---|
[911dbe4] | 73 | self._data = None |
---|
| 74 | self.set_data(data, scale) |
---|
[3b8efec] | 75 | self.lowerq = lowerq |
---|
| 76 | self.upperq = upperq |
---|
[5878a9ea] | 77 | self.background = self.compute_background() |
---|
[a2db1ab] | 78 | self._transform_thread = None |
---|
[911dbe4] | 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 |
---|
[5878a9ea] | 109 | elif upperq is None and self.upperq is None: return 0 |
---|
[911dbe4] | 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 |
---|
[3b8efec] | 115 | |
---|
| 116 | def compute_extrapolation(self): |
---|
[6970e51] | 117 | """ |
---|
| 118 | Extrapolate and interpolate scattering data |
---|
| 119 | |
---|
| 120 | :return: The extrapolated data |
---|
| 121 | """ |
---|
[3b8efec] | 122 | q = self._data.x |
---|
| 123 | iq = self._data.y |
---|
| 124 | |
---|
[ff11b21] | 125 | params, s2 = self._fit_data(q, iq) |
---|
[3b8efec] | 126 | qs = np.arange(0, q[-1]*100, (q[1]-q[0])) |
---|
| 127 | iqs = s2(qs) |
---|
| 128 | |
---|
| 129 | extrapolation = Data1D(qs, iqs) |
---|
| 130 | |
---|
[412c509] | 131 | return params, extrapolation, s2 |
---|
[3b8efec] | 132 | |
---|
[412c509] | 133 | def compute_transform(self, extrapolation, trans_type, extrap_fn=None, |
---|
| 134 | background=None, completefn=None, updatefn=None): |
---|
[3b8efec] | 135 | """ |
---|
| 136 | Transform an extrapolated scattering curve into a correlation function. |
---|
[6970e51] | 137 | |
---|
| 138 | :param extrapolation: The extrapolated data |
---|
| 139 | :param background: The background value (if not provided, previously |
---|
| 140 | calculated value will be used) |
---|
[412c509] | 141 | :param extrap_fn: A callable function representing the extraoplated data |
---|
[a2db1ab] | 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 |
---|
[6970e51] | 146 | :return: The transformed data |
---|
[3b8efec] | 147 | """ |
---|
[a2db1ab] | 148 | if self._transform_thread is not None: |
---|
| 149 | if self._transform_thread.isrunning(): return |
---|
| 150 | |
---|
[911dbe4] | 151 | if background is None: background = self.background |
---|
[3b8efec] | 152 | |
---|
[d03228e] | 153 | if trans_type == 'fourier': |
---|
| 154 | self._transform_thread = FourierThread(self._data, extrapolation, |
---|
[412c509] | 155 | background, extrap_fn=extrap_fn, completefn=completefn, |
---|
| 156 | updatefn=updatefn) |
---|
[d03228e] | 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 | |
---|
[a2db1ab] | 165 | self._transform_thread.queue() |
---|
[3b8efec] | 166 | |
---|
[a2db1ab] | 167 | def transform_isrunning(self): |
---|
| 168 | if self._transform_thread is None: return False |
---|
| 169 | return self._transform_thread.isrunning() |
---|
[3b8efec] | 170 | |
---|
[a2db1ab] | 171 | def stop_transform(self): |
---|
| 172 | if self._transform_thread.isrunning(): |
---|
| 173 | self._transform_thread.stop() |
---|
[3b8efec] | 174 | |
---|
[033c14c] | 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)/m # 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 |
---|
[eb320682] | 222 | return { 'max': self._round_sig_figs(x[maxs[0]], 6) } |
---|
[033c14c] | 223 | dtr = x[mask[0]] # Beginning of Linear Section |
---|
| 224 | d0 = x[mask[-1]] # End of Linear Section |
---|
| 225 | GammaMax = y[mask[-1]] |
---|
[8bdc103] | 226 | A = np.abs(GammaMin/GammaMax) # Normalized depth of minimum |
---|
[033c14c] | 227 | |
---|
| 228 | params = { |
---|
| 229 | 'max': x[maxs[0]], |
---|
| 230 | 'dtr': dtr, |
---|
| 231 | 'Lc': Lc, |
---|
| 232 | 'd0': d0, |
---|
| 233 | 'A': A, |
---|
[a684c64] | 234 | 'fill': Lc/x[maxs[0]] |
---|
[033c14c] | 235 | } |
---|
| 236 | |
---|
| 237 | return params |
---|
| 238 | |
---|
| 239 | |
---|
[911dbe4] | 240 | def _porod(self, q, K, sigma, bg): |
---|
[3b8efec] | 241 | """Equation for the Porod region of the data""" |
---|
[911dbe4] | 242 | return bg + (K*q**(-4))*np.exp(-q**2*sigma**2) |
---|
[3b8efec] | 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 | |
---|
[911dbe4] | 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, |
---|
[c97d76a] | 252 | q, iq*q**2, bounds=([-np.inf, 0, -np.inf], [np.inf, np.inf, np.inf]))[0] |
---|
[911dbe4] | 253 | k, sigma, bg = fitp |
---|
| 254 | return k, sigma, bg |
---|
[3b8efec] | 255 | |
---|
| 256 | def _fit_data(self, q, iq): |
---|
[6970e51] | 257 | """ |
---|
| 258 | Given a data set, extrapolate out to large q with Porod and |
---|
| 259 | to q=0 with Guinier |
---|
| 260 | """ |
---|
[3b8efec] | 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 |
---|
[9454a27] | 265 | k, sigma, _ = self._fit_porod(q[mask], iq[mask]) |
---|
| 266 | bg = self.background |
---|
[3b8efec] | 267 | |
---|
| 268 | # Smooths between the best-fit porod function and the data to produce a |
---|
| 269 | # better fitting curve |
---|
| 270 | data = interp1d(q, iq) |
---|
[77e9ac6] | 271 | s1 = self._Interpolator(data, |
---|
[911dbe4] | 272 | lambda x: self._porod(x, k, sigma, bg), self.upperq[0], q[-1]) |
---|
[3b8efec] | 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 | |
---|
[ff11b21] | 283 | params = {'A': g[1], 'B': g[0], 'K': k, 'sigma': sigma} |
---|
| 284 | |
---|
| 285 | return params, s2 |
---|