source: sasmodels/sasmodels/core.py @ 07646b6

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 07646b6 was 07646b6, checked in by Paul Kienzle <pkienzle@…>, 12 months ago

Merge branch 'cuda-test' into beta_approx

  • Property mode set to 100644
File size: 14.7 KB
Line 
1"""
2Core model handling routines.
3"""
4from __future__ import print_function
5
6__all__ = [
7    "list_models", "load_model", "load_model_info",
8    "build_model", "precompile_dlls",
9    ]
10
11import os
12from os.path import basename, join as joinpath
13from glob import glob
14import re
15
16import numpy as np # type: ignore
17
18from . import generate
19from . import modelinfo
20from . import product
21from . import mixture
22from . import kernelpy
23from . import kernelcuda
24from . import kernelcl
25from . import kerneldll
26from . import custom
27
28# pylint: disable=unused-import
29try:
30    from typing import List, Union, Optional, Any
31    from .kernel import KernelModel
32    from .modelinfo import ModelInfo
33except ImportError:
34    pass
35# pylint: enable=unused-import
36
37CUSTOM_MODEL_PATH = os.environ.get('SAS_MODELPATH', "")
38if CUSTOM_MODEL_PATH == "":
39    CUSTOM_MODEL_PATH = joinpath(os.path.expanduser("~"), ".sasmodels", "custom_models")
40    #if not os.path.isdir(CUSTOM_MODEL_PATH):
41    #    os.makedirs(CUSTOM_MODEL_PATH)
42
43# TODO: refactor composite model support
44# The current load_model_info/build_model does not reuse existing model
45# definitions when loading a composite model, instead reloading and
46# rebuilding the kernel for each component model in the expression.  This
47# is fine in a scripting environment where the model is built when the script
48# starts and is thrown away when the script ends, but may not be the best
49# solution in a long-lived application.  This affects the following functions:
50#
51#    load_model
52#    load_model_info
53#    build_model
54
55KINDS = ("all", "py", "c", "double", "single", "opencl", "1d", "2d",
56         "nonmagnetic", "magnetic")
57def list_models(kind=None):
58    # type: (str) -> List[str]
59    """
60    Return the list of available models on the model path.
61
62    *kind* can be one of the following:
63
64        * all: all models
65        * py: python models only
66        * c: compiled models only
67        * single: models which support single precision
68        * double: models which require double precision
69        * opencl: controls if OpenCL is supperessed
70        * 1d: models which are 1D only, or 2D using abs(q)
71        * 2d: models which can be 2D
72        * magnetic: models with an sld
73        * nommagnetic: models without an sld
74
75    For multiple conditions, combine with plus.  For example, *c+single+2d*
76    would return all oriented models implemented in C which can be computed
77    accurately with single precision arithmetic.
78    """
79    if kind and any(k not in KINDS for k in kind.split('+')):
80        raise ValueError("kind not in " + ", ".join(KINDS))
81    files = sorted(glob(joinpath(generate.MODEL_PATH, "[a-zA-Z]*.py")))
82    available_models = [basename(f)[:-3] for f in files]
83    if kind and '+' in kind:
84        all_kinds = kind.split('+')
85        condition = lambda name: all(_matches(name, k) for k in all_kinds)
86    else:
87        condition = lambda name: _matches(name, kind)
88    selected = [name for name in available_models if condition(name)]
89
90    return selected
91
92def _matches(name, kind):
93    if kind is None or kind == "all":
94        return True
95    info = load_model_info(name)
96    pars = info.parameters.kernel_parameters
97    if kind == "py" and callable(info.Iq):
98        return True
99    elif kind == "c" and not callable(info.Iq):
100        return True
101    elif kind == "double" and not info.single:
102        return True
103    elif kind == "single" and info.single:
104        return True
105    elif kind == "opencl" and info.opencl:
106        return True
107    elif kind == "2d" and any(p.type == 'orientation' for p in pars):
108        return True
109    elif kind == "1d" and all(p.type != 'orientation' for p in pars):
110        return True
111    elif kind == "magnetic" and any(p.type == 'sld' for p in pars):
112        return True
113    elif kind == "nonmagnetic" and any(p.type != 'sld' for p in pars):
114        return True
115    return False
116
117def load_model(model_name, dtype=None, platform='ocl'):
118    # type: (str, str, str) -> KernelModel
119    """
120    Load model info and build model.
121
122    *model_name* is the name of the model, or perhaps a model expression
123    such as sphere*hardsphere or sphere+cylinder.
124
125    *dtype* and *platform* are given by :func:`build_model`.
126    """
127    return build_model(load_model_info(model_name),
128                       dtype=dtype, platform=platform)
129
130def load_model_info(model_string):
131    # type: (str) -> modelinfo.ModelInfo
132    """
133    Load a model definition given the model name.
134
135    *model_string* is the name of the model, or perhaps a model expression
136    such as sphere*cylinder or sphere+cylinder. Use '@' for a structure
137    factor product, e.g. sphere@hardsphere. Custom models can be specified by
138    prefixing the model name with 'custom.', e.g. 'custom.MyModel+sphere'.
139
140    This returns a handle to the module defining the model.  This can be
141    used with functions in generate to build the docs or extract model info.
142    """
143    if "+" in model_string:
144        parts = [load_model_info(part)
145                 for part in model_string.split("+")]
146        return mixture.make_mixture_info(parts, operation='+')
147    elif "*" in model_string:
148        parts = [load_model_info(part)
149                 for part in model_string.split("*")]
150        return mixture.make_mixture_info(parts, operation='*')
151    elif "@" in model_string:
152        p_info, q_info = [load_model_info(part)
153                          for part in model_string.split("@")]
154        return product.make_product_info(p_info, q_info)
155    # We are now dealing with a pure model
156    elif "custom." in model_string:
157        pattern = "custom.([A-Za-z0-9_-]+)"
158        result = re.match(pattern, model_string)
159        if result is None:
160            raise ValueError("Model name in invalid format: " + model_string)
161        model_name = result.group(1)
162        # Use ModelName to find the path to the custom model file
163        model_path = joinpath(CUSTOM_MODEL_PATH, model_name + ".py")
164        if not os.path.isfile(model_path):
165            raise ValueError("The model file {} doesn't exist".format(model_path))
166        kernel_module = custom.load_custom_kernel_module(model_path)
167        return modelinfo.make_model_info(kernel_module)
168    kernel_module = generate.load_kernel_module(model_string)
169    return modelinfo.make_model_info(kernel_module)
170
171
172def build_model(model_info, dtype=None, platform="ocl"):
173    # type: (modelinfo.ModelInfo, str, str) -> KernelModel
174    """
175    Prepare the model for the default execution platform.
176
177    This will return an OpenCL model, a DLL model or a python model depending
178    on the model and the computing platform.
179
180    *model_info* is the model definition structure returned from
181    :func:`load_model_info`.
182
183    *dtype* indicates whether the model should use single or double precision
184    for the calculation.  Choices are 'single', 'double', 'quad', 'half',
185    or 'fast'.  If *dtype* ends with '!', then force the use of the DLL rather
186    than OpenCL for the calculation.
187
188    *platform* should be "dll" to force the dll to be used for C models,
189    otherwise it uses the default "ocl".
190    """
191    composition = model_info.composition
192    if composition is not None:
193        composition_type, parts = composition
194        models = [build_model(p, dtype=dtype, platform=platform) for p in parts]
195        if composition_type == 'mixture':
196            return mixture.MixtureModel(model_info, models)
197        elif composition_type == 'product':
198            P, S = models
199            return product.ProductModel(model_info, P, S)
200        else:
201            raise ValueError('unknown mixture type %s'%composition_type)
202
203    # If it is a python model, return it immediately
204    if callable(model_info.Iq):
205        return kernelpy.PyModel(model_info)
206
207    numpy_dtype, fast, platform = parse_dtype(model_info, dtype, platform)
208    source = generate.make_source(model_info)
209    if platform == "dll":
210        #print("building dll", numpy_dtype)
211        return kerneldll.load_dll(source['dll'], model_info, numpy_dtype)
212    elif platform == "cuda":
213        return kernelcuda.GpuModel(source, model_info, numpy_dtype, fast=fast)
214    else:
215        #print("building ocl", numpy_dtype)
216        return kernelcl.GpuModel(source, model_info, numpy_dtype, fast=fast)
217
218def precompile_dlls(path, dtype="double"):
219    # type: (str, str) -> List[str]
220    """
221    Precompile the dlls for all builtin models, returning a list of dll paths.
222
223    *path* is the directory in which to save the dlls.  It will be created if
224    it does not already exist.
225
226    This can be used when build the windows distribution of sasmodels
227    which may be missing the OpenCL driver and the dll compiler.
228    """
229    numpy_dtype = np.dtype(dtype)
230    if not os.path.exists(path):
231        os.makedirs(path)
232    compiled_dlls = []
233    for model_name in list_models():
234        model_info = load_model_info(model_name)
235        if not callable(model_info.Iq):
236            source = generate.make_source(model_info)['dll']
237            old_path = kerneldll.SAS_DLL_PATH
238            try:
239                kerneldll.SAS_DLL_PATH = path
240                dll = kerneldll.make_dll(source, model_info, dtype=numpy_dtype)
241            finally:
242                kerneldll.SAS_DLL_PATH = old_path
243            compiled_dlls.append(dll)
244    return compiled_dlls
245
246def parse_dtype(model_info, dtype=None, platform=None):
247    # type: (ModelInfo, str, str) -> (np.dtype, bool, str)
248    """
249    Interpret dtype string, returning np.dtype, fast flag and platform.
250
251    Possible types include 'half', 'single', 'double' and 'quad'.  If the
252    type is 'fast', then this is equivalent to dtype 'single' but using
253    fast native functions rather than those with the precision level
254    guaranteed by the OpenCL standard.  'default' will choose the appropriate
255    default for the model and platform.
256
257    Platform preference can be specfied ("ocl", "cuda", "dll"), with the
258    default being OpenCL or CUDA if available, otherwise DLL.  If the dtype
259    name ends with '!' then platform is forced to be DLL rather than GPU.
260    The default platform is set by the environment variable SAS_OPENCL,
261    SAS_OPENCL=driver:device for OpenCL, SAS_OPENCL=cuda:device for CUDA
262    or SAS_OPENCL=none for DLL.
263
264    This routine ignores the preferences within the model definition.  This
265    is by design.  It allows us to test models in single precision even when
266    we have flagged them as requiring double precision so we can easily check
267    the performance on different platforms without having to change the model
268    definition.
269    """
270    # Assign default platform, overriding ocl with dll if OpenCL is unavailable
271    # If opencl=False OpenCL is switched off
272    if platform is None:
273        platform = "ocl"
274
275    # Check if type indicates dll regardless of which platform is given
276    if dtype is not None and dtype.endswith('!'):
277        platform = "dll"
278        dtype = dtype[:-1]
279
280    # Make sure model allows opencl/gpu
281    if not model_info.opencl:
282        platform = "dll"
283
284    # Make sure opencl is available, or fallback to cuda then to dll
285    if platform == "ocl" and not kernelcl.use_opencl():
286        platform = "cuda" if kernelcuda.use_cuda() else "dll"
287
288    # Convert special type names "half", "fast", and "quad"
289    fast = (dtype == "fast")
290    if fast:
291        dtype = "single"
292    elif dtype == "quad":
293        dtype = "longdouble"
294    elif dtype == "half":
295        dtype = "float16"
296
297    # Convert dtype string to numpy dtype.  Use single precision for GPU
298    # if model allows it, otherwise use double precision.
299    if dtype is None or dtype == "default":
300        numpy_dtype = (generate.F32 if model_info.single and platform in ("ocl", "cuda")
301                       else generate.F64)
302    else:
303        numpy_dtype = np.dtype(dtype)
304
305    # Make sure that the type is supported by GPU, otherwise use dll
306    if platform == "ocl":
307        env = kernelcl.environment()
308    elif platform == "cuda":
309        env = kernelcuda.environment()
310    else:
311        env = None
312    if env is not None and not env.has_type(numpy_dtype):
313        platform = "dll"
314        if dtype is None:
315            numpy_dtype = generate.F64
316
317    return numpy_dtype, fast, platform
318
319def list_models_main():
320    # type: () -> None
321    """
322    Run list_models as a main program.  See :func:`list_models` for the
323    kinds of models that can be requested on the command line.
324    """
325    import sys
326    kind = sys.argv[1] if len(sys.argv) > 1 else "all"
327    print("\n".join(list_models(kind)))
328
329def test_composite_order():
330    def test_models(fst, snd):
331        """Confirm that two models produce the same parameters"""
332        fst = load_model(fst)
333        snd = load_model(snd)
334        # Un-disambiguate parameter names so that we can check if the same
335        # parameters are in a pair of composite models. Since each parameter in
336        # the mixture model is tagged as e.g., A_sld, we ought to use a
337        # regex subsitution s/^[A-Z]+_/_/, but removing all uppercase letters
338        # is good enough.
339        fst = [[x for x in p.name if x == x.lower()] for p in fst.info.parameters.kernel_parameters]
340        snd = [[x for x in p.name if x == x.lower()] for p in snd.info.parameters.kernel_parameters]
341        assert sorted(fst) == sorted(snd), "{} != {}".format(fst, snd)
342
343    def build_test(first, second):
344        test = lambda description: test_models(first, second)
345        description = first + " vs. " + second
346        return test, description
347
348    yield build_test(
349        "cylinder+sphere",
350        "sphere+cylinder")
351    yield build_test(
352        "cylinder*sphere",
353        "sphere*cylinder")
354    yield build_test(
355        "cylinder@hardsphere*sphere",
356        "sphere*cylinder@hardsphere")
357    yield build_test(
358        "barbell+sphere*cylinder@hardsphere",
359        "sphere*cylinder@hardsphere+barbell")
360    yield build_test(
361        "barbell+cylinder@hardsphere*sphere",
362        "cylinder@hardsphere*sphere+barbell")
363    yield build_test(
364        "barbell+sphere*cylinder@hardsphere",
365        "barbell+cylinder@hardsphere*sphere")
366    yield build_test(
367        "sphere*cylinder@hardsphere+barbell",
368        "cylinder@hardsphere*sphere+barbell")
369    yield build_test(
370        "barbell+sphere*cylinder@hardsphere",
371        "cylinder@hardsphere*sphere+barbell")
372    yield build_test(
373        "barbell+cylinder@hardsphere*sphere",
374        "sphere*cylinder@hardsphere+barbell")
375
376def test_composite():
377    # type: () -> None
378    """Check that model load works"""
379    #Test the the model produces the parameters that we would expect
380    model = load_model("cylinder@hardsphere*sphere")
381    actual = [p.name for p in model.info.parameters.kernel_parameters]
382    target = ("sld sld_solvent radius length theta phi"
383              " radius_effective volfraction "
384              " structure_factor_mode radius_effective_mode"
385              " A_sld A_sld_solvent A_radius").split()
386    assert target == actual, "%s != %s"%(target, actual)
387
388if __name__ == "__main__":
389    list_models_main()
Note: See TracBrowser for help on using the repository browser.