Send data mask from theory calculator into plotter

example/fit.py

@@ -13,6 +13,9 @@
 set_beam_stop(radial_data, 0.00669, outer=0.025)
 set_top(radial_data, -.0185)
 
+# Set data limits on fit
+#radial_data.qmin, radial_data.qmax = 0.01, 0.02
+
 tan_data = load_data('DEC07266.DAT')
 set_beam_stop(tan_data, 0.00669, outer=0.025)
 set_top(tan_data, -.0185)

example/simul_fit.py

@@ -1,13 +1,15 @@
 import numpy as np
 from bumps.names import FitProblem, FreeVariables
 
+import sasdata
+
 from sasmodels.bumps_model import Experiment, Model
 from sasmodels.core import load_model
 from sasmodels.data import load_data
 
 # latex data, same sample usans and sans
 # particles radius ~2300, uniform dispersity
-datasets = load_data('latex_smeared.xml', index='all')
+datasets = load_data(str(sasdata.data_path / '1d_data' / 'latex_smeared.xml'), index='all')
 #[print(data) for data in datasets]
 
 # A single sphere model to share between the datasets.  We will use
@@ -57,7 +59,9 @@
 #    model1.background.range(0, 2)
 #    model2.background.range(0, 2)
 
+# Set data fitting limits
+#datasets[1].qmin, datasets[1].qmax = 1e-4, 4e-3
+
 # Setup the experiments, sharing the same model across all datasets.
 M = [Experiment(data=data, model=model, name=data.run[0]) for data in datasets]
-
 problem = FitProblem(M, freevars=free)

sasmodels/TwoYukawa/CalTYSk.py

@@ -1,8 +1,8 @@
 import numpy as np
-from numpy import pi, mean
+from numpy import mean, pi
 
-from .Ecoefficient import TYCoeff
 from .CalcRealRoot import CalRealRoot
+from .Ecoefficient import TYCoeff
 from .TInvFourier import TInvFourier
 
 # Supplied Q vector must go out to at least this value to calculate g(r).

sasmodels/TwoYukawa/CalcRealRoot.py

@@ -3,6 +3,7 @@
 from .Ecoefficient import TYCoeff
 from .Epoly import make_poly
 
+
 def CalRealRoot(coeff: TYCoeff):
 
     Ecoefficient = coeff.Ecoefficient

sasmodels/TwoYukawa/Ecoefficient.py

@@ -1,7 +1,7 @@
 from typing import Callable, Tuple
 
 import numpy as np
-from numpy import exp, pi, cos, sin, cosh
+from numpy import cos, cosh, exp, pi, sin
 from numpy.typing import NDArray
 
 # CalCoeff.m
@@ -194,17 +194,17 @@ def ABC12(d1, d2):
               Ccd2_22*Cdd1_12*d1*d2**2 -
               Ccd1_21*Cdd2_22*d1*d2**2 -
               Ccd2_22*d2*k1 + Ccd1_21*d1*k2)))/
-        ((d1*((
+        (d1*(
               (-Ccd1_21)*Ccd2_12*d2 +
-                Ccd1_11*Ccd2_22*d2)))))
+                Ccd1_11*Ccd2_22*d2)))
         C2 = ((-((Ccd2_12*d2*
                 (((-Cd1_1)*d1 - Cdd1_11*d1**2 -
                     Cdd1_12*d1*d2 + k1)) -
               Ccd1_11*d1*
                 (((-Cd2_2)*d2 - Cdd2_12*d1*d2 -
                     Cdd2_22*d2**2 + k2))))) /
-        (((-Ccd1_21)*Ccd2_12*d1*d2 +
-            Ccd1_11*Ccd2_22*d1*d2)))
+        ((-Ccd1_21)*Ccd2_12*d1*d2 +
+            Ccd1_11*Ccd2_22*d1*d2))
         return A, B, C1, C2
     self.ABC12 = ABC12
 
@@ -370,10 +370,10 @@ def tau_d21(s):
 
     E1d02 = 12*c1F01*phi*sigma_d01(z1) - 12*c1F01*exp(-z1)*phi*tau_d01(z1)
 
-    E1d11 = (((12*c1F10*phi*sigma_d01(z1) + \
+    E1d11 = (12*c1F10*phi*sigma_d01(z1) + \
           12*c1F01*phi*sigma_d10(z1) - \
           12*c1F10*exp(-z1)*phi*tau_d01(z1) - \
-          12*c1F01*exp(-z1)*phi*tau_d10(z1))))
+          12*c1F01*exp(-z1)*phi*tau_d10(z1))
     E1d12 = (-c1F01 + 12*c1F11*phi*sigma_d01(z1) +
              12*c1F01*phi*sigma_d11(z1) -
              12*c1F11*exp(-z1)*phi*tau_d01(z1) -

sasmodels/TwoYukawa/TFourier.py

@@ -1,6 +1,7 @@
-from numpy import exp, pi, arange, linspace, abs, ceil, log2, interp
+from numpy import abs, arange, ceil, exp, interp, linspace, log2, pi
 from numpy.fft import fft
 
+
 def TFourier(x, y, deltaX):
     """
     Compute the Fourier transform of a function y(x) with sampling interval deltaX.

sasmodels/TwoYukawa/TInvFourier.py

@@ -1,6 +1,7 @@
-from numpy import exp, ceil, log2, pi, arange, interp
+from numpy import arange, ceil, exp, interp, log2, pi
 from numpy.fft import fft
 
+
 def TInvFourier(x, y, deltaX):
     """
     Inverse Fourier transform implementation

sasmodels/bumps_model.py

@@ -266,7 +266,7 @@ def plot(self, view=None):
         data, theory, resid = self._data, self.theory(), self.residuals()
         # TODO: hack to display oriented usans 2-D pattern
         Iq_calc = self.Iq_calc if isinstance(self.Iq_calc, tuple) else None
-        plot_theory(data, theory, resid, view, Iq_calc=Iq_calc)
+        plot_theory(data, theory, resid, view, Iq_calc=Iq_calc, index=self.index)
 
     def simulate_data(self, noise=None):
         # type: (float) -> None

sasmodels/data.py

@@ -40,10 +40,11 @@
 
 # pylint: disable=unused-import
 try:
-    from typing import Callable, Optional, Union
+    from typing import Optional, Union
     Data = Union["Data1D", "Data2D", "SesansData"]
     OptArray = Optional[np.ndarray]
     OptLimits = Optional[tuple[float, float]]
+    OptIndex = np.ndarray | slice | None
     OptString = Optional[str]
 except ImportError:
     pass
@@ -413,8 +414,8 @@ def empty_data2D(qx, qy=None, resolution=0.0):
     return data
 
 
-def plot_data(data, view=None, limits=None):
-    # type: (Data, str, OptLimits) -> None
+def plot_data(data, view=None, limits=None, index=None):
+    # type: (Data, str, OptLimits, OptIndex) -> None
     """
     Plot data loaded by the sasview loader.
 
@@ -431,14 +432,14 @@ def plot_data(data, view=None, limits=None):
     if hasattr(data, 'isSesans') and data.isSesans:
         _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)
+        _plot_result2D(data, None, None, view, use_data=True, limits=limits, index=index)
     else:
-        _plot_result1D(data, None, None, view, use_data=True, limits=limits)
+        _plot_result1D(data, None, None, view, use_data=True, limits=limits, index=index)
 
 
 def plot_theory(data, theory, resid=None, view=None, use_data=True,
-                limits=None, Iq_calc=None):
-    # type: (Data, OptArray, OptArray, OptString, bool, OptLimits, OptArray) -> None
+                limits=None, Iq_calc=None, index=None):
+    # type: (Data, OptArray, OptArray, OptString, bool, OptLimits, OptArray, OptIndex) -> None
     """
     Plot theory calculation.
 
@@ -461,10 +462,10 @@ def plot_theory(data, theory, resid=None, view=None, use_data=True,
     if hasattr(data, 'isSesans') and data.isSesans:
         _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)
+        _plot_result2D(data, theory, resid, view, use_data, limits=limits, index=index)
     else:
         _plot_result1D(data, theory, resid, view, use_data,
-                       limits=limits, Iq_calc=Iq_calc)
+                       limits=limits, Iq_calc=Iq_calc, index=index)
 
 
 def protect(func):
@@ -490,7 +491,7 @@ def wrapper(*args, **kw):
 
 @protect
 def _plot_result1D(data, theory, resid, view, use_data,
-                   limits=None, Iq_calc=None):
+                   limits=None, Iq_calc=None, index=None):
     # type: (Data1D, OptArray, OptArray, str, bool, OptLimits, OptArray) -> None
     """
     Plot the data and residuals for 1D data.
@@ -515,6 +516,9 @@ def _plot_result1D(data, theory, resid, view, use_data,
 
     scale = 1e8 * data.x**4 if view == 'q4' else 1.0
 
+    if index is None:
+        index = slice(None)
+
     if use_data or use_theory:
         if num_plots > 1:
             plt.subplot(1, num_plots, 1)
@@ -538,8 +542,8 @@ def _plot_result1D(data, theory, resid, view, use_data,
             # Note: masks merge, so any masked theory points will stay masked,
             # and the data mask will be added to it.
             #mtheory = masked_array(theory, data.mask.copy())
-            theory_x = data.x[data.mask == 0]
-            theory_scale = scale if np.isscalar(scale) else scale[data.mask == 0]
+            theory_x = data.x[index]
+            theory_scale = scale if np.isscalar(scale) else scale[index]
             mtheory = masked_array(theory)
             mtheory[~np.isfinite(mtheory)] = masked
             if view == 'log':
@@ -573,7 +577,7 @@ def _plot_result1D(data, theory, resid, view, use_data,
         #plt.axis('equal')
 
     if use_resid:
-        theory_x = data.x[data.mask == 0]
+        theory_x = data.x[index]
         mresid = masked_array(resid)
         mresid[~np.isfinite(mresid)] = masked
         some_present = (mresid.count() > 0)
@@ -641,7 +645,7 @@ def _plot_result_sesans(data, theory, resid, view, use_data, limits=None):
 
 
 @protect
-def _plot_result2D(data, theory, resid, view, use_data, limits=None):
+def _plot_result2D(data, theory, resid, view, use_data, limits=None, index=None):
     # type: (Data2D, OptArray, OptArray, str, bool, OptLimits) -> None
     """
     Plot the data and residuals for 2D data.
@@ -653,12 +657,15 @@ def _plot_result2D(data, theory, resid, view, use_data, limits=None):
     use_theory = theory is not None
     use_resid = resid is not None
     num_plots = use_data + use_theory + use_resid
+    mask = ~data.mask
+    if index is None:
+        index = mask
 
     # Put theory and data on a common colormap scale
     vmin, vmax = np.inf, -np.inf
     target = None  # type: OptArray
     if use_data:
-        target = data.data[~data.mask]
+        target = data.data[mask] # full data
         datamin = target[target > 0].min() if view == 'log' else target.min()
         datamax = target.max()
         vmin = min(vmin, datamin)
@@ -677,7 +684,7 @@ def _plot_result2D(data, theory, resid, view, use_data, limits=None):
     if use_data:
         if num_plots > 1:
             plt.subplot(1, num_plots, 1)
-        _plot_2d_signal(data, target, view=view, vmin=vmin, vmax=vmax)
+        _plot_2d_signal(data, target, view=view, vmin=vmin, vmax=vmax, index=mask)
         plt.title('data')
         h = plt.colorbar()
         h.set_label('$I(q)$')
@@ -686,7 +693,7 @@ def _plot_result2D(data, theory, resid, view, use_data, limits=None):
     if use_theory:
         if num_plots > 1:
             plt.subplot(1, num_plots, use_data+1)
-        _plot_2d_signal(data, theory, view=view, vmin=vmin, vmax=vmax)
+        _plot_2d_signal(data, theory, view=view, vmin=vmin, vmax=vmax, index=index)
         plt.title('theory')
         h = plt.colorbar()
         h.set_label(r'$\log_{10}I(q)$' if view == 'log'
@@ -697,14 +704,14 @@ def _plot_result2D(data, theory, resid, view, use_data, limits=None):
     if use_resid:
         if num_plots > 1:
             plt.subplot(1, num_plots, use_data+use_theory+1)
-        _plot_2d_signal(data, resid, view='linear')
+        _plot_2d_signal(data, resid, view='linear', index=index)
         plt.title('residuals')
         h = plt.colorbar()
         h.set_label(r'$\Delta I(q)$')
 
 
 @protect
-def _plot_2d_signal(data, signal, vmin=None, vmax=None, view=None):
+def _plot_2d_signal(data, signal, vmin=None, vmax=None, view=None, index=None):
     # type: (Data2D, np.ndarray, Optional[float], Optional[float], str) -> tuple[float, float]
     """
     Plot the target value for the data.  This could be the data itself,
@@ -719,8 +726,8 @@ def _plot_2d_signal(data, signal, vmin=None, vmax=None, view=None):
         view = 'log'
 
     image = np.zeros_like(data.qx_data)
-    image[~data.mask] = signal
-    valid = np.isfinite(image) & ~data.mask
+    image[index] = signal
+    valid = np.isfinite(image) & index
     if view == 'log':
         valid &= image > 0
         if vmin is None:

sasmodels/models/two_yukawa.py

@@ -94,7 +94,7 @@
 from numpy import inf
 
 # TODO: pep8 says packages and modules should not use camel case
-from sasmodels.TwoYukawa.CalTYSk import CalTYSk, K_MIN, Z_MIN, Z_MIN_DIFF
+from sasmodels.TwoYukawa.CalTYSk import K_MIN, Z_MIN, Z_MIN_DIFF, CalTYSk
 
 # If you want a customized version of two_yukawa as a plugin (for example,
 # because you want to use the high precision polynomial root solver from mpmath)
@@ -104,9 +104,9 @@
 #     https://github.com/SasView/sasmodels/tree/master/sasmodels/models/two_yukawa.py
 #     https://github.com/SasView/sasmodels/tree/master/sasmodels/TwoYukawa
 if 0:
+    import importlib.util
     import sys
     from pathlib import Path
-    import importlib.util
 
     # Remove existing TwoYukawa from sys.modules to force a reload
     remove = [modname for modname in sys.modules if modname.startswith('TwoYukawa.') or modname == 'TwoYukawa']
@@ -121,7 +121,7 @@
     sys.modules['TwoYukawa'] = module
 
     # Override sasmodels library symbols with the local symbols.
-    from TwoYukawa.CalTYSk import CalTYSk, K_MIN, Z_MIN, Z_MIN_DIFF
+    from TwoYukawa.CalTYSk import K_MIN, Z_MIN, Z_MIN_DIFF, CalTYSk
 
 name = "two_yukawa"
 title = "User model for two Yukawa structure factor (S(q))"