source: sasmodels/sasmodels/data.py @ 630156b

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

sascomp: improve data file handling; add -help as alias to -html; default to -mono

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