-                      r5d80bbf
+                      r7cf2cfd
 how far the polydispersity integral extends.
-A variety of helper functions are provided:
-    :func:`load_data` loads a sasview data file.
-    :func:`empty_data1D` creates an empty dataset, which is useful for plotting
-    a theory function before the data is measured.
-    :func:`empty_data2D` creates an empty 2D dataset.
-    :func:`set_beam_stop` masks the beam stop from the data.
-    :func:`set_half` selects the right or left half of the data, which can
-    be useful for shear measurements which have not been properly corrected
-    for path length and reflections.
-    :func:`set_top` cuts the top part off the data.
-    :func:`plot_data` plots the data file.
-    :func:`plot_theory` plots a calculated result from the model.
 """
 import datetime
 import warnings
-import traceback
 import numpy as np
+from bumps.names import Parameter
 from . import sesans
+from .resolution import Perfect1D, Pinhole1D, Slit1D
+from .resolution2d import Pinhole2D
+# CRUFT python 2.6
+if not hasattr(datetime.timedelta, 'total_seconds'):
+    def delay(dt):
+        """Return number date-time delta as number seconds"""
+        return dt.days * 86400 + dt.seconds + 1e-6 * dt.microseconds
+else:
+    def delay(dt):
+        """Return number date-time delta as number seconds"""
+        return dt.total_seconds()
+from . import weights
+from .data import plot_theory
+from .direct_model import DataMixin
 # CRUFT: old style bumps wrapper which doesn't separate data and model
 …
-def tic():
-    """
-    Timer function.
-    Use "toc=tic()" to start the clock and "toc()" to measure
-    a time interval.
-    """
-    then = datetime.datetime.now()
-    return lambda: delay(datetime.datetime.now() - then)
-def load_data(filename):
-    """
-    Load data using a sasview loader.
-    """
-    from sas.dataloader.loader import Loader
-    loader = Loader()
-    data = loader.load(filename)
-    if data is None:
-        raise IOError("Data %r could not be loaded" % filename)
-    return data
-def plot_data(data, view='log'):
-    """
-    Plot data loaded by the sasview loader.
-    """
-    if hasattr(data, 'qx_data'):
-        _plot_2d_signal(data, data.data, view=view)
-    else:
-        # Note: kind of weird using the _plot_result1D to plot just the
-        # data, but it handles the masking and graph markup already, so
-        # do not repeat.
-        _plot_result1D(data, None, None, view)
-def plot_theory(data, theory, view='log'):
-    if hasattr(data, 'qx_data'):
-        _plot_2d_signal(data, theory, view=view)
-    else:
-        _plot_result1D(data, theory, None, view, include_data=False)
-def empty_data1D(q, resolution=0.05):
-    """
-    Create empty 1D data using the given *q* as the x value.
-    *resolution* dq/q defaults to 5%.
-    """
-    from sas.dataloader.data_info import Data1D
-    Iq = 100 * np.ones_like(q)
-    dIq = np.sqrt(Iq)
-    data = Data1D(q, Iq, dx=resolution * q, dy=dIq)
-    data.filename = "fake data"
-    data.qmin, data.qmax = q.min(), q.max()
-    data.mask = np.zeros(len(Iq), dtype='bool')
-    return data
-def empty_data2D(qx, qy=None, resolution=0.05):
-    """
-    Create empty 2D data using the given mesh.
-    If *qy* is missing, create a square mesh with *qy=qx*.
-    *resolution* dq/q defaults to 5%.
-    """
-    from sas.dataloader.data_info import Data2D, Detector
-    if qy is None:
-        qy = qx
-    Qx, Qy = np.meshgrid(qx, qy)
-    Qx, Qy = Qx.flatten(), Qy.flatten()
-    Iq = 100 * np.ones_like(Qx)
-    dIq = np.sqrt(Iq)
-    mask = np.ones(len(Iq), dtype='bool')
-    data = Data2D()
-    data.filename = "fake data"
-    data.qx_data = Qx
-    data.qy_data = Qy
-    data.data = Iq
-    data.err_data = dIq
-    data.mask = mask
-    data.qmin = 1e-16
-    data.qmax = np.inf
-    # 5% dQ/Q resolution
-    if resolution != 0:
-        # https://www.ncnr.nist.gov/staff/hammouda/distance_learning/chapter_15.pdf
-        # Should have an additional constant which depends on distances and
-        # radii of the aperture, pixel dimensions and wavelength spread
-        # Instead, assume radial dQ/Q is constant, and perpendicular matches
-        # radial (which instead it should be inverse).
-        Q = np.sqrt(Qx**2 + Qy**2)
-        data.dqx_data = resolution * Q
-        data.dqy_data = resolution * Q
-    detector = Detector()
-    detector.pixel_size.x = 5 # mm
-    detector.pixel_size.y = 5 # mm
-    detector.distance = 4 # m
-    data.detector.append(detector)
-    data.xbins = qx
-    data.ybins = qy
-    data.source.wavelength = 5 # angstroms
-    data.source.wavelength_unit = "A"
-    data.Q_unit = "1/A"
-    data.I_unit = "1/cm"
-    data.q_data = np.sqrt(Qx ** 2 + Qy ** 2)
-    data.xaxis("Q_x", "A^{-1}")
-    data.yaxis("Q_y", "A^{-1}")
-    data.zaxis("Intensity", r"\text{cm}^{-1}")
-    return data
-def set_beam_stop(data, radius, outer=None):
-    """
-    Add a beam stop of the given *radius*.  If *outer*, make an annulus.
-    """
-    from sas.dataloader.manipulations import Ringcut
-    if hasattr(data, 'qx_data'):
-        data.mask = Ringcut(0, radius)(data)
-        if outer is not None:
-            data.mask += Ringcut(outer, np.inf)(data)
-    else:
-        data.mask = (data.x >= radius)
-        if outer is not None:
-            data.mask &= (data.x < outer)
-def set_half(data, half):
-    """
-    Select half of the data, either "right" or "left".
-    """
-    from sas.dataloader.manipulations import Boxcut
-    if half == 'right':
-        data.mask += \
-            Boxcut(x_min=-np.inf, x_max=0.0, y_min=-np.inf, y_max=np.inf)(data)
-    if half == 'left':
-        data.mask += \
-            Boxcut(x_min=0.0, x_max=np.inf, y_min=-np.inf, y_max=np.inf)(data)
-def set_top(data, cutoff):
-    """
-    Chop the top off the data, above *cutoff*.
-    """
-    from sas.dataloader.manipulations import Boxcut
-    data.mask += \
-        Boxcut(x_min=-np.inf, x_max=np.inf, y_min=-np.inf, y_max=cutoff)(data)
-def protect(fn):
-    def wrapper(*args, **kw):
-        try:
-            return fn(*args, **kw)
-        except:
-            traceback.print_exc()
-    return wrapper
-@protect
-def _plot_result1D(data, theory, resid, view, include_data=True):
-    """
-    Plot the data and residuals for 1D data.
-    """
-    import matplotlib.pyplot as plt
-    from numpy.ma import masked_array, masked
-    #print "not a number",sum(np.isnan(data.y))
-    #data.y[data.y<0.05] = 0.5
-    mdata = masked_array(data.y, data.mask)
-    mdata[~np.isfinite(mdata)] = masked
-    if view is 'log':
-        mdata[mdata <= 0] = masked
-    scale = data.x**4 if view == 'q4' else 1.0
-    if resid is not None:
-        plt.subplot(121)
-    positive = False
-    if include_data:
-        plt.errorbar(data.x, scale*mdata, yerr=data.dy, fmt='.')
-        positive = positive or (mdata>0).any()
-    if theory is not None:
-        mtheory = masked_array(theory, mdata.mask)
-        plt.plot(data.x, scale*mtheory, '-', hold=True)
-        positive = positive or (mtheory>0).any()
-    plt.xscale(view)
-    plt.yscale('linear' if view == 'q4' or not positive else view)
-    plt.xlabel('Q')
-    plt.ylabel('I(Q)')
-    if resid is not None:
-        mresid = masked_array(resid, mdata.mask)
-        plt.subplot(122)
-        plt.plot(data.x, mresid, 'x')
-        plt.ylabel('residuals')
-        plt.xlabel('Q')
-        plt.xscale(view)
-# pylint: disable=unused-argument
-@protect
-def _plot_sesans(data, theory, resid, view):
-    import matplotlib.pyplot as plt
-    plt.subplot(121)
-    plt.errorbar(data.x, data.y, yerr=data.dy)
-    plt.plot(data.x, theory, '-', hold=True)
-    plt.xlabel('spin echo length (nm)')
-    plt.ylabel('polarization (P/P0)')
-    plt.subplot(122)
-    plt.plot(data.x, resid, 'x')
-    plt.xlabel('spin echo length (nm)')
-    plt.ylabel('residuals (P/P0)')
-@protect
-def _plot_result2D(data, theory, resid, view):
-    """
-    Plot the data and residuals for 2D data.
-    """
-    import matplotlib.pyplot as plt
-    target = data.data[~data.mask]
-    if view == 'log':
-        vmin = min(target[target>0].min(), theory[theory>0].min())
-        vmax = max(target.max(), theory.max())
-    else:
-        vmin = min(target.min(), theory.min())
-        vmax = max(target.max(), theory.max())
-    #print vmin, vmax
-    plt.subplot(131)
-    _plot_2d_signal(data, target, view=view, vmin=vmin, vmax=vmax)
-    plt.title('data')
-    plt.colorbar()
-    plt.subplot(132)
-    _plot_2d_signal(data, theory, view=view, vmin=vmin, vmax=vmax)
-    plt.title('theory')
-    plt.colorbar()
-    plt.subplot(133)
-    _plot_2d_signal(data, resid, view='linear')
-    plt.title('residuals')
-    plt.colorbar()
-@protect
-def _plot_2d_signal(data, signal, vmin=None, vmax=None, view='log'):
-    """
-    Plot the target value for the data.  This could be the data itself,
-    the theory calculation, or the residuals.
-    *scale* can be 'log' for log scale data, or 'linear'.
-    """
-    import matplotlib.pyplot as plt
-    from numpy.ma import masked_array
-    image = np.zeros_like(data.qx_data)
-    image[~data.mask] = signal
-    valid = np.isfinite(image)
-    if view == 'log':
-        valid[valid] = (image[valid] > 0)
-        image[valid] = np.log10(image[valid])
-    elif view == 'q4':
-        image[valid] *= (data.qx_data[valid]**2+data.qy_data[valid]**2)**2
-    image[~valid | data.mask] = 0
-    #plottable = Iq
-    plottable = masked_array(image, ~valid | data.mask)
-    xmin, xmax = min(data.qx_data), max(data.qx_data)
-    ymin, ymax = min(data.qy_data), max(data.qy_data)
-    # TODO: fix vmin, vmax so it is shared for theory/resid
-    vmin = vmax = None
-    try:
-        if vmin is None: vmin = image[valid & ~data.mask].min()
-        if vmax is None: vmax = image[valid & ~data.mask].max()
-    except:
-        vmin, vmax = 0, 1
-    #print vmin,vmax
-    plt.imshow(plottable.reshape(128, 128),
-               interpolation='nearest', aspect=1, origin='upper',
-               extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax)
 class Model(object):
+    def __init__(self, kernel, **kw):
+        from bumps.names import Parameter
+        self.kernel = kernel
+        partype = kernel.info['partype']
+    def __init__(self, model, **kw):
+        self._sasmodel = model
+        partype = model.info['partype']
         pars = []
         for p in kernel.info['parameters']:
+        for p in model.info['parameters']:
             name, default, limits = p[0], p[2], p[3]
             value = kw.pop(name, default)
 …
         return dict((k, getattr(self, k)) for k in self._parameter_names)
+class Experiment(object):
+class Experiment(DataMixin):
     """
     Return a bumps wrapper for a SAS model.
 …
         # remember inputs so we can inspect from outside
-        self.data = data
         self.model = model
         self.cutoff = cutoff
+        if hasattr(data, 'lam'):
+            self.data_type = 'sesans'
+        elif hasattr(data, 'qx_data'):
+            self.data_type = 'Iqxy'
+        else:
+            self.data_type = 'Iq'
+        # interpret data
+        partype = model.kernel.info['partype']
+        if self.data_type == 'sesans':
+            q = sesans.make_q(data.sample.zacceptance, data.Rmax)
+            self.index = slice(None, None)
+            self.Iq = data.y
+            self.dIq = data.dy
+            #self._theory = np.zeros_like(q)
+            q_vectors = [q]
+        elif self.data_type == 'Iqxy':
+            q = np.sqrt(data.qx_data**2 + data.qy_data**2)
+            qmin = getattr(data, 'qmin', 1e-16)
+            qmax = getattr(data, 'qmax', np.inf)
+            accuracy = getattr(data, 'accuracy', 'Low')
+            self.index = (~data.mask) & (~np.isnan(data.data)) \
+                         & (q >= qmin) & (q <= qmax)
+            self.Iq = data.data[self.index]
+            self.dIq = data.err_data[self.index]
+            self.resolution = Pinhole2D(data=data, index=self.index,
+                                        nsigma=3.0, accuracy=accuracy)
+            #self._theory = np.zeros_like(self.Iq)
+            if not partype['orientation'] and not partype['magnetic']:
+                raise ValueError("not 2D without orientation or magnetic parameters")
+                #qx,qy = self.resolution.q_calc
+                #q_vectors = [np.sqrt(qx**2 + qy**2)]
+            else:
+                q_vectors = self.resolution.q_calc
+        elif self.data_type == 'Iq':
+            self.index = (data.x >= data.qmin) & (data.x <= data.qmax) & ~np.isnan(data.y)
+            self.Iq = data.y[self.index]
+            self.dIq = data.dy[self.index]
+            if getattr(data, 'dx', None) is not None:
+                q, dq = data.x[self.index], data.dx[self.index]
+                if (dq>0).any():
+                    self.resolution = Pinhole1D(q, dq)
+                else:
+                    self.resolution = Perfect1D(q)
+            elif (getattr(data, 'dxl', None) is not None and
+                  getattr(data, 'dxw', None) is not None):
+                q = data.x[self.index]
+                width = data.dxh[self.index]  # Note: dx
+                self.resolution = Slit1D(data.x[self.index],
+                                         width=data.dxh[self.index],
+                                         height=data.dxw[self.index])
+            else:
+                self.resolution = Perfect1D(data.x[self.index])
+            #self._theory = np.zeros_like(self.Iq)
+            q_vectors = [self.resolution.q_calc]
+        else:
+            raise ValueError("Unknown data type") # never gets here
+        # Remember function inputs so we can delay loading the function and
+        # so we can save/restore state
+        self._fn_inputs = [v for v in q_vectors]
+        self._fn = None
+        self._interpret_data(data, model._sasmodel)
         self.update()
 …
     def theory(self):
         if 'theory' not in self._cache:
+            pars = dict((k, v.value) for k,v in self.model.parameters().items())
+            self._cache['theory'] = self._calc_theory(pars, cutoff=self.cutoff)
+            """
             if self._fn is None:
                 q_input = self.model.kernel.make_input(self._fn_inputs)
+                q_input = self.model.kernel.make_input(self._kernel_inputs)
                 self._fn = self.model.kernel(q_input)
 …
                 result = sesans.hankel(self.data.x, self.data.lam * 1e-9,
                                        self.data.sample.thickness / 10,
                                        self._fn_inputs[0], Iq_calc)
+                                       self._kernel_inputs[0], Iq_calc)
                 self._cache['theory'] = result
             else:
                 Iq = self.resolution.apply(Iq_calc)
                 self._cache['theory'] = Iq
+            """
         return self._cache['theory']
 …
         Plot the data and residuals.
         """
+        data, theory, resid = self.data, self.theory(), self.residuals()
+        if self.data_type == 'Iq':
+            _plot_result1D(data, theory, resid, view)
+        elif self.data_type == 'Iqxy':
+            _plot_result2D(data, theory, resid, view)
+        elif self.data_type == 'sesans':
+            _plot_sesans(data, theory, resid, view)
+        else:
+            raise ValueError("Unknown data type")
+        data, theory, resid = self._data, self.theory(), self.residuals()
+        plot_theory(data, theory, resid, view)
     def simulate_data(self, noise=None):
+        theory = self.theory()
+        if noise is not None:
+            self.dIq = theory*noise*0.01
+        dy = self.dIq
+        y = theory + np.random.randn(*dy.shape) * dy
+        self.Iq = y
+        if self.data_type == 'Iq':
+            self.data.dy[self.index] = dy
+            self.data.y[self.index] = y
+        elif self.data_type == 'Iqxy':
+            self.data.data[self.index] = y
+        elif self.data_type == 'sesans':
+            self.data.y[self.index] = y
+        else:
+            raise ValueError("Unknown model")
+        Iq = self.theory()
+        self._set_data(Iq, noise)
     def save(self, basename):
         pass
     def _get_weights(self, par):
+    def remove_get_weights(self, name):
         """
         Get parameter dispersion weights
         """
+        from . import weights
+        relative = self.model.kernel.info['partype']['pd-rel']
+        limits = self.model.kernel.info['limits']
+        info = self.model.kernel.info
+        relative = name in info['partype']['pd-rel']
+        limits = info['limits'][name]
         disperser, value, npts, width, nsigma = [
             getattr(self.model, par + ext)
+            getattr(self.model, name + ext)
             for ext in ('_pd_type', '', '_pd_n', '_pd', '_pd_nsigma')]
         value, weight = weights.get_weights(
             disperser, int(npts.value), width.value, nsigma.value,
             value.value, limits[par], par in relative)
+            value.value, limits, relative)
         return value, weight / np.sum(weight)
 …
         # Can't pickle gpu functions, so instead make them lazy
         state = self.__dict__.copy()
         state['_fn'] = None
+        state['_kernel'] = None
         return state
 …
         # pylint: disable=attribute-defined-outside-init
         self.__dict__ = state
-def demo():
-    data = load_data('DEC07086.DAT')
-    set_beam_stop(data, 0.004)
-    plot_data(data)
-    import matplotlib.pyplot as plt; plt.show()
-if __name__ == "__main__":
-    demo()

Note: See TracChangeset for help on using the changeset viewer.

SasView

Changeset 7cf2cfd in sasmodels for sasmodels/bumps_model.py

Legend:

sasmodels/bumps_model.py

Download in other formats: