source: sasmodels/sasmodels/data.py @ d18582e

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since d18582e was d18582e, checked in by Paul Kienzle <pkienzle@…>, 7 years ago

default to double precision if single=False is set in model file

  • Property mode set to 100644
File size: 17.1 KB
Line 
1"""
2SAS data representations.
3
4Plotting functions for data sets:
5
6    :func:`plot_data` plots the data file.
7
8    :func:`plot_theory` plots a calculated result from the model.
9
10Wrappers for the sasview data loader and data manipulations:
11
12    :func:`load_data` loads a sasview data file.
13
14    :func:`set_beam_stop` masks the beam stop from the data.
15
16    :func:`set_half` selects the right or left half of the data, which can
17    be useful for shear measurements which have not been properly corrected
18    for path length and reflections.
19
20    :func:`set_top` cuts the top part off the data.
21
22
23Empty data sets for evaluating models without data:
24
25    :func:`empty_data1D` creates an empty dataset, which is useful for plotting
26    a theory function before the data is measured.
27
28    :func:`empty_data2D` creates an empty 2D dataset.
29
30Note that the empty datasets use a minimal representation of the SasView
31objects so that models can be run without SasView on the path.  You could
32also use these for your own data loader.
33
34"""
35import traceback
36
37import numpy as np
38
39def load_data(filename):
40    """
41    Load data using a sasview loader.
42    """
43    from sas.dataloader.loader import Loader
44    loader = Loader()
45    data = loader.load(filename)
46    if data is None:
47        raise IOError("Data %r could not be loaded" % filename)
48    return data
49
50
51def set_beam_stop(data, radius, outer=None):
52    """
53    Add a beam stop of the given *radius*.  If *outer*, make an annulus.
54    """
55    from sas.dataloader.manipulations import Ringcut
56    if hasattr(data, 'qx_data'):
57        data.mask = Ringcut(0, radius)(data)
58        if outer is not None:
59            data.mask += Ringcut(outer, np.inf)(data)
60    else:
61        data.mask = (data.x < radius)
62        if outer is not None:
63            data.mask |= (data.x >= outer)
64
65
66def set_half(data, half):
67    """
68    Select half of the data, either "right" or "left".
69    """
70    from sas.dataloader.manipulations import Boxcut
71    if half == 'right':
72        data.mask += \
73            Boxcut(x_min=-np.inf, x_max=0.0, y_min=-np.inf, y_max=np.inf)(data)
74    if half == 'left':
75        data.mask += \
76            Boxcut(x_min=0.0, x_max=np.inf, y_min=-np.inf, y_max=np.inf)(data)
77
78
79def set_top(data, cutoff):
80    """
81    Chop the top off the data, above *cutoff*.
82    """
83    from sas.dataloader.manipulations import Boxcut
84    data.mask += \
85        Boxcut(x_min=-np.inf, x_max=np.inf, y_min=-np.inf, y_max=cutoff)(data)
86
87
88class Data1D(object):
89    """
90    1D data object.
91
92    Note that this definition matches the attributes from sasview, with
93    some generic 1D data vectors and some SAS specific definitions.  Some
94    refactoring to allow consistent naming conventions between 1D, 2D and
95    SESANS data would be helpful.
96
97    **Attributes**
98
99    *x*, *dx*: $q$ vector and gaussian resolution
100
101    *y*, *dy*: $I(q)$ vector and measurement uncertainty
102
103    *mask*: values to include in plotting/analysis
104
105    *dxl*: slit widths for slit smeared data, with *dx* ignored
106
107    *qmin*, *qmax*: range of $q$ values in *x*
108
109    *filename*: label for the data line
110
111    *_xaxis*, *_xunit*: label and units for the *x* axis
112
113    *_yaxis*, *_yunit*: label and units for the *y* axis
114    """
115    def __init__(self, x=None, y=None, dx=None, dy=None):
116        self.x, self.y, self.dx, self.dy = x, y, dx, dy
117        self.dxl = None
118        self.filename = None
119        self.qmin = x.min() if x is not None else np.NaN
120        self.qmax = x.max() if x is not None else np.NaN
121        # TODO: why is 1D mask False and 2D mask True?
122        self.mask = (np.isnan(y) if y is not None
123                     else np.zeros_like(x, 'b') if x is not None
124                     else None)
125        self._xaxis, self._xunit = "x", ""
126        self._yaxis, self._yunit = "y", ""
127
128    def xaxis(self, label, unit):
129        """
130        set the x axis label and unit
131        """
132        self._xaxis = label
133        self._xunit = unit
134
135    def yaxis(self, label, unit):
136        """
137        set the y axis label and unit
138        """
139        self._yaxis = label
140        self._yunit = unit
141
142
143
144class Data2D(object):
145    """
146    2D data object.
147
148    Note that this definition matches the attributes from sasview. Some
149    refactoring to allow consistent naming conventions between 1D, 2D and
150    SESANS data would be helpful.
151
152    **Attributes**
153
154    *qx_data*, *dqx_data*: $q_x$ matrix and gaussian resolution
155
156    *qy_data*, *dqy_data*: $q_y$ matrix and gaussian resolution
157
158    *data*, *err_data*: $I(q)$ matrix and measurement uncertainty
159
160    *mask*: values to exclude from plotting/analysis
161
162    *qmin*, *qmax*: range of $q$ values in *x*
163
164    *filename*: label for the data line
165
166    *_xaxis*, *_xunit*: label and units for the *x* axis
167
168    *_yaxis*, *_yunit*: label and units for the *y* axis
169
170    *_zaxis*, *_zunit*: label and units for the *y* axis
171
172    *Q_unit*, *I_unit*: units for Q and intensity
173
174    *x_bins*, *y_bins*: grid steps in *x* and *y* directions
175    """
176    def __init__(self, x=None, y=None, z=None, dx=None, dy=None, dz=None):
177        self.qx_data, self.dqx_data = x, dx
178        self.qy_data, self.dqy_data = y, dy
179        self.data, self.err_data = z, dz
180        self.mask = (~np.isnan(z) if z is not None
181                     else np.ones_like(x) if x is not None
182                     else None)
183        self.q_data = np.sqrt(x**2 + y**2)
184        self.qmin = 1e-16
185        self.qmax = np.inf
186        self.detector = []
187        self.source = Source()
188        self.Q_unit = "1/A"
189        self.I_unit = "1/cm"
190        self.xaxis("Q_x", "1/A")
191        self.yaxis("Q_y", "1/A")
192        self.zaxis("Intensity", "1/cm")
193        self._xaxis, self._xunit = "x", ""
194        self._yaxis, self._yunit = "y", ""
195        self._zaxis, self._zunit = "z", ""
196        self.x_bins, self.y_bins = None, None
197
198    def xaxis(self, label, unit):
199        """
200        set the x axis label and unit
201        """
202        self._xaxis = label
203        self._xunit = unit
204
205    def yaxis(self, label, unit):
206        """
207        set the y axis label and unit
208        """
209        self._yaxis = label
210        self._yunit = unit
211
212    def zaxis(self, label, unit):
213        """
214        set the y axis label and unit
215        """
216        self._zaxis = label
217        self._zunit = unit
218
219
220class Vector(object):
221    """
222    3-space vector of *x*, *y*, *z*
223    """
224    def __init__(self, x=None, y=None, z=None):
225        self.x, self.y, self.z = x, y, z
226
227class Detector(object):
228    """
229    Detector attributes.
230    """
231    def __init__(self, pixel_size=(None, None), distance=None):
232        self.pixel_size = Vector(*pixel_size)
233        self.distance = distance
234
235class Source(object):
236    """
237    Beam attributes.
238    """
239    def __init__(self):
240        self.wavelength = np.NaN
241        self.wavelength_unit = "A"
242
243
244def empty_data1D(q, resolution=0.0):
245    """
246    Create empty 1D data using the given *q* as the x value.
247
248    *resolution* dq/q defaults to 5%.
249    """
250
251    #Iq = 100 * np.ones_like(q)
252    #dIq = np.sqrt(Iq)
253    Iq, dIq = None, None
254    q = np.asarray(q)
255    data = Data1D(q, Iq, dx=resolution * q, dy=dIq)
256    data.filename = "fake data"
257    return data
258
259
260def empty_data2D(qx, qy=None, resolution=0.0):
261    """
262    Create empty 2D data using the given mesh.
263
264    If *qy* is missing, create a square mesh with *qy=qx*.
265
266    *resolution* dq/q defaults to 5%.
267    """
268    if qy is None:
269        qy = qx
270    qx, qy = np.asarray(qx), np.asarray(qy)
271    # 5% dQ/Q resolution
272    Qx, Qy = np.meshgrid(qx, qy)
273    Qx, Qy = Qx.flatten(), Qy.flatten()
274    Iq = 100 * np.ones_like(Qx)
275    dIq = np.sqrt(Iq)
276    if resolution != 0:
277        # https://www.ncnr.nist.gov/staff/hammouda/distance_learning/chapter_15.pdf
278        # Should have an additional constant which depends on distances and
279        # radii of the aperture, pixel dimensions and wavelength spread
280        # Instead, assume radial dQ/Q is constant, and perpendicular matches
281        # radial (which instead it should be inverse).
282        Q = np.sqrt(Qx**2 + Qy**2)
283        dqx = resolution * Q
284        dqy = resolution * Q
285    else:
286        dqx = dqy = None
287
288    data = Data2D(x=Qx, y=Qy, z=Iq, dx=dqx, dy=dqy, dz=dIq)
289    data.x_bins = qx
290    data.y_bins = qy
291    data.filename = "fake data"
292
293    # pixel_size in mm, distance in m
294    detector = Detector(pixel_size=(5, 5), distance=4)
295    data.detector.append(detector)
296    data.source.wavelength = 5 # angstroms
297    data.source.wavelength_unit = "A"
298    return data
299
300
301def plot_data(data, view='log', limits=None):
302    """
303    Plot data loaded by the sasview loader.
304
305    *data* is a sasview data object, either 1D, 2D or SESANS.
306
307    *view* is log or linear.
308
309    *limits* sets the intensity limits on the plot; if None then the limits
310    are inferred from the data.
311    """
312    # Note: kind of weird using the plot result functions to plot just the
313    # data, but they already handle the masking and graph markup already, so
314    # do not repeat.
315    if hasattr(data, 'lam'):
316        _plot_result_sesans(data, None, None, use_data=True, limits=limits)
317    elif hasattr(data, 'qx_data'):
318        _plot_result2D(data, None, None, view, use_data=True, limits=limits)
319    else:
320        _plot_result1D(data, None, None, view, use_data=True, limits=limits)
321
322
323def plot_theory(data, theory, resid=None, view='log',
324                use_data=True, limits=None):
325    """
326    Plot theory calculation.
327
328    *data* is needed to define the graph properties such as labels and
329    units, and to define the data mask.
330
331    *theory* is a matrix of the same shape as the data.
332
333    *view* is log or linear
334
335    *use_data* is True if the data should be plotted as well as the theory.
336
337    *limits* sets the intensity limits on the plot; if None then the limits
338    are inferred from the data.
339    """
340    if hasattr(data, 'lam'):
341        _plot_result_sesans(data, theory, resid, use_data=True, limits=limits)
342    elif hasattr(data, 'qx_data'):
343        _plot_result2D(data, theory, resid, view, use_data, limits=limits)
344    else:
345        _plot_result1D(data, theory, resid, view, use_data, limits=limits)
346
347
348def protect(fn):
349    """
350    Decorator to wrap calls in an exception trapper which prints the
351    exception and continues.  Keyboard interrupts are ignored.
352    """
353    def wrapper(*args, **kw):
354        """
355        Trap and print errors from function.
356        """
357        try:
358            return fn(*args, **kw)
359        except KeyboardInterrupt:
360            raise
361        except:
362            traceback.print_exc()
363
364    return wrapper
365
366
367@protect
368def _plot_result1D(data, theory, resid, view, use_data, limits=None):
369    """
370    Plot the data and residuals for 1D data.
371    """
372    import matplotlib.pyplot as plt
373    from numpy.ma import masked_array, masked
374
375    use_data = use_data and data.y is not None
376    use_theory = theory is not None
377    use_resid = resid is not None
378    num_plots = (use_data or use_theory) + use_resid
379
380    scale = data.x**4 if view == 'q4' else 1.0
381
382    if use_data or use_theory:
383        #print(vmin, vmax)
384        all_positive = True
385        some_present = False
386        if use_data:
387            mdata = masked_array(data.y, data.mask.copy())
388            mdata[~np.isfinite(mdata)] = masked
389            if view is 'log':
390                mdata[mdata <= 0] = masked
391            plt.errorbar(data.x/10, scale*mdata, yerr=data.dy, fmt='.')
392            all_positive = all_positive and (mdata > 0).all()
393            some_present = some_present or (mdata.count() > 0)
394
395
396        if use_theory:
397            mtheory = masked_array(theory, data.mask.copy())
398            mtheory[~np.isfinite(mtheory)] = masked
399            if view is 'log':
400                mtheory[mtheory <= 0] = masked
401            plt.plot(data.x/10, scale*mtheory, '-', hold=True)
402            all_positive = all_positive and (mtheory > 0).all()
403            some_present = some_present or (mtheory.count() > 0)
404
405        if limits is not None:
406            plt.ylim(*limits)
407
408        if num_plots > 1:
409            plt.subplot(1, num_plots, 1)
410        plt.xscale('linear' if not some_present else view)
411        plt.yscale('linear'
412                   if view == 'q4' or not some_present or not all_positive
413                   else view)
414        plt.xlabel("$q$/nm$^{-1}$")
415        plt.ylabel('$I(q)$')
416
417    if use_resid:
418        mresid = masked_array(resid, data.mask.copy())
419        mresid[~np.isfinite(mresid)] = masked
420        some_present = (mresid.count() > 0)
421
422        if num_plots > 1:
423            plt.subplot(1, num_plots, (use_data or use_theory) + 1)
424        plt.plot(data.x/10, mresid, '-')
425        plt.xlabel("$q$/nm$^{-1}$")
426        plt.ylabel('residuals')
427        plt.xscale('linear' if not some_present else view)
428
429
430@protect
431def _plot_result_sesans(data, theory, resid, use_data, limits=None):
432    """
433    Plot SESANS results.
434    """
435    import matplotlib.pyplot as plt
436    use_data = use_data and data.y is not None
437    use_theory = theory is not None
438    use_resid = resid is not None
439    num_plots = (use_data or use_theory) + use_resid
440
441    if use_data or use_theory:
442        if num_plots > 1:
443            plt.subplot(1, num_plots, 1)
444        if use_data:
445            plt.errorbar(data.x, data.y, yerr=data.dy)
446        if theory is not None:
447            plt.plot(data.x, theory, '-', hold=True)
448        if limits is not None:
449            plt.ylim(*limits)
450        plt.xlabel('spin echo length (nm)')
451        plt.ylabel('polarization (P/P0)')
452
453    if resid is not None:
454        if num_plots > 1:
455            plt.subplot(1, num_plots, (use_data or use_theory) + 1)
456        plt.plot(data.x, resid, 'x')
457        plt.xlabel('spin echo length (nm)')
458        plt.ylabel('residuals (P/P0)')
459
460
461@protect
462def _plot_result2D(data, theory, resid, view, use_data, limits=None):
463    """
464    Plot the data and residuals for 2D data.
465    """
466    import matplotlib.pyplot as plt
467    use_data = use_data and data.data is not None
468    use_theory = theory is not None
469    use_resid = resid is not None
470    num_plots = use_data + use_theory + use_resid
471
472    # Put theory and data on a common colormap scale
473    vmin, vmax = np.inf, -np.inf
474    if use_data:
475        target = data.data[~data.mask]
476        datamin = target[target > 0].min() if view == 'log' else target.min()
477        datamax = target.max()
478        vmin = min(vmin, datamin)
479        vmax = max(vmax, datamax)
480    if use_theory:
481        theorymin = theory[theory > 0].min() if view == 'log' else theory.min()
482        theorymax = theory.max()
483        vmin = min(vmin, theorymin)
484        vmax = max(vmax, theorymax)
485
486    # Override data limits from the caller
487    if limits is not None:
488        vmin, vmax = limits
489
490    # Plot data
491    if use_data:
492        if num_plots > 1:
493            plt.subplot(1, num_plots, 1)
494        _plot_2d_signal(data, target, view=view, vmin=vmin, vmax=vmax)
495        plt.title('data')
496        h = plt.colorbar()
497        h.set_label('$I(q)$')
498
499    # plot theory
500    if use_theory:
501        if num_plots > 1:
502            plt.subplot(1, num_plots, use_data+1)
503        _plot_2d_signal(data, theory, view=view, vmin=vmin, vmax=vmax)
504        plt.title('theory')
505        h = plt.colorbar()
506        h.set_label(r'$\log_{10}I(q)$' if view == 'log'
507                    else r'$q^4 I(q)$' if view == 'q4'
508                    else '$I(q)$')
509
510    # plot resid
511    if use_resid:
512        if num_plots > 1:
513            plt.subplot(1, num_plots, use_data+use_theory+1)
514        _plot_2d_signal(data, resid, view='linear')
515        plt.title('residuals')
516        h = plt.colorbar()
517        h.set_label(r'$\Delta I(q)$')
518
519
520@protect
521def _plot_2d_signal(data, signal, vmin=None, vmax=None, view='log'):
522    """
523    Plot the target value for the data.  This could be the data itself,
524    the theory calculation, or the residuals.
525
526    *scale* can be 'log' for log scale data, or 'linear'.
527    """
528    import matplotlib.pyplot as plt
529    from numpy.ma import masked_array
530
531    image = np.zeros_like(data.qx_data)
532    image[~data.mask] = signal
533    valid = np.isfinite(image)
534    if view == 'log':
535        valid[valid] = (image[valid] > 0)
536        if vmin is None: vmin = image[valid & ~data.mask].min()
537        if vmax is None: vmax = image[valid & ~data.mask].max()
538        image[valid] = np.log10(image[valid])
539    elif view == 'q4':
540        image[valid] *= (data.qx_data[valid]**2+data.qy_data[valid]**2)**2
541        if vmin is None: vmin = image[valid & ~data.mask].min()
542        if vmax is None: vmax = image[valid & ~data.mask].max()
543    else:
544        if vmin is None: vmin = image[valid & ~data.mask].min()
545        if vmax is None: vmax = image[valid & ~data.mask].max()
546
547    image[~valid | data.mask] = 0
548    #plottable = Iq
549    plottable = masked_array(image, ~valid | data.mask)
550    xmin, xmax = min(data.qx_data)/10, max(data.qx_data)/10
551    ymin, ymax = min(data.qy_data)/10, max(data.qy_data)/10
552    if view == 'log':
553        vmin, vmax = np.log10(vmin), np.log10(vmax)
554    plt.imshow(plottable.reshape(len(data.x_bins), len(data.y_bins)),
555               interpolation='nearest', aspect=1, origin='upper',
556               extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax)
557    plt.xlabel("$q_x$/nm$^{-1}$")
558    plt.ylabel("$q_y$/nm$^{-1}$")
559    return vmin, vmax
560
561def demo():
562    """
563    Load and plot a SAS dataset.
564    """
565    data = load_data('DEC07086.DAT')
566    set_beam_stop(data, 0.004)
567    plot_data(data)
568    import matplotlib.pyplot as plt; plt.show()
569
570
571if __name__ == "__main__":
572    demo()
Note: See TracBrowser for help on using the repository browser.