Merge pull request #465 from SasView/fix-sesans-simulations

Add -sesans option to sascomp simulator

sasmodels/compare.py

@@ -42,7 +42,7 @@
 from . import kerneldll
 from . import kernelcl
 from . import kernelcuda
-from .data import plot_theory, empty_data1D, empty_data2D, load_data
+from .data import plot_theory, empty_data1D, empty_data2D, empty_sesans, load_data
 from .direct_model import DirectModel, get_mesh
 from .generate import FLOAT_RE, set_integration_size
 
@@ -80,9 +80,9 @@ def __call__(self, **par: float) -> np.ndarray: ...
     -noise=0 sets the measurement error dI/I
     -res=0 sets the resolution width dQ/Q if calculating with resolution
     -lowq*/-midq/-highq/-exq use q values up to 0.05, 0.2, 1.0, 10.0
-    -q=min:max alternative specification of qrange
+    -q=min:max alternative specification of qrange; sesans uses 1/qmax:1/qmin
     -nq=128 sets the number of Q points in the data set
-    -1d*/-2d computes 1d or 2d data
+    -1d*/-2d/-sesans computes 1d, 2d or sesans data
     -zero indicates that q=0 should be included
 
     === model parameters ===
@@ -685,7 +685,7 @@ def make_data(opts):
     Generate an empty dataset, used with the model to set Q points
     and resolution.
 
-    *opts* contains the options, with 'qmax', 'nq', 'res',
+    *opts* contains the options, with 'qmax', 'nq', 'sesans', 'res',
     'accuracy', 'is2d' and 'view' parsed from the command line.
     """
     qmin, qmax, nq, res = opts['qmin'], opts['qmax'], opts['nq'], opts['res']
@@ -695,6 +695,13 @@ def make_data(opts):
         data.accuracy = opts['accuracy']
         set_beam_stop(data, qmin)
         index = ~data.mask
+    elif opts['sesans']:
+        if opts['view'] == 'log':
+            z = np.logspace(-math.log10(qmax), -math.log10(qmin), nq)
+        else:
+            z = np.linspace(1/qmax, 1/qmin, nq)
+        data = empty_sesans(z)
+        index = slice(None, None)
     else:
         if opts['view'] == 'log' and not opts['zero']:
             q = np.logspace(math.log10(qmin), math.log10(qmax), nq)
@@ -761,6 +768,9 @@ def compare(opts, limits=None, maxdim=None):
 
     *maxdim* **DEPRECATED** Use opts['maxdim'] instead.
     """
+    if (opts['sesans'] and limits is None) or opts['sets'] > 1:
+        limits = (-np.inf, np.inf)
+
     # CRUFT: remove maxdim parameter
     if maxdim is not None:
         opts['maxdim'] = maxdim
@@ -1016,7 +1026,7 @@ def plot_models(opts, result, limits=None, setnum=0):
     # Data generation
     'data=', 'noise=', 'res=', 'nq=', 'q=',
     'lowq', 'midq', 'highq', 'exq', 'zero',
-    '2d', '1d',
+    '2d', '1d', 'sesans',
 
     # Parameter set
     'preset', 'random', 'random=', 'sets=',
@@ -1176,6 +1186,7 @@ def parse_opts(argv):
         'qmin'      : None,
         'qmax'      : 0.05,
         'nq'        : 128,
+        'sesans'    : False,
         'res'       : '0.0',
         'noise'     : 0.0,
         'accuracy'  : 'Low',
@@ -1213,6 +1224,7 @@ def parse_opts(argv):
         elif arg == '-highq':   opts['qmax'] = 1.0
         elif arg == '-midq':    opts['qmax'] = 0.2
         elif arg == '-lowq':    opts['qmax'] = 0.05
+        elif arg == '-sesans':  opts['sesans'] = True
         elif arg == '-zero':    opts['zero'] = True
         elif arg.startswith('-nq='):       opts['nq'] = int(arg[4:])
         elif arg.startswith('-q='):

sasmodels/data.py

@@ -33,6 +33,7 @@
 
 """
 import traceback
+from functools import wraps
 
 import numpy as np  # type: ignore
 from numpy import sqrt, sin, cos, pi
@@ -43,6 +44,7 @@
     Data = Union["Data1D", "Data2D", "SesansData"]
     OptArray = Optional[np.ndarray]
     OptLimits = Optional[Tuple[float, float]]
+    OptString = Optional[str]
 except ImportError:
     pass
 # pylint: enable=unused-import
@@ -64,7 +66,9 @@ def load_data(filename, index=0):
     if not isinstance(datasets, list):
         datasets = [datasets]
     for data in datasets:
-        if hasattr(data, 'x'):
+        if getattr(data, 'isSesans', False):
+            pass
+        elif hasattr(data, 'x'):
             data.qmin, data.qmax = data.x.min(), data.x.max()
             data.mask = (np.isnan(data.y) if data.y is not None
                          else np.zeros_like(data.x, dtype='bool'))
@@ -113,6 +117,11 @@ def set_top(data, cutoff):
         Boxcut(x_min=-np.inf, x_max=np.inf, y_min=-np.inf, y_max=cutoff)(data)
 
 
+class Source:
+    ...
+class Sample:
+    ...
+
 class Data1D(object):
     """
     1D data object.
@@ -182,6 +191,8 @@ class SesansData(Data1D):
     def __init__(self, **kw):
         Data1D.__init__(self, **kw)
         self.lam = None  # type: OptArray
+        self.xaxis("SE length", "A")
+        self.yaxis("log(P)/(t L^2)", "1/A^2 1/cm")
 
 class Data2D(object):
     """
@@ -298,6 +309,25 @@ def __init__(self):
         # type: () -> None
         pass
 
+def empty_sesans(z, wavelength=None, zacceptance=None):
+    data = SesansData(x=z, y=None, dx=None, dy=None)
+    data.filename = "fake data"
+    DEFAULT_WAVELENGTH = 5
+    if wavelength is None:
+        wavelength = DEFAULT_WAVELENGTH
+    if np.isscalar(wavelength):
+        wavelength = np.full_like(z, wavelength)
+    if zacceptance is None:
+        zacceptance = (np.pi/2, 'radians')
+    source = Source()
+    source.wavelength = wavelength
+    source.wavelength_unit = "A"
+    sample = Sample()
+    sample.zacceptance = zacceptance
+    data.source = source
+    data.sample = sample
+    return data
+
 def empty_data1D(q, resolution=0.0, L=0., dL=0.):
     # type: (np.ndarray, float, float, float) -> Data1D
     r"""
@@ -379,7 +409,7 @@ def plot_data(data, view=None, limits=None):
 
     *data* is a sasview data object, either 1D, 2D or SESANS.
 
-    *view* is log or linear.
+    *view* is log, linear or normed.
 
     *limits* sets the intensity limits on the plot; if None then the limits
     are inferred from the data.
@@ -388,7 +418,7 @@ def plot_data(data, view=None, limits=None):
     # data, but they already handle the masking and graph markup already, so
     # do not repeat.
     if hasattr(data, 'isSesans') and data.isSesans:
-        _plot_result_sesans(data, None, None, use_data=True, limits=limits)
+        _plot_result_sesans(data, None, None, view, use_data=True, limits=limits)
     elif hasattr(data, 'qx_data') and not getattr(data, 'radial', False):
         _plot_result2D(data, None, None, view, use_data=True, limits=limits)
     else:
@@ -397,7 +427,7 @@ def plot_data(data, view=None, limits=None):
 
 def plot_theory(data, theory, resid=None, view=None, use_data=True,
                 limits=None, Iq_calc=None):
-    # type: (Data, OptArray, OptArray, str, bool, OptLimits, OptArray) -> None
+    # type: (Data, OptArray, OptArray, OptString, bool, OptLimits, OptArray) -> None
     """
     Plot theory calculation.
 
@@ -406,17 +436,19 @@ def plot_theory(data, theory, resid=None, view=None, use_data=True,
 
     *theory* is a matrix of the same shape as the data.
 
-    *view* is log or linear
+    *view* is log, linear or normed
 
     *use_data* is True if the data should be plotted as well as the theory.
 
     *limits* sets the intensity limits on the plot; if None then the limits
-    are inferred from the data.
+    are inferred from the data. If (-inf, inf) then use auto limits.
 
     *Iq_calc* is the raw theory values without resolution smearing
     """
+    if limits is not None and np.isinf(limits[0]):
+        limits = None
     if hasattr(data, 'isSesans') and data.isSesans:
-        _plot_result_sesans(data, theory, resid, use_data=True, limits=limits)
+        _plot_result_sesans(data, theory, resid, view, use_data=True, limits=limits)
     elif hasattr(data, 'qx_data') and not getattr(data, 'radial', False):
         _plot_result2D(data, theory, resid, view, use_data, limits=limits)
     else:
@@ -430,14 +462,17 @@ def protect(func):
     Decorator to wrap calls in an exception trapper which prints the
     exception and continues.  Keyboard interrupts are ignored.
     """
+    @wraps(func)
     def wrapper(*args, **kw):
         """
         Trap and print errors from function.
         """
         try:
             return func(*args, **kw)
-        except Exception:
-            traceback.print_exc()
+        except Exception as exc:
+            print("Traceback (most recent call last):")
+            print("".join(traceback.format_list(traceback.extract_stack(limit=4)[:-1]))[:-1])
+            print(f"{exc.__class__.__name__}: {exc}")
 
     return wrapper
 
@@ -543,8 +578,8 @@ def _plot_result1D(data, theory, resid, view, use_data,
 
 
 @protect
-def _plot_result_sesans(data, theory, resid, use_data, limits=None):
-    # type: (SesansData, OptArray, OptArray, bool, OptLimits) -> None
+def _plot_result_sesans(data, theory, resid, view, use_data, limits=None):
+    # type: (SesansData, OptArray, OptArray, OptString, bool, OptLimits) -> None
     """
     Plot SESANS results.
     """
@@ -554,6 +589,12 @@ def _plot_result_sesans(data, theory, resid, use_data, limits=None):
     use_resid = resid is not None
     num_plots = (use_data or use_theory) + use_resid
 
+    normed = (view == "normed")
+    #normed = True
+    offset, scale = 0, 1
+    if normed and theory is not None:
+        offset, scale = theory[-1], theory[0] - theory[-1]
+
     if use_data or use_theory:
         is_tof = data.lam is not None and (data.lam != data.lam[0]).any()
         if num_plots > 1:
@@ -563,28 +604,29 @@ def _plot_result_sesans(data, theory, resid, use_data, limits=None):
                 plt.errorbar(data.x, np.log(data.y)/(data.lam*data.lam),
                              yerr=data.dy/data.y/(data.lam*data.lam))
             else:
-                plt.errorbar(data.x, data.y, yerr=data.dy)
+                #plt.errorbar(data.x, data.y, yerr=data.dy)
+                plt.errorbar(data.x, (data.y-offset)/scale, yerr=data.dy/scale)
         if theory is not None:
             if is_tof:
                 plt.plot(data.x, np.log(theory)/(data.lam*data.lam), '-')
             else:
-                plt.plot(data.x, theory, '-')
+                #plt.plot(data.x, theory, '-')
+                plt.plot(data.x, (theory-offset)/scale, '-')
         if limits is not None:
             plt.ylim(*limits)
 
         plt.xlabel('spin echo length ({})'.format(data._xunit))
-        if is_tof:
-            plt.ylabel(r'(Log (P/P$_0$))/$\lambda^2$')
-        else:
-            plt.ylabel('polarization (P/P0)')
+        plt.ylabel(r'$\log(P)/(t\lambda^2) (\mathrm{A}^{-2}\mathrm{cm}^{-1})$')
+        plt.xscale('log')
 
 
     if resid is not None:
         if num_plots > 1:
             plt.subplot(1, num_plots, (use_data or use_theory) + 1)
         plt.plot(data.x, resid, 'x')
         plt.xlabel('spin echo length ({})'.format(data._xunit))
-        plt.ylabel('residuals (P/P0)')
+        plt.ylabel('polarization residuals')
+        plt.xscale('log')
 
 
 @protect

sasmodels/direct_model.py

@@ -159,19 +159,19 @@ def _vol_pars(model_info, values):
 
 
 def _make_sesans_transform(data):
-    from sas.sascalc.data_util.nxsunit import Converter
-
     # Pre-compute the Hankel matrix (H)
-    SElength = Converter(data._xunit)(data.x, "A")
-
-    theta_max = Converter("radians")(data.sample.zacceptance)[0]
-    q_max = 2 * np.pi / np.max(data.source.wavelength) * np.sin(theta_max)
-    zaccept = Converter("1/A")(q_max, "1/" + data.source.wavelength_unit)
+    SElength, SEunits = data.x, data._xunit
+    wavelength, wunits = data.source.wavelength, data.source.wavelength_unit
+    theta_max, theta_units = data.sample.zacceptance
+    if SEunits != "A" or wunits != "A" or theta_units != "radians":
+        from sas.sascalc.data_util.nxsunit import Converter
+        SElength = Converter("A")(SElength, units=SEunits)
+        wavelength = Converter("A")(wavelength, units=wunits)
+        theta_max = Converter("radian")(theta_max, units=theta_units)
 
     Rmax = 10000000
-    hankel = sesans.SesansTransform(data.x, SElength,
-                                    data.source.wavelength,
-                                    zaccept, Rmax)
+    zaccept = 2 * np.pi / np.max(wavelength) * np.sin(theta_max)
+    hankel = sesans.SesansTransform(data.x, SElength, wavelength, zaccept, Rmax)
     return hankel
 
 
@@ -204,7 +204,7 @@ def _interpret_data(self, data: Data, model: KernelModel) -> None:
         self._model = model
 
         # interpret data
-        if hasattr(data, 'isSesans') and data.isSesans:
+        if getattr(data, 'isSesans', False):
             self.data_type = 'sesans'
         elif hasattr(data, 'qx_data'):
             self.data_type = 'Iqxy'
@@ -327,7 +327,9 @@ def _calc_theory(self, pars, cutoff=0.0):
 
         # Need to pull background out of resolution for multiple scattering
         default_background = self._model.info.parameters.common_parameters[1].default
-        background = pars.get('background', default_background)
+        background = (
+            pars.get('background', default_background)
+            if self.data_type != 'sesans' else 0.)
         pars = pars.copy()
         pars['background'] = 0.