Skip to content
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions example/fit.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
8 changes: 6 additions & 2 deletions example/simul_fit.py
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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)
2 changes: 1 addition & 1 deletion sasmodels/bumps_model.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
49 changes: 28 additions & 21 deletions sasmodels/data.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -418,8 +419,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.

Expand All @@ -436,14 +437,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.

Expand All @@ -466,10 +467,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):
Expand All @@ -495,7 +496,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.
Expand All @@ -520,6 +521,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)
Expand All @@ -543,8 +547,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':
Expand Down Expand Up @@ -578,7 +582,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)
Expand Down Expand Up @@ -646,7 +650,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.
Expand All @@ -658,12 +662,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)
Expand All @@ -682,7 +689,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)$')
Expand All @@ -691,7 +698,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'
Expand All @@ -702,14 +709,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,
Expand All @@ -724,8 +731,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:
Expand Down