Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Allow pure python Fq for beta approximation. #570

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Allow pure python Fq for beta approximation. #570

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
8822358
Allow pure python Fq for beta approximation.
Jun 23, 2023
e8935f2
Add flag for whether model is compiled or not.
Jun 23, 2023
ad9fb9a
Merge branch 'master' into pr-python-fq
Jan 3, 2024
aa03630
Merge branch 'master' into pr-python-fq
pkienzle Nov 26, 2024
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
2 changes: 1 addition & 1 deletion sasmodels/compare.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -551,7 +551,7 @@ def constrain_pars(model_info, pars):
name = name.split('*')[0]

# Suppress magnetism for python models (not yet implemented)
if callable(model_info.Iq):
if not model_info.compiled:
pars.update(suppress_magnetism(pars))

if name == 'barbell':
Expand Down
6 changes: 3 additions & 3 deletions sasmodels/core.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -96,7 +96,7 @@ def _matches(name, kind):
info = load_model_info(name)
pars = info.parameters.kernel_parameters
# TODO: may be adding Fq to the list at some point
is_pure_py = callable(info.Iq)
is_pure_py = not info.compiled
if kind == "py":
return is_pure_py
elif kind == "c":
Expand Down Expand Up @@ -327,7 +327,7 @@ def build_model(model_info, dtype=None, platform="ocl"):
raise ValueError('unknown mixture type %s'%composition_type)

# If it is a python model, return it immediately
if callable(model_info.Iq):
if not model_info.compiled:
from . import kernelpy
return kernelpy.PyModel(model_info)

Expand Down Expand Up @@ -364,7 +364,7 @@ def precompile_dlls(path, dtype="double"):
compiled_dlls = []
for model_name in list_models():
model_info = load_model_info(model_name)
if not callable(model_info.Iq):
if model_info.compiled:
source = generate.make_source(model_info)['dll']
old_path = kerneldll.SAS_DLL_PATH
try:
Expand Down
3 changes: 1 addition & 2 deletions sasmodels/generate.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -935,9 +935,8 @@ def make_source(model_info):
Uses source files found in the given search path. Returns None if this
is a pure python model, with no C source components.
"""
if callable(model_info.Iq):
if not model_info.compiled:
raise ValueError("can't compile python model")
#return None

# TODO: need something other than volume to indicate dispersion parameters
# No volume normalization despite having a volume parameter.
Expand Down
37 changes: 34 additions & 3 deletions sasmodels/kernelpy.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -161,10 +161,17 @@ def __init__(self, model_info, q_input):
form = model_info.Iqxy
qx, qy = q_input.q[:, 0], q_input.q[:, 1]
self._form = lambda: form(qx, qy, *kernel_args)
self._nq_copies = 1
else:
form = model_info.Iq
q = q_input.q
self._form = lambda: form(q, *kernel_args)
if model_info.have_Fq:
form = model_info.Fq
self._form = lambda: _pack_fq(*form(q, *kernel_args))
self._nq_copies = 2
else:
form = model_info.Iq
self._form = lambda: form(q, *kernel_args)
self._nq_copies = 1

# Generate a closure which calls the form_volume if it exists.
self._volume_args = volume_args
Expand All @@ -188,7 +195,7 @@ def _call_kernel(self, call_details, values, cutoff, magnetic, radius_effective_
else (lambda: self._radius(radius_effective_mode)))
self.result = _loops(
self._parameter_vector, self._form, self._volume, radius,
self.q_input.nq, call_details, values, cutoff)
self.q_input.nq*self._nq_copies, call_details, values, cutoff)

def release(self):
# type: () -> None
Expand All @@ -198,6 +205,13 @@ def release(self):
self.q_input.release()
self.q_input = None

def _pack_fq(fq, fqsq):
result = np.empty(2*fq.shape[0], dtype=fq.dtype)
# PAK: For reasons unknown, the original author (me) decided to pack the
# values into the returned array as (f^2, f) pairs.
result[0::2] = fqsq
result[1::2] = fq
return result

def _loops(parameters, form, form_volume, form_radius, nq, call_details,
values, cutoff):
Expand Down Expand Up @@ -291,6 +305,7 @@ def _create_default_functions(model_info):
"""
# Note: Must call create_vector_Iq before create_vector_Iqxy.
_create_vector_Iq(model_info)
_create_vector_Fq(model_info)
_create_vector_Iqxy(model_info)


Expand All @@ -309,6 +324,22 @@ def vector_Iq(q, *args):
vector_Iq.vectorized = True
model_info.Iq = vector_Iq

def _create_vector_Fq(model_info):
"""
Define Fq as a vector function if it exists.
"""
# Note that this is doing slightly too much work since we are composing
# two vector results which are then interlaced via _pack_fq above.
Fq = model_info.Fq
if callable(Fq) and not getattr(Fq, 'vectorized', False):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I take it you're using callable here because model_info.compiled isn't exactly what you want. If it is, consider being consistent with the other changes.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

info.compiled is for info.Fq and info.Iq jointly. Only one of them will be defined, and it may be defined using C source or using a callable python function.

It occurs to me that we could use some numba transforms here so that the code runs as fast as a C model on the CPU (or faster since loop parallelization can run over multiple cores). Python models get much easier to write when we don't have to vectorize across if statements.

def vector_Fq(q, *args):
"""
Vectorized 1D kernel returning Fq, Fq**2
"""
fq, fqsq = zip(*(Fq(qi, *args) for qi in q))
return np.array(fq), np.array(fqsq)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not certain, but slicing a numpy array might be considerably faster than zip(*)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Surprisingly, the zip version is faster:

In [15]: %timeit x,y = [np.array(v) for v in zip(*(((i, i) for i in np.arange(1000))))]
242 μs ± 7.15 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

In [16]: %timeit x,y = (a:=np.array(list((i, i) for i in np.arange(1000))))[:,0],a[:,1]
389 μs ± 4.75 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

But in any case, this cost is small compared to evaluating the function at each q.

vector_Fq.vectorized = True
model_info.Fq = vector_Fq

def _create_vector_Iqxy(model_info):
"""
Expand Down
4 changes: 2 additions & 2 deletions sasmodels/model_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -111,10 +111,10 @@ def _add_model_to_suite(loaders, suite, model_info):
#print('found tests in', model_name)
#print('------')

# if ispy then use the dll loader to call pykernel
# if is python then use the dll loader to call pykernel
# don't try to call cl kernel since it will not be
# available in some environmentes.
is_py = callable(model_info.Iq)
is_py = not model_info.compiled

# Some OpenCL drivers seem to be flaky, and are not producing the
# expected result. Since we don't have known test values yet for
Expand Down
54 changes: 29 additions & 25 deletions sasmodels/modelinfo.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -862,6 +862,7 @@ def isstr(x):
return isinstance(x, str)


# Note: adding 'Fq' even though we can't yet use it to define C code in python.
#: Set of variables defined in the model that might contain C code
C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc',
'form_volume', 'shell_volume', 'c_code', 'valid']
Expand All @@ -882,8 +883,7 @@ def _find_source_lines(model_info, kernel_module):
"""
# Only need line numbers if we are creating a C module and the C symbols
# are defined.
if (callable(model_info.Iq)
or not any(hasattr(model_info, s) for s in C_SYMBOLS)):
if not model_info.compiled:
return

# load the module source if we can
Expand Down Expand Up @@ -943,9 +943,7 @@ def make_model_info(kernel_module):
info.docs = kernel_module.__doc__
info.category = getattr(kernel_module, 'category', None)
info.structure_factor = getattr(kernel_module, 'structure_factor', False)
# TODO: find Fq by inspection
info.radius_effective_modes = getattr(kernel_module, 'radius_effective_modes', None)
info.have_Fq = getattr(kernel_module, 'have_Fq', False)
info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y'])
# Note: custom.load_custom_kernel_module assumes the C sources are defined
# by this attribute.
Expand All @@ -958,15 +956,19 @@ def make_model_info(kernel_module):
info.form_volume = getattr(kernel_module, 'form_volume', None) # type: ignore
info.shell_volume = getattr(kernel_module, 'shell_volume', None) # type: ignore
info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore
info.Fq = getattr(kernel_module, 'Fq', None) # type: ignore
# TODO: We should be able to find Fq in C code by inspection.
info.have_Fq = getattr(kernel_module, 'have_Fq', (info.Fq is not None))
info.Iqxy = getattr(kernel_module, 'Iqxy', None) # type: ignore
info.Iqac = getattr(kernel_module, 'Iqac', None) # type: ignore
info.Iqabc = getattr(kernel_module, 'Iqabc', None) # type: ignore
info.Imagnetic = getattr(kernel_module, 'Imagnetic', None) # type: ignore
info.profile = getattr(kernel_module, 'profile', None) # type: ignore
info.sesans = getattr(kernel_module, 'sesans', None) # type: ignore
# Default single and opencl to True for C models. Python models have callable Iq.
info.opencl = getattr(kernel_module, 'opencl', not callable(info.Iq))
info.single = getattr(kernel_module, 'single', not callable(info.Iq))
info.compiled = not callable(info.Iq) and not callable(info.Fq)
info.opencl = getattr(kernel_module, 'opencl', info.compiled)
info.single = getattr(kernel_module, 'single', info.compiled)
info.random = getattr(kernel_module, 'random', None)
info.hidden = getattr(kernel_module, 'hidden', None) # type: ignore

Expand All @@ -975,7 +977,7 @@ def make_model_info(kernel_module):
if control is not None:
parameters[control].is_control = True

if callable(info.Iq) and parameters.has_2d:
if not info.compiled and parameters.has_2d:
raise ValueError("oriented python models not supported")

# CRUFT: support old-style ER() for effective radius
Expand All @@ -989,7 +991,7 @@ def make_model_info(kernel_module):
# so just issue a warning if we see ER in a C model.
ER = getattr(kernel_module, 'ER', None)
if ER is not None:
if callable(info.Iq) and info.radius_effective is None:
if not info.compiled and info.radius_effective is None:
info.radius_effective_modes = ['ER']
info.radius_effective = lambda mode, *args: ER(*args)
# TODO: uncomment the following for the sasview 4.3 release
Expand Down Expand Up @@ -1092,13 +1094,12 @@ class ModelInfo(object):
#: the model cannot be run in opencl (e.g., because the model passes
#: functions by reference), then set this to false.
opencl = None # type: bool
#: True if the model is compiled with C or OpenCL
compiled = None # type: bool
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It might make sense to make this a dataclass and avoid having variables initialised to values outside their type

#: True if the model is a structure factor used to model the interaction
#: between form factor models. This will default to False if it is not
#: provided in the file.
structure_factor = None # type: bool
#: True if the model defines an Fq function with signature
#: ``void Fq(double q, double *F1, double *F2, ...)``
have_Fq = False
#: List of options for computing the effective radius of the shape,
#: or None if the model is not usable as a form factor model.
radius_effective_modes = None # type: List[str]
Expand Down Expand Up @@ -1136,20 +1137,16 @@ class ModelInfo(object):
#: monodisperse approximation for non-dilute solutions, P@S. The first
#: argument is the integer effective radius mode, with default 0.
radius_effective = None # type: Union[None, Callable[[int, np.ndarray], float]]
#: Returns *I(q, a, b, ...)* for parameters *a*, *b*, etc. defined
#: by the parameter table. *Iq* can be defined as a python function, or
#: as a C function. If it is defined in C, then set *Iq* to the body of
#: the C function, including the return statement. This function takes
#: values for *q* and each of the parameters as separate *double* values
#: (which may be converted to float or long double by sasmodels). All
#: source code files listed in :attr:`source` will be loaded before the
#: *Iq* function is defined. If *Iq* is not present, then sources should
#: define *static double Iq(double q, double a, double b, ...)* which
#: will return *I(q, a, b, ...)*. Multiplicity parameters are sent as
#: pointers to doubles. Constants in floating point expressions should
#: include the decimal point. See :mod:`.generate` for more details. If
#: *have_Fq* is True, then Iq should return an interleaved array of
#: $[\sum F(q_1), \sum F^2(q_1), \ldots, \sum F(q_n), \sum F^2(q_n)]$.
#: Returns *I(q, a, b, ...)* for parameters *a*, *b*, etc. defined by the
#: parameter table. Multiplicity parameters such as the number of shells are
#: sent as floating point values. If the function can operate with a vector
#: of *q* values (that is, you aren't doing simple comparisons such as
#: *q > 0*), then set *Iq.vectorized = True* to make your code run faster.
#: You can also set *Iq* to a string containing the body of a C function.
#: For example, ``Iq = return (a*q + b)*q + c;`` will generate a quadratic.
#: All source code files listed in :attr:`source` will be loaded before the
#: C function is defined. Constants in floating point expressions should
#: include the decimal point. See :mod:`.generate` for more details.
Iq = None # type: Union[None, str, Callable[[...], np.ndarray]]
#: Returns *I(qx, qy, a, b, ...)*. The interface follows :attr:`Iq`.
Iqxy = None # type: Union[None, str, Callable[[...], np.ndarray]]
Expand All @@ -1159,6 +1156,13 @@ class ModelInfo(object):
Iqabc = None # type: Union[None, str, Callable[[...], np.ndarray]]
#: Returns *I(qx, qy, a, b, ...)*. The interface follows :attr:`Iq`.
Imagnetic = None # type: Union[None, str, Callable[[...], np.ndarray]]
#: Returns *F(q, a, b, ...), F^2(q, a, b, c, ...)*. Note: you cannot assign
#: a C source code body to *Fq*.
Fq = None # type: Union[None, Callable[[...], Tuple[np.ndarray,np.ndarray]]]
#: True if the model defines an Fq function in C with signature
#: ``void Fq(double q, double *F1, double *F2, ...)``
#: or in python as ``Fq(q, ...) -> (fq, fq^2)``.
have_Fq = False
#: Returns a model profile curve *x, y*. If *profile* is defined, this
#: curve will appear in response to the *Show* button in SasView. Use
#: :attr:`profile_axes` to set the axis labels. Note that *y* values
Expand Down
32 changes: 19 additions & 13 deletions sasmodels/models/_spherepy.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -79,23 +79,29 @@ def radius_effective(mode, radius):
"""Calculate R_eff for sphere"""
return radius if mode else 0.

def Iq(q, sld, sld_solvent, radius):
"""Calculate I(q) for sphere"""
# Variants for testing purposes: toggle True/False
vectorized = True # Whether to call a q vector or loop over q scalars.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If they were for prototyping, consider removing them, if they're for testing, they and their dependents should be settable after loading the file.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

prototyping. It could be moved to the explore directory.

have_Fq = True # Whether to use Fq or Iq for evaluation.
def Fq(q, sld, sld_solvent, radius):
"""Calculate F(q), F^2(q) for sphere"""
#print "q",q
#print "sld,r",sld,sld_solvent,radius
qr = q * radius
sn, cn = sin(qr), cos(qr)
## The natural expression for the Bessel function is the following:
## bes = 3 * (sn-qr*cn)/qr**3 if qr>0 else 1
## however, to support vector q values we need to handle the conditional
## as a vector, which we do by first evaluating the full expression
## everywhere, then fixing it up where it is broken. We should probably
## set numpy to ignore the 0/0 error before we do though...
bes = 3 * (sn - qr * cn) / qr ** 3 # may be 0/0 but we fix that next line
bes[qr == 0] = 1
fq = bes * (sld - sld_solvent) * form_volume(radius)
return 1.0e-4 * fq ** 2
Iq.vectorized = True # Iq accepts an array of q values
if vectorized:
with np.errstate(all='ignore'):
bes = 3 * (sn - qr * cn) / qr ** 3 # may be 0/0 but we fix that next line
bes[qr == 0] = 1
else:
bes = 3 * (sn-qr*cn)/qr**3 if qr != 0 else 1
fq = bes * (1e-2 * (sld - sld_solvent) * form_volume(radius))
return fq, fq**2
Fq.vectorized = vectorized # Fq accepts an array of q value

def Iq(q, sld, sld_solvent, radius):
"""Calculate I(q) for sphere"""
return Fq(q, sld, sld_solvent, radius)[1]
Iq.vectorized = vectorized # Iq accepts an array of q value

def sesans(z, sld, sld_solvent, radius):
"""
Expand Down