| 1 | """ |
|---|
| 2 | Core model handling routines. |
|---|
| 3 | """ |
|---|
| 4 | from __future__ import print_function |
|---|
| 5 | |
|---|
| 6 | __all__ = [ |
|---|
| 7 | "list_models", "load_model", "load_model_info", |
|---|
| 8 | "build_model", "precompile_dlls", |
|---|
| 9 | ] |
|---|
| 10 | |
|---|
| 11 | import os |
|---|
| 12 | from os.path import basename, join as joinpath |
|---|
| 13 | from glob import glob |
|---|
| 14 | import re |
|---|
| 15 | |
|---|
| 16 | import numpy as np # type: ignore |
|---|
| 17 | |
|---|
| 18 | from . import generate |
|---|
| 19 | from . import modelinfo |
|---|
| 20 | from . import product |
|---|
| 21 | from . import mixture |
|---|
| 22 | from . import kernelpy |
|---|
| 23 | from . import kernelcuda |
|---|
| 24 | from . import kernelcl |
|---|
| 25 | from . import kerneldll |
|---|
| 26 | from . import custom |
|---|
| 27 | |
|---|
| 28 | # pylint: disable=unused-import |
|---|
| 29 | try: |
|---|
| 30 | from typing import List, Union, Optional, Any |
|---|
| 31 | from .kernel import KernelModel |
|---|
| 32 | from .modelinfo import ModelInfo |
|---|
| 33 | except ImportError: |
|---|
| 34 | pass |
|---|
| 35 | # pylint: enable=unused-import |
|---|
| 36 | |
|---|
| 37 | CUSTOM_MODEL_PATH = os.environ.get('SAS_MODELPATH', "") |
|---|
| 38 | if 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 | |
|---|
| 55 | KINDS = ("all", "py", "c", "double", "single", "opencl", "1d", "2d", |
|---|
| 56 | "nonmagnetic", "magnetic") |
|---|
| 57 | def 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 | |
|---|
| 92 | def _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 | |
|---|
| 117 | def 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 | |
|---|
| 130 | def 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 | |
|---|
| 172 | def 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 | |
|---|
| 218 | def 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 | |
|---|
| 246 | def 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 | |
|---|
| 319 | def 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 | |
|---|
| 329 | def 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 | |
|---|
| 376 | def 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 | |
|---|
| 388 | if __name__ == "__main__": |
|---|
| 389 | list_models_main() |
|---|
