source: sasmodels/sasmodels/kerneldll.py @ 01c8d9e

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 01c8d9e was 01c8d9e, checked in by Suczewski <ges3@…>, 13 months ago

beta approximation, first pass

  • Property mode set to 100644
File size: 16.1 KB
Line 
1r"""
2DLL driver for C kernels
3
4If the environment variable *SAS_OPENMP* is set, then sasmodels
5will attempt to compile with OpenMP flags so that the model can use all
6available kernels.  This may or may not be available on your compiler
7toolchain.  Depending on operating system and environment.
8
9Windows does not have provide a compiler with the operating system.
10Instead, we assume that TinyCC is installed and available.  This can
11be done with a simple pip command if it is not already available::
12
13    pip install tinycc
14
15If Microsoft Visual C++ is available (because VCINSTALLDIR is
16defined in the environment), then that will be used instead.
17Microsoft Visual C++ for Python is available from Microsoft:
18
19    `<http://www.microsoft.com/en-us/download/details.aspx?id=44266>`_
20
21If neither compiler is available, sasmodels will check for *MinGW*,
22the GNU compiler toolchain. This available in packages such as Anaconda
23and PythonXY, or available stand alone. This toolchain has had
24difficulties on some systems, and may or may not work for you.
25
26You can control which compiler to use by setting SAS_COMPILER in the
27environment:
28
29  - tinycc (Windows): use the TinyCC compiler shipped with SasView
30  - msvc (Windows): use the Microsoft Visual C++ compiler
31  - mingw (Windows): use the MinGW GNU cc compiler
32  - unix (Linux): use the system cc compiler.
33  - unix (Mac): use the clang compiler. You will need XCode installed, and
34    the XCode command line tools. Mac comes with OpenCL drivers, so generally
35    this will not be needed.
36
37Both *msvc* and *mingw* require that the compiler is available on your path.
38For *msvc*, this can done by running vcvarsall.bat in a windows terminal.
39Install locations are system dependent, such as:
40
41    C:\Program Files (x86)\Common Files\Microsoft\Visual
42    C++ for Python\9.0\vcvarsall.bat
43
44or maybe
45
46    C:\Users\yourname\AppData\Local\Programs\Common\Microsoft\Visual
47    C++ for Python\9.0\vcvarsall.bat
48
49OpenMP for *msvc* requires the Microsoft vcomp90.dll library, which doesn't
50seem to be included with the compiler, nor does there appear to be a public
51download location.  There may be one on your machine already in a location
52such as:
53
54    C:\Windows\winsxs\x86_microsoft.vc90.openmp*\vcomp90.dll
55
56If you copy this to somewhere on your path, such as the python directory or
57the install directory for this application, then OpenMP should be supported.
58
59For full control of the compiler, define a function
60*compile_command(source,output)* which takes the name of the source file
61and the name of the output file and returns a compile command that can be
62evaluated in the shell.  For even more control, replace the entire
63*compile(source,output)* function.
64
65The global attribute *ALLOW_SINGLE_PRECISION_DLLS* should be set to *False* if
66you wish to prevent single precision floating point evaluation for the compiled
67models, otherwise set it defaults to *True*.
68"""
69from __future__ import print_function
70
71import sys
72import os
73from os.path import join as joinpath, splitext
74import subprocess
75import tempfile
76import ctypes as ct  # type: ignore
77import _ctypes as _ct
78import logging
79
80import numpy as np  # type: ignore
81
82try:
83    import tinycc
84except ImportError:
85    tinycc = None
86
87from . import generate
88from .kernel import KernelModel, Kernel
89from .kernelpy import PyInput
90from .exception import annotate_exception
91from .generate import F16, F32, F64
92
93# pylint: disable=unused-import
94try:
95    from typing import Tuple, Callable, Any
96    from .modelinfo import ModelInfo
97    from .details import CallDetails
98except ImportError:
99    pass
100# pylint: enable=unused-import
101
102if "SAS_COMPILER" in os.environ:
103    COMPILER = os.environ["SAS_COMPILER"]
104elif os.name == 'nt':
105    if tinycc is not None:
106        COMPILER = "tinycc"
107    elif "VCINSTALLDIR" in os.environ:
108        # If vcvarsall.bat has been called, then VCINSTALLDIR is in the environment
109        # and we can use the MSVC compiler.  Otherwise, if tinycc is available
110        # the use it.  Otherwise, hope that mingw is available.
111        COMPILER = "msvc"
112    else:
113        COMPILER = "mingw"
114else:
115    COMPILER = "unix"
116
117ARCH = "" if ct.sizeof(ct.c_void_p) > 4 else "x86"  # 4 byte pointers on x86
118if COMPILER == "unix":
119    # Generic unix compile
120    # On mac users will need the X code command line tools installed
121    #COMPILE = "gcc-mp-4.7 -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm -lgomp"
122    CC = "cc -shared -fPIC -std=c99 -O2 -Wall".split()
123    # add openmp support if not running on a mac
124    if sys.platform != "darwin":
125        # OpenMP seems to be broken on gcc 5.4.0 (ubuntu 16.04.9)
126        # Shut it off for all unix until we can investigate.
127        #CC.append("-fopenmp")
128        pass
129    def compile_command(source, output):
130        """unix compiler command"""
131        return CC + [source, "-o", output, "-lm"]
132elif COMPILER == "msvc":
133    # Call vcvarsall.bat before compiling to set path, headers, libs, etc.
134    # MSVC compiler is available, so use it.  OpenMP requires a copy of
135    # vcomp90.dll on the path.  One may be found here:
136    #       C:/Windows/winsxs/x86_microsoft.vc90.openmp*/vcomp90.dll
137    # Copy this to the python directory and uncomment the OpenMP COMPILE
138    # TODO: remove intermediate OBJ file created in the directory
139    # TODO: maybe don't use randomized name for the c file
140    # TODO: maybe ask distutils to find MSVC
141    CC = "cl /nologo /Ox /MD /W3 /GS- /DNDEBUG".split()
142    if "SAS_OPENMP" in os.environ:
143        CC.append("/openmp")
144    LN = "/link /DLL /INCREMENTAL:NO /MANIFEST".split()
145    def compile_command(source, output):
146        """MSVC compiler command"""
147        return CC + ["/Tp%s"%source] + LN + ["/OUT:%s"%output]
148elif COMPILER == "tinycc":
149    # TinyCC compiler.
150    CC = [tinycc.TCC] + "-shared -rdynamic -Wall".split()
151    def compile_command(source, output):
152        """tinycc compiler command"""
153        return CC + [source, "-o", output]
154elif COMPILER == "mingw":
155    # MinGW compiler.
156    CC = "gcc -shared -std=c99 -O2 -Wall".split()
157    if "SAS_OPENMP" in os.environ:
158        CC.append("-fopenmp")
159    def compile_command(source, output):
160        """mingw compiler command"""
161        return CC + [source, "-o", output, "-lm"]
162
163# Assume the default location of module DLLs is in .sasmodels/compiled_models.
164DLL_PATH = os.path.join(os.path.expanduser("~"), ".sasmodels", "compiled_models")
165
166ALLOW_SINGLE_PRECISION_DLLS = True
167
168def compile(source, output):
169    # type: (str, str) -> None
170    """
171    Compile *source* producing *output*.
172
173    Raises RuntimeError if the compile failed or the output wasn't produced.
174    """
175    command = compile_command(source=source, output=output)
176    command_str = " ".join('"%s"'%p if ' ' in p else p for p in command)
177    logging.info(command_str)
178    try:
179        # need shell=True on windows to keep console box from popping up
180        shell = (os.name == 'nt')
181        subprocess.check_output(command, shell=shell, stderr=subprocess.STDOUT)
182    except subprocess.CalledProcessError as exc:
183        raise RuntimeError("compile failed.\n%s\n%s"%(command_str, exc.output))
184    if not os.path.exists(output):
185        raise RuntimeError("compile failed.  File is in %r"%source)
186
187def dll_name(model_info, dtype):
188    # type: (ModelInfo, np.dtype) ->  str
189    """
190    Name of the dll containing the model.  This is the base file name without
191    any path or extension, with a form such as 'sas_sphere32'.
192    """
193    bits = 8*dtype.itemsize
194    basename = "sas%d_%s"%(bits, model_info.id)
195    basename += ARCH + ".so"
196
197    # Hack to find precompiled dlls
198    path = joinpath(generate.DATA_PATH, '..', 'compiled_models', basename)
199    if os.path.exists(path):
200        return path
201
202    return joinpath(DLL_PATH, basename)
203
204
205def dll_path(model_info, dtype):
206    # type: (ModelInfo, np.dtype) -> str
207    """
208    Complete path to the dll for the model.  Note that the dll may not
209    exist yet if it hasn't been compiled.
210    """
211    return os.path.join(DLL_PATH, dll_name(model_info, dtype))
212
213
214def make_dll(source, model_info, dtype=F64):
215    # type: (str, ModelInfo, np.dtype) -> str
216    """
217    Returns the path to the compiled model defined by *kernel_module*.
218
219    If the model has not been compiled, or if the source file(s) are newer
220    than the dll, then *make_dll* will compile the model before returning.
221    This routine does not load the resulting dll.
222
223    *dtype* is a numpy floating point precision specifier indicating whether
224    the model should be single, double or long double precision.  The default
225    is double precision, *np.dtype('d')*.
226
227    Set *sasmodels.ALLOW_SINGLE_PRECISION_DLLS* to False if single precision
228    models are not allowed as DLLs.
229
230    Set *sasmodels.kerneldll.DLL_PATH* to the compiled dll output path.
231    The default is in ~/.sasmodels/compiled_models.
232    """
233    if dtype == F16:
234        raise ValueError("16 bit floats not supported")
235    if dtype == F32 and not ALLOW_SINGLE_PRECISION_DLLS:
236        dtype = F64  # Force 64-bit dll
237    # Note: dtype may be F128 for long double precision
238
239    dll = dll_path(model_info, dtype)
240
241    if not os.path.exists(dll):
242        need_recompile = True
243    else:
244        dll_time = os.path.getmtime(dll)
245        newest_source = generate.dll_timestamp(model_info)
246        need_recompile = dll_time < newest_source
247    if need_recompile:
248        # Make sure the DLL path exists
249        if not os.path.exists(DLL_PATH):
250            os.makedirs(DLL_PATH)
251        basename = splitext(os.path.basename(dll))[0] + "_"
252        system_fd, filename = tempfile.mkstemp(suffix=".c", prefix=basename)
253        source = generate.convert_type(source, dtype)
254        with os.fdopen(system_fd, "w") as file_handle:
255            file_handle.write(source)
256        compile(source=filename, output=dll)
257        # comment the following to keep the generated c file
258        # Note: if there is a syntax error then compile raises an error
259        # and the source file will not be deleted.
260        #os.unlink(filename)
261        print("saving compiled file in %r"%filename)
262    return dll
263
264
265def load_dll(source, model_info, dtype=F64):
266    # type: (str, ModelInfo, np.dtype) -> "DllModel"
267    """
268    Create and load a dll corresponding to the source, info pair returned
269    from :func:`sasmodels.generate.make` compiled for the target precision.
270
271    See :func:`make_dll` for details on controlling the dll path and the
272    allowed floating point precision.
273    """
274    filename = make_dll(source, model_info, dtype=dtype)
275    return DllModel(filename, model_info, dtype=dtype)
276
277
278class DllModel(KernelModel):
279    """
280    ctypes wrapper for a single model.
281
282    *source* and *model_info* are the model source and interface as returned
283    from :func:`gen.make`.
284
285    *dtype* is the desired model precision.  Any numpy dtype for single
286    or double precision floats will do, such as 'f', 'float32' or 'single'
287    for single and 'd', 'float64' or 'double' for double.  Double precision
288    is an optional extension which may not be available on all devices.
289
290    Call :meth:`release` when done with the kernel.
291    """
292    def __init__(self, dllpath, model_info, dtype=generate.F32):
293        # type: (str, ModelInfo, np.dtype) -> None
294        self.info = model_info
295        self.dllpath = dllpath
296        self._dll = None  # type: ct.CDLL
297        self._kernels = None # type: List[Callable, Callable]
298        self.dtype = np.dtype(dtype)
299
300    def _load_dll(self):
301        # type: () -> None
302        try:
303            self._dll = ct.CDLL(self.dllpath)
304        except:
305            annotate_exception("while loading "+self.dllpath)
306            raise
307
308        float_type = (ct.c_float if self.dtype == generate.F32
309                      else ct.c_double if self.dtype == generate.F64
310                      else ct.c_longdouble)
311
312        # int, int, int, int*, double*, double*, double*, double*, double
313        argtypes = [ct.c_int32]*3 + [ct.c_void_p]*4 + [float_type]
314        names = [generate.kernel_name(self.info, variant)
315                 for variant in ("Iq", "Iqxy", "Imagnetic")]
316        self._kernels = [self._dll[name] for name in names]
317        for k in self._kernels:
318            k.argtypes = argtypes
319
320    def __getstate__(self):
321        # type: () -> Tuple[ModelInfo, str]
322        return self.info, self.dllpath
323
324    def __setstate__(self, state):
325        # type: (Tuple[ModelInfo, str]) -> None
326        self.info, self.dllpath = state
327        self._dll = None
328
329    def make_kernel(self, q_vectors):
330        # type: (List[np.ndarray]) -> DllKernel
331        q_input = PyInput(q_vectors, self.dtype)
332        # Note: pickle not supported for DllKernel
333        if self._dll is None:
334            self._load_dll()
335        is_2d = len(q_vectors) == 2
336        kernel = self._kernels[1:3] if is_2d else [self._kernels[0]]*2
337        return DllKernel(kernel, self.info, q_input)
338
339    def release(self):
340        # type: () -> None
341        """
342        Release any resources associated with the model.
343        """
344        dll_handle = self._dll._handle
345        if os.name == 'nt':
346            ct.windll.kernel32.FreeLibrary(dll_handle)
347        else:
348            _ct.dlclose(dll_handle)
349        del self._dll
350        self._dll = None
351
352class DllKernel(Kernel):
353    """
354    Callable SAS kernel.
355
356    *kernel* is the c function to call.
357
358    *model_info* is the module information
359
360    *q_input* is the DllInput q vectors at which the kernel should be
361    evaluated.
362
363    The resulting call method takes the *pars*, a list of values for
364    the fixed parameters to the kernel, and *pd_pars*, a list of (value, weight)
365    vectors for the polydisperse parameters.  *cutoff* determines the
366    integration limits: any points with combined weight less than *cutoff*
367    will not be calculated.
368
369    Call :meth:`release` when done with the kernel instance.
370    """
371    def __init__(self, kernel, model_info, q_input):
372        # type: (Callable[[], np.ndarray], ModelInfo, PyInput) -> None
373        #,model_info,q_input)
374        self.kernel = kernel
375        self.info = model_info
376        self.q_input = q_input
377        self.dtype = q_input.dtype
378        self.dim = '2d' if q_input.is_2d else '1d'
379        self.result = np.empty(2*q_input.nq+2, q_input.dtype)
380        self.real = (np.float32 if self.q_input.dtype == generate.F32
381                     else np.float64 if self.q_input.dtype == generate.F64
382                     else np.float128)
383
384    def Iq(self, call_details, values, cutoff, magnetic):
385        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray
386        self._call_kernel(call_details, values, cutoff, magnetic)
387        #print("returned",self.q_input.q, self.result)
388        pd_norm = self.result[self.q_input.nq]
389        scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0)
390        background = values[1]
391        #print("scale",scale,background)
392        return scale*self.result[:self.q_input.nq] + background
393    __call__ = Iq
394
395    def beta(self, call_details, values, cutoff, magnetic):
396        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray
397        self._call_kernel(call_details, values, cutoff, magnetic)
398        w_norm = self.result[2*self.q_input.nq + 1]
399        pd_norm = self.result[self.q_input.nq]
400        if w_norm == 0.:
401            w_norm = 1.
402        F2 = self.result[:self.q_input.nq]/w_norm
403        F1 = self.result[self.q_input.nq+1:2*self.q_input.nq+1]/w_norm
404        volume_avg = pd_norm/w_norm
405        return F1, F2, volume_avg
406
407    def _call_kernel(self, call_details, values, cutoff, magnetic):
408        kernel = self.kernel[1 if magnetic else 0]
409        args = [
410            self.q_input.nq, # nq
411            None, # pd_start
412            None, # pd_stop pd_stride[MAX_PD]
413            call_details.buffer.ctypes.data, # problem
414            values.ctypes.data,  # pars
415            self.q_input.q.ctypes.data, # q
416            self.result.ctypes.data,   # results
417            self.real(cutoff), # cutoff
418        ]
419        #print(self.beta_result.ctypes.data)
420#        print("Calling DLL line 397")
421#        print("values", values)
422        #call_details.show(values)
423        step = 100
424        for start in range(0, call_details.num_eval, step):
425            stop = min(start + step, call_details.num_eval)
426            args[1:3] = [start, stop]
427            kernel(*args) # type: ignore
428
429    def release(self):
430        # type: () -> None
431        """
432        Release any resources associated with the kernel.
433        """
434        self.q_input.release()
Note: See TracBrowser for help on using the repository browser.