Send data mask from theory calculator into plotter

Author: Paul Kienzle

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,14 @@
 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 +58,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/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

@@ -413,7 +413,7 @@ def empty_data2D(qx, qy=None, resolution=0.0):
     return data
 
 
-def plot_data(data, view=None, limits=None):
+def plot_data(data, view=None, limits=None, index=None):
     # type: (Data, str, OptLimits) -> None
     """
     Plot data loaded by the sasview loader.
@@ -431,13 +431,13 @@ 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):
+                limits=None, Iq_calc=None, index=None):
     # type: (Data, OptArray, OptArray, OptString, bool, OptLimits, OptArray) -> None
     """
     Plot theory calculation.
@@ -461,10 +461,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 +490,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 +515,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 +541,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 +576,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 +644,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 +656,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 +683,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 +692,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 +703,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 +725,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: