1 | """ |
---|
2 | Core model handling routines. |
---|
3 | """ |
---|
4 | #__all__ = [ |
---|
5 | # "list_models", "load_model_info", "precompile_dll", |
---|
6 | # "build_model", "make_kernel", "call_kernel", "call_ER_VR", |
---|
7 | # ] |
---|
8 | |
---|
9 | import os |
---|
10 | from os.path import basename, dirname, join as joinpath, splitext |
---|
11 | from glob import glob |
---|
12 | |
---|
13 | import numpy as np |
---|
14 | |
---|
15 | from . import models |
---|
16 | from . import weights |
---|
17 | from . import generate |
---|
18 | # TODO: remove circular references between product and core |
---|
19 | # product uses call_ER/call_VR, core uses make_product_info/ProductModel |
---|
20 | #from . import product |
---|
21 | from . import mixture |
---|
22 | from . import kernelpy |
---|
23 | from . import kerneldll |
---|
24 | try: |
---|
25 | from . import kernelcl |
---|
26 | HAVE_OPENCL = True |
---|
27 | except: |
---|
28 | HAVE_OPENCL = False |
---|
29 | |
---|
30 | try: |
---|
31 | np.meshgrid([]) |
---|
32 | meshgrid = np.meshgrid |
---|
33 | except ValueError: |
---|
34 | # CRUFT: np.meshgrid requires multiple vectors |
---|
35 | def meshgrid(*args): |
---|
36 | if len(args) > 1: |
---|
37 | return np.meshgrid(*args) |
---|
38 | else: |
---|
39 | return [np.asarray(v) for v in args] |
---|
40 | |
---|
41 | # TODO: refactor composite model support |
---|
42 | # The current load_model_info/build_model does not reuse existing model |
---|
43 | # definitions when loading a composite model, instead reloading and |
---|
44 | # rebuilding the kernel for each component model in the expression. This |
---|
45 | # is fine in a scripting environment where the model is built when the script |
---|
46 | # starts and is thrown away when the script ends, but may not be the best |
---|
47 | # solution in a long-lived application. This affects the following functions: |
---|
48 | # |
---|
49 | # load_model |
---|
50 | # load_model_info |
---|
51 | # build_model |
---|
52 | |
---|
53 | def list_models(): |
---|
54 | """ |
---|
55 | Return the list of available models on the model path. |
---|
56 | """ |
---|
57 | root = dirname(__file__) |
---|
58 | available_models = [] |
---|
59 | files = sorted(glob(joinpath(root, 'models', "[a-zA-Z]*.py"))) |
---|
60 | if not files and os.name=='nt': |
---|
61 | # Look inside library.zip on windows, |
---|
62 | # being careful as the library stores only .pyc. |
---|
63 | import zipfile |
---|
64 | location = root[:root.find('.zip')+4] |
---|
65 | zf = zipfile.ZipFile(location,'r') |
---|
66 | for filepath in zf.namelist(): |
---|
67 | # directory structure in library.zip uses "/" |
---|
68 | models_loc = "sasmodels/models" |
---|
69 | if models_loc in filepath: |
---|
70 | base = basename(filepath)[:-4] |
---|
71 | # Careful with non-models. |
---|
72 | if base[0] != "_": |
---|
73 | available_models.append(base) |
---|
74 | else: |
---|
75 | available_models = [basename(f)[:-3] for f in files] |
---|
76 | return available_models |
---|
77 | |
---|
78 | def isstr(s): |
---|
79 | """ |
---|
80 | Return True if *s* is a string-like object. |
---|
81 | """ |
---|
82 | try: s + '' |
---|
83 | except: return False |
---|
84 | return True |
---|
85 | |
---|
86 | def load_model(model_name, **kw): |
---|
87 | """ |
---|
88 | Load model info and build model. |
---|
89 | """ |
---|
90 | return build_model(load_model_info(model_name), **kw) |
---|
91 | |
---|
92 | |
---|
93 | def load_model_info(model_name): |
---|
94 | """ |
---|
95 | Load a model definition given the model name. |
---|
96 | |
---|
97 | This returns a handle to the module defining the model. This can be |
---|
98 | used with functions in generate to build the docs or extract model info. |
---|
99 | """ |
---|
100 | parts = model_name.split('+') |
---|
101 | if len(parts) > 1: |
---|
102 | model_info_list = [load_model_info(p) for p in parts] |
---|
103 | return mixture.make_mixture_info(model_info_list) |
---|
104 | |
---|
105 | parts = model_name.split('*') |
---|
106 | if len(parts) > 1: |
---|
107 | from . import product |
---|
108 | # Note: currently have circular reference |
---|
109 | if len(parts) > 2: |
---|
110 | raise ValueError("use P*S to apply structure factor S to model P") |
---|
111 | P_info, Q_info = [load_model_info(p) for p in parts] |
---|
112 | return product.make_product_info(P_info, Q_info) |
---|
113 | |
---|
114 | kernel_module = generate.load_kernel_module(model_name) |
---|
115 | return generate.make_model_info(kernel_module) |
---|
116 | |
---|
117 | |
---|
118 | def build_model(model_info, dtype=None, platform="ocl"): |
---|
119 | """ |
---|
120 | Prepare the model for the default execution platform. |
---|
121 | |
---|
122 | This will return an OpenCL model, a DLL model or a python model depending |
---|
123 | on the model and the computing platform. |
---|
124 | |
---|
125 | *model_info* is the model definition structure returned from |
---|
126 | :func:`load_model_info`. |
---|
127 | |
---|
128 | *dtype* indicates whether the model should use single or double precision |
---|
129 | for the calculation. Any valid numpy single or double precision identifier |
---|
130 | is valid, such as 'single', 'f', 'f32', or np.float32 for single, or |
---|
131 | 'double', 'd', 'f64' and np.float64 for double. If *None*, then use |
---|
132 | 'single' unless the model defines single=False. |
---|
133 | |
---|
134 | *platform* should be "dll" to force the dll to be used for C models, |
---|
135 | otherwise it uses the default "ocl". |
---|
136 | """ |
---|
137 | composition = model_info.get('composition', None) |
---|
138 | if composition is not None: |
---|
139 | composition_type, parts = composition |
---|
140 | models = [build_model(p, dtype=dtype, platform=platform) for p in parts] |
---|
141 | if composition_type == 'mixture': |
---|
142 | return mixture.MixtureModel(model_info, models) |
---|
143 | elif composition_type == 'product': |
---|
144 | from . import product |
---|
145 | P, S = models |
---|
146 | return product.ProductModel(model_info, P, S) |
---|
147 | else: |
---|
148 | raise ValueError('unknown mixture type %s'%composition_type) |
---|
149 | |
---|
150 | ## for debugging: |
---|
151 | ## 1. uncomment open().write so that the source will be saved next time |
---|
152 | ## 2. run "python -m sasmodels.direct_model $MODELNAME" to save the source |
---|
153 | ## 3. recomment the open.write() and uncomment open().read() |
---|
154 | ## 4. rerun "python -m sasmodels.direct_model $MODELNAME" |
---|
155 | ## 5. uncomment open().read() so that source will be regenerated from model |
---|
156 | # open(model_info['name']+'.c','w').write(source) |
---|
157 | # source = open(model_info['name']+'.cl','r').read() |
---|
158 | if callable(model_info.get('Iq', None)): |
---|
159 | return kernelpy.PyModel(model_info) |
---|
160 | source = generate.make_source(model_info) |
---|
161 | default_dtype = 'single' if model_info['single'] else 'double' |
---|
162 | ocl_dtype = default_dtype if dtype is None else dtype |
---|
163 | dll_dtype = 'double' if dtype is None else dtype |
---|
164 | if (platform == "dll" |
---|
165 | or not HAVE_OPENCL |
---|
166 | or not kernelcl.environment().has_type(ocl_dtype)): |
---|
167 | return kerneldll.load_dll(source, model_info, dll_dtype) |
---|
168 | else: |
---|
169 | return kernelcl.GpuModel(source, model_info, ocl_dtype) |
---|
170 | |
---|
171 | def precompile_dlls(path, dtype="double"): |
---|
172 | """ |
---|
173 | Precompile the dlls for all builtin models, returning a list of dll paths. |
---|
174 | |
---|
175 | *path* is the directory in which to save the dlls. It will be created if |
---|
176 | it does not already exist. |
---|
177 | |
---|
178 | This can be used when build the windows distribution of sasmodels |
---|
179 | which may be missing the OpenCL driver and the dll compiler. |
---|
180 | """ |
---|
181 | if not os.path.exists(path): |
---|
182 | os.makedirs(path) |
---|
183 | compiled_dlls = [] |
---|
184 | for model_name in list_models(): |
---|
185 | model_info = load_model_info(model_name) |
---|
186 | source = generate.make_source(model_info) |
---|
187 | if source: |
---|
188 | old_path = kerneldll.DLL_PATH |
---|
189 | try: |
---|
190 | kerneldll.DLL_PATH = path |
---|
191 | dll = kerneldll.make_dll(source, model_info, dtype=dtype) |
---|
192 | finally: |
---|
193 | kerneldll.DLL_PATH = old_path |
---|
194 | compiled_dlls.append(dll) |
---|
195 | return compiled_dlls |
---|
196 | |
---|
197 | def get_weights(model_info, pars, name): |
---|
198 | """ |
---|
199 | Generate the distribution for parameter *name* given the parameter values |
---|
200 | in *pars*. |
---|
201 | |
---|
202 | Uses "name", "name_pd", "name_pd_type", "name_pd_n", "name_pd_sigma" |
---|
203 | from the *pars* dictionary for parameter value and parameter dispersion. |
---|
204 | """ |
---|
205 | relative = name in model_info['partype']['pd-rel'] |
---|
206 | limits = model_info['limits'][name] |
---|
207 | disperser = pars.get(name+'_pd_type', 'gaussian') |
---|
208 | value = pars.get(name, model_info['defaults'][name]) |
---|
209 | npts = pars.get(name+'_pd_n', 0) |
---|
210 | width = pars.get(name+'_pd', 0.0) |
---|
211 | nsigma = pars.get(name+'_pd_nsigma', 3.0) |
---|
212 | value, weight = weights.get_weights( |
---|
213 | disperser, npts, width, nsigma, value, limits, relative) |
---|
214 | return value, weight / np.sum(weight) |
---|
215 | |
---|
216 | def dispersion_mesh(pars): |
---|
217 | """ |
---|
218 | Create a mesh grid of dispersion parameters and weights. |
---|
219 | |
---|
220 | Returns [p1,p2,...],w where pj is a vector of values for parameter j |
---|
221 | and w is a vector containing the products for weights for each |
---|
222 | parameter set in the vector. |
---|
223 | """ |
---|
224 | value, weight = zip(*pars) |
---|
225 | value = [v.flatten() for v in meshgrid(*value)] |
---|
226 | weight = np.vstack([v.flatten() for v in meshgrid(*weight)]) |
---|
227 | weight = np.prod(weight, axis=0) |
---|
228 | return value, weight |
---|
229 | |
---|
230 | def call_kernel(kernel, pars, cutoff=0, mono=False): |
---|
231 | """ |
---|
232 | Call *kernel* returned from *model.make_kernel* with parameters *pars*. |
---|
233 | |
---|
234 | *cutoff* is the limiting value for the product of dispersion weights used |
---|
235 | to perform the multidimensional dispersion calculation more quickly at a |
---|
236 | slight cost to accuracy. The default value of *cutoff=0* integrates over |
---|
237 | the entire dispersion cube. Using *cutoff=1e-5* can be 50% faster, but |
---|
238 | with an error of about 1%, which is usually less than the measurement |
---|
239 | uncertainty. |
---|
240 | |
---|
241 | *mono* is True if polydispersity should be set to none on all parameters. |
---|
242 | """ |
---|
243 | fixed_pars = [pars.get(name, kernel.info['defaults'][name]) |
---|
244 | for name in kernel.fixed_pars] |
---|
245 | if mono: |
---|
246 | pd_pars = [( np.array([pars[name]]), np.array([1.0]) ) |
---|
247 | for name in kernel.pd_pars] |
---|
248 | else: |
---|
249 | pd_pars = [get_weights(kernel.info, pars, name) for name in kernel.pd_pars] |
---|
250 | return kernel(fixed_pars, pd_pars, cutoff=cutoff) |
---|
251 | |
---|
252 | def call_ER_VR(model_info, vol_pars): |
---|
253 | """ |
---|
254 | Return effect radius and volume ratio for the model. |
---|
255 | |
---|
256 | *info* is either *kernel.info* for *kernel=make_kernel(model,q)* |
---|
257 | or *model.info*. |
---|
258 | |
---|
259 | *pars* are the parameters as expected by :func:`call_kernel`. |
---|
260 | """ |
---|
261 | ER = model_info.get('ER', None) |
---|
262 | VR = model_info.get('VR', None) |
---|
263 | value, weight = dispersion_mesh(vol_pars) |
---|
264 | |
---|
265 | individual_radii = ER(*value) if ER else 1.0 |
---|
266 | whole, part = VR(*value) if VR else (1.0, 1.0) |
---|
267 | |
---|
268 | effect_radius = np.sum(weight*individual_radii) / np.sum(weight) |
---|
269 | volume_ratio = np.sum(weight*part)/np.sum(weight*whole) |
---|
270 | return effect_radius, volume_ratio |
---|
271 | |
---|
272 | |
---|
273 | def call_ER(model_info, values): |
---|
274 | """ |
---|
275 | Call the model ER function using *values*. *model_info* is either |
---|
276 | *model.info* if you have a loaded model, or *kernel.info* if you |
---|
277 | have a model kernel prepared for evaluation. |
---|
278 | """ |
---|
279 | ER = model_info.get('ER', None) |
---|
280 | if ER is None: |
---|
281 | return 1.0 |
---|
282 | else: |
---|
283 | vol_pars = [get_weights(model_info, values, name) |
---|
284 | for name in model_info['partype']['volume']] |
---|
285 | value, weight = dispersion_mesh(vol_pars) |
---|
286 | individual_radii = ER(*value) |
---|
287 | #print(values[0].shape, weights.shape, fv.shape) |
---|
288 | return np.sum(weight*individual_radii) / np.sum(weight) |
---|
289 | |
---|
290 | def call_VR(model_info, values): |
---|
291 | """ |
---|
292 | Call the model VR function using *pars*. |
---|
293 | *info* is either *model.info* if you have a loaded model, or *kernel.info* |
---|
294 | if you have a model kernel prepared for evaluation. |
---|
295 | """ |
---|
296 | VR = model_info.get('VR', None) |
---|
297 | if VR is None: |
---|
298 | return 1.0 |
---|
299 | else: |
---|
300 | vol_pars = [get_weights(model_info, values, name) |
---|
301 | for name in model_info['partype']['volume']] |
---|
302 | value, weight = dispersion_mesh(vol_pars) |
---|
303 | whole, part = VR(*value) |
---|
304 | return np.sum(weight*part)/np.sum(weight*whole) |
---|
305 | |
---|
306 | # TODO: remove call_ER, call_VR |
---|
307 | |
---|