[aa4946b] | 1 | """ |
---|
| 2 | Core model handling routines. |
---|
| 3 | """ |
---|
| 4 | __all__ = ["list_models", "load_model_definition", "precompile_dll", |
---|
| 5 | "load_model", "make_kernel", "call_kernel", "call_ER", "call_VR" ] |
---|
[f734e7d] | 6 | |
---|
| 7 | from os.path import basename, dirname, join as joinpath |
---|
| 8 | from glob import glob |
---|
| 9 | |
---|
| 10 | import numpy as np |
---|
| 11 | |
---|
| 12 | from . import models |
---|
| 13 | from . import weights |
---|
[aa4946b] | 14 | from . import generate |
---|
[f734e7d] | 15 | |
---|
[aa4946b] | 16 | from . import kernelpy |
---|
| 17 | from . import kerneldll |
---|
[f734e7d] | 18 | try: |
---|
[aa4946b] | 19 | from . import kernelcl |
---|
| 20 | HAVE_OPENCL = True |
---|
[3c56da87] | 21 | except: |
---|
[aa4946b] | 22 | HAVE_OPENCL = False |
---|
| 23 | |
---|
[f734e7d] | 24 | |
---|
| 25 | def list_models(): |
---|
[aa4946b] | 26 | """ |
---|
| 27 | Return the list of available models on the model path. |
---|
| 28 | """ |
---|
[f734e7d] | 29 | root = dirname(__file__) |
---|
| 30 | files = sorted(glob(joinpath(root, 'models', "[a-zA-Z]*.py"))) |
---|
| 31 | available_models = [basename(f)[:-3] for f in files] |
---|
| 32 | return available_models |
---|
| 33 | |
---|
[aa4946b] | 34 | |
---|
[f734e7d] | 35 | def load_model_definition(model_name): |
---|
[aa4946b] | 36 | """ |
---|
| 37 | Load a model definition given the model name. |
---|
| 38 | """ |
---|
[f734e7d] | 39 | __import__('sasmodels.models.'+model_name) |
---|
| 40 | model_definition = getattr(models, model_name, None) |
---|
| 41 | return model_definition |
---|
| 42 | |
---|
[aa4946b] | 43 | |
---|
| 44 | def precompile_dll(model_name, dtype="double"): |
---|
| 45 | """ |
---|
| 46 | Precompile the dll for a model. |
---|
| 47 | |
---|
| 48 | Returns the path to the compiled model. |
---|
| 49 | |
---|
| 50 | This can be used when build the windows distribution of sasmodels |
---|
| 51 | (which may be missing the OpenCL driver and the dll compiler), or |
---|
| 52 | otherwise sharing models with windows users who do not have a compiler. |
---|
| 53 | |
---|
| 54 | See :func:`sasmodels.kerneldll.make_dll` for details on controlling the |
---|
| 55 | dll path and the allowed floating point precision. |
---|
| 56 | """ |
---|
| 57 | model_definition = load_model_definition(model_name) |
---|
| 58 | source, info = generate.make(model_definition) |
---|
| 59 | return kerneldll.make_dll(source, info, dtype=dtype) |
---|
| 60 | |
---|
| 61 | |
---|
[bcd3aa3] | 62 | def isstr(s): |
---|
| 63 | try: s + '' |
---|
| 64 | except: return False |
---|
| 65 | return True |
---|
| 66 | |
---|
[aa4946b] | 67 | def load_model(model_definition, dtype="single", platform="ocl"): |
---|
| 68 | """ |
---|
| 69 | Prepare the model for the default execution platform. |
---|
| 70 | |
---|
| 71 | This will return an OpenCL model, a DLL model or a python model depending |
---|
| 72 | on the model and the computing platform. |
---|
| 73 | |
---|
[bcd3aa3] | 74 | *model_definition* is the python module which defines the model. If the |
---|
| 75 | model name is given instead, then :func:`load_model_definition` will be |
---|
| 76 | called with the model name. |
---|
| 77 | |
---|
[aa4946b] | 78 | *dtype* indicates whether the model should use single or double precision |
---|
| 79 | for the calculation. Any valid numpy single or double precision identifier |
---|
| 80 | is valid, such as 'single', 'f', 'f32', or np.float32 for single, or |
---|
| 81 | 'double', 'd', 'f64' and np.float64 for double. |
---|
| 82 | |
---|
| 83 | *platform* should be "dll" to force the dll to be used for C models, |
---|
| 84 | otherwise it uses the default "ocl". |
---|
| 85 | """ |
---|
[bcd3aa3] | 86 | if isstr(model_definition): |
---|
| 87 | model_definition = load_model_definition(model_definition) |
---|
[aa4946b] | 88 | source, info = generate.make(model_definition) |
---|
| 89 | if callable(info.get('Iq', None)): |
---|
| 90 | return kernelpy.PyModel(info) |
---|
| 91 | |
---|
| 92 | ## for debugging: |
---|
| 93 | ## 1. uncomment open().write so that the source will be saved next time |
---|
| 94 | ## 2. run "python -m sasmodels.direct_model $MODELNAME" to save the source |
---|
| 95 | ## 3. recomment the open.write() and uncomment open().read() |
---|
| 96 | ## 4. rerun "python -m sasmodels.direct_model $MODELNAME" |
---|
| 97 | ## 5. uncomment open().read() so that source will be regenerated from model |
---|
| 98 | # open(info['name']+'.c','w').write(source) |
---|
| 99 | # source = open(info['name']+'.cl','r').read() |
---|
| 100 | |
---|
| 101 | dtype = np.dtype(dtype) |
---|
| 102 | if (platform=="dll" |
---|
| 103 | or not HAVE_OPENCL |
---|
| 104 | or (dtype == np.float64 and not kernelcl.environment().has_double)): |
---|
| 105 | return kerneldll.load_dll(source, info, dtype) |
---|
| 106 | else: |
---|
| 107 | return kernelcl.GpuModel(source, info, dtype) |
---|
| 108 | |
---|
[f734e7d] | 109 | def make_kernel(model, q_vectors): |
---|
| 110 | """ |
---|
| 111 | Return a computation kernel from the model definition and the q input. |
---|
| 112 | """ |
---|
[3c56da87] | 113 | model_input = model.make_input(q_vectors) |
---|
| 114 | return model(model_input) |
---|
[f734e7d] | 115 | |
---|
[3c56da87] | 116 | def get_weights(info, pars, name): |
---|
[f734e7d] | 117 | """ |
---|
| 118 | Generate the distribution for parameter *name* given the parameter values |
---|
| 119 | in *pars*. |
---|
| 120 | |
---|
[aa4946b] | 121 | Uses "name", "name_pd", "name_pd_type", "name_pd_n", "name_pd_sigma" |
---|
| 122 | from the *pars* dictionary for parameter value and parameter dispersion. |
---|
[f734e7d] | 123 | """ |
---|
[3c56da87] | 124 | relative = name in info['partype']['pd-rel'] |
---|
[7cf2cfd] | 125 | limits = info['limits'][name] |
---|
[f734e7d] | 126 | disperser = pars.get(name+'_pd_type', 'gaussian') |
---|
[3c56da87] | 127 | value = pars.get(name, info['defaults'][name]) |
---|
[f734e7d] | 128 | npts = pars.get(name+'_pd_n', 0) |
---|
| 129 | width = pars.get(name+'_pd', 0.0) |
---|
| 130 | nsigma = pars.get(name+'_pd_nsigma', 3.0) |
---|
[3c56da87] | 131 | value,weight = weights.get_weights( |
---|
[7cf2cfd] | 132 | disperser, npts, width, nsigma, value, limits, relative) |
---|
| 133 | return value, weight / np.sum(weight) |
---|
[f734e7d] | 134 | |
---|
[9890053] | 135 | def dispersion_mesh(pars): |
---|
| 136 | """ |
---|
| 137 | Create a mesh grid of dispersion parameters and weights. |
---|
| 138 | |
---|
| 139 | Returns [p1,p2,...],w where pj is a vector of values for parameter j |
---|
| 140 | and w is a vector containing the products for weights for each |
---|
| 141 | parameter set in the vector. |
---|
| 142 | """ |
---|
[3c56da87] | 143 | value, weight = zip(*pars) |
---|
| 144 | if len(value) > 1: |
---|
| 145 | value = [v.flatten() for v in np.meshgrid(*value)] |
---|
| 146 | weight = np.vstack([v.flatten() for v in np.meshgrid(*weight)]) |
---|
| 147 | weight = np.prod(weight, axis=0) |
---|
| 148 | return value, weight |
---|
[9890053] | 149 | |
---|
[aa4946b] | 150 | def call_kernel(kernel, pars, cutoff=0): |
---|
| 151 | """ |
---|
| 152 | Call *kernel* returned from :func:`make_kernel` with parameters *pars*. |
---|
| 153 | |
---|
| 154 | *cutoff* is the limiting value for the product of dispersion weights used |
---|
| 155 | to perform the multidimensional dispersion calculation more quickly at a |
---|
| 156 | slight cost to accuracy. The default value of *cutoff=0* integrates over |
---|
| 157 | the entire dispersion cube. Using *cutoff=1e-5* can be 50% faster, but |
---|
| 158 | with an error of about 1%, which is usually less than the measurement |
---|
| 159 | uncertainty. |
---|
| 160 | """ |
---|
[f734e7d] | 161 | fixed_pars = [pars.get(name, kernel.info['defaults'][name]) |
---|
| 162 | for name in kernel.fixed_pars] |
---|
[3c56da87] | 163 | pd_pars = [get_weights(kernel.info, pars, name) for name in kernel.pd_pars] |
---|
[f734e7d] | 164 | return kernel(fixed_pars, pd_pars, cutoff=cutoff) |
---|
| 165 | |
---|
[3c56da87] | 166 | def call_ER(info, pars): |
---|
[aa4946b] | 167 | """ |
---|
| 168 | Call the model ER function using *pars*. |
---|
| 169 | |
---|
| 170 | *info* is either *model.info* if you have a loaded model, or *kernel.info* |
---|
| 171 | if you have a model kernel prepared for evaluation. |
---|
| 172 | """ |
---|
[3c56da87] | 173 | ER = info.get('ER', None) |
---|
[9890053] | 174 | if ER is None: |
---|
| 175 | return 1.0 |
---|
| 176 | else: |
---|
[3c56da87] | 177 | vol_pars = [get_weights(info, pars, name) |
---|
| 178 | for name in info['partype']['volume']] |
---|
| 179 | value, weight = dispersion_mesh(vol_pars) |
---|
| 180 | individual_radii = ER(*value) |
---|
[9404dd3] | 181 | #print(values[0].shape, weights.shape, fv.shape) |
---|
[3c56da87] | 182 | return np.sum(weight*individual_radii) / np.sum(weight) |
---|
[9890053] | 183 | |
---|
[3c56da87] | 184 | def call_VR(info, pars): |
---|
[aa4946b] | 185 | """ |
---|
| 186 | Call the model VR function using *pars*. |
---|
| 187 | |
---|
| 188 | *info* is either *model.info* if you have a loaded model, or *kernel.info* |
---|
| 189 | if you have a model kernel prepared for evaluation. |
---|
| 190 | """ |
---|
[3c56da87] | 191 | VR = info.get('VR', None) |
---|
[9890053] | 192 | if VR is None: |
---|
| 193 | return 1.0 |
---|
| 194 | else: |
---|
[3c56da87] | 195 | vol_pars = [get_weights(info, pars, name) |
---|
| 196 | for name in info['partype']['volume']] |
---|
| 197 | value, weight = dispersion_mesh(vol_pars) |
---|
| 198 | whole,part = VR(*value) |
---|
| 199 | return np.sum(weight*part)/np.sum(weight*whole) |
---|
[9890053] | 200 | |
---|