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