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 | |
---|
57 | def __init__(self, data, lowerq, upperq, background=0, scale=1): |
---|
58 | """ |
---|
59 | Initialize the class. |
---|
60 | |
---|
61 | :param data: Data of the type DataLoader.Data1D |
---|
62 | :param background: Background value. Will be subtracted from the data |
---|
63 | before processing |
---|
64 | :param scale: Scaling factor for I(q) |
---|
65 | """ |
---|
66 | print "Before: {}".format(data.y[0]) |
---|
67 | self._data = self._get_data(data, background, scale) |
---|
68 | print "After: {}".format(self._data.y[0]) |
---|
69 | self.lowerq = lowerq |
---|
70 | self.upperq = upperq |
---|
71 | |
---|
72 | def compute_extrapolation(self): |
---|
73 | q = self._data.x |
---|
74 | iq = self._data.y |
---|
75 | |
---|
76 | s2 = self._fit_data(q, iq) |
---|
77 | qs = np.arange(0, q[-1]*100, (q[1]-q[0])) |
---|
78 | iqs = s2(qs) |
---|
79 | |
---|
80 | extrapolation = Data1D(qs, iqs) |
---|
81 | |
---|
82 | return extrapolation |
---|
83 | |
---|
84 | def compute_transform(self, extrapolation): |
---|
85 | """ |
---|
86 | Transform an extrapolated scattering curve into a correlation function. |
---|
87 | """ |
---|
88 | qs = extrapolation.x |
---|
89 | iqs = extrapolation.y |
---|
90 | |
---|
91 | gamma = dct(iqs*qs**2) |
---|
92 | gamma = gamma / gamma.max() |
---|
93 | xs = np.pi*np.arange(len(qs),dtype=np.float32)/(q[1]-q[0])/len(qs) |
---|
94 | |
---|
95 | transform = Data1D(xs, gamma) |
---|
96 | |
---|
97 | return transform |
---|
98 | |
---|
99 | def _porod(self, q, K, sigma): |
---|
100 | """Equation for the Porod region of the data""" |
---|
101 | return (K*q**(-4))*np.exp(-q**2*sigma**2) |
---|
102 | |
---|
103 | def _fit_guinier(self, q, iq): |
---|
104 | """Fit the Guinier region of the curve""" |
---|
105 | A = np.vstack([q**2, np.ones(q.shape)]).T |
---|
106 | return lstsq(A, np.log(iq)) |
---|
107 | |
---|
108 | def _get_data(self, data, background, scale): |
---|
109 | """ |
---|
110 | Prepares the data for analysis |
---|
111 | |
---|
112 | :return: new_data = data * scale - background |
---|
113 | """ |
---|
114 | # Only process data of the class Data1D |
---|
115 | if not issubclass(data.__class__, Data1D): |
---|
116 | raise ValueError, "Data must be of the type DataLoader.Data1D" |
---|
117 | |
---|
118 | # Prepare the data |
---|
119 | new_data = (scale * data) |
---|
120 | new_data.y -= background |
---|
121 | |
---|
122 | # Check the vector lengths are equal |
---|
123 | assert len(new_data.x) == len(new_data.y) |
---|
124 | |
---|
125 | # Ensure the errors are set correctly |
---|
126 | if new_data.dy is None or len(new_data.x) != len(new_data.dy) or \ |
---|
127 | (min(new_data.dy) == 0 and max(new_data.dy) == 0): |
---|
128 | new_data.dy = np.ones(len(new_data.x)) |
---|
129 | |
---|
130 | return new_data |
---|
131 | |
---|
132 | def _fit_data(self, q, iq): |
---|
133 | """Given a data set, extrapolate out to large q with Porod |
---|
134 | and to q=0 with Guinier""" |
---|
135 | mask = np.logical_and(q > self.upperq[0], q < self.upperq[1]) |
---|
136 | |
---|
137 | # Returns an array where the 1st and 2nd elements are the values of k |
---|
138 | # and sigma for the best-fit Porod function |
---|
139 | fitp = curve_fit(lambda q, k, sig: self._porod(q, k, sig)*q**2, |
---|
140 | q[mask], iq[mask]*q[mask]**2)[0] |
---|
141 | |
---|
142 | # Smooths between the best-fit porod function and the data to produce a |
---|
143 | # better fitting curve |
---|
144 | data = interp1d(q, iq) |
---|
145 | s1 = self._Interpolator(data, lambda x: self._porod(x, fitp[0], fitp[1]), |
---|
146 | self.upperq[0], q[-1]) |
---|
147 | |
---|
148 | mask = np.logical_and(q < self.lowerq, 0 < q) |
---|
149 | |
---|
150 | # Returns parameters for the best-fit Guinier function |
---|
151 | g = self._fit_guinier(q[mask], iq[mask])[0] |
---|
152 | |
---|
153 | # Smooths between the best-fit Guinier function and the Porod curve |
---|
154 | s2 = self._Interpolator((lambda x: (np.exp(g[1]+g[0]*x**2))), s1, q[0], |
---|
155 | self.lowerq) |
---|
156 | |
---|
157 | return s2 |
---|