Changes in sasmodels/product.py [99658f6:b171acd] in sasmodels
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/product.py
r99658f6 rb171acd 13 13 from __future__ import print_function, division 14 14 15 from collections import OrderedDict16 17 15 from copy import copy 18 16 import numpy as np # type: ignore 19 17 20 from .modelinfo import ParameterTable, ModelInfo , Parameter, parse_parameter18 from .modelinfo import ParameterTable, ModelInfo 21 19 from .kernel import KernelModel, Kernel 22 20 from .details import make_details, dispersion_mesh … … 24 22 # pylint: disable=unused-import 25 23 try: 26 from typing import Tuple , Callable, Union24 from typing import Tuple 27 25 except ImportError: 28 26 pass … … 37 35 #] 38 36 39 STRUCTURE_MODE_ID = "structure_factor_mode" 40 RADIUS_MODE_ID = "radius_effective_mode" 41 RADIUS_ID = "radius_effective" 42 VOLFRAC_ID = "volfraction" 43 def make_extra_pars(p_info): 44 pars = [] 45 if p_info.have_Fq: 46 par = parse_parameter( 47 STRUCTURE_MODE_ID, 48 "", 49 0, 50 [["P*S","P*(1+beta*(S-1))"]], 51 "", 52 "Structure factor calculation") 53 pars.append(par) 54 if p_info.effective_radius_type is not None: 55 par = parse_parameter( 56 RADIUS_MODE_ID, 57 "", 58 1, 59 [["unconstrained"] + p_info.effective_radius_type], 60 "", 61 "Effective radius calculation") 62 pars.append(par) 63 return pars 64 37 ER_ID = "radius_effective" 38 VF_ID = "volfraction" 39 40 # TODO: core_shell_sphere model has suppressed the volume ratio calculation 41 # revert it after making VR and ER available at run time as constraints. 65 42 def make_product_info(p_info, s_info): 66 43 # type: (ModelInfo, ModelInfo) -> ModelInfo … … 73 50 # have any magnetic parameters 74 51 if not len(s_info.parameters.kernel_parameters) >= 2: 75 raise TypeError("S needs {} and {} as its first parameters".format( RADIUS_ID, VOLFRAC_ID))76 if not s_info.parameters.kernel_parameters[0].id == RADIUS_ID:77 raise TypeError("S needs {} as first parameter".format( RADIUS_ID))78 if not s_info.parameters.kernel_parameters[1].id == V OLFRAC_ID:79 raise TypeError("S needs {} as second parameter".format(V OLFRAC_ID))52 raise TypeError("S needs {} and {} as its first parameters".format(ER_ID, VF_ID)) 53 if not s_info.parameters.kernel_parameters[0].id == ER_ID: 54 raise TypeError("S needs {} as first parameter".format(ER_ID)) 55 if not s_info.parameters.kernel_parameters[1].id == VF_ID: 56 raise TypeError("S needs {} as second parameter".format(VF_ID)) 80 57 if not s_info.parameters.magnetism_index == []: 81 58 raise TypeError("S should not have SLD parameters") … … 83 60 s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters 84 61 85 # Create list of parameters for the combined model. If there 62 # Create list of parameters for the combined model. Skip the first 63 # parameter of S, which we verified above is effective radius. If there 86 64 # are any names in P that overlap with those in S, modify the name in S 87 65 # to distinguish it. 88 66 p_set = set(p.id for p in p_pars.kernel_parameters) 89 67 s_list = [(_tag_parameter(par) if par.id in p_set else par) 90 for par in s_pars.kernel_parameters ]68 for par in s_pars.kernel_parameters[1:]] 91 69 # Check if still a collision after renaming. This could happen if for 92 70 # example S has volfrac and P has both volfrac and volfrac_S. … … 94 72 raise TypeError("name collision: P has P.name and P.name_S while S has S.name") 95 73 96 # make sure effective radius is not a polydisperse parameter in product97 s_list[0] = copy(s_list[0])98 s_list[0].polydisperse = False99 100 74 translate_name = dict((old.id, new.id) for old, new 101 in zip(s_pars.kernel_parameters , s_list))102 combined_pars = p_pars.kernel_parameters + s_list + make_extra_pars(p_info)75 in zip(s_pars.kernel_parameters[1:], s_list)) 76 combined_pars = p_pars.kernel_parameters + s_list 103 77 parameters = ParameterTable(combined_pars) 104 78 parameters.max_pd = p_pars.max_pd + s_pars.max_pd 105 79 def random(): 106 80 combined_pars = p_info.random() 107 s_names = set(par.id for par in s_pars.kernel_parameters )81 s_names = set(par.id for par in s_pars.kernel_parameters[1:]) 108 82 combined_pars.update((translate_name[k], v) 109 83 for k, v in s_info.random().items() … … 126 100 #model_info.tests = [] 127 101 #model_info.source = [] 102 # Iq, Iqxy, form_volume, ER, VR and sesans 128 103 # Remember the component info blocks so we can build the model 129 104 model_info.composition = ('product', [p_info, s_info]) 130 model_info.control = p_info.control131 105 model_info.hidden = p_info.hidden 132 106 if getattr(p_info, 'profile', None) is not None: … … 167 141 return par 168 142 169 def _intermediates(170 F1, # type: np.ndarray171 F2, # type: np.ndarray172 S, # type: np.ndarray173 scale, # type: float174 effective_radius, # type: float175 beta_mode, # type: bool176 ):177 # type: (...) -> OrderedDict[str, Union[np.ndarray, float]]178 """179 Returns intermediate results for beta approximation-enabled product.180 The result may be an array or a float.181 """182 if beta_mode:183 # TODO: 1. include calculated Q vector184 # TODO: 2. consider implications if there are intermediate results in P(Q)185 parts = OrderedDict((186 ("P(Q)", scale*F2),187 ("S(Q)", S),188 ("beta(Q)", F1**2 / F2),189 ("S_eff(Q)", 1 + (F1**2 / F2)*(S-1)),190 ("effective_radius", effective_radius),191 # ("I(Q)", scale*(F2 + (F1**2)*(S-1)) + bg),192 ))193 else:194 parts = OrderedDict((195 ("P(Q)", scale*F2),196 ("S(Q)", S),197 ("effective_radius", effective_radius),198 ))199 return parts200 201 143 class ProductModel(KernelModel): 202 144 def __init__(self, model_info, P, S): … … 208 150 #: Structure factor modelling interaction between particles. 209 151 self.S = S 210 211 152 #: Model precision. This is not really relevant, since it is the 212 153 #: individual P and S models that control the effective dtype, … … 226 167 # in opencl; or both in opencl, but one in single precision and the 227 168 # other in double precision). 228 229 169 p_kernel = self.P.make_kernel(q_vectors) 230 170 s_kernel = self.S.make_kernel(q_vectors) … … 251 191 def __call__(self, call_details, values, cutoff, magnetic): 252 192 # type: (CallDetails, np.ndarray, float, bool) -> np.ndarray 253 254 193 p_info, s_info = self.info.composition[1] 255 p_npars = p_info.parameters.npars256 p_length = call_details.length[:p_npars]257 p_offset = call_details.offset[:p_npars]258 s_npars = s_info.parameters.npars259 s_length = call_details.length[p_npars:p_npars+s_npars]260 s_offset = call_details.offset[p_npars:p_npars+s_npars]261 262 # Beta mode parameter is the first parameter after P and S parameters263 have_beta_mode = p_info.have_Fq264 beta_mode_offset = 2+p_npars+s_npars265 beta_mode = (values[beta_mode_offset] > 0) if have_beta_mode else False266 if beta_mode and self.p_kernel.dim== '2d':267 raise NotImplementedError("beta not yet supported for 2D")268 269 # R_eff type parameter is the second parameter after P and S parameters270 # unless the model doesn't support beta mode, in which case it is first271 have_radius_type = p_info.effective_radius_type is not None272 radius_type_offset = 2+p_npars+s_npars + (1 if have_beta_mode else 0)273 radius_type = int(values[radius_type_offset]) if have_radius_type else 0274 275 # Retrieve the volume fraction, which is the second of the276 # 'S' parameters in the parameter list, or 2+np in 0-origin,277 # as well as the scale and background.278 volfrac = values[3+p_npars]279 scale, background = values[0], values[1]280 194 281 195 # if there are magnetic parameters, they will only be on the … … 292 206 293 207 # Construct the calling parameters for P. 208 p_npars = p_info.parameters.npars 209 p_length = call_details.length[:p_npars] 210 p_offset = call_details.offset[:p_npars] 294 211 p_details = make_details(p_info, p_length, p_offset, nweights) 295 p_values = [ 296 [1., 0.], # scale=1, background=0, 297 values[2:2+p_npars], 298 magnetism, 299 weights] 212 # Set p scale to the volume fraction in s, which is the first of the 213 # 'S' parameters in the parameter list, or 2+np in 0-origin. 214 volfrac = values[2+p_npars] 215 p_values = [[volfrac, 0.0], values[2:2+p_npars], magnetism, weights] 300 216 spacer = (32 - sum(len(v) for v in p_values)%32)%32 301 217 p_values.append([0.]*spacer) 302 218 p_values = np.hstack(p_values).astype(self.p_kernel.dtype) 303 219 220 # Call ER and VR for P since these are needed for S. 221 p_er, p_vr = calc_er_vr(p_info, p_details, p_values) 222 s_vr = (volfrac/p_vr if p_vr != 0. else volfrac) 223 #print("volfrac:%g p_er:%g p_vr:%g s_vr:%g"%(volfrac,p_er,p_vr,s_vr)) 224 304 225 # Construct the calling parameters for S. 305 if radius_type > 0: 306 # If R_eff comes from form factor, make sure it is monodisperse. 307 # weight is set to 1 later, after the value array is created 308 s_length[0] = 1 309 s_details = make_details(s_info, s_length, s_offset, nweights) 226 # The effective radius is not in the combined parameter list, so 227 # the number of 'S' parameters is one less than expected. The 228 # computed effective radius needs to be added into the weights 229 # vector, especially since it is a polydisperse parameter in the 230 # stand-alone structure factor models. We will added it at the 231 # end so the remaining offsets don't need to change. 232 s_npars = s_info.parameters.npars-1 233 s_length = call_details.length[p_npars:p_npars+s_npars] 234 s_offset = call_details.offset[p_npars:p_npars+s_npars] 235 s_length = np.hstack((1, s_length)) 236 s_offset = np.hstack((nweights, s_offset)) 237 s_details = make_details(s_info, s_length, s_offset, nweights+1) 238 v, w = weights[:nweights], weights[nweights:] 310 239 s_values = [ 311 [1., 0.], # scale=1, background=0, 312 values[2+p_npars:2+p_npars+s_npars], 313 weights, 240 # scale=1, background=0, radius_effective=p_er, volfraction=s_vr 241 [1., 0., p_er, s_vr], 242 # structure factor parameters start after scale, background and 243 # all the form factor parameters. Skip the volfraction parameter 244 # as well, since it is computed elsewhere, and go to the end of the 245 # parameter list. 246 values[2+p_npars+1:2+p_npars+s_npars], 247 # no magnetism parameters to include for S 248 # add er into the (value, weights) pairs 249 v, [p_er], w, [1.0] 314 250 ] 315 251 spacer = (32 - sum(len(v) for v in s_values)%32)%32 … … 317 253 s_values = np.hstack(s_values).astype(self.s_kernel.dtype) 318 254 319 # Call the form factor kernel to compute <F> and <F^2>. 320 # If the model doesn't support Fq the returned <F> will be None. 321 F1, F2, effective_radius, shell_volume, volume_ratio = self.p_kernel.Fq( 322 p_details, p_values, cutoff, magnetic, radius_type) 323 324 # Call the structure factor kernel to compute S. 325 # Plug R_eff from the form factor into structure factor parameters 326 # and scale volume fraction by form:shell volume ratio. These changes 327 # needs to be both in the initial value slot as well as the 328 # polydispersity distribution slot in the values array due to 329 # implementation details in kernel_iq.c. 330 #print("R_eff=%d:%g, volfrac=%g, volume ratio=%g"%(radius_type, effective_radius, volfrac, volume_ratio)) 331 if radius_type > 0: 332 # set the value to the model R_eff and set the weight to 1 333 s_values[2] = s_values[2+s_npars+s_offset[0]] = effective_radius 334 s_values[2+s_npars+s_offset[0]+nweights] = 1.0 335 s_values[3] = s_values[2+s_npars+s_offset[1]] = volfrac*volume_ratio 336 S = self.s_kernel.Iq(s_details, s_values, cutoff, False) 337 338 # Determine overall scale factor. Hollow shapes are weighted by 339 # shell_volume, so that is needed for volume normalization. For 340 # solid shapes we can use shell_volume as well since it is equal 341 # to form volume. 342 combined_scale = scale*volfrac/shell_volume 343 344 # Combine form factor and structure factor 345 #print("beta", beta_mode, F1, F2, S) 346 PS = F2 + F1**2*(S-1) if beta_mode else F2*S 347 final_result = combined_scale*PS + background 348 349 # Capture intermediate values so user can see them. These are 350 # returned as a lazy evaluator since they are only needed in the 351 # GUI, and not for each evaluation during a fit. 352 # TODO: return the results structure with the final results 353 # That way the model calcs are idempotent. Further, we can 354 # generalize intermediates to various other model types if we put it 355 # kernel calling interface. Could do this as an "optional" 356 # return value in the caller, though in that case we could return 357 # the results directly rather than through a lazy evaluator. 358 self.results = lambda: _intermediates( 359 F1, F2, S, combined_scale, effective_radius, beta_mode) 360 361 return final_result 255 # Call the kernels 256 p_result = self.p_kernel(p_details, p_values, cutoff, magnetic) 257 s_result = self.s_kernel(s_details, s_values, cutoff, False) 258 259 #print("p_npars",p_npars,s_npars,p_er,s_vr,values[2+p_npars+1:2+p_npars+s_npars]) 260 #call_details.show(values) 261 #print("values", values) 262 #p_details.show(p_values) 263 #print("=>", p_result) 264 #s_details.show(s_values) 265 #print("=>", s_result) 266 267 # remember the parts for plotting later 268 self.results = [p_result, s_result] 269 270 #import pylab as plt 271 #plt.subplot(211); plt.loglog(self.p_kernel.q_input.q, p_result, '-') 272 #plt.subplot(212); plt.loglog(self.s_kernel.q_input.q, s_result, '-') 273 #plt.figure() 274 275 return values[0]*(p_result*s_result) + values[1] 362 276 363 277 def release(self): … … 365 279 self.p_kernel.release() 366 280 self.s_kernel.release() 281 282 283 def calc_er_vr(model_info, call_details, values): 284 # type: (ModelInfo, ParameterSet) -> Tuple[float, float] 285 286 if model_info.ER is None and model_info.VR is None: 287 return 1.0, 1.0 288 289 nvalues = model_info.parameters.nvalues 290 value = values[nvalues:nvalues + call_details.num_weights] 291 weight = values[nvalues + call_details.num_weights: nvalues + 2*call_details.num_weights] 292 npars = model_info.parameters.npars 293 # Note: changed from pairs ([v], [w]) to triples (p, [v], [w]), but the 294 # dispersion mesh code doesn't actually care about the nominal parameter 295 # value, p, so set it to None. 296 pairs = [(None, value[offset:offset+length], weight[offset:offset+length]) 297 for p, offset, length 298 in zip(model_info.parameters.call_parameters[2:2+npars], 299 call_details.offset, 300 call_details.length) 301 if p.type == 'volume'] 302 value, weight = dispersion_mesh(model_info, pairs) 303 304 if model_info.ER is not None: 305 individual_radii = model_info.ER(*value) 306 radius_effective = np.sum(weight*individual_radii) / np.sum(weight) 307 else: 308 radius_effective = 1.0 309 310 if model_info.VR is not None: 311 whole, part = model_info.VR(*value) 312 volume_ratio = np.sum(weight*part)/np.sum(weight*whole) 313 else: 314 volume_ratio = 1.0 315 316 return radius_effective, volume_ratio
Note: See TracChangeset
for help on using the changeset viewer.