source: sasmodels/sasmodels/data.py @ ba7302a

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since ba7302a was ba7302a, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

Merge branch 'master' into flexible-cylinder

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