From 9a3a6638aa638a64e07e4c7c405010104e372a70 Mon Sep 17 00:00:00 2001
From: Paul Kienzle <pkienzle@nist.gov>
Date: Thu, 30 Mar 2023 17:20:57 -0400
Subject: [PATCH 01/18] Simplified interface for computing I(q,dq)

---
 sasmodels/data.py         |  18 ++--
 sasmodels/direct_model.py | 167 ++++++++++++++++++++++++++++++--------
 sasmodels/resolution.py   |  24 ++++--
 sasmodels/sesans.py       |  24 +++++-
 4 files changed, 180 insertions(+), 53 deletions(-)

diff --git a/sasmodels/data.py b/sasmodels/data.py
index 9766a957e..b3c60c343 100644
--- a/sasmodels/data.py
+++ b/sasmodels/data.py
@@ -134,6 +134,9 @@ class Source:
 class Sample:
     ...
 
+def _as_numpy(data):
+    return None if data is None else np.asarray(data)
+
 class Data1D(object):
     """
     1D data object.
@@ -163,11 +166,12 @@ class Data1D(object):
     """
     def __init__(self, x=None, y=None, dx=None, dy=None):
         # type: (OptArray, OptArray, OptArray, OptArray) -> None
-        self.x, self.y, self.dx, self.dy = x, y, dx, dy
+        self.x, self.dx = _as_numpy(x), _as_numpy(dx)
+        self.y, self.cy = _as_numpy(y), _as_numpy(dy)
         self.dxl = None
         self.filename = None
-        self.qmin = x.min() if x is not None else np.NaN
-        self.qmax = x.max() if x is not None else np.NaN
+        self.qmin = self.x.min() if self.x is not None else np.NaN
+        self.qmax = self.x.max() if self.x is not None else np.NaN
         # TODO: why is 1D mask False and 2D mask True?
         self.mask = (np.isnan(y) if y is not None
                      else np.zeros_like(x, 'b') if x is not None
@@ -240,13 +244,13 @@ class Data2D(object):
     """
     def __init__(self, x=None, y=None, z=None, dx=None, dy=None, dz=None):
         # type: (OptArray, OptArray, OptArray, OptArray, OptArray, OptArray) -> None
-        self.qx_data, self.dqx_data = x, dx
-        self.qy_data, self.dqy_data = y, dy
-        self.data, self.err_data = z, dz
+        self.qx_data, self.dqx_data = _as_numpy(x), _as_numpy(dx)
+        self.qy_data, self.dqy_data = _as_numpy(y), _as_numpy(dy)
+        self.data, self.err_data = _as_numpy(z), _as_numpy(dz)
         self.mask = (np.isnan(z) if z is not None
                      else np.zeros_like(x, dtype='bool') if x is not None
                      else None)
-        self.q_data = np.sqrt(x**2 + y**2)
+        self.q_data = np.sqrt(self.qx_data**2 + self.qy_data**2)
         self.qmin = 1e-16
         self.qmax = np.inf
         self.detector = []
diff --git a/sasmodels/direct_model.py b/sasmodels/direct_model.py
index 5b9b23478..4f3ea1a05 100644
--- a/sasmodels/direct_model.py
+++ b/sasmodels/direct_model.py
@@ -118,12 +118,15 @@ def get_mesh(model_info, values, dim='1d', mono=False):
         active = lambda name: True
 
     #print("in get_mesh: pars:",[p.id for p in parameters.call_parameters])
-    mesh = [_get_par_weights(p, values, active(p.name))
+    values = values.copy()
+    mesh = [_pop_par_weights(p, values, active(p.name))
             for p in parameters.call_parameters]
+    if values:
+        raise TypeError(f"Unused parameters in call: {', '.join(values.keys())}")
     return mesh
 
 
-def _get_par_weights(parameter, values, active=True):
+def _pop_par_weights(parameter, values, active=True):
     # type: (Parameter, Dict[str, float], bool) -> Tuple[float, np.ndarray, np.ndarray]
     """
     Generate the distribution for parameter *name* given the parameter values
@@ -132,21 +135,24 @@ def _get_par_weights(parameter, values, active=True):
     Uses "name", "name_pd", "name_pd_type", "name_pd_n", "name_pd_sigma"
     from the *pars* dictionary for parameter value and parameter dispersion.
     """
-    value = float(values.get(parameter.name, parameter.default))
-    npts = values.get(parameter.name+'_pd_n', 0)
-    width = values.get(parameter.name+'_pd', 0.0)
-    relative = parameter.relative_pd
-    if npts == 0 or width == 0.0 or not active:
-        # Note: orientation parameters have the viewing angle as the parameter
-        # value and the jitter in the distribution, so be sure to set the
-        # empty pd for orientation parameters to 0.
-        pd = [value if relative or not parameter.polydisperse else 0.0], [1.0]
+    value = float(values.pop(parameter.name, parameter.default))
+    if parameter.polydisperse:
+        npts = values.pop(parameter.name+'_pd_n', 0)
+        width = values.pop(parameter.name+'_pd', 0.0)
+        nsigma = values.pop(parameter.name+'_pd_nsigma', 3.0)
+        distribution = values.pop(parameter.name+'_pd_type', 'gaussian')
+        relative = parameter.relative_pd
+        if npts == 0 or width == 0.0 or not active:
+            # Note: orientation parameters have the viewing angle as the parameter
+            # value and the jitter in the distribution, so be sure to set the
+            # empty pd for orientation parameters to 0.
+            pd = [value if relative else 0.0], [1.0]
+        else:
+            limits = parameter.limits
+            pd = weights.get_weights(distribution, npts, width, nsigma,
+                                    value, limits, relative)
     else:
-        limits = parameter.limits
-        disperser = values.get(parameter.name+'_pd_type', 'gaussian')
-        nsigma = values.get(parameter.name+'_pd_nsigma', 3.0)
-        pd = weights.get_weights(disperser, npts, width, nsigma,
-                                 value, limits, relative)
+        pd = [value], [1.0]
     return value, pd[0], pd[1]
 
 
@@ -261,9 +267,10 @@ def _interpret_data(self, data: Data, model: KernelModel) -> None:
                     res = resolution.Perfect1D(q)
             elif (getattr(data, 'dxl', None) is not None
                   and getattr(data, 'dxw', None) is not None):
-                res = resolution.Slit1D(data.x[index],
-                                        qx_width=data.dxl[index],
-                                        qy_width=data.dxw[index])
+                res = resolution.Slit1D(
+                    data.x[index],
+                    q_length=data.dxl[index],
+                    q_width=data.dxw[index])
             else:
                 res = resolution.Perfect1D(data.x[index])
         elif self.data_type == 'Iq-oriented':
@@ -456,6 +463,65 @@ def test_reparameterize():
     except Exception:
         pass
 
+def _direct_calculate(model, data, pars):
+    from .core import load_model_info, build_model
+    model_info = load_model_info(model)
+    kernel = build_model(model_info)
+    calculator = DirectModel(data, kernel)
+    return calculator(**pars)
+
+def Iq(model, q, dq=None, ql=None, qw=None, **pars):
+    """
+    Compute I(q) for *model*. Resolution is *dq* for pinhole or *ql* and *qw*
+    for slit geometry.
+
+    Model is the name of a builtin or custom model, or a model expression, such
+    as sphere+sphere for a mixture of spheres of different radii, or
+    sphere@hardsphere for concentrated solutions where the dilute approximation
+    no longer applies.
+
+    Use additional keywords for model parameters, tagged with *_pd*, *_pd_n*,
+    *_pd_nsigma*, *_pd_type* to set polydispersity parameters, or *_M0*,
+    *_mphi*, *_mtheta* for magnetic parameters.
+
+    This is not intended for use when the same I(q) is evaluated many times
+    with different parameter values. For that you should set up the model
+    with `model = build_model(load_model_info(model_name))`, set up a data
+    object to define q values and resolution, then use
+    `calculator = DirectModel(data, model)` to set up a calculator, or
+    `problem = bumps.FitProblem(sasmodels.bumps_model.Experiment(data, model))`
+    to define a fit problem for uses with the bumps optimizer. Data can be
+    loaded using the `sasdata` package, or use one of the empty data generators
+    from `sasmodels.data`.
+
+    Models are cached. Custom models will not be reloaded even if the
+    underlying files have changed. If you are using this in a long running
+    application then you will need to call
+    `sasmodels.direct_model._model_cache.clear()` to reset the cache and force
+    custom model reload.
+    """
+    from .data import Data1D, _as_numpy
+    data = Data1D(x=q, dx=dq)
+    data.dxl, data.dxw = _as_numpy(ql), _as_numpy(qw)
+    return _direct_calculate(model, data, pars)
+
+def Iqxy(model, qx, qy, dqx=None, dqy=None, **pars):
+    """
+    Compute I(qx, qy) for *model*. Resolution is *dqx* and *dqy*.
+    See :func:`Iq` for details on model and parameters.
+    """
+    from .data import Data2D
+    data = Data2D(x=qx, y=qy, dx=dqx, dy=dqy)
+    return _direct_calculate(model, data, pars)
+
+def Gxi(model, xi, **pars):
+    """
+    Compute SESANS correlation G' = G(xi) - G(0) for *model*.
+    See :func:`Iq` for details on model and parameters.
+    """
+    from .data import empty_sesans
+    data = empty_sesans(z=xi)
+    return _direct_calculate(model, data, pars)
 
 def main():
     # type: () -> None
@@ -463,36 +529,71 @@ def main():
     Program to evaluate a particular model at a set of q values.
     """
     import sys
-    from .data import empty_data1D, empty_data2D
-    from .core import load_model_info, build_model
 
     if len(sys.argv) < 3:
         print("usage: python -m sasmodels.direct_model modelname (q|qx,qy) par=val ...")
         sys.exit(1)
-    model_name = sys.argv[1]
+    model = sys.argv[1]
     call = sys.argv[2].upper()
+    pars = dict((k, (float(v) if not k.endswith("_pd_type") else v))
+                for pair in sys.argv[3:]
+                for k, v in [pair.split('=')])
     try:
         values = [float(v) for v in call.split(',')]
     except ValueError:
         values = []
     if len(values) == 1:
         q, = values
-        data = empty_data1D([q])
+        dq = dqw = dql = None
+        #dq = [q*0.05] # 5% pinhole resolution
+        #dqw, dql = [q*0.05], [1.0] # 5% horizontal slit resolution
+        print(Iq(model, [q], dq=dq, qw=dqw, ql=dql, **pars)[0])
+        #print(Gxi(model, [q], **pars)[0])
     elif len(values) == 2:
         qx, qy = values
-        data = empty_data2D([qx], [qy])
+        dq = None
+        #dq = [0.005] # 5% pinhole resolution at q = 0.1
+        print(Iqxy(model, [qx], [qy], dqx=dq, dqy=dq, **pars)[0])
     else:
         print("use q or qx,qy")
         sys.exit(1)
 
-    model_info = load_model_info(model_name)
-    model = build_model(model_info)
-    calculator = DirectModel(data, model)
-    pars = dict((k, (float(v) if not k.endswith("_pd_type") else v))
-                for pair in sys.argv[3:]
-                for k, v in [pair.split('=')])
-    Iq = calculator(**pars)
-    print(Iq[0])
+def test_simple_interface():
+    def near(value, target):
+        """Close enough in single precision"""
+        #print(f"value: {value}, target: {target}")
+        return np.allclose(value, target, rtol=1e-6, atol=0, equal_nan=True)
+    # Note: target values taken from running main() on parameters.
+    # Resolution was 5% dq/q.
+    pars = dict(radius=200)
+    # simple sphere in 1D (perfect, pinhole, slit)
+    assert near(Iq('sphere', [0.1], **pars), [0.6200146273894904])
+    assert near(Iq('sphere', [0.1], dq=[0.005], **pars), [2.3019224683980215])
+    assert near(Iq('sphere', [0.1], qw=[0.005], ql=[1.0], **pars), [0.3673431784535172])
+    # simple sphere in 2D (perfect, pinhole)
+    assert near(Iqxy('sphere', [0.1], [0.1], **pars), [1.1781532874802199])
+    assert near(Iqxy('sphere', [0.1], [0.1], dqx=[0.005], dqy=[0.005], **pars), 
+        [0.8177780778578667])
+    # sesans
+    assert near(Gxi('sphere', [100], **pars), [-0.19146959126623486])
+    # Check that single point sesans matches value in an array
+    xi = np.logspace(1, 3, 100)
+    y = Gxi('sphere', xi, **pars)
+    for k in (0, len(xi)//5, len(xi)//2, len(xi)-1):
+        ysingle = Gxi('sphere', [xi[k]], **pars)[0]
+        print(f"SESANS point check {k}: xi={xi[k]:.1f} single={ysingle:.4f} vector={y[k]:.4f}")
+        assert abs((ysingle-y[k])/y[k]) < 0.1, "SESANS point value not matching vector value within 10%"
+    # magnetic 2D
+    pars = dict(radius=200, sld_M0=3, sld_mtheta=30)
+    assert near(Iqxy('sphere', [0.1], [0.1], **pars), [1.5577852226925908])
+    # polydisperse 1D
+    pars = dict(
+        radius=200, radius_pd=0.1, radius_pd_n=15, radius_pd_nsigma=2.5,
+        radius_pd_type="uniform")
+    assert near(Iq('sphere', [0.1], **pars), [2.703169824954617])
 
 if __name__ == "__main__":
-    main()
+    import logging
+    logging.disable(logging.ERROR)
+    #main()
+    test_simple_interface()
diff --git a/sasmodels/resolution.py b/sasmodels/resolution.py
index 27d7b255d..e4ea7e6e1 100644
--- a/sasmodels/resolution.py
+++ b/sasmodels/resolution.py
@@ -451,12 +451,14 @@ def linear_extrapolation(q, q_min, q_max):
     """
     q = np.sort(q)
     if q_min + 2*MINIMUM_RESOLUTION < q[0]:
-        n_low = int(np.ceil((q[0]-q_min) / (q[1]-q[0]))) if q[1] > q[0] else 15
+        delta = q[1] - q[0] if len(q) > 1 else 0
+        n_low = int(np.ceil((q[0]-q_min) / delta)) if delta > 0 else 15
         q_low = np.linspace(q_min, q[0], n_low+1)[:-1]
     else:
         q_low = []
     if q_max - 2*MINIMUM_RESOLUTION > q[-1]:
-        n_high = int(np.ceil((q_max-q[-1]) / (q[-1]-q[-2]))) if q[-1] > q[-2] else 15
+        delta = q[-1] - q[-2] if len(q) > 1 else 0
+        n_high = int(np.ceil((q_max-q[-1]) / delta)) if delta > 0 else 15
         q_high = np.linspace(q[-1], q_max, n_high+1)[1:]
     else:
         q_high = []
@@ -496,21 +498,25 @@ def geometric_extrapolation(q, q_min, q_max, points_per_decade=None):
          n_\text{extend} = (n-1) (\log q_\text{max} - \log q_n)
             / (\log q_n - \log q_1)
     """
-    q = np.sort(q)
+    DEFAULT_POINTS_PER_DECADE = 10
+    data_min, data_max = q.min(), q.max()
     if points_per_decade is None:
-        log_delta_q = (len(q) - 1) / (log(q[-1]) - log(q[0]))
+        if data_max > data_min:
+            log_delta_q = (len(q) - 1) / (log(data_max) - log(data_min))
+        else:
+            log_delta_q = log(10.) / DEFAULT_POINTS_PER_DECADE
     else:
         log_delta_q = log(10.) / points_per_decade
-    if q_min < q[0]:
+    if q_min < data_min:
         if q_min < 0:
-            q_min = q[0]*MINIMUM_ABSOLUTE_Q
+            q_min = data_min*MINIMUM_ABSOLUTE_Q
         n_low = int(np.ceil(log_delta_q * (log(q[0])-log(q_min))))
         q_low = np.logspace(log10(q_min), log10(q[0]), n_low+1)[:-1]
     else:
         q_low = []
-    if q_max > q[-1]:
-        n_high = int(np.ceil(log_delta_q * (log(q_max)-log(q[-1]))))
-        q_high = np.logspace(log10(q[-1]), log10(q_max), n_high+1)[1:]
+    if q_max > data_max:
+        n_high = int(np.ceil(log_delta_q * (log(q_max)-log(data_max))))
+        q_high = np.logspace(log10(data_max), log10(q_max), n_high+1)[1:]
     else:
         q_high = []
     return np.concatenate([q_low, q, q_high])
diff --git a/sasmodels/sesans.py b/sasmodels/sesans.py
index 8eb65d2d2..b4bc6caba 100644
--- a/sasmodels/sesans.py
+++ b/sasmodels/sesans.py
@@ -60,10 +60,18 @@ def apply(self, Iq):
     def _set_hankel(self, SElength, lam, zaccept, Rmax):
         # type: (np.ndarray, float, float, float) -> None
         SElength = np.asarray(SElength)
-        q_max = 2*pi / (SElength[1] - SElength[0])
-        q_min = 0.1 * 2*pi / (np.size(SElength) * SElength[-1])
+        if len(SElength) == 1:
+            # TODO: Do we care that this fails for xi = 0?
+            q_min, q_max = 0.01 * 2*pi/SElength[-1], 10*2*pi / SElength[0]
+        else:
+            # TODO: Why does q_min depend on the number of correlation lengths?
+            # TODO: Why does q_max depend on the correlation step size?
+            q_min = 0.1 * 2*pi / (np.size(SElength) * SElength[-1])
+            q_max = 2*pi / (SElength[1] - SElength[0])
+        #print("Hankel xi, Qmin, Qmax", SElength[0], q_min, q_max, len(SElength))
         q = np.exp(np.arange(np.log(q_min), np.log(q_max),
                              np.log(self.log_spacing)))
+        #print(q)
 
         dq = np.diff(q)
         dq = np.insert(dq, 0, dq[0])
@@ -75,9 +83,17 @@ def _set_hankel(self, SElength, lam, zaccept, Rmax):
         H *= (dq * q / (2*pi)).reshape((-1, 1))
 
         reptheta = np.outer(q, lam/(2*pi))
-        np.arcsin(reptheta, out=reptheta)
-        mask = reptheta > zaccept
+        # Note: Using inplace update with reptheta => arcsin(reptheta).
+        # When q L / 2 pi > 1 that means wavelength is too large to
+        # reach that q value at any angle. These should produce theta = NaN
+        # without any warnings.
+        with np.errstate(invalid='ignore'):
+            np.arcsin(reptheta, out=reptheta)
+        # Reverse the condition to protect against NaN. We can't use
+        # theta > zaccept since all comparisons with NaN return False.
+        mask = ~(reptheta <= zaccept)
         H[mask] = 0
+        #print("number of masked points", np.sum(mask))
 
         self.q_calc = q
         self._H, self._H0 = H, H0

From 47c80b571cc7a11b5460b99f67a15f1602ba42fb Mon Sep 17 00:00:00 2001
From: Paul Kienzle <pkienzle@nist.gov>
Date: Thu, 30 Mar 2023 17:39:34 -0400
Subject: [PATCH 02/18] Restore q sorting for geometric extrapolation.

---
 sasmodels/resolution.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/sasmodels/resolution.py b/sasmodels/resolution.py
index e4ea7e6e1..757e8d79c 100644
--- a/sasmodels/resolution.py
+++ b/sasmodels/resolution.py
@@ -499,7 +499,8 @@ def geometric_extrapolation(q, q_min, q_max, points_per_decade=None):
             / (\log q_n - \log q_1)
     """
     DEFAULT_POINTS_PER_DECADE = 10
-    data_min, data_max = q.min(), q.max()
+    q = np.sort(q)
+    data_min, data_max = q[0], q[-1]
     if points_per_decade is None:
         if data_max > data_min:
             log_delta_q = (len(q) - 1) / (log(data_max) - log(data_min))

From f1761b6bb17982d74697a1e3e39ee5ba039122bd Mon Sep 17 00:00:00 2001
From: Paul Kienzle <pkienzle@nist.gov>
Date: Thu, 11 May 2023 15:13:21 -0400
Subject: [PATCH 03/18] Accept None or zero for unlimited slit resolution

---
 sasmodels/data.py         |  2 +-
 sasmodels/direct_model.py | 16 +++-------------
 sasmodels/resolution.py   | 17 ++++++++++++-----
 3 files changed, 16 insertions(+), 19 deletions(-)

diff --git a/sasmodels/data.py b/sasmodels/data.py
index b3c60c343..c1d1105d9 100644
--- a/sasmodels/data.py
+++ b/sasmodels/data.py
@@ -167,7 +167,7 @@ class Data1D(object):
     def __init__(self, x=None, y=None, dx=None, dy=None):
         # type: (OptArray, OptArray, OptArray, OptArray) -> None
         self.x, self.dx = _as_numpy(x), _as_numpy(dx)
-        self.y, self.cy = _as_numpy(y), _as_numpy(dy)
+        self.y, self.dy = _as_numpy(y), _as_numpy(dy)
         self.dxl = None
         self.filename = None
         self.qmin = self.x.min() if self.x is not None else np.NaN
diff --git a/sasmodels/direct_model.py b/sasmodels/direct_model.py
index 4f3ea1a05..5804650ac 100644
--- a/sasmodels/direct_model.py
+++ b/sasmodels/direct_model.py
@@ -156,16 +156,6 @@ def _pop_par_weights(parameter, values, active=True):
     return value, pd[0], pd[1]
 
 
-def _vol_pars(model_info, values):
-    # type: (ModelInfo, ParameterSet) -> Tuple[np.ndarray, np.ndarray]
-    vol_pars = [_get_par_weights(p, values)
-                for p in model_info.parameters.call_parameters
-                if p.type == 'volume']
-    #import pylab; pylab.plot(vol_pars[0][0],vol_pars[0][1]); pylab.show()
-    dispersity, weight = dispersion_mesh(model_info, vol_pars)
-    return dispersity, weight
-
-
 def _make_sesans_transform(data):
     # Pre-compute the Hankel matrix (H)
     SElength, SEunits = data.x, data._xunit
@@ -473,7 +463,7 @@ def _direct_calculate(model, data, pars):
 def Iq(model, q, dq=None, ql=None, qw=None, **pars):
     """
     Compute I(q) for *model*. Resolution is *dq* for pinhole or *ql* and *qw*
-    for slit geometry.
+    for slit geometry. Use 0 or None for infinite slits.
 
     Model is the name of a builtin or custom model, or a model expression, such
     as sphere+sphere for a mixture of spheres of different radii, or
@@ -595,5 +585,5 @@ def near(value, target):
 if __name__ == "__main__":
     import logging
     logging.disable(logging.ERROR)
-    #main()
-    test_simple_interface()
+    main()
+    #test_simple_interface()
diff --git a/sasmodels/resolution.py b/sasmodels/resolution.py
index 757e8d79c..df001bfe6 100644
--- a/sasmodels/resolution.py
+++ b/sasmodels/resolution.py
@@ -124,25 +124,31 @@ class Slit1D(Resolution):
 
     """
 
-    def __init__(self, q, q_length, q_width=0., q_calc=None):
+    def __init__(self, q, q_length=None, q_width=None, q_calc=None):
         # Remember what width/dqy was used even though we won't need them
         # after the weight matrix is constructed
         self.q_length, self.q_width = q_length, q_width
 
         # Allow independent resolution on each point even though it is not
         # needed in practice.
-        if np.isscalar(q_width):
+        if q_width is None:
+            q_width = np.zeros(len(q))
+        elif np.isscalar(q_width):
             q_width = np.ones(len(q))*q_width
         else:
             q_width = np.asarray(q_width)
-        if np.isscalar(q_length):
+        if q_length is None:
+            q_length = np.zeros(len(q))
+        elif np.isscalar(q_length):
             q_length = np.ones(len(q))*q_length
         else:
             q_length = np.asarray(q_length)
 
         self.q = q.flatten()
-        self.q_calc = slit_extend_q(q, q_width, q_length) \
+        self.q_calc = (
+            slit_extend_q(q, q_width, q_length)
             if q_calc is None else np.sort(q_calc)
+        )
 
         # Protect against models which are not defined for very low q.  Limit
         # the smallest q value evaluated (in absolute) to 0.02*min
@@ -150,8 +156,9 @@ def __init__(self, q, q_length, q_width=0., q_calc=None):
         self.q_calc = self.q_calc[abs(self.q_calc) >= cutoff]
 
         # Build weight matrix from calculated q values
-        self.weight_matrix = \
+        self.weight_matrix = (
             slit_resolution(self.q_calc, self.q, q_length, q_width)
+        )
         self.q_calc = abs(self.q_calc)
 
     def apply(self, theory):

From 67770f9bf5ca81b01fd1e699320eae1184be676f Mon Sep 17 00:00:00 2001
From: Paul Kienzle <pkienzle@nist.gov>
Date: Wed, 7 Jun 2023 14:36:31 -0400
Subject: [PATCH 04/18] Allow None or scalar for dxl and dxw in simplified
 interface.

---
 sasmodels/direct_model.py | 20 +++++++++++++-------
 1 file changed, 13 insertions(+), 7 deletions(-)

diff --git a/sasmodels/direct_model.py b/sasmodels/direct_model.py
index 5804650ac..e172de6b9 100644
--- a/sasmodels/direct_model.py
+++ b/sasmodels/direct_model.py
@@ -256,11 +256,11 @@ def _interpret_data(self, data: Data, model: KernelModel) -> None:
                 else:
                     res = resolution.Perfect1D(q)
             elif (getattr(data, 'dxl', None) is not None
-                  and getattr(data, 'dxw', None) is not None):
+                  or getattr(data, 'dxw', None) is not None):
                 res = resolution.Slit1D(
                     data.x[index],
-                    q_length=data.dxl[index],
-                    q_width=data.dxw[index])
+                    q_length=None if data.dxl is None else data.dxl[index],
+                    q_width=None if data.dxw is None else data.dxw[index])
             else:
                 res = resolution.Perfect1D(data.x[index])
         elif self.data_type == 'Iq-oriented':
@@ -275,9 +275,10 @@ def _interpret_data(self, data: Data, model: KernelModel) -> None:
                     or getattr(data, 'dxw', None) is None):
                 raise ValueError("oriented sample with 1D data needs slit resolution")
 
-            res = resolution2d.Slit2D(data.x[index],
-                                      qx_width=data.dxw[index],
-                                      qy_width=data.dxl[index])
+            res = resolution2d.Slit2D(
+                data.x[index],
+                qx_width=data.dxw[index],
+                qy_width=data.dxl[index])
         else:
             raise ValueError("Unknown data type") # never gets here
 
@@ -492,7 +493,12 @@ def Iq(model, q, dq=None, ql=None, qw=None, **pars):
     """
     from .data import Data1D, _as_numpy
     data = Data1D(x=q, dx=dq)
-    data.dxl, data.dxw = _as_numpy(ql), _as_numpy(qw)
+    def broadcast(v):
+        return (
+            None if v is None
+            else np.full(len(q), v) if np.isscalar(v)
+            else _as_numpy(v))
+    data.dxl, data.dxw = broadcast(ql), broadcast(qw)
     return _direct_calculate(model, data, pars)
 
 def Iqxy(model, qx, qy, dqx=None, dqy=None, **pars):

From 05f5ca2085bf82e526f2d16238381551e8a8820d Mon Sep 17 00:00:00 2001
From: tsole0 <locke.apogee@gmail.com>
Date: Thu, 3 Aug 2023 16:39:42 -0400
Subject: [PATCH 05/18] fixes docs so math will display correctly

---
 sasmodels/models/lamellar_hg.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/sasmodels/models/lamellar_hg.py b/sasmodels/models/lamellar_hg.py
index ebdd183ba..887673135 100644
--- a/sasmodels/models/lamellar_hg.py
+++ b/sasmodels/models/lamellar_hg.py
@@ -29,7 +29,7 @@
 and $\Delta\rho_T$ is tail contrast (*sld* $-$ *sld_solvent*).
 
 The total thickness of the lamellar sheet is
-a_H + \delta_T + \delta_T + \delta_H$. Note that in a non aqueous solvent
+$a_H + \delta_T + \delta_T + \delta_H$. Note that in a non aqueous solvent
 the chemical "head" group may be the "Tail region" and vice-versa.
 
 The 2D scattering intensity is calculated in the same way as 1D, where

From b7cfdbaf437d50e8e3c4504c125b9ad91d1dae84 Mon Sep 17 00:00:00 2001
From: lucas-wilkins <lucas.wilkins@stfc.ac.uk>
Date: Fri, 4 Aug 2023 21:02:10 +0100
Subject: [PATCH 06/18] Disable name check

---
 sasmodels/direct_model.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/sasmodels/direct_model.py b/sasmodels/direct_model.py
index e172de6b9..b106878bf 100644
--- a/sasmodels/direct_model.py
+++ b/sasmodels/direct_model.py
@@ -121,8 +121,11 @@ def get_mesh(model_info, values, dim='1d', mono=False):
     values = values.copy()
     mesh = [_pop_par_weights(p, values, active(p.name))
             for p in parameters.call_parameters]
-    if values:
-        raise TypeError(f"Unused parameters in call: {', '.join(values.keys())}")
+
+    # TODO: This check has been removed because it currently prevents DirectModel from functioning correctly
+    # if values:
+    #     raise TypeError(f"Unused parameters in call: {', '.join(values.keys())}")
+
     return mesh
 
 

From 40076125e539c65e17447b98e2976abeffff4cba Mon Sep 17 00:00:00 2001
From: Paul Kienzle <pkienzle@nist.gov>
Date: Tue, 8 Aug 2023 11:18:23 -0400
Subject: [PATCH 07/18] Re-enable parameter name checking in kernel call after
 fixing jitter to use correct parameter names

---
 sasmodels/direct_model.py | 25 +++++++++++++---------
 sasmodels/jitter.py       | 44 ++++++++++++++++++++++-----------------
 2 files changed, 40 insertions(+), 29 deletions(-)

diff --git a/sasmodels/direct_model.py b/sasmodels/direct_model.py
index b106878bf..0bff6600f 100644
--- a/sasmodels/direct_model.py
+++ b/sasmodels/direct_model.py
@@ -122,9 +122,8 @@ def get_mesh(model_info, values, dim='1d', mono=False):
     mesh = [_pop_par_weights(p, values, active(p.name))
             for p in parameters.call_parameters]
 
-    # TODO: This check has been removed because it currently prevents DirectModel from functioning correctly
-    # if values:
-    #     raise TypeError(f"Unused parameters in call: {', '.join(values.keys())}")
+    if values:
+        raise TypeError(f"Unused parameters in call: {', '.join(values.keys())}")
 
     return mesh
 
@@ -564,16 +563,17 @@ def near(value, target):
         return np.allclose(value, target, rtol=1e-6, atol=0, equal_nan=True)
     # Note: target values taken from running main() on parameters.
     # Resolution was 5% dq/q.
-    pars = dict(radius=200)
+    pars = dict(radius=200, background=0)  # default background=1e-3, scale=1
     # simple sphere in 1D (perfect, pinhole, slit)
-    assert near(Iq('sphere', [0.1], **pars), [0.6200146273894904])
-    assert near(Iq('sphere', [0.1], dq=[0.005], **pars), [2.3019224683980215])
-    assert near(Iq('sphere', [0.1], qw=[0.005], ql=[1.0], **pars), [0.3673431784535172])
+    perfect_target = 0.6190146273894904
+    assert near(Iq('sphere', [0.1], **pars), [perfect_target])
+    assert near(Iq('sphere', [0.1], dq=[0.005], **pars), [2.3009224683980215])
+    assert near(Iq('sphere', [0.1], qw=[0.005], ql=[1.0], **pars), [0.3663431784535172])
     # simple sphere in 2D (perfect, pinhole)
-    assert near(Iqxy('sphere', [0.1], [0.1], **pars), [1.1781532874802199])
+    assert near(Iqxy('sphere', [0.1], [0.1], **pars), [1.1771532874802199])
     assert near(Iqxy('sphere', [0.1], [0.1], dqx=[0.005], dqy=[0.005], **pars), 
-        [0.8177780778578667])
-    # sesans
+        [0.8167780778578667])
+    # sesans (no background or scale)
     assert near(Gxi('sphere', [100], **pars), [-0.19146959126623486])
     # Check that single point sesans matches value in an array
     xi = np.logspace(1, 3, 100)
@@ -590,6 +590,11 @@ def near(value, target):
         radius=200, radius_pd=0.1, radius_pd_n=15, radius_pd_nsigma=2.5,
         radius_pd_type="uniform")
     assert near(Iq('sphere', [0.1], **pars), [2.703169824954617])
+    # background and scale
+    background, scale = 1e-4, 0.1
+    pars = dict(radius=200, background=background, scale=scale)
+    assert near(Iq('sphere', [0.1], **pars), [perfect_target*scale + background])
+
 
 if __name__ == "__main__":
     import logging
diff --git a/sasmodels/jitter.py b/sasmodels/jitter.py
index cdaf64ba1..b8429c0d4 100755
--- a/sasmodels/jitter.py
+++ b/sasmodels/jitter.py
@@ -881,7 +881,7 @@ def build_model(model_name, n=150, qmax=0.5, **pars):
 
     # stuff the values for non-orientation parameters into the calculator
     calculator.pars = pars.copy()
-    calculator.pars.setdefault('backgound', 1e-3)
+    calculator.pars.setdefault('background', 1e-3)
 
     # fix the data limits so that we can see if the pattern fades
     # under rotation or angular dispersion
@@ -908,47 +908,53 @@ def select_calculator(model_name, n=150, size=(10, 40, 100)):
     a, b, c = size
     d_factor = 0.06  # for paracrystal models
     if model_name == 'sphere':
-        calculator = build_model('sphere', n=n, radius=c)
+        calculator = build_model(
+            'sphere', n=n, radius=c)
         a = b = c
     elif model_name == 'sc_paracrystal':
         a = b = c
         dnn = c
         radius = 0.5*c
-        calculator = build_model('sc_paracrystal', n=n, dnn=dnn,
-                                 d_factor=d_factor, radius=(1-d_factor)*radius,
-                                 background=0)
+        calculator = build_model(
+            'sc_paracrystal', n=n,
+            dnn=dnn, d_factor=d_factor, radius=(1-d_factor)*radius,
+            background=0)
     elif model_name == 'fcc_paracrystal':
         a = b = c
         # nearest neigbour distance dnn should be 2 radius, but I think the
         # model uses lattice spacing rather than dnn in its calculations
         dnn = 0.5*c
         radius = sqrt(2)/4 * c
-        calculator = build_model('fcc_paracrystal', n=n, dnn=dnn,
-                                 d_factor=d_factor, radius=(1-d_factor)*radius,
-                                 background=0)
+        calculator = build_model(
+            'fcc_paracrystal', n=n,
+            dnn=dnn, d_factor=d_factor, radius=(1-d_factor)*radius,
+            background=0)
     elif model_name == 'bcc_paracrystal':
         a = b = c
         # nearest neigbour distance dnn should be 2 radius, but I think the
         # model uses lattice spacing rather than dnn in its calculations
         dnn = 0.5*c
         radius = sqrt(3)/2 * c
-        calculator = build_model('bcc_paracrystal', n=n, dnn=dnn,
-                                 d_factor=d_factor, radius=(1-d_factor)*radius,
-                                 background=0)
+        calculator = build_model(
+            'bcc_paracrystal', n=n,
+            dnn=dnn, d_factor=d_factor, radius=(1-d_factor)*radius,
+            background=0)
     elif model_name == 'cylinder':
-        calculator = build_model('cylinder', n=n, qmax=0.3, radius=b, length=c)
+        calculator = build_model(
+            'cylinder', n=n, qmax=0.3, radius=b, length=c)
         a = b
     elif model_name == 'ellipsoid':
-        calculator = build_model('ellipsoid', n=n, qmax=1.0,
-                                 radius_polar=c, radius_equatorial=b)
+        calculator = build_model(
+            'ellipsoid', n=n, qmax=1.0,
+            radius_polar=c, radius_equatorial=b)
         a = b
     elif model_name == 'triaxial_ellipsoid':
-        calculator = build_model('triaxial_ellipsoid', n=n, qmax=0.5,
-                                 radius_equat_minor=a,
-                                 radius_equat_major=b,
-                                 radius_polar=c)
+        calculator = build_model(
+            'triaxial_ellipsoid', n=n, qmax=0.5,
+            radius_equat_minor=a, radius_equat_major=b, radius_polar=c)
     elif model_name == 'parallelepiped':
-        calculator = build_model('parallelepiped', n=n, a=a, b=b, c=c)
+        calculator = build_model(
+            'parallelepiped', n=n, length_a=a, length_b=b, length_c=c)
     else:
         raise ValueError("unknown model %s"%model_name)
 

From 2643975430834da931ca6aecb8298b80e3171ad7 Mon Sep 17 00:00:00 2001
From: smk78 <smk78@NDLT1227.isis.cclrc.ac.uk>
Date: Sat, 16 Sep 2023 17:06:35 +0100
Subject: [PATCH 08/18] Fix doc build warning (not in toctree)

---
 doc/guide/theory.rst | 1 +
 1 file changed, 1 insertion(+)

diff --git a/doc/guide/theory.rst b/doc/guide/theory.rst
index 0bd3f371c..507867231 100644
--- a/doc/guide/theory.rst
+++ b/doc/guide/theory.rst
@@ -1,3 +1,4 @@
+.. currentmodule:: sasmodels
 .. theory.rst
 
 .. Much of the following text was scraped from fitting_sq.py

From e240c781aec93aa861ecb84441150c297f0159bb Mon Sep 17 00:00:00 2001
From: Wojciech Potrzebowski <wojciech.potrzebowski@ess.eu>
Date: Sat, 16 Sep 2023 17:46:22 +0100
Subject: [PATCH 09/18] Storing generated docs

---
 .github/workflows/test.yml | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 6e7a07ef7..f8fe1f15a 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -59,3 +59,12 @@ jobs:
       if: ${{ matrix.os == 'ubuntu-latest' }}
       run: |
         make -j 4 -C doc SPHINXOPTS="-W --keep-going -n" html
+
+    - name: Publish samodels docs
+      if: ${{ matrix.installer }}
+      uses: actions/upload-artifact@v3
+      with:
+        name: sasmodels-docs-${{ matrix.os }}-${{ matrix.python-version }}
+        path: |
+          docs/html
+        if-no-files-found: error

From bb5c83325e54a1bd8d002cb25fe2e0ea5ad4bf77 Mon Sep 17 00:00:00 2001
From: Wojciech Potrzebowski <wojciech.potrzebowski@ess.eu>
Date: Sat, 16 Sep 2023 18:06:07 +0100
Subject: [PATCH 10/18] Storing generated sasmodels docs

---
 .github/workflows/test.yml | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index f8fe1f15a..86ab19961 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -59,12 +59,13 @@ jobs:
       if: ${{ matrix.os == 'ubuntu-latest' }}
       run: |
         make -j 4 -C doc SPHINXOPTS="-W --keep-going -n" html
+        tar zcf sasmodels-docs.tar.gz html
 
     - name: Publish samodels docs
-      if: ${{ matrix.installer }}
+      if: ${{ matrix.os == 'ubuntu-latest' }}
       uses: actions/upload-artifact@v3
       with:
         name: sasmodels-docs-${{ matrix.os }}-${{ matrix.python-version }}
         path: |
-          docs/html
+          sasmodels_docs.tar.gz
         if-no-files-found: error

From 39ac6f2536f9aecd6b9282f6e0a999263dbef279 Mon Sep 17 00:00:00 2001
From: Wojciech Potrzebowski <wojciech.potrzebowski@ess.eu>
Date: Sat, 16 Sep 2023 18:26:38 +0100
Subject: [PATCH 11/18] Fixing path

---
 .github/workflows/test.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 86ab19961..1be93dec7 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -59,7 +59,7 @@ jobs:
       if: ${{ matrix.os == 'ubuntu-latest' }}
       run: |
         make -j 4 -C doc SPHINXOPTS="-W --keep-going -n" html
-        tar zcf sasmodels-docs.tar.gz html
+        tar zcf sasmodels-docs.tar.gz _build/html
 
     - name: Publish samodels docs
       if: ${{ matrix.os == 'ubuntu-latest' }}

From 1d29862b85ce7596601e91999f3ad453dc8a22ab Mon Sep 17 00:00:00 2001
From: Wojciech Potrzebowski <wojciech.potrzebowski@ess.eu>
Date: Sat, 16 Sep 2023 21:12:47 +0100
Subject: [PATCH 12/18] fixing path again

---
 .github/workflows/test.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 1be93dec7..57d8f45fe 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -59,7 +59,7 @@ jobs:
       if: ${{ matrix.os == 'ubuntu-latest' }}
       run: |
         make -j 4 -C doc SPHINXOPTS="-W --keep-going -n" html
-        tar zcf sasmodels-docs.tar.gz _build/html
+        tar zcf sasmodels-docs.tar.gz doc/_build/html
 
     - name: Publish samodels docs
       if: ${{ matrix.os == 'ubuntu-latest' }}

From 1b32f2b2e2a60e5c3aa73030489642b004ab8d83 Mon Sep 17 00:00:00 2001
From: Wojciech Potrzebowski <wojciech.potrzebowski@ess.eu>
Date: Sat, 16 Sep 2023 21:41:41 +0100
Subject: [PATCH 13/18] fixing path for upload

---
 .github/workflows/test.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 57d8f45fe..58788e2c3 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -67,5 +67,5 @@ jobs:
       with:
         name: sasmodels-docs-${{ matrix.os }}-${{ matrix.python-version }}
         path: |
-          sasmodels_docs.tar.gz
+          sasmodels/sasmodels_docs.tar.gz
         if-no-files-found: error

From a463e89ec8fbcb3e6e1aa81b4e3cf3ba4bb54b0c Mon Sep 17 00:00:00 2001
From: Wojciech Potrzebowski <wojciech.potrzebowski@ess.eu>
Date: Sun, 17 Sep 2023 08:25:23 +0100
Subject: [PATCH 14/18] fixing path for upload

---
 .github/workflows/test.yml | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 58788e2c3..aa3e2f435 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -60,6 +60,7 @@ jobs:
       run: |
         make -j 4 -C doc SPHINXOPTS="-W --keep-going -n" html
         tar zcf sasmodels-docs.tar.gz doc/_build/html
+        ls -l sasmodels-docs.tar.gz
 
     - name: Publish samodels docs
       if: ${{ matrix.os == 'ubuntu-latest' }}
@@ -67,5 +68,5 @@ jobs:
       with:
         name: sasmodels-docs-${{ matrix.os }}-${{ matrix.python-version }}
         path: |
-          sasmodels/sasmodels_docs.tar.gz
+          sasmodels/sasmodels/sasmodels_docs.tar.gz
         if-no-files-found: error

From 02af1f205bfd8eed06e79b85d3c7492fcab388bc Mon Sep 17 00:00:00 2001
From: Wojciech Potrzebowski <wojciech.potrzebowski@ess.eu>
Date: Sun, 17 Sep 2023 08:41:08 +0100
Subject: [PATCH 15/18] fixing path for upload

---
 .github/workflows/test.yml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index aa3e2f435..59b8ba2de 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -60,6 +60,7 @@ jobs:
       run: |
         make -j 4 -C doc SPHINXOPTS="-W --keep-going -n" html
         tar zcf sasmodels-docs.tar.gz doc/_build/html
+        pwd
         ls -l sasmodels-docs.tar.gz
 
     - name: Publish samodels docs

From 2f1abd8cd0e7bf0a127c35b816fa6c817363cb6c Mon Sep 17 00:00:00 2001
From: Wojciech Potrzebowski <wojciech.potrzebowski@ess.eu>
Date: Sun, 17 Sep 2023 09:34:18 +0100
Subject: [PATCH 16/18] fixing file name

---
 .github/workflows/test.yml | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 59b8ba2de..07c571f90 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -60,8 +60,6 @@ jobs:
       run: |
         make -j 4 -C doc SPHINXOPTS="-W --keep-going -n" html
         tar zcf sasmodels-docs.tar.gz doc/_build/html
-        pwd
-        ls -l sasmodels-docs.tar.gz
 
     - name: Publish samodels docs
       if: ${{ matrix.os == 'ubuntu-latest' }}
@@ -69,5 +67,5 @@ jobs:
       with:
         name: sasmodels-docs-${{ matrix.os }}-${{ matrix.python-version }}
         path: |
-          sasmodels/sasmodels/sasmodels_docs.tar.gz
+          sasmodels-docs.tar.gz
         if-no-files-found: error

From d4d4b4ac47c56ab2252cd56f7fe2697601622069 Mon Sep 17 00:00:00 2001
From: Wojciech Potrzebowski <wojciech.potrzebowski@ess.eu>
Date: Sun, 17 Sep 2023 10:07:50 +0100
Subject: [PATCH 17/18] fixing file name

---
 .github/workflows/test.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 07c571f90..f350a475d 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -62,7 +62,7 @@ jobs:
         tar zcf sasmodels-docs.tar.gz doc/_build/html
 
     - name: Publish samodels docs
-      if: ${{ matrix.os == 'ubuntu-latest' }}
+      if: ${{ matrix.os == 'ubuntu-latest' && matrix.python-version == '3.10'}}
       uses: actions/upload-artifact@v3
       with:
         name: sasmodels-docs-${{ matrix.os }}-${{ matrix.python-version }}

From dbfb021567aa06ae2c43a5190210f2a8db2ef973 Mon Sep 17 00:00:00 2001
From: Wojciech Potrzebowski <Wojciech.Potrzebowski@esss.se>
Date: Wed, 20 Sep 2023 12:20:11 +0100
Subject: [PATCH 18/18] Update test.yml

Skipping tar compression
---
 .github/workflows/test.yml | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index f350a475d..ff0441e5a 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -59,7 +59,6 @@ jobs:
       if: ${{ matrix.os == 'ubuntu-latest' }}
       run: |
         make -j 4 -C doc SPHINXOPTS="-W --keep-going -n" html
-        tar zcf sasmodels-docs.tar.gz doc/_build/html
 
     - name: Publish samodels docs
       if: ${{ matrix.os == 'ubuntu-latest' && matrix.python-version == '3.10'}}
@@ -67,5 +66,5 @@ jobs:
       with:
         name: sasmodels-docs-${{ matrix.os }}-${{ matrix.python-version }}
         path: |
-          sasmodels-docs.tar.gz
+          doc/_build/html
         if-no-files-found: error