[8d79a87] | 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 |
---|
| 12 | |
---|
| 13 | class CorfuncCalculator(object): |
---|
| 14 | |
---|
| 15 | # Helper class |
---|
| 16 | class _Struct: |
---|
| 17 | def __init__(self, **entries): |
---|
| 18 | self.__dict__.update(entries) |
---|
| 19 | |
---|
| 20 | class _Interpolator(object): |
---|
| 21 | """ |
---|
| 22 | Interpolates between curve f and curve g over the range start:stop and |
---|
| 23 | caches the result of the function when it's called |
---|
| 24 | |
---|
| 25 | :param f: The first curve to interpolate |
---|
| 26 | :param g: The second curve to interpolate |
---|
| 27 | :param start: The value at which to start the interpolation |
---|
| 28 | :param stop: The value at which to stop the interpolation |
---|
| 29 | """ |
---|
| 30 | def __init__(self, f, g, start, stop): |
---|
| 31 | self.f = f |
---|
| 32 | self.g = g |
---|
| 33 | self.start = start |
---|
| 34 | self.stop = stop |
---|
| 35 | self._lastx = [] |
---|
| 36 | self._lasty = [] |
---|
| 37 | |
---|
| 38 | def __call__(self, x): |
---|
| 39 | if self._lastx == [] or x.tolist() != self._lastx.tolist(): |
---|
| 40 | self._lasty = self._smoothed_function(x) |
---|
| 41 | self._lastx = x |
---|
| 42 | return self._lasty |
---|
| 43 | |
---|
| 44 | def _smoothed_function(self,x): |
---|
| 45 | ys = np.zeros(x.shape) |
---|
| 46 | ys[x <= self.start] = self.f(x[x <= self.start]) |
---|
| 47 | ys[x >= self.stop] = self.g(x[x >= self.stop]) |
---|
| 48 | with warnings.catch_warnings(): |
---|
| 49 | # Ignore divide by zero error |
---|
| 50 | warnings.simplefilter('ignore') |
---|
| 51 | h = 1/(1+(x-self.stop)**2/(self.start-x)**2) |
---|
| 52 | mask = np.logical_and(x > self.start, x < self.stop) |
---|
| 53 | ys[mask] = h[mask]*self.g(x[mask])+(1-h[mask])*self.f(x[mask]) |
---|
| 54 | return ys |
---|
| 55 | |
---|
| 56 | |
---|
[bd9f10e] | 57 | def __init__(self, data=None, lowerq=None, upperq=None, scale=1): |
---|
[8d79a87] | 58 | """ |
---|
| 59 | Initialize the class. |
---|
| 60 | |
---|
| 61 | :param data: Data of the type DataLoader.Data1D |
---|
[bd9f10e] | 62 | :param lowerq: The Q value to use as the boundary for |
---|
| 63 | Guinier extrapolation |
---|
| 64 | :param upperq: A tuple of the form (lower, upper). |
---|
| 65 | Values between lower and upper will be used for Porod extrapolation |
---|
[8d79a87] | 66 | :param scale: Scaling factor for I(q) |
---|
| 67 | """ |
---|
[bd9f10e] | 68 | self._data = None |
---|
| 69 | self.set_data(data, scale) |
---|
[8d79a87] | 70 | self.lowerq = lowerq |
---|
| 71 | self.upperq = upperq |
---|
[bd9f10e] | 72 | self.background = 0 |
---|
| 73 | |
---|
| 74 | def set_data(self, data, scale=1): |
---|
| 75 | """ |
---|
| 76 | Prepares the data for analysis |
---|
| 77 | |
---|
| 78 | :return: new_data = data * scale - background |
---|
| 79 | """ |
---|
| 80 | if data is None: |
---|
| 81 | return |
---|
| 82 | # Only process data of the class Data1D |
---|
| 83 | if not issubclass(data.__class__, Data1D): |
---|
| 84 | raise ValueError, "Data must be of the type DataLoader.Data1D" |
---|
| 85 | |
---|
| 86 | # Prepare the data |
---|
| 87 | new_data = Data1D(x=data.x, y=data.y) |
---|
| 88 | new_data *= scale |
---|
| 89 | |
---|
| 90 | # Ensure the errors are set correctly |
---|
| 91 | if new_data.dy is None or len(new_data.x) != len(new_data.dy) or \ |
---|
| 92 | (min(new_data.dy) == 0 and max(new_data.dy) == 0): |
---|
| 93 | new_data.dy = np.ones(len(new_data.x)) |
---|
| 94 | |
---|
| 95 | self._data = new_data |
---|
| 96 | |
---|
| 97 | def compute_background(self, upperq=None): |
---|
| 98 | """ |
---|
| 99 | Compute the background level from the Porod region of the data |
---|
| 100 | """ |
---|
| 101 | if self._data is None: return 0 |
---|
| 102 | elif upperq is None and self.upperq is not None: upperq = self.upperq |
---|
| 103 | elif upperq == 0: return 0 |
---|
| 104 | q = self._data.x |
---|
| 105 | mask = np.logical_and(q > upperq[0], q < upperq[1]) |
---|
| 106 | _, _, bg = self._fit_porod(q[mask], self._data.y[mask]) |
---|
| 107 | |
---|
| 108 | return bg |
---|
[8d79a87] | 109 | |
---|
| 110 | def compute_extrapolation(self): |
---|
[fd5cfaf] | 111 | """ |
---|
| 112 | Extrapolate and interpolate scattering data |
---|
| 113 | |
---|
| 114 | :return: The extrapolated data |
---|
| 115 | """ |
---|
[8d79a87] | 116 | q = self._data.x |
---|
| 117 | iq = self._data.y |
---|
| 118 | |
---|
| 119 | s2 = self._fit_data(q, iq) |
---|
| 120 | qs = np.arange(0, q[-1]*100, (q[1]-q[0])) |
---|
| 121 | iqs = s2(qs) |
---|
| 122 | |
---|
| 123 | extrapolation = Data1D(qs, iqs) |
---|
| 124 | |
---|
| 125 | return extrapolation |
---|
| 126 | |
---|
[bd9f10e] | 127 | def compute_transform(self, extrapolation, background=None): |
---|
[8d79a87] | 128 | """ |
---|
| 129 | Transform an extrapolated scattering curve into a correlation function. |
---|
[fd5cfaf] | 130 | |
---|
| 131 | :param extrapolation: The extrapolated data |
---|
| 132 | :param background: The background value (if not provided, previously |
---|
| 133 | calculated value will be used) |
---|
| 134 | :return: The transformed data |
---|
[8d79a87] | 135 | """ |
---|
| 136 | qs = extrapolation.x |
---|
| 137 | iqs = extrapolation.y |
---|
[bd9f10e] | 138 | q = self._data.x |
---|
| 139 | if background is None: background = self.background |
---|
[8d79a87] | 140 | |
---|
[bd9f10e] | 141 | gamma = dct((iqs-background)*qs**2) |
---|
[8d79a87] | 142 | gamma = gamma / gamma.max() |
---|
| 143 | xs = np.pi*np.arange(len(qs),dtype=np.float32)/(q[1]-q[0])/len(qs) |
---|
| 144 | |
---|
| 145 | transform = Data1D(xs, gamma) |
---|
| 146 | |
---|
| 147 | return transform |
---|
| 148 | |
---|
[bd9f10e] | 149 | def _porod(self, q, K, sigma, bg): |
---|
[8d79a87] | 150 | """Equation for the Porod region of the data""" |
---|
[bd9f10e] | 151 | return bg + (K*q**(-4))*np.exp(-q**2*sigma**2) |
---|
[8d79a87] | 152 | |
---|
| 153 | def _fit_guinier(self, q, iq): |
---|
| 154 | """Fit the Guinier region of the curve""" |
---|
| 155 | A = np.vstack([q**2, np.ones(q.shape)]).T |
---|
| 156 | return lstsq(A, np.log(iq)) |
---|
| 157 | |
---|
[bd9f10e] | 158 | def _fit_porod(self, q, iq): |
---|
| 159 | """Fit the Porod region of the curve""" |
---|
| 160 | fitp = curve_fit(lambda q, k, sig, bg: self._porod(q, k, sig, bg)*q**2, |
---|
| 161 | q, iq*q**2)[0] |
---|
| 162 | k, sigma, bg = fitp |
---|
| 163 | return k, sigma, bg |
---|
[8d79a87] | 164 | |
---|
| 165 | def _fit_data(self, q, iq): |
---|
[fd5cfaf] | 166 | """ |
---|
| 167 | Given a data set, extrapolate out to large q with Porod and |
---|
| 168 | to q=0 with Guinier |
---|
| 169 | """ |
---|
[8d79a87] | 170 | mask = np.logical_and(q > self.upperq[0], q < self.upperq[1]) |
---|
| 171 | |
---|
| 172 | # Returns an array where the 1st and 2nd elements are the values of k |
---|
| 173 | # and sigma for the best-fit Porod function |
---|
[01a26f2] | 174 | k, sigma, _ = self._fit_porod(q[mask], iq[mask]) |
---|
| 175 | bg = self.background |
---|
[8d79a87] | 176 | |
---|
| 177 | # Smooths between the best-fit porod function and the data to produce a |
---|
| 178 | # better fitting curve |
---|
| 179 | data = interp1d(q, iq) |
---|
[a9b7a8f] | 180 | s1 = self._Interpolator(data, |
---|
[bd9f10e] | 181 | lambda x: self._porod(x, k, sigma, bg), self.upperq[0], q[-1]) |
---|
[8d79a87] | 182 | |
---|
| 183 | mask = np.logical_and(q < self.lowerq, 0 < q) |
---|
| 184 | |
---|
| 185 | # Returns parameters for the best-fit Guinier function |
---|
| 186 | g = self._fit_guinier(q[mask], iq[mask])[0] |
---|
| 187 | |
---|
| 188 | # Smooths between the best-fit Guinier function and the Porod curve |
---|
| 189 | s2 = self._Interpolator((lambda x: (np.exp(g[1]+g[0]*x**2))), s1, q[0], |
---|
| 190 | self.lowerq) |
---|
| 191 | |
---|
| 192 | return s2 |
---|