source: sasmodels/sasmodels/data.py @ e65c3ba

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

lint

  • Property mode set to 100644
File size: 21.8 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') and not getattr(data, 'radial', False):
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') and not getattr(data, 'radial', False):
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    if getattr(data, 'radial', False):
429        data.x = data.q_data
430        data.y = data.data
431
432    use_data = use_data and data.y is not None
433    use_theory = theory is not None
434    use_resid = resid is not None
435    use_calc = use_theory and Iq_calc is not None
436    num_plots = (use_data or use_theory) + use_calc + use_resid
437    non_positive_x = (data.x <= 0.0).any()
438
439    scale = data.x**4 if view == 'q4' else 1.0
440    xscale = yscale = 'linear' if view == 'linear' else 'log'
441
442    if use_data or use_theory:
443        if num_plots > 1:
444            plt.subplot(1, num_plots, 1)
445
446        #print(vmin, vmax)
447        all_positive = True
448        some_present = False
449        if use_data:
450            mdata = masked_array(data.y, data.mask.copy())
451            mdata[~np.isfinite(mdata)] = masked
452            if view is 'log':
453                mdata[mdata <= 0] = masked
454            plt.errorbar(data.x, scale*mdata, yerr=data.dy, fmt='.')
455            all_positive = all_positive and (mdata > 0).all()
456            some_present = some_present or (mdata.count() > 0)
457
458
459        if use_theory:
460            # Note: masks merge, so any masked theory points will stay masked,
461            # and the data mask will be added to it.
462            mtheory = masked_array(theory, data.mask.copy())
463            mtheory[~np.isfinite(mtheory)] = masked
464            if view is 'log':
465                mtheory[mtheory <= 0] = masked
466            plt.plot(data.x, scale*mtheory, '-')
467            all_positive = all_positive and (mtheory > 0).all()
468            some_present = some_present or (mtheory.count() > 0)
469
470        if limits is not None:
471            plt.ylim(*limits)
472
473
474        xscale = ('linear' if not some_present or non_positive_x
475                  else view if view is not None
476                  else 'log')
477        yscale = ('linear'
478                  if view == 'q4' or not some_present or not all_positive
479                  else view if view is not None
480                  else 'log')
481        plt.xscale(xscale)
482        plt.xlabel("$q$/A$^{-1}$")
483        plt.yscale(yscale)
484        plt.ylabel('$I(q)$')
485        title = ("data and model" if use_theory and use_data
486                 else "data" if use_data
487                 else "model")
488        plt.title(title)
489
490    if use_calc:
491        # Only have use_calc if have use_theory
492        plt.subplot(1, num_plots, 2)
493        qx, qy, Iqxy = Iq_calc
494        plt.pcolormesh(qx, qy[qy > 0], np.log10(Iqxy[qy > 0, :]))
495        plt.xlabel("$q_x$/A$^{-1}$")
496        plt.xlabel("$q_y$/A$^{-1}$")
497        plt.xscale('log')
498        plt.yscale('log')
499        #plt.axis('equal')
500
501    if use_resid:
502        mresid = masked_array(resid, data.mask.copy())
503        mresid[~np.isfinite(mresid)] = masked
504        some_present = (mresid.count() > 0)
505
506        if num_plots > 1:
507            plt.subplot(1, num_plots, use_calc + 2)
508        plt.plot(data.x, mresid, '.')
509        plt.xlabel("$q$/A$^{-1}$")
510        plt.ylabel('residuals')
511        plt.title('(model - Iq)/dIq')
512        plt.xscale(xscale)
513        plt.yscale('linear')
514
515
516@protect
517def _plot_result_sesans(data, theory, resid, use_data, limits=None):
518    # type: (SesansData, Optional[np.ndarray], Optional[np.ndarray], bool, Optional[Tuple[float, float]]) -> None
519    """
520    Plot SESANS results.
521    """
522    import matplotlib.pyplot as plt  # type: ignore
523    use_data = use_data and data.y is not None
524    use_theory = theory is not None
525    use_resid = resid is not None
526    num_plots = (use_data or use_theory) + use_resid
527
528    if use_data or use_theory:
529        is_tof = (data.lam != data.lam[0]).any()
530        if num_plots > 1:
531            plt.subplot(1, num_plots, 1)
532        if use_data:
533            if is_tof:
534                plt.errorbar(data.x, np.log(data.y)/(data.lam*data.lam),
535                             yerr=data.dy/data.y/(data.lam*data.lam))
536            else:
537                plt.errorbar(data.x, data.y, yerr=data.dy)
538        if theory is not None:
539            if is_tof:
540                plt.plot(data.x, np.log(theory)/(data.lam*data.lam), '-')
541            else:
542                plt.plot(data.x, theory, '-')
543        if limits is not None:
544            plt.ylim(*limits)
545
546        plt.xlabel('spin echo length ({})'.format(data._xunit))
547        if is_tof:
548            plt.ylabel(r'(Log (P/P$_0$))/$\lambda^2$')
549        else:
550            plt.ylabel('polarization (P/P0)')
551
552
553    if resid is not None:
554        if num_plots > 1:
555            plt.subplot(1, num_plots, (use_data or use_theory) + 1)
556        plt.plot(data.x, resid, 'x')
557        plt.xlabel('spin echo length ({})'.format(data._xunit))
558        plt.ylabel('residuals (P/P0)')
559
560
561@protect
562def _plot_result2D(data, theory, resid, view, use_data, limits=None):
563    # type: (Data2D, Optional[np.ndarray], Optional[np.ndarray], str, bool, Optional[Tuple[float,float]]) -> None
564    """
565    Plot the data and residuals for 2D data.
566    """
567    import matplotlib.pyplot as plt  # type: ignore
568    use_data = use_data and data.data is not None
569    use_theory = theory is not None
570    use_resid = resid is not None
571    num_plots = use_data + use_theory + use_resid
572
573    # Put theory and data on a common colormap scale
574    vmin, vmax = np.inf, -np.inf
575    target = None # type: Optional[np.ndarray]
576    if use_data:
577        target = data.data[~data.mask]
578        datamin = target[target > 0].min() if view == 'log' else target.min()
579        datamax = target.max()
580        vmin = min(vmin, datamin)
581        vmax = max(vmax, datamax)
582    if use_theory:
583        theorymin = theory[theory > 0].min() if view == 'log' else theory.min()
584        theorymax = theory.max()
585        vmin = min(vmin, theorymin)
586        vmax = max(vmax, theorymax)
587
588    # Override data limits from the caller
589    if limits is not None:
590        vmin, vmax = limits
591
592    # Plot data
593    if use_data:
594        if num_plots > 1:
595            plt.subplot(1, num_plots, 1)
596        _plot_2d_signal(data, target, view=view, vmin=vmin, vmax=vmax)
597        plt.title('data')
598        h = plt.colorbar()
599        h.set_label('$I(q)$')
600
601    # plot theory
602    if use_theory:
603        if num_plots > 1:
604            plt.subplot(1, num_plots, use_data+1)
605        _plot_2d_signal(data, theory, view=view, vmin=vmin, vmax=vmax)
606        plt.title('theory')
607        h = plt.colorbar()
608        h.set_label(r'$\log_{10}I(q)$' if view == 'log'
609                    else r'$q^4 I(q)$' if view == 'q4'
610                    else '$I(q)$')
611
612    # plot resid
613    if use_resid:
614        if num_plots > 1:
615            plt.subplot(1, num_plots, use_data+use_theory+1)
616        _plot_2d_signal(data, resid, view='linear')
617        plt.title('residuals')
618        h = plt.colorbar()
619        h.set_label(r'$\Delta I(q)$')
620
621
622@protect
623def _plot_2d_signal(data, signal, vmin=None, vmax=None, view='log'):
624    # type: (Data2D, np.ndarray, Optional[float], Optional[float], str) -> Tuple[float, float]
625    """
626    Plot the target value for the data.  This could be the data itself,
627    the theory calculation, or the residuals.
628
629    *scale* can be 'log' for log scale data, or 'linear'.
630    """
631    import matplotlib.pyplot as plt  # type: ignore
632    from numpy.ma import masked_array  # type: ignore
633
634    image = np.zeros_like(data.qx_data)
635    image[~data.mask] = signal
636    valid = np.isfinite(image)
637    if view == 'log':
638        valid[valid] = (image[valid] > 0)
639        if vmin is None: vmin = image[valid & ~data.mask].min()
640        if vmax is None: vmax = image[valid & ~data.mask].max()
641        image[valid] = np.log10(image[valid])
642    elif view == 'q4':
643        image[valid] *= (data.qx_data[valid]**2+data.qy_data[valid]**2)**2
644        if vmin is None: vmin = image[valid & ~data.mask].min()
645        if vmax is None: vmax = image[valid & ~data.mask].max()
646    else:
647        if vmin is None: vmin = image[valid & ~data.mask].min()
648        if vmax is None: vmax = image[valid & ~data.mask].max()
649
650    image[~valid | data.mask] = 0
651    #plottable = Iq
652    plottable = masked_array(image, ~valid | data.mask)
653    # Divide range by 10 to convert from angstroms to nanometers
654    xmin, xmax = min(data.qx_data), max(data.qx_data)
655    ymin, ymax = min(data.qy_data), max(data.qy_data)
656    if view == 'log':
657        vmin_scaled, vmax_scaled= np.log10(vmin), np.log10(vmax)
658    else:
659        vmin_scaled, vmax_scaled = vmin, vmax
660    plt.imshow(plottable.reshape(len(data.x_bins), len(data.y_bins)),
661               interpolation='nearest', aspect=1, origin='lower',
662               extent=[xmin, xmax, ymin, ymax],
663               vmin=vmin_scaled, vmax=vmax_scaled)
664    plt.xlabel("$q_x$/A$^{-1}$")
665    plt.ylabel("$q_y$/A$^{-1}$")
666    return vmin, vmax
667
668def demo():
669    # type: () -> None
670    """
671    Load and plot a SAS dataset.
672    """
673    data = load_data('DEC07086.DAT')
674    set_beam_stop(data, 0.004)
675    plot_data(data)
676    import matplotlib.pyplot as plt  # type: ignore
677    plt.show()
678
679
680if __name__ == "__main__":
681    demo()
Note: See TracBrowser for help on using the repository browser.