# Changeset 89dba62 in sasmodels

Ignore:
Timestamp:
Mar 29, 2019 2:22:20 PM (7 weeks ago)
Branches:
ticket-1257-vesicle-product
Children:
93cac17
Parents:
98c045a (diff), de032da (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'ticket_822_v5_unit_tests' into ticket-1257-vesicle-product

Location:
sasmodels
Files:
7 edited
1 moved

Unmodified
Removed
• ## sasmodels/conversion_table.py

 rdb1d9d5 "TwoPowerLawModel", ], "unified_power_rg": [ "unified_power_Rg": [ "UnifiedPowerRg", dict(((field_new+str(index), field_old+str(index))
• ## sasmodels/convert.py

 rdb1d9d5 if pars['volfraction_a'] > 0.5: pars['volfraction_a'] = 1.0 - pars['volfraction_a'] elif name == 'unified_power_rg': elif name == 'unified_power_Rg': pars['level'] = int(pars['level'])
• ## sasmodels/kernel_iq.c

 ra34b811 static void set_spin_weights(double in_spin, double out_spin, double weight[6]) { double norm; in_spin = clip(in_spin, 0.0, 1.0); out_spin = clip(out_spin, 0.0, 1.0); // However, since the weights are applied to the final intensity and // are not interned inside the I(q) function, we want the full // weight and not the square root.  Any function using // set_spin_weights as part of calculating an amplitude will need to // manually take that square root, but there is currently no such // function. weight[0] = (1.0-in_spin) * (1.0-out_spin); // dd weight[1] = (1.0-in_spin) * out_spin;       // du weight[2] = in_spin * (1.0-out_spin);       // ud weight[3] = in_spin * out_spin;             // uu // weight and not the square root.  Anyway no function will ever use // set_spin_weights as part of calculating an amplitude, as the weights are // related to polarisation efficiency of the instrument. The weights serve to // construct various magnet scattering cross sections, which are linear combinations // of the spin-resolved cross sections. The polarisation efficiency e_in and e_out // are parameters ranging from 0.5 (unpolarised) beam to 1 (perfect optics). // For in_spin or out_spin <0.5 one assumes a CS, where the spin is reversed/flipped // with respect to the initial supermirror polariser. The actual polarisation efficiency // in this case is however e_in/out = 1-in/out_spin. if (out_spin < 0.5){norm=1-out_spin;} else{norm=out_spin;} // The norm is needed to make sure that the scattering cross sections are //correctly weighted, such that the sum of spin-resolved measurements adds up to // the unpolarised or half-polarised scattering cross section. No intensity weighting // needed on the incoming polariser side (assuming that a user), has normalised // to the incoming flux with polariser in for SANSPOl and unpolarised beam, respectively. weight[0] = (1.0-in_spin) * (1.0-out_spin) / norm; // dd weight[1] = (1.0-in_spin) * out_spin / norm;       // du weight[2] = in_spin * (1.0-out_spin) / norm;       // ud weight[3] = in_spin * out_spin / norm;             // uu weight[4] = weight[1]; // du.imag weight[5] = weight[2]; // ud.imag switch (xs) { default: // keep compiler happy; condition ensures xs in [0,1,2,3] case 0: // uu => sld - D M_perpx case 0: // dd => sld - D M_perpx return sld - px*perp; case 1: // ud.real => -D M_perpy case 1: // du.real => -D M_perpy return py*perp; case 2: // du.real => -D M_perpy case 2: // ud.real => -D M_perpy return py*perp; case 3: // dd => sld + D M_perpx case 3: // uu => sld + D M_perpx return sld + px*perp; }
• ## sasmodels/models/spherical_sld.py

 ra34b811 sub-shell is described by a line function, with *n_steps* sub-shells per interface. The form factor is normalized by the total volume of the sphere. .. note:: *n_shells* must be an integer. *n_steps* must be an ODD integer. Interface shapes are as follows: 3 \rho_\text{solvent} V(r_N) \Big[ \frac{\sin(qr_N) - qr_N \cos(qr_N)} {qr_N^3} \Big] Here we assumed that the SLDs of the core and solvent are constant in $r$. \end{align*} We assume $\rho_{\text{inter}_j} (r)$ is approximately linear within the sub-shell $j$. when $P(Q) * S(Q)$ is applied. References ---------- Authorship and Verification ---------------------------- --------------------------- * **Author:** Jae-Hie Cho **Date:** Nov 1, 2010 * **Last Modified by:** Paul Kienzle **Date:** Dec 20, 2016 * **Last Reviewed by:** Paul Butler **Date:** September 8, 2018 * **Last Reviewed by:** Steve King **Date:** March 29, 2019 * **Source added by :** Steve King **Date:** March 25, 2019 """ name = "spherical_sld" title = "Sperical SLD intensity calculation" title = "Spherical SLD intensity calculation" description = """ I(q) = # pylint: disable=bad-whitespace, line-too-long #            ["name", "units", default, [lower, upper], "type", "description"], parameters = [["n_shells",             "",           1,      [1, 10],        "volume", "number of shells"], parameters = [["n_shells",             "",           1,      [1, 10],        "volume", "number of shells (must be integer)"], ["sld_solvent",          "1e-6/Ang^2", 1.0,    [-inf, inf],    "sld", "solvent sld"], ["sld[n_shells]",        "1e-6/Ang^2", 4.06,   [-inf, inf],    "sld", "sld of the shell"],
• ## sasmodels/models/unified_power_Rg.py

 r7a5f8af ------ unified_power_rg.py _ unified_power_Rg.py _ Authorship and Verification category = "shape-independent" name = "unified_power_rg" name = "unified_power_Rg" title = "Unified Power Rg" description = """
• ## sasmodels/mixture.py

 rb297ba9 combined_pars.append(p) parameters = ParameterTable(combined_pars) # Allow for the scenario in which each component has all its PD parameters # active simultaneously.  details.make_details() will throw an error if # too many are used from any one component. parameters.max_pd = sum(part.parameters.max_pd for part in parts)
• ## sasmodels/modelinfo.py

 ra34b811 processed.append(parse_parameter(*p)) partable = ParameterTable(processed) partable.check_angles() partable.check_angles(strict=True) return partable # properties, such as default=0.0 for structure factor backgrounds. self.common_parameters = [Parameter(*p) for p in COMMON_PARAMETERS] self.kernel_parameters = parameters self._set_vector_lengths() self.pd_2d = set(p.name for p in self.call_parameters if p.polydisperse) # Final checks self.check_duplicates() self.check_angles() def set_zero_background(self): """ self.defaults = self._get_defaults() def check_angles(self): def check_angles(self, strict=False): """ Check that orientation angles are theta, phi and possibly psi. *strict* should be True when checking a parameter table defined in a model file, but False when checking from mixture models, etc., where the parameters aren't being passed to a calculator directly. """ theta = phi = psi = -1 if p.type != 'orientation': raise TypeError("psi must be an orientation parameter") elif p.type == 'orientation': elif p.type == 'orientation' and strict: raise TypeError("only theta, phi and psi can be orientation parameters") if theta >= 0 and phi >= 0: if psi >= 0 and psi != phi+1: raise TypeError("psi must follow phi") # TODO: Why must theta/phi/psi be at the end?  Consistency only? if (psi >= 0 and psi != last_par) or (psi < 0 and phi != last_par): raise TypeError("orientation parameters must appear at the " "end of the parameter table") if strict: raise TypeError("orientation parameters must appear at the " "end of the parameter table") elif theta >= 0 or phi >= 0 or psi >= 0: raise TypeError("oriented shapes must have both theta and phi and maybe psi") def check_duplicates(self): """ Check for duplicate parameter names """ checked, dups = set(), set() for p in self.call_parameters: if p.id in checked: dups.add(p.id) else: checked.add(p.id) if dups: raise TypeError("Duplicate parameters: {}" .format(", ".join(sorted(dups)))) def __getitem__(self, key):
• ## sasmodels/product.py

 r8795b6f To use it, first load form factor P and structure factor S, then create *make_product_info(P, S)*. The P@S models is somewhat complicated because there are many special parameters that need to be handled in particular ways.  Much of the code is used to figure out what special parameters we have, where to find them in the P@S model inputs and how to distribute them to the underlying P and S model calculators. The parameter packet received by the P@S is a :class:details.CallDetails structure along with a data vector. The CallDetails structure indicates which parameters are polydisperse, the length of the distribution, and where to find it in the data vector.  The distributions are ordered from longest to shortest, with length 1 distributions filling out the distribution set.  That way the kernel calculator doesn't have to check if it needs another nesting level since it is always there.  The data vector consists of a list of target values for the parameters, followed by a concatenation of the distribution values, and then followed by a concatenation of the distribution weights. Given the combined details and data for P@S, we must decompose them in to details for P and details for S separately, which unfortunately requires intimate knowledge of the data structures and tricky code. The special parameters are: * *scale* and *background*: First two parameters of the value list in each of P, S and P@S. When decomposing P@S parameters, ignore *scale* and *background*, instead using 1 and 0 for the first two slots of both P and S. After calling P and S individually, the results are combined as :code:volfraction*scale*P*S + background.  The *scale* and *background* do not show up in the polydispersity structure so they are easy to handle. * *volfraction*: Always the first parameter of S, but it may also be in P. If it is in P, then *P.volfraction* is used in the combined P@S list, and *S.volfraction* is elided, otherwise *S.volfraction* is used. If we are using *volfraction* from P we can treat it like all the other P parameters when calling P, but when calling S we need to insert the *P.volfraction* into data vector for S and assign a slot of length 1 in the distribution. Because we are using the original layout of the distribution vectors from P@S, but copying it into private data vectors for S and P, we are free to "borrow" a P slots to store the missing *S.volfraction* distribution.  We use the *P.volfraction* slot itself but any slot will work. For hollow shapes, *volfraction* represents the volume fraction of material but S needs the volume fraction enclosed by the shape. The answer is to scale the user specified volume fraction by the form:shell ratio computed from the average form volume and average shell volume returned from P. Use the original *volfraction* divided by *shell_volume* to compute the number density, and scale P@S by that to get absolute scaling on the final *I(q)*. The *scale* for P@S should therefore usually be one. * *radius_effective*: Always the second parameter of S and always part of P@S, but never in P. The value may be calculated using *P.radius_effective()* or it may be set to the *radius_effective* value in P@S, depending on *radius_effective_mode*.  If part of S, the value may be polydisperse. If calculated by P, then it will be the weighted average of effective radii computed for the polydisperse shape parameters. * *structure_factor_mode* If P@S supports beta approximation (i.e., if it has the *Fq* function that returns and ), then *structure_factor_mode* will be added to the P@S parameters right after the S parameters.  This mode may be 0 for the monodisperse approximation or 1 for the beta approximation.  We will add more values here as we implemented more complicated operations, but for now P and S must be computed separately.  If beta, then we return *I = scale volfrac/volume ( + ^2 (S-1)) + background*. If not beta then return *I = scale/volume P S + background* .  In both cases, return the appropriate immediate values. * *radius_effective_mode* If P defines the *radius_effective* function (and therefore *P.info.radius_effective_modes* is a list of effective radius modes), then *radius_effective_mode* will be the final parameter in P@S.  Mode will be zero if *radius_effective* is defined by the user using the S parameter; any other value and the *radius_effective* parameter will be filled in from the value computed in P.  In the latter case, the polydispersity information for *S.radius_effective* will need to be suppressed, with pd length set to 1, the first value set to the effective radius and the first weight set to 1.  Do this after composing the S data vector so the inputs are left untouched. * *regular parameters* The regular P parameters form a block of length *P.info.npars* at the start of the data vector (after scale and background).  These will be followed by *S.effective_radius*, and *S.volfraction* (if *P.volfraction* is absent), and then the regular S parameters.  The P and S blocks can be copied as a group into the respective P and S data vectors. We can copy the distribution value and weight vectors untouched to both the P and S data vectors since they are referenced by offset and length. We can update the radius_effective slots in the P data vector with *P.radius_effective()* if needed. * *magnetic parameters* For each P parameter that is an SLD there will be a set of three magnetic parameters tacked on to P@S after the regular P and S and after the special *structure_factor_mode* and *radius_effective_mode*.  These can be copied as a group after the regular P parameters.  There won't be any magnetic S parameters. """ from __future__ import print_function, division if not s_info.parameters.magnetism_index == []: raise TypeError("S should not have SLD parameters") if RADIUS_ID in p_info.parameters: raise TypeError("P should not have {}".format(RADIUS_ID)) p_id, p_name, p_pars = p_info.id, p_info.name, p_info.parameters s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters # Create list of parameters for the combined model.  If there # are any names in P that overlap with those in S, modify the name in S # to distinguish it. p_has_volfrac = VOLFRAC_ID in p_info.parameters # Create list of parameters for the combined model.  If a name in # S overlaps a name in P, tag the S parameter name to distinguish it. # If the tagged name also collides it will be caught by the parameter # table builder.  Similarly if any special names are abused.  Need the # pairs to create the translation table for random model generation. p_set = set(p.id for p in p_pars.kernel_parameters) s_list = [(_tag_parameter(par) if par.id in p_set else par) for par in s_pars.kernel_parameters] # Check if still a collision after renaming.  This could happen if for # example S has volfrac and P has both volfrac and volfrac_S. if any(p.id in p_set for p in s_list): raise TypeError("name collision: P has P.name and P.name_S while S has S.name") # make sure effective radius is not a polydisperse parameter in product s_list[0] = copy(s_list[0]) s_list[0].polydisperse = False translate_name = dict((old.id, new.id) for old, new in zip(s_pars.kernel_parameters, s_list)) s_pairs = [(par, (_tag_parameter(par) if par.id in p_set else par)) for par in s_pars.kernel_parameters # Note: exclude volfraction from s_list if volfraction in p if par.id != VOLFRAC_ID or not p_has_volfrac] s_list = [pair[0] for pair in s_pairs] # Build combined parameter table combined_pars = p_pars.kernel_parameters + s_list + make_extra_pars(p_info) parameters = ParameterTable(combined_pars) parameters.max_pd = p_pars.max_pd + s_pars.max_pd # Allow for the scenario in which each component has all its PD parameters # active simultaneously.  details.make_details() will throw an error if # too many are used from any one component. parameters.Pumax_pd = p_pars.max_pd + s_pars.max_pd # TODO: does user-defined polydisperse S.radius_effective make sense? # make sure effective radius is not a polydisperse parameter in product #s_list[0] = copy(s_list[0]) #s_list[0].polydisperse = False s_translate = {old.id: new.id for old, new in s_pairs} def random(): """Random set of model parameters for product model""" combined_pars = p_info.random() s_names = set(par.id for par in s_pars.kernel_parameters) combined_pars.update((translate_name[k], v) combined_pars.update((s_translate[k], v) for k, v in s_info.random().items() if k in s_names) if k in s_translate) return combined_pars def _intermediates( F1,               # type: np.ndarray F2,               # type: np.ndarray F,                # type: np.ndarray Fsq,              # type: np.ndarray S,                # type: np.ndarray scale,            # type: float # TODO: 2. consider implications if there are intermediate results in P(Q) parts = OrderedDict(( ("P(Q)", scale*F2), ("P(Q)", scale*Fsq), ("S(Q)", S), ("beta(Q)", F1**2 / F2), ("S_eff(Q)", 1 + (F1**2 / F2)*(S-1)), ("beta(Q)", F**2 / Fsq), ("S_eff(Q)", 1 + (F**2 / Fsq)*(S-1)), ("effective_radius", radius_effective), # ("I(Q)", scale*(F2 + (F1**2)*(S-1)) + bg), # ("I(Q)", scale*(Fsq + (F**2)*(S-1)) + bg), )) else: parts = OrderedDict(( ("P(Q)", scale*F2), ("P(Q)", scale*Fsq), ("S(Q)", S), ("effective_radius", radius_effective), self.results = []  # type: List[np.ndarray] # Find index of volfraction parameter in parameter list for k, p in enumerate(model_info.parameters.call_parameters): if p.id == VOLFRAC_ID: self._volfrac_index = k break else: raise RuntimeError("no %s parameter in %s"%(VOLFRAC_ID, self)) p_info, s_info = self.info.composition[1] p_npars = p_info.parameters.npars s_npars = s_info.parameters.npars have_beta_mode = p_info.have_Fq have_er_mode = p_info.radius_effective_modes is not None volfrac_in_p = self._volfrac_index < p_npars + 2  # scale & background # Slices into the details length/offset structure for P@S. # Made complicated by the possibly missing volfraction in S. self._p_detail_slice = slice(0, p_npars) self._s_detail_slice = slice(p_npars, p_npars+s_npars-volfrac_in_p) self._volfrac_in_p = volfrac_in_p # P block from data vector, without scale and background first_p = 2 last_p = p_npars + 2 self._p_value_slice = slice(first_p, last_p) # radius_effective is the first parameter in S from the data vector. self._er_index = last_p # S block from data vector, without scale, background, volfrac or er. first_s = last_p + 2 - volfrac_in_p last_s = first_s + s_npars - 2 self._s_value_slice = slice(first_s, last_s) # S distribution block in S data vector starts after all S values self._s_dist_slice = slice(2 + s_npars, None) # structure_factor_mode is the first parameter after P and S.  Skip # 2 for scale and background, and subtract 1 in case there is no # volfraction in S. self._beta_mode_index = last_s if have_beta_mode else 0 # radius_effective_mode is the second parameter after P and S # unless structure_factor_mode isn't available, in which case it # is first. self._er_mode_index = last_s + have_beta_mode if have_er_mode else 0 # Magnetic parameters are after everything else.  If they exist, # they will only be for form factor P, not structure factor S. first_mag = last_s + have_beta_mode + have_er_mode mag_pars = 3*p_info.parameters.nmagnetic last_mag = first_mag + (mag_pars + 3 if mag_pars else 0) self._magentic_slice = slice(first_mag, last_mag) def Iq(self, call_details, values, cutoff, magnetic): # type: (CallDetails, np.ndarray, float, bool) -> np.ndarray p_info, s_info = self.info.composition[1] p_npars = p_info.parameters.npars p_length = call_details.length[:p_npars] p_offset = call_details.offset[:p_npars] s_npars = s_info.parameters.npars s_length = call_details.length[p_npars:p_npars+s_npars] s_offset = call_details.offset[p_npars:p_npars+s_npars] # Beta mode parameter is the first parameter after P and S parameters have_beta_mode = p_info.have_Fq beta_mode_offset = 2+p_npars+s_npars beta_mode = (values[beta_mode_offset] > 0) if have_beta_mode else False if beta_mode and self.p_kernel.dim == '2d': raise NotImplementedError("beta not yet supported for 2D") # R_eff type parameter is the second parameter after P and S parameters # unless the model doesn't support beta mode, in which case it is first have_radius_type = p_info.radius_effective_modes is not None #print(p_npars,s_npars) radius_type_offset = 2+p_npars+s_npars + (1 if have_beta_mode else 0) #print(values[radius_type_offset]) radius_type = int(values[radius_type_offset]) if have_radius_type else 0 # Retrieve the volume fraction, which is the second of the # 'S' parameters in the parameter list, or 2+np in 0-origin, # as well as the scale and background. volfrac = values[3+p_npars] # Retrieve values from the data vector scale, background = values[0], values[1] # if there are magnetic parameters, they will only be on the # form factor P, not the structure factor S. nmagnetic = len(self.info.parameters.magnetism_index) if nmagnetic: spin_index = self.info.parameters.npars + 2 magnetism = values[spin_index: spin_index+3+3*nmagnetic] else: magnetism = [] volfrac = values[self._volfrac_index] er_mode = (int(values[self._er_mode_index]) if self._er_mode_index > 0 else 0) beta_mode = (values[self._beta_mode_index] > 0 if self._beta_mode_index > 0 else False) nvalues = self.info.parameters.nvalues nweights = call_details.num_weights weights = values[nvalues:nvalues + 2*nweights] # Can't do 2d and beta_mode just yet if beta_mode and self.p_kernel.dim == '2d': raise NotImplementedError("beta not yet supported for 2D") # Construct the calling parameters for P. p_length = call_details.length[self._p_detail_slice] p_offset = call_details.offset[self._p_detail_slice] p_details = make_details(p_info, p_length, p_offset, nweights) p_values = [ [1., 0.], # scale=1, background=0, values[2:2+p_npars], magnetism, values[self._p_value_slice], values[self._magentic_slice], weights] spacer = (32 - sum(len(v) for v in p_values)%32)%32 p_values = np.hstack(p_values).astype(self.p_kernel.dtype) # Call the form factor kernel to compute and . # If the model doesn't support Fq the returned will be None. F, Fsq, radius_effective, shell_volume, volume_ratio \ = self.p_kernel.Fq(p_details, p_values, cutoff, magnetic, er_mode) # TODO: async call to the GPU # Construct the calling parameters for S. if radius_type > 0: # If R_eff comes from form factor, make sure it is monodisperse. # weight is set to 1 later, after the value array is created s_length = call_details.length[self._s_detail_slice] s_offset = call_details.offset[self._s_detail_slice] if self._volfrac_in_p: # Volfrac is in P and missing from S so insert a slot for it.  Say # the distribution is length 1 and use the slot for volfraction # from the P distribution. s_length = np.insert(s_length, 1, 1) s_offset = np.insert(s_offset, 1, p_offset[self._volfrac_index - 2]) if er_mode > 0: # If effective_radius comes from P, make sure it is monodisperse. # Weight is set to 1 later, after the value array is created s_length[0] = 1 s_details = make_details(s_info, s_length, s_offset, nweights) s_values = [ [1., 0.], # scale=1, background=0, values[2+p_npars:2+p_npars+s_npars], [1., # scale=1 0., # background=0, values[self._er_index], # S.radius_effective; may be replaced by P 0.], # volfraction; will be replaced by volfrac * volume_ratio # followed by S parameters after effective_radius and volfraction values[self._s_value_slice], weights, ] s_values = np.hstack(s_values).astype(self.s_kernel.dtype) # Call the form factor kernel to compute and . # If the model doesn't support Fq the returned will be None. F1, F2, radius_effective, shell_volume, volume_ratio = self.p_kernel.Fq( p_details, p_values, cutoff, magnetic, radius_type) # Call the structure factor kernel to compute S. # Plug R_eff from the form factor into structure factor parameters # and scale volume fraction by form:shell volume ratio. These changes #print("R_eff=%d:%g, volfrac=%g, volume ratio=%g" #      % (radius_type, radius_effective, volfrac, volume_ratio)) if radius_type > 0: s_dist = s_values[self._s_dist_slice] if er_mode > 0: # set the value to the model R_eff and set the weight to 1 s_values[2] = s_values[2+s_npars+s_offset[0]] = radius_effective s_values[2+s_npars+s_offset[0]+nweights] = 1.0 s_values[3] = s_values[2+s_npars+s_offset[1]] = volfrac*volume_ratio s_values[2] = s_dist[s_offset[0]] = radius_effective s_dist[s_offset[0]+nweights] = 1.0 s_values[3] = s_dist[s_offset[1]] = volfrac*volume_ratio s_dist[s_offset[1]+nweights] = 1.0 # Call the structure factor kernel to compute S. S = self.s_kernel.Iq(s_details, s_values, cutoff, False) #print("P", Fsq[:10]) #print("S", S[:10]) #print(radius_effective, volfrac*volume_ratio) # Combine form factor and structure factor #print("beta", beta_mode, F, Fsq, S) PS = Fsq + F**2*(S-1) if beta_mode else Fsq*S # Determine overall scale factor. Hollow shapes are weighted by # shell_volume, so that is needed for volume normalization.  For # solid shapes we can use shell_volume as well since it is equal # to form volume. combined_scale = scale*volfrac/shell_volume # Combine form factor and structure factor #print("beta", beta_mode, F1, F2, S) PS = F2 + F1**2*(S-1) if beta_mode else F2*S final_result = combined_scale*PS + background # shell_volume, so that is needed for number density estimation. # For solid shapes we can use shell_volume as well since it is # equal to form volume.  If P already has a volfraction parameter, # then assume that it is already on absolute scale, and don't # include volfrac in the combined_scale. combined_scale = scale*(volfrac if not self._volfrac_in_p else 1.0) final_result = combined_scale/shell_volume*PS + background # Capture intermediate values so user can see them.  These are # the results directly rather than through a lazy evaluator. self.results = lambda: _intermediates( F1, F2, S, combined_scale, radius_effective, beta_mode) F, Fsq, S, combined_scale, radius_effective, beta_mode) return final_result
Note: See TracChangeset for help on using the changeset viewer.