From 9d70932f0a036297ae0371e58a6ab4edd06a3541 Mon Sep 17 00:00:00 2001 From: Paul Kienzle Date: Wed, 28 Jul 2021 14:37:27 -0400 Subject: [PATCH 1/3] Add -sesans option to sascomp simulator --- sasmodels/compare.py | 22 ++++++++--- sasmodels/data.py | 78 ++++++++++++++++++++++++++++++--------- sasmodels/direct_model.py | 15 ++++---- 3 files changed, 84 insertions(+), 31 deletions(-) diff --git a/sasmodels/compare.py b/sasmodels/compare.py index d6f84d2fd..d7f4a1457 100755 --- a/sasmodels/compare.py +++ b/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 @@ -75,9 +75,9 @@ -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 === @@ -680,7 +680,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'] @@ -690,6 +690,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) @@ -754,6 +761,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 @@ -1009,7 +1019,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=', @@ -1169,6 +1179,7 @@ def parse_opts(argv): 'qmin' : None, 'qmax' : 0.05, 'nq' : 128, + 'sesans' : False, 'res' : '0.0', 'noise' : 0.0, 'accuracy' : 'Low', @@ -1206,6 +1217,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='): diff --git a/sasmodels/data.py b/sasmodels/data.py index a1427f4b4..73b7a1425 100644 --- a/sasmodels/data.py +++ b/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 = (90., 'degrees') + 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,20 +604,20 @@ 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: @@ -584,7 +625,8 @@ def _plot_result_sesans(data, theory, resid, use_data, limits=None): 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 diff --git a/sasmodels/direct_model.py b/sasmodels/direct_model.py index beaa555a5..e3bb3be4d 100644 --- a/sasmodels/direct_model.py +++ b/sasmodels/direct_model.py @@ -166,15 +166,14 @@ def _make_sesans_transform(data): # 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) + wavelength, wunits = data.source.wavelength, data.source.wavelength_unit + zacceptance, zunits = data.sample.zacceptance + wavelength = Converter("A")(wavelength, units=wunits) + theta_max = Converter("radian")(zacceptance, units=zunits) + zaccept = 2 * np.pi / np.max(wavelength) * np.sin(theta_max) Rmax = 10000000 - hankel = sesans.SesansTransform(data.x, SElength, - data.source.wavelength, - zaccept, Rmax) + hankel = sesans.SesansTransform(data.x, SElength, wavelength, zaccept, Rmax) return hankel @@ -207,7 +206,7 @@ def _interpret_data(self, data, model): 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' From f1fb5cecb476b8f6b391aee2891afa5d11e3965b Mon Sep 17 00:00:00 2001 From: Paul Kienzle Date: Thu, 29 Jul 2021 17:35:47 -0400 Subject: [PATCH 2/3] ignore background parameter in sesans calculation --- sasmodels/direct_model.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/sasmodels/direct_model.py b/sasmodels/direct_model.py index e3bb3be4d..cb88fcd4b 100644 --- a/sasmodels/direct_model.py +++ b/sasmodels/direct_model.py @@ -329,7 +329,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. From cb130d7b1c2df8680a828ce110e987f5c5537afc Mon Sep 17 00:00:00 2001 From: Paul Kienzle Date: Thu, 21 Apr 2022 12:57:17 -0400 Subject: [PATCH 3/3] delay loading nxsutil until it is needed --- sasmodels/data.py | 2 +- sasmodels/direct_model.py | 15 ++++++++------- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/sasmodels/data.py b/sasmodels/data.py index 73b7a1425..49fba06b5 100644 --- a/sasmodels/data.py +++ b/sasmodels/data.py @@ -318,7 +318,7 @@ def empty_sesans(z, wavelength=None, zacceptance=None): if np.isscalar(wavelength): wavelength = np.full_like(z, wavelength) if zacceptance is None: - zacceptance = (90., 'degrees') + zacceptance = (np.pi/2, 'radians') source = Source() source.wavelength = wavelength source.wavelength_unit = "A" diff --git a/sasmodels/direct_model.py b/sasmodels/direct_model.py index cb88fcd4b..b56cb2416 100644 --- a/sasmodels/direct_model.py +++ b/sasmodels/direct_model.py @@ -162,17 +162,18 @@ 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") + SElength, SEunits = data.x, data._xunit wavelength, wunits = data.source.wavelength, data.source.wavelength_unit - zacceptance, zunits = data.sample.zacceptance - wavelength = Converter("A")(wavelength, units=wunits) - theta_max = Converter("radian")(zacceptance, units=zunits) - zaccept = 2 * np.pi / np.max(wavelength) * np.sin(theta_max) + 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 + zaccept = 2 * np.pi / np.max(wavelength) * np.sin(theta_max) hankel = sesans.SesansTransform(data.x, SElength, wavelength, zaccept, Rmax) return hankel