Changeset e44432d in sasmodels for sasmodels/product.py
- Timestamp:
- Oct 19, 2018 3:46:26 PM (5 years ago)
- Branches:
- master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- be43e39
- Parents:
- 2586ab72
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/product.py
r2773c66 re44432d 168 168 return par 169 169 170 def _intermediates(P, S, effective_radius): 171 # type: (np.ndarray, np.ndarray, float) -> OrderedDict[str, np.ndarray] 172 """ 173 Returns intermediate results for standard product (P(Q)*S(Q)) 174 """ 175 return OrderedDict(( 176 ("P(Q)", P), 177 ("S(Q)", S), 178 ("effective_radius", effective_radius), 179 )) 180 181 def _intermediates_beta(F1, # type: np.ndarray 182 F2, # type: np.ndarray 183 S, # type: np.ndarray 184 scale, # type: np.ndarray 185 bg, # type: np.ndarray 186 effective_radius # type: float 187 ): 170 def _intermediates( 171 F1, # type: np.ndarray 172 F2, # type: np.ndarray 173 S, # type: np.ndarray 174 scale, # type: float 175 effective_radius, # type: float 176 beta_mode, # type: bool 177 ): 188 178 # type: (...) -> OrderedDict[str, Union[np.ndarray, float]] 189 179 """ 190 Returns intermediate results for beta approximation-enabled product. The result may be an array or a float. 191 """ 192 # TODO: 1. include calculated Q vector 193 # TODO: 2. consider implications if there are intermediate results in P(Q) 194 return OrderedDict(( 195 ("P(Q)", scale*F2), 196 ("S(Q)", S), 197 ("beta(Q)", F1**2 / F2), 198 ("S_eff(Q)", 1 + (F1**2 / F2)*(S-1)), 199 ("effective_radius", effective_radius), 200 # ("I(Q)", scale*(F2 + (F1**2)*(S-1)) + bg), 201 )) 180 Returns intermediate results for beta approximation-enabled product. 181 The result may be an array or a float. 182 """ 183 if beta_mode: 184 # TODO: 1. include calculated Q vector 185 # TODO: 2. consider implications if there are intermediate results in P(Q) 186 parts = OrderedDict(( 187 ("P(Q)", scale*F2), 188 ("S(Q)", S), 189 ("beta(Q)", F1**2 / F2), 190 ("S_eff(Q)", 1 + (F1**2 / F2)*(S-1)), 191 ("effective_radius", effective_radius), 192 # ("I(Q)", scale*(F2 + (F1**2)*(S-1)) + bg), 193 )) 194 else: 195 parts = OrderedDict(( 196 ("P(Q)", scale*F2), 197 ("S(Q)", S), 198 ("effective_radius", effective_radius), 199 )) 200 return parts 202 201 203 202 class ProductModel(KernelModel): … … 253 252 def __call__(self, call_details, values, cutoff, magnetic): 254 253 # type: (CallDetails, np.ndarray, float, bool) -> np.ndarray 254 255 255 p_info, s_info = self.info.composition[1] 256 p_npars = p_info.parameters.npars 257 p_length = call_details.length[:p_npars] 258 p_offset = call_details.offset[:p_npars] 259 s_npars = s_info.parameters.npars 260 s_length = call_details.length[p_npars:p_npars+s_npars] 261 s_offset = call_details.offset[p_npars:p_npars+s_npars] 262 263 # Beta mode parameter is the first parameter after P and S parameters 264 have_beta_mode = p_info.have_Fq 265 beta_mode_offset = 2+p_npars+s_npars 266 beta_mode = (values[beta_mode_offset] > 0) if have_beta_mode else False 267 if beta_mode and self.p_kernel.dim== '2d': 268 raise NotImplementedError("beta not yet supported for 2D") 269 270 # R_eff type parameter is the second parameter after P and S parameters 271 # unless the model doesn't support beta mode, in which case it is first 272 have_radius_type = p_info.effective_radius_type is not None 273 radius_type_offset = 2+p_npars+s_npars + (1 if have_beta_mode else 0) 274 radius_type = int(values[radius_type_offset]) if have_radius_type else 0 275 276 # Retrieve the volume fraction, which is the second of the 277 # 'S' parameters in the parameter list, or 2+np in 0-origin, 278 # as well as the scale and background. 279 volfrac = values[3+p_npars] 280 scale, background = values[0], values[1] 256 281 257 282 # if there are magnetic parameters, they will only be on the … … 268 293 269 294 # Construct the calling parameters for P. 270 p_npars = p_info.parameters.npars271 p_length = call_details.length[:p_npars]272 p_offset = call_details.offset[:p_npars]273 295 p_details = make_details(p_info, p_length, p_offset, nweights) 274 275 # Set p scale to the volume fraction in s, which is the first of the276 # 'S' parameters in the parameter list, or 2+np in 0-origin.277 volfrac = values[2+p_npars]278 p_values = [[volfrac, 0.0], values[2:2+p_npars], magnetism,weights]296 p_values = [ 297 [1., 0.], # scale=1, background=0, 298 values[2:2+p_npars], 299 magnetism, 300 weights] 279 301 spacer = (32 - sum(len(v) for v in p_values)%32)%32 280 302 p_values.append([0.]*spacer) … … 282 304 283 305 # Construct the calling parameters for S. 284 s_npars = s_info.parameters.npars 285 s_length = call_details.length[p_npars:p_npars+s_npars] 286 s_offset = call_details.offset[p_npars:p_npars+s_npars] 287 s_length = np.hstack((1, s_length)) 288 s_offset = np.hstack((nweights, s_offset)) 289 s_details = make_details(s_info, s_length, s_offset, nweights+1) 306 if radius_type > 0: 307 # If R_eff comes from form factor, make sure it is monodisperse. 308 # weight is set to 1 later, after the value array is created 309 s_length[0] = 1 310 s_details = make_details(s_info, s_length, s_offset, nweights) 290 311 s_values = [ 291 # scale=1, background=0, 292 [1., 0.], 312 [1., 0.], # scale=1, background=0, 293 313 values[2+p_npars:2+p_npars+s_npars], 294 314 weights, … … 298 318 s_values = np.hstack(s_values).astype(self.s_kernel.dtype) 299 319 300 # beta mode is the first parameter after the structure factor pars 301 extra_offset = 2+p_npars+s_npars 302 if p_info.have_Fq: 303 beta_mode = values[extra_offset] 304 extra_offset += 1 305 else: 306 beta_mode = 0 307 if p_info.effective_radius_type is not None: 308 effective_radius_type = int(values[extra_offset]) 309 extra_offset += 1 310 else: 311 effective_radius_type = 0 312 313 # Call the kernels 314 scale, background = values[0], values[1] 315 if beta_mode: 316 F1, F2, volume_avg, effective_radius = self.p_kernel.beta( 317 p_details, p_values, cutoff, magnetic, effective_radius_type) 318 if effective_radius_type > 0: 319 # Plug R_eff from p into S model (initial value and pd value) 320 s_values[2] = s_values[2+s_npars+s_offset[0]] = effective_radius 321 s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False) 322 combined_scale = scale*volfrac/volume_avg 323 324 self.results = lambda: _intermediates_beta(F1, F2, s_result, volfrac/volume_avg, background, effective_radius) 325 final_result = combined_scale*(F2 + (F1**2)*(s_result - 1)) + background 326 327 else: 328 p_result, effective_radius = self.p_kernel.Pq_Reff( 329 p_details, p_values, cutoff, magnetic, effective_radius_type) 330 if effective_radius_type > 0: 331 # Plug R_eff from p into S model (initial value and pd value) 332 s_values[2] = s_values[2+s_npars+s_offset[0]] = effective_radius 333 s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False) 334 # remember the parts for plotting later 335 self.results = lambda: _intermediates(p_result, s_result, effective_radius) 336 final_result = scale*(p_result*s_result) + background 337 338 #call_details.show(values) 339 #print("values", values) 340 #p_details.show(p_values) 341 #print("=>", p_result) 342 #s_details.show(s_values) 343 #print("=>", s_result) 344 #import pylab as plt 345 #plt.subplot(211); plt.loglog(self.p_kernel.q_input.q, p_result, '-') 346 #plt.subplot(212); plt.loglog(self.s_kernel.q_input.q, s_result, '-') 347 #plt.figure() 320 # Call the form factor kernel to compute <F> and <F^2>. 321 # If the model doesn't support Fq the returned <F> will be None. 322 F1, F2, effective_radius, shell_volume, volume_ratio = self.p_kernel.Fq( 323 p_details, p_values, cutoff, magnetic, radius_type) 324 325 # Call the structure factor kernel to compute S. 326 # Plug R_eff from the form factor into structure factor parameters 327 # and scale volume fraction by form:shell volume ratio. These changes 328 # needs to be both in the initial value slot as well as the 329 # polydispersity distribution slot in the values array due to 330 # implementation details in kernel_iq.c. 331 #print("R_eff=%d:%g, volfrac=%g, volume ratio=%g"%(radius_type, effective_radius, volfrac, volume_ratio)) 332 if radius_type > 0: 333 # set the value to the model ER and set the weight to 1 334 s_values[2] = s_values[2+s_npars+s_offset[0]] = effective_radius 335 s_values[2+s_npars+s_offset[0]+nweights] = 1.0 336 s_values[3] = s_values[2+s_npars+s_offset[1]] = volfrac*volume_ratio 337 S = self.s_kernel.Iq(s_details, s_values, cutoff, False) 338 339 # Determine overall scale factor. Hollow shapes are weighted by 340 # shell_volume, so that is needed for volume normalization. For 341 # solid shapes we can use shell_volume as well since it is equal 342 # to form volume. 343 combined_scale = scale*volfrac/shell_volume 344 345 # Combine form factor and structure factor 346 #print("beta", beta_mode, F1, F2, S) 347 PS = F2 + F1**2*(S-1) if beta_mode else F2*S 348 final_result = combined_scale*PS + background 349 350 # Capture intermediate values so user can see them. These are 351 # returned as a lazy evaluator since they are only needed in the 352 # GUI, and not for each evaluation during a fit. 353 # TODO: return the results structure with the final results 354 # That way the model calcs are idempotent. Further, we can 355 # generalize intermediates to various other model types if we put it 356 # kernel calling interface. Could do this as an "optional" 357 # return value in the caller, though in that case we could return 358 # the results directly rather than through a lazy evaluator. 359 self.results = lambda: _intermediates( 360 F1, F2, S, combined_scale, effective_radius, beta_mode) 361 348 362 return final_result 349 363
Note: See TracChangeset
for help on using the changeset viewer.