source: sasmodels/sasmodels/data.py @ 74b0495

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

update oriented sesans and simultaneous fit examples

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