Changeset 1198f90 in sasmodels

Timestamp:

May 17, 2018 8:07:02 PM (6 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

df0d2ca

Parents:

33969b6

Message:

limit pinhole resolution integral to ± 2.5 dq

Location:

sasmodels

Files:

• compare.py (modified) (12 diffs)
• jitter.py (modified) (1 diff)
• resolution.py (modified) (6 diffs)

sasmodels/compare.py

@@ -645 +645 @@
 
 def make_data(opts):
-    # type: (Dict[str, Any]) -> Tuple[Data, np.ndarray]
+    # type: (Dict[str, Any], float) -> Tuple[Data, np.ndarray]
     """
     Generate an empty dataset, used with the model to set Q points
@@ -786 +786 @@
     base_n, comp_n = opts['count']
     base_pars, comp_pars = opts['pars']
-    data = opts['data']
+    base_data, comp_data = opts['data']
 
     comparison = comp is not None
@@ -800 +800 @@
             print("%s t=%.2f ms, intensity=%.0f"
                   % (base.engine, base_time, base_value.sum()))
-        _show_invalid(data, base_value)
+        _show_invalid(base_data, base_value)
     except ImportError:
         traceback.print_exc()
@@ -812 +812 @@
                 print("%s t=%.2f ms, intensity=%.0f"
                       % (comp.engine, comp_time, comp_value.sum()))
-            _show_invalid(data, comp_value)
+            _show_invalid(base_data, comp_value)
         except ImportError:
             traceback.print_exc()
@@ -866 +866 @@
     have_base, have_comp = (base_value is not None), (comp_value is not None)
     base, comp = opts['engines']
-    data = opts['data']
+    base_data, comp_data = opts['data']
     use_data = (opts['datafile'] is not None) and (have_base ^ have_comp)
 
@@ -884 +884 @@
         if have_comp:
             plt.subplot(131)
-        plot_theory(data, base_value, view=view, use_data=use_data, limits=limits)
+        plot_theory(base_data, base_value, view=view, use_data=use_data, limits=limits)
         plt.title("%s t=%.2f ms"%(base.engine, base_time))
         #cbar_title = "log I"
@@ -891 +891 @@
             plt.subplot(132)
         if not opts['is2d'] and have_base:
-            plot_theory(data, base_value, view=view, use_data=use_data, limits=limits)
-        plot_theory(data, comp_value, view=view, use_data=use_data, limits=limits)
+            plot_theory(comp_data, base_value, view=view, use_data=use_data, limits=limits)
+        plot_theory(comp_data, comp_value, view=view, use_data=use_data, limits=limits)
         plt.title("%s t=%.2f ms"%(comp.engine, comp_time))
         #cbar_title = "log I"
@@ -908 +908 @@
             err[err > cutoff] = cutoff
         #err,errstr = base/comp,"ratio"
-        plot_theory(data, None, resid=err, view=errview, use_data=use_data)
+        # Note: base_data only since base and comp have same q values (though
+        # perhaps different resolution), and we are plotting the difference
+        # at each q
+        plot_theory(base_data, None, resid=err, view=errview, use_data=use_data)
         plt.xscale('log' if view == 'log' and not opts['is2d'] else 'linear')
         plt.legend(['P%d'%(k+1) for k in range(setnum+1)], loc='best')
@@ -1075 +1078 @@
         'qmax'      : 0.05,
         'nq'        : 128,
-        'res'       : 0.0,
+        'res'       : '0.0',
         'noise'     : 0.0,
         'accuracy'  : 'Low',
@@ -1115 +1118 @@
         elif arg.startswith('-q='):
             opts['qmin'], opts['qmax'] = [float(v) for v in arg[3:].split(':')]
-        elif arg.startswith('-res='):      opts['res'] = float(arg[5:])
+        elif arg.startswith('-res='):      opts['res'] = arg[5:]
         elif arg.startswith('-noise='):    opts['noise'] = float(arg[7:])
         elif arg.startswith('-sets='):     opts['sets'] = int(arg[6:])
@@ -1173 +1176 @@
     if opts['qmin'] is None:
         opts['qmin'] = 0.001*opts['qmax']
-    if opts['datafile'] is not None:
-        data = load_data(os.path.expanduser(opts['datafile']))
-    else:
-        data, _ = make_data(opts)
 
     comparison = any(PAR_SPLIT in v for v in values)
@@ -1216 +1215 @@
         opts['cutoff'] = [float(opts['cutoff'])]*2
 
-    base = make_engine(model_info[0], data, opts['engine'][0],
+    if PAR_SPLIT in opts['res']:
+        opts['res'] = [float(k) for k in opts['res'].split(PAR_SPLIT, 2)]
+        comparison = True
+    else:
+        opts['res'] = [float(opts['res'])]*2
+
+    if opts['datafile'] is not None:
+        data = load_data(os.path.expanduser(opts['datafile']))
+    else:
+        # Hack around the fact that make_data doesn't take a pair of resolutions
+        res = opts['res']
+        opts['res'] = res[0]
+        data0, _ = make_data(opts)
+        if res[0] != res[1]:
+            opts['res'] = res[1]
+            data1, _ = make_data(opts)
+        else:
+            data1 = data0
+        opts['res'] = res
+        data = data0, data1
+
+    base = make_engine(model_info[0], data[0], opts['engine'][0],
                        opts['cutoff'][0], opts['ngauss'][0])
     if comparison:
-        comp = make_engine(model_info[1], data, opts['engine'][1],
+        comp = make_engine(model_info[1], data[1], opts['engine'][1],
                            opts['cutoff'][1], opts['ngauss'][1])
     else:

sasmodels/jitter.py

@@ -774 +774 @@
         # set small jitter as 0 if multiple pd dims
         dims = sum(v > 0 for v in jitter)
-        limit = [0, 0, 0.5, 5][dims]
+        limit = [0, 0.5, 5][dims]
         jitter = [0 if v < limit else v for v in jitter]
         axes.cla()

sasmodels/resolution.py

@@ -20 +20 @@
 MINIMUM_RESOLUTION = 1e-8
 MINIMUM_ABSOLUTE_Q = 0.02  # relative to the minimum q in the data
+PINHOLE_N_SIGMA = 2.5 # From: Barker & Pedersen 1995 JAC
 
 class Resolution(object):
@@ -65 +66 @@
     *q_calc* is the list of points to calculate, or None if this should
     be estimated from the *q* and *q_width*.
-    """
-    def __init__(self, q, q_width, q_calc=None, nsigma=3):
+
+    *nsigma* is the width of the resolution function.  Should be 2.5.
+    See :func:`pinhole_resolution` for details.
+    """
+    def __init__(self, q, q_width, q_calc=None, nsigma=PINHOLE_N_SIGMA):
         #*min_step* is the minimum point spacing to use when computing the
         #underlying model.  It should be on the order of
@@ -88 +92 @@
         # Build weight matrix from calculated q values
         self.weight_matrix = pinhole_resolution(
-            self.q_calc, self.q, np.maximum(q_width, MINIMUM_RESOLUTION))
+            self.q_calc, self.q, np.maximum(q_width, MINIMUM_RESOLUTION),
+            nsigma=nsigma)
         self.q_calc = abs(self.q_calc)
 
@@ -154 +159 @@
 
 
-def pinhole_resolution(q_calc, q, q_width):
-    """
+def pinhole_resolution(q_calc, q, q_width, nsigma=PINHOLE_N_SIGMA):
+    r"""
     Compute the convolution matrix *W* for pinhole resolution 1-D data.
 
@@ -162 +167 @@
     *W*, the resolution smearing can be computed using *dot(W,q)*.
 
+    Note that resolution is limited to $\pm 2.5 \sigma$.[1]  The true resolution
+    function is a broadened triangle, and does not extend over the entire
+    range $(-\infty, +\infty)$.  It is important to impose this limitation
+    since some models fall so steeply that the weighted value in gaussian
+    tails would otherwise dominate the integral.
+
     *q_calc* must be increasing.  *q_width* must be greater than zero.
+
+    [1] Barker, J. G., and J. S. Pedersen. 1995. Instrumental Smearing Effects
+    in Radially Symmetric Small-Angle Neutron Scattering by Numerical and
+    Analytical Methods. Journal of Applied Crystallography 28 (2): 105--14.
+    https://doi.org/10.1107/S0021889894010095.
     """
     # The current algorithm is a midpoint rectangle rule.  In the test case,
@@ -170 +186 @@
     cdf = erf((edges[:, None] - q[None, :]) / (sqrt(2.0)*q_width)[None, :])
     weights = cdf[1:] - cdf[:-1]
+    # Limit q range to +/- 2.5 sigma
+    weights[q_calc[:, None] < (q - nsigma*q_width)[None, :]] = 0.
+    weights[q_calc[:, None] > (q + nsigma*q_width)[None, :]] = 0.
     weights /= np.sum(weights, axis=0)[None, :]
     return weights