Changeset 89dba62 in sasmodels


Ignore:
Timestamp:
Mar 29, 2019 2:22:20 PM (6 years ago)
Author:
Paul Kienzle <pkienzle@…>
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

Legend:

Unmodified
Added
Removed

  • sasmodels/conversion_table.py

    rdb1d9d5 ref476e6  
    854854            "TwoPowerLawModel", 
    855855        ], 
    856         "unified_power_rg": [ 
     856        "unified_power_Rg": [ 
    857857            "UnifiedPowerRg", 
    858858            dict(((field_new+str(index), field_old+str(index)) 

  • sasmodels/convert.py

    rdb1d9d5 ref476e6  
    606606            if pars['volfraction_a'] > 0.5: 
    607607                pars['volfraction_a'] = 1.0 - pars['volfraction_a'] 
    608         elif name == 'unified_power_rg': 
     608        elif name == 'unified_power_Rg': 
    609609            pars['level'] = int(pars['level']) 
    610610 

  • sasmodels/kernel_iq.c

    ra34b811 rde032da  
    8585static void set_spin_weights(double in_spin, double out_spin, double weight[6]) 
    8686{ 
     87 
     88  double norm; 
    8789  in_spin = clip(in_spin, 0.0, 1.0); 
    8890  out_spin = clip(out_spin, 0.0, 1.0); 
    … …  
    9496  // However, since the weights are applied to the final intensity and 
    9597  // are not interned inside the I(q) function, we want the full 
    96   // weight and not the square root.  Any function using 
    97   // set_spin_weights as part of calculating an amplitude will need to 
    98   // manually take that square root, but there is currently no such 
    99   // function. 
    100   weight[0] = (1.0-in_spin) * (1.0-out_spin); // dd 
    101   weight[1] = (1.0-in_spin) * out_spin;       // du 
    102   weight[2] = in_spin * (1.0-out_spin);       // ud 
    103   weight[3] = in_spin * out_spin;             // uu 
     98  // weight and not the square root.  Anyway no function will ever use 
     99  // set_spin_weights as part of calculating an amplitude, as the weights are 
     100  // related to polarisation efficiency of the instrument. The weights serve to 
     101  // construct various magnet scattering cross sections, which are linear combinations 
     102  // of the spin-resolved cross sections. The polarisation efficiency e_in and e_out 
     103  // are parameters ranging from 0.5 (unpolarised) beam to 1 (perfect optics). 
     104  // For in_spin or out_spin <0.5 one assumes a CS, where the spin is reversed/flipped 
     105  // with respect to the initial supermirror polariser. The actual polarisation efficiency 
     106  // in this case is however e_in/out = 1-in/out_spin. 
     107 
     108  if (out_spin < 0.5){norm=1-out_spin;} 
     109  else{norm=out_spin;} 
     110 
     111 
     112// The norm is needed to make sure that the scattering cross sections are 
     113//correctly weighted, such that the sum of spin-resolved measurements adds up to 
     114// the unpolarised or half-polarised scattering cross section. No intensity weighting 
     115// needed on the incoming polariser side (assuming that a user), has normalised 
     116// to the incoming flux with polariser in for SANSPOl and unpolarised beam, respectively. 
     117 
     118 
     119  weight[0] = (1.0-in_spin) * (1.0-out_spin) / norm; // dd 
     120  weight[1] = (1.0-in_spin) * out_spin / norm;       // du 
     121  weight[2] = in_spin * (1.0-out_spin) / norm;       // ud 
     122  weight[3] = in_spin * out_spin / norm;             // uu 
    104123  weight[4] = weight[1]; // du.imag 
    105124  weight[5] = weight[2]; // ud.imag 
    … …  
    119138    switch (xs) { 
    120139      default: // keep compiler happy; condition ensures xs in [0,1,2,3] 
    121       case 0: // uu => sld - D M_perpx 
     140      case 0: // dd => sld - D M_perpx 
    122141          return sld - px*perp; 
    123       case 1: // ud.real => -D M_perpy 
     142      case 1: // du.real => -D M_perpy 
    124143          return py*perp; 
    125       case 2: // du.real => -D M_perpy 
     144      case 2: // ud.real => -D M_perpy 
    126145          return py*perp; 
    127       case 3: // dd => sld + D M_perpx 
     146      case 3: // uu => sld + D M_perpx 
    128147          return sld + px*perp; 
    129148    } 

  • sasmodels/models/spherical_sld.py

    ra34b811 r627b68b  
    1818sub-shell is described by a line function, with *n_steps* sub-shells per 
    1919interface. The form factor is normalized by the total volume of the sphere. 
     20 
     21.. note:: 
     22 
     23   *n_shells* must be an integer. *n_steps* must be an ODD integer. 
    2024 
    2125Interface shapes are as follows: 
    … …  
    7377    3 \rho_\text{solvent} V(r_N) 
    7478    \Big[ \frac{\sin(qr_N) - qr_N \cos(qr_N)} {qr_N^3} \Big] 
    75  
    7679 
    7780Here we assumed that the SLDs of the core and solvent are constant in $r$. 
    … …  
    156159    \end{align*} 
    157160 
    158  
    159161We assume $\rho_{\text{inter}_j} (r)$ is approximately linear 
    160162within the sub-shell $j$. 
    … …  
    179181    when $P(Q) * S(Q)$ is applied. 
    180182 
    181  
    182183References 
    183184---------- 
    … …  
    194195 
    195196Authorship and Verification 
    196 ---------------------------- 
     197--------------------------- 
    197198 
    198199* **Author:** Jae-Hie Cho **Date:** Nov 1, 2010 
    199200* **Last Modified by:** Paul Kienzle **Date:** Dec 20, 2016 
    200 * **Last Reviewed by:** Paul Butler **Date:** September 8, 2018 
     201* **Last Reviewed by:** Steve King **Date:** March 29, 2019 
    201202* **Source added by :** Steve King **Date:** March 25, 2019 
    202203""" 
    … …  
    207208 
    208209name = "spherical_sld" 
    209 title = "Sperical SLD intensity calculation" 
     210title = "Spherical SLD intensity calculation" 
    210211description = """ 
    211212            I(q) = 
    … …  
    219220# pylint: disable=bad-whitespace, line-too-long 
    220221#            ["name", "units", default, [lower, upper], "type", "description"], 
    221 parameters = [["n_shells",             "",           1,      [1, 10],        "volume", "number of shells"], 
     222parameters = [["n_shells",             "",           1,      [1, 10],        "volume", "number of shells (must be integer)"], 
    222223              ["sld_solvent",          "1e-6/Ang^2", 1.0,    [-inf, inf],    "sld", "solvent sld"], 
    223224              ["sld[n_shells]",        "1e-6/Ang^2", 4.06,   [-inf, inf],    "sld", "sld of the shell"], 

  • sasmodels/models/unified_power_Rg.py

    r7a5f8af ref476e6  
    7070------ 
    7171 
    72 `unified_power_rg.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/unified_power_rg.py>`_ 
     72`unified_power_Rg.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/unified_power_Rg.py>`_ 
    7373 
    7474Authorship and Verification 
    … …  
    8888 
    8989category = "shape-independent" 
    90 name = "unified_power_rg" 
     90name = "unified_power_Rg" 
    9191title = "Unified Power Rg" 
    9292description = """ 

  • sasmodels/mixture.py

    rb297ba9 rb2f0e5d  
    117117            combined_pars.append(p) 
    118118    parameters = ParameterTable(combined_pars) 
     119    # Allow for the scenario in which each component has all its PD parameters 
     120    # active simultaneously.  details.make_details() will throw an error if 
     121    # too many are used from any one component. 
    119122    parameters.max_pd = sum(part.parameters.max_pd for part in parts) 
    120123 

  • sasmodels/modelinfo.py

    ra34b811 r98c045a  
    7070        processed.append(parse_parameter(*p)) 
    7171    partable = ParameterTable(processed) 
    72     partable.check_angles() 
     72    partable.check_angles(strict=True) 
    7373    return partable 
    7474 
    … …  
    446446        # properties, such as default=0.0 for structure factor backgrounds. 
    447447        self.common_parameters = [Parameter(*p) for p in COMMON_PARAMETERS] 
    448  
    449448        self.kernel_parameters = parameters 
    450449        self._set_vector_lengths() 
    … …  
    495494        self.pd_2d = set(p.name for p in self.call_parameters if p.polydisperse) 
    496495 
     496        # Final checks 
     497        self.check_duplicates() 
     498        self.check_angles() 
     499 
    497500    def set_zero_background(self): 
    498501        """ 
    … …  
    506509        self.defaults = self._get_defaults() 
    507510 
    508     def check_angles(self): 
     511    def check_angles(self, strict=False): 
    509512        """ 
    510513        Check that orientation angles are theta, phi and possibly psi. 
     514 
     515        *strict* should be True when checking a parameter table defined 
     516        in a model file, but False when checking from mixture models, etc., 
     517        where the parameters aren't being passed to a calculator directly. 
    511518        """ 
    512519        theta = phi = psi = -1 
    … …  
    524531                if p.type != 'orientation': 
    525532                    raise TypeError("psi must be an orientation parameter") 
    526             elif p.type == 'orientation': 
     533            elif p.type == 'orientation' and strict: 
    527534                raise TypeError("only theta, phi and psi can be orientation parameters") 
    528535        if theta >= 0 and phi >= 0: 
    … …  
    532539            if psi >= 0 and psi != phi+1: 
    533540                raise TypeError("psi must follow phi") 
     541            # TODO: Why must theta/phi/psi be at the end?  Consistency only? 
    534542            if (psi >= 0 and psi != last_par) or (psi < 0 and phi != last_par): 
    535                 raise TypeError("orientation parameters must appear at the " 
    536                                 "end of the parameter table") 
     543                if strict: 
     544                    raise TypeError("orientation parameters must appear at the " 
     545                                    "end of the parameter table") 
    537546        elif theta >= 0 or phi >= 0 or psi >= 0: 
    538547            raise TypeError("oriented shapes must have both theta and phi and maybe psi") 
     548 
     549    def check_duplicates(self): 
     550        """ 
     551        Check for duplicate parameter names 
     552        """ 
     553        checked, dups = set(), set() 
     554        for p in self.call_parameters: 
     555            if p.id in checked: 
     556                dups.add(p.id) 
     557            else: 
     558                checked.add(p.id) 
     559        if dups: 
     560            raise TypeError("Duplicate parameters: {}" 
     561                            .format(", ".join(sorted(dups)))) 
    539562 
    540563    def __getitem__(self, key): 

  • sasmodels/product.py

    r8795b6f rb2f0e5d  
    1010To use it, first load form factor P and structure factor S, then create 
    1111*make_product_info(P, S)*. 
     12 
     13The P@S models is somewhat complicated because there are many special 
     14parameters that need to be handled in particular ways.  Much of the 
     15code is used to figure out what special parameters we have, where to 
     16find them in the P@S model inputs and how to distribute them to the underlying 
     17P and S model calculators. 
     18 
     19The parameter packet received by the P@S is a :class:`details.CallDetails` 
     20structure along with a data vector. The CallDetails structure indicates which 
     21parameters are polydisperse, the length of the distribution, and where to 
     22find it in the data vector.  The distributions are ordered from longest to 
     23shortest, with length 1 distributions filling out the distribution set.  That 
     24way the kernel calculator doesn't have to check if it needs another nesting 
     25level since it is always there.  The data vector consists of a list of target 
     26values for the parameters, followed by a concatenation of the distribution 
     27values, and then followed by a concatenation of the distribution weights. 
     28Given the combined details and data for P@S, we must decompose them in to 
     29details for P and details for S separately, which unfortunately requires 
     30intimate knowledge of the data structures and tricky code. 
     31 
     32The special parameters are: 
     33 
     34* *scale* and *background*: 
     35    First two parameters of the value list in each of P, S and P@S. 
     36    When decomposing P@S parameters, ignore *scale* and *background*, 
     37    instead using 1 and 0 for the first two slots of both P and S. 
     38    After calling P and S individually, the results are combined as 
     39    :code:`volfraction*scale*P*S + background`.  The *scale* and 
     40    *background* do not show up in the polydispersity structure so 
     41    they are easy to handle. 
     42 
     43* *volfraction*: 
     44    Always the first parameter of S, but it may also be in P. If it is in P, 
     45    then *P.volfraction* is used in the combined P@S list, and 
     46    *S.volfraction* is elided, otherwise *S.volfraction* is used. If we 
     47    are using *volfraction* from P we can treat it like all the other P 
     48    parameters when calling P, but when calling S we need to insert the 
     49    *P.volfraction* into data vector for S and assign a slot of length 1 
     50    in the distribution. Because we are using the original layout of the 
     51    distribution vectors from P@S, but copying it into private data 
     52    vectors for S and P, we are free to "borrow" a P slots to store the 
     53    missing *S.volfraction* distribution.  We use the *P.volfraction* 
     54    slot itself but any slot will work. 
     55 
     56    For hollow shapes, *volfraction* represents the volume fraction of 
     57    material but S needs the volume fraction enclosed by the shape. The 
     58    answer is to scale the user specified volume fraction by the form:shell 
     59    ratio computed from the average form volume and average shell volume 
     60    returned from P. Use the original *volfraction* divided by *shell_volume* 
     61    to compute the number density, and scale P@S by that to get absolute 
     62    scaling on the final *I(q)*. The *scale* for P@S should therefore usually 
     63    be one. 
     64 
     65* *radius_effective*: 
     66    Always the second parameter of S and always part of P@S, but never in P. 
     67    The value may be calculated using *P.radius_effective()* or it 
     68    may be set to the *radius_effective* value in P@S, depending on 
     69    *radius_effective_mode*.  If part of S, the value may be polydisperse. 
     70    If calculated by P, then it will be the weighted average of effective 
     71    radii computed for the polydisperse shape parameters. 
     72 
     73* *structure_factor_mode* 
     74    If P@S supports beta approximation (i.e., if it has the *Fq* function that 
     75    returns <FF*> and <F><F*>), then *structure_factor_mode* will be added 
     76    to the P@S parameters right after the S parameters.  This mode may be 0 
     77    for the monodisperse approximation or 1 for the beta approximation.  We 
     78    will add more values here as we implemented more complicated operations, 
     79    but for now P and S must be computed separately.  If beta, then we 
     80    return *I = scale volfrac/volume ( <FF> + <F>^2 (S-1)) + background*. 
     81    If not beta then return *I = scale/volume P S + background* .  In both 
     82    cases, return the appropriate immediate values. 
     83 
     84* *radius_effective_mode* 
     85    If P defines the *radius_effective* function (and therefore 
     86    *P.info.radius_effective_modes* is a list of effective radius modes), 
     87    then *radius_effective_mode* will be the final parameter in P@S.  Mode 
     88    will be zero if *radius_effective* is defined by the user using the S 
     89    parameter; any other value and the *radius_effective* parameter will be 
     90    filled in from the value computed in P.  In the latter case, the 
     91    polydispersity information for *S.radius_effective* will need to be 
     92    suppressed, with pd length set to 1, the first value set to the 
     93    effective radius and the first weight set to 1.  Do this after composing 
     94    the S data vector so the inputs are left untouched. 
     95 
     96* *regular parameters* 
     97    The regular P parameters form a block of length *P.info.npars* at the 
     98    start of the data vector (after scale and background).  These will be 
     99    followed by *S.effective_radius*, and *S.volfraction* (if *P.volfraction* 
     100    is absent), and then the regular S parameters.  The P and S blocks can 
     101    be copied as a group into the respective P and S data vectors. 
     102    We can copy the distribution value and weight vectors untouched to both 
     103    the P and S data vectors since they are referenced by offset and length. 
     104    We can update the radius_effective slots in the P data vector with 
     105    *P.radius_effective()* if needed. 
     106 
     107* *magnetic parameters* 
     108    For each P parameter that is an SLD there will be a set of three magnetic 
     109    parameters tacked on to P@S after the regular P and S and after the 
     110    special *structure_factor_mode* and *radius_effective_mode*.  These 
     111    can be copied as a group after the regular P parameters.  There won't 
     112    be any magnetic S parameters. 
     113 
    12114""" 
    13115from __future__ import print_function, division 
    … …  
    84186    if not s_info.parameters.magnetism_index == []: 
    85187        raise TypeError("S should not have SLD parameters") 
     188    if RADIUS_ID in p_info.parameters: 
     189        raise TypeError("P should not have {}".format(RADIUS_ID)) 
    86190    p_id, p_name, p_pars = p_info.id, p_info.name, p_info.parameters 
    87191    s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters 
    88  
    89     # Create list of parameters for the combined model.  If there 
    90     # are any names in P that overlap with those in S, modify the name in S 
    91     # to distinguish it. 
     192    p_has_volfrac = VOLFRAC_ID in p_info.parameters 
     193 
     194    # Create list of parameters for the combined model.  If a name in 
     195    # S overlaps a name in P, tag the S parameter name to distinguish it. 
     196    # If the tagged name also collides it will be caught by the parameter 
     197    # table builder.  Similarly if any special names are abused.  Need the 
     198    # pairs to create the translation table for random model generation. 
    92199    p_set = set(p.id for p in p_pars.kernel_parameters) 
    93     s_list = [(_tag_parameter(par) if par.id in p_set else par) 
    94               for par in s_pars.kernel_parameters] 
    95     # Check if still a collision after renaming.  This could happen if for 
    96     # example S has volfrac and P has both volfrac and volfrac_S. 
    97     if any(p.id in p_set for p in s_list): 
    98         raise TypeError("name collision: P has P.name and P.name_S while S has S.name") 
    99  
    100     # make sure effective radius is not a polydisperse parameter in product 
    101     s_list[0] = copy(s_list[0]) 
    102     s_list[0].polydisperse = False 
    103  
    104     translate_name = dict((old.id, new.id) for old, new 
    105                           in zip(s_pars.kernel_parameters, s_list)) 
     200    s_pairs = [(par, (_tag_parameter(par) if par.id in p_set else par)) 
     201               for par in s_pars.kernel_parameters 
     202               # Note: exclude volfraction from s_list if volfraction in p 
     203               if par.id != VOLFRAC_ID or not p_has_volfrac] 
     204    s_list = [pair[0] for pair in s_pairs] 
     205 
     206    # Build combined parameter table 
    106207    combined_pars = p_pars.kernel_parameters + s_list + make_extra_pars(p_info) 
    107208    parameters = ParameterTable(combined_pars) 
    108     parameters.max_pd = p_pars.max_pd + s_pars.max_pd 
     209    # Allow for the scenario in which each component has all its PD parameters 
     210    # active simultaneously.  details.make_details() will throw an error if 
     211    # too many are used from any one component. 
     212    parameters.Pumax_pd = p_pars.max_pd + s_pars.max_pd 
     213 
     214    # TODO: does user-defined polydisperse S.radius_effective make sense? 
     215    # make sure effective radius is not a polydisperse parameter in product 
     216    #s_list[0] = copy(s_list[0]) 
     217    #s_list[0].polydisperse = False 
     218 
     219    s_translate = {old.id: new.id for old, new in s_pairs} 
    109220    def random(): 
    110221        """Random set of model parameters for product model""" 
    111222        combined_pars = p_info.random() 
    112         s_names = set(par.id for par in s_pars.kernel_parameters) 
    113         combined_pars.update((translate_name[k], v) 
     223        combined_pars.update((s_translate[k], v) 
    114224                             for k, v in s_info.random().items() 
    115                              if k in s_names) 
     225                             if k in s_translate) 
    116226        return combined_pars 
    117227 
    … …  
    173283 
    174284def _intermediates( 
    175         F1,               # type: np.ndarray 
    176         F2,               # type: np.ndarray 
     285        F,                # type: np.ndarray 
     286        Fsq,              # type: np.ndarray 
    177287        S,                # type: np.ndarray 
    178288        scale,            # type: float 
    … …  
    189299        # TODO: 2. consider implications if there are intermediate results in P(Q) 
    190300        parts = OrderedDict(( 
    191             ("P(Q)", scale*F2), 
     301            ("P(Q)", scale*Fsq), 
    192302            ("S(Q)", S), 
    193             ("beta(Q)", F1**2 / F2), 
    194             ("S_eff(Q)", 1 + (F1**2 / F2)*(S-1)), 
     303            ("beta(Q)", F**2 / Fsq), 
     304            ("S_eff(Q)", 1 + (F**2 / Fsq)*(S-1)), 
    195305            ("effective_radius", radius_effective), 
    196             # ("I(Q)", scale*(F2 + (F1**2)*(S-1)) + bg), 
     306            # ("I(Q)", scale*(Fsq + (F**2)*(S-1)) + bg), 
    197307        )) 
    198308    else: 
    199309        parts = OrderedDict(( 
    200             ("P(Q)", scale*F2), 
     310            ("P(Q)", scale*Fsq), 
    201311            ("S(Q)", S), 
    202312            ("effective_radius", radius_effective), 
    … …  
    261371        self.results = []  # type: List[np.ndarray] 
    262372 
     373        # Find index of volfraction parameter in parameter list 
     374        for k, p in enumerate(model_info.parameters.call_parameters): 
     375            if p.id == VOLFRAC_ID: 
     376                self._volfrac_index = k 
     377                break 
     378        else: 
     379            raise RuntimeError("no %s parameter in %s"%(VOLFRAC_ID, self)) 
     380 
     381        p_info, s_info = self.info.composition[1] 
     382        p_npars = p_info.parameters.npars 
     383        s_npars = s_info.parameters.npars 
     384 
     385        have_beta_mode = p_info.have_Fq 
     386        have_er_mode = p_info.radius_effective_modes is not None 
     387        volfrac_in_p = self._volfrac_index < p_npars + 2  # scale & background 
     388 
     389        # Slices into the details length/offset structure for P@S. 
     390        # Made complicated by the possibly missing volfraction in S. 
     391        self._p_detail_slice = slice(0, p_npars) 
     392        self._s_detail_slice = slice(p_npars, p_npars+s_npars-volfrac_in_p) 
     393        self._volfrac_in_p = volfrac_in_p 
     394 
     395        # P block from data vector, without scale and background 
     396        first_p = 2 
     397        last_p = p_npars + 2 
     398        self._p_value_slice = slice(first_p, last_p) 
     399 
     400        # radius_effective is the first parameter in S from the data vector. 
     401        self._er_index = last_p 
     402 
     403        # S block from data vector, without scale, background, volfrac or er. 
     404        first_s = last_p + 2 - volfrac_in_p 
     405        last_s = first_s + s_npars - 2 
     406        self._s_value_slice = slice(first_s, last_s) 
     407 
     408        # S distribution block in S data vector starts after all S values 
     409        self._s_dist_slice = slice(2 + s_npars, None) 
     410 
     411        # structure_factor_mode is the first parameter after P and S.  Skip 
     412        # 2 for scale and background, and subtract 1 in case there is no 
     413        # volfraction in S. 
     414        self._beta_mode_index = last_s if have_beta_mode else 0 
     415 
     416        # radius_effective_mode is the second parameter after P and S 
     417        # unless structure_factor_mode isn't available, in which case it 
     418        # is first. 
     419        self._er_mode_index = last_s + have_beta_mode if have_er_mode else 0 
     420 
     421        # Magnetic parameters are after everything else.  If they exist, 
     422        # they will only be for form factor P, not structure factor S. 
     423        first_mag = last_s + have_beta_mode + have_er_mode 
     424        mag_pars = 3*p_info.parameters.nmagnetic 
     425        last_mag = first_mag + (mag_pars + 3 if mag_pars else 0) 
     426        self._magentic_slice = slice(first_mag, last_mag) 
     427 
    263428    def Iq(self, call_details, values, cutoff, magnetic): 
    264429        # type: (CallDetails, np.ndarray, float, bool) -> np.ndarray 
    265  
    266430        p_info, s_info = self.info.composition[1] 
    267         p_npars = p_info.parameters.npars 
    268         p_length = call_details.length[:p_npars] 
    269         p_offset = call_details.offset[:p_npars] 
    270         s_npars = s_info.parameters.npars 
    271         s_length = call_details.length[p_npars:p_npars+s_npars] 
    272         s_offset = call_details.offset[p_npars:p_npars+s_npars] 
    273  
    274         # Beta mode parameter is the first parameter after P and S parameters 
    275         have_beta_mode = p_info.have_Fq 
    276         beta_mode_offset = 2+p_npars+s_npars 
    277         beta_mode = (values[beta_mode_offset] > 0) if have_beta_mode else False 
    278         if beta_mode and self.p_kernel.dim == '2d': 
    279             raise NotImplementedError("beta not yet supported for 2D") 
    280  
    281         # R_eff type parameter is the second parameter after P and S parameters 
    282         # unless the model doesn't support beta mode, in which case it is first 
    283         have_radius_type = p_info.radius_effective_modes is not None 
    284         #print(p_npars,s_npars) 
    285         radius_type_offset = 2+p_npars+s_npars + (1 if have_beta_mode else 0) 
    286         #print(values[radius_type_offset]) 
    287         radius_type = int(values[radius_type_offset]) if have_radius_type else 0 
    288  
    289         # Retrieve the volume fraction, which is the second of the 
    290         # 'S' parameters in the parameter list, or 2+np in 0-origin, 
    291         # as well as the scale and background. 
    292         volfrac = values[3+p_npars] 
     431 
     432        # Retrieve values from the data vector 
    293433        scale, background = values[0], values[1] 
    294  
    295         # if there are magnetic parameters, they will only be on the 
    296         # form factor P, not the structure factor S. 
    297         nmagnetic = len(self.info.parameters.magnetism_index) 
    298         if nmagnetic: 
    299             spin_index = self.info.parameters.npars + 2 
    300             magnetism = values[spin_index: spin_index+3+3*nmagnetic] 
    301         else: 
    302             magnetism = [] 
     434        volfrac = values[self._volfrac_index] 
     435        er_mode = (int(values[self._er_mode_index]) 
     436                   if self._er_mode_index > 0 else 0) 
     437        beta_mode = (values[self._beta_mode_index] > 0 
     438                     if self._beta_mode_index > 0 else False) 
     439 
    303440        nvalues = self.info.parameters.nvalues 
    304441        nweights = call_details.num_weights 
    305442        weights = values[nvalues:nvalues + 2*nweights] 
    306443 
     444        # Can't do 2d and beta_mode just yet 
     445        if beta_mode and self.p_kernel.dim == '2d': 
     446            raise NotImplementedError("beta not yet supported for 2D") 
     447 
    307448        # Construct the calling parameters for P. 
     449        p_length = call_details.length[self._p_detail_slice] 
     450        p_offset = call_details.offset[self._p_detail_slice] 
    308451        p_details = make_details(p_info, p_length, p_offset, nweights) 
    309452        p_values = [ 
    310453            [1., 0.], # scale=1, background=0, 
    311             values[2:2+p_npars], 
    312             magnetism, 
     454            values[self._p_value_slice], 
     455            values[self._magentic_slice], 
    313456            weights] 
    314457        spacer = (32 - sum(len(v) for v in p_values)%32)%32 
    … …  
    316459        p_values = np.hstack(p_values).astype(self.p_kernel.dtype) 
    317460 
     461        # Call the form factor kernel to compute <F> and <F^2>. 
     462        # If the model doesn't support Fq the returned <F> will be None. 
     463        F, Fsq, radius_effective, shell_volume, volume_ratio \ 
     464            = self.p_kernel.Fq(p_details, p_values, cutoff, magnetic, er_mode) 
     465 
     466        # TODO: async call to the GPU 
     467 
    318468        # Construct the calling parameters for S. 
    319         if radius_type > 0: 
    320             # If R_eff comes from form factor, make sure it is monodisperse. 
    321             # weight is set to 1 later, after the value array is created 
     469        s_length = call_details.length[self._s_detail_slice] 
     470        s_offset = call_details.offset[self._s_detail_slice] 
     471        if self._volfrac_in_p: 
     472            # Volfrac is in P and missing from S so insert a slot for it.  Say 
     473            # the distribution is length 1 and use the slot for volfraction 
     474            # from the P distribution. 
     475            s_length = np.insert(s_length, 1, 1) 
     476            s_offset = np.insert(s_offset, 1, p_offset[self._volfrac_index - 2]) 
     477        if er_mode > 0: 
     478            # If effective_radius comes from P, make sure it is monodisperse. 
     479            # Weight is set to 1 later, after the value array is created 
    322480            s_length[0] = 1 
    323481        s_details = make_details(s_info, s_length, s_offset, nweights) 
    324482        s_values = [ 
    325             [1., 0.], # scale=1, background=0, 
    326             values[2+p_npars:2+p_npars+s_npars], 
     483            [1., # scale=1 
     484             0., # background=0, 
     485             values[self._er_index], # S.radius_effective; may be replaced by P 
     486             0.], # volfraction; will be replaced by volfrac * volume_ratio 
     487            # followed by S parameters after effective_radius and volfraction 
     488            values[self._s_value_slice], 
    327489            weights, 
    328490        ] 
    … …  
    331493        s_values = np.hstack(s_values).astype(self.s_kernel.dtype) 
    332494 
    333         # Call the form factor kernel to compute <F> and <F^2>. 
    334         # If the model doesn't support Fq the returned <F> will be None. 
    335         F1, F2, radius_effective, shell_volume, volume_ratio = self.p_kernel.Fq( 
    336             p_details, p_values, cutoff, magnetic, radius_type) 
    337  
    338         # Call the structure factor kernel to compute S. 
    339495        # Plug R_eff from the form factor into structure factor parameters 
    340496        # and scale volume fraction by form:shell volume ratio. These changes 
    … …  
    344500        #print("R_eff=%d:%g, volfrac=%g, volume ratio=%g" 
    345501        #      % (radius_type, radius_effective, volfrac, volume_ratio)) 
    346         if radius_type > 0: 
     502        s_dist = s_values[self._s_dist_slice] 
     503        if er_mode > 0: 
    347504            # set the value to the model R_eff and set the weight to 1 
    348             s_values[2] = s_values[2+s_npars+s_offset[0]] = radius_effective 
    349             s_values[2+s_npars+s_offset[0]+nweights] = 1.0 
    350         s_values[3] = s_values[2+s_npars+s_offset[1]] = volfrac*volume_ratio 
     505            s_values[2] = s_dist[s_offset[0]] = radius_effective 
     506            s_dist[s_offset[0]+nweights] = 1.0 
     507        s_values[3] = s_dist[s_offset[1]] = volfrac*volume_ratio 
     508        s_dist[s_offset[1]+nweights] = 1.0 
     509 
     510        # Call the structure factor kernel to compute S. 
    351511        S = self.s_kernel.Iq(s_details, s_values, cutoff, False) 
     512        #print("P", Fsq[:10]) 
     513        #print("S", S[:10]) 
     514        #print(radius_effective, volfrac*volume_ratio) 
     515 
     516        # Combine form factor and structure factor 
     517        #print("beta", beta_mode, F, Fsq, S) 
     518        PS = Fsq + F**2*(S-1) if beta_mode else Fsq*S 
    352519 
    353520        # Determine overall scale factor. Hollow shapes are weighted by 
    354         # shell_volume, so that is needed for volume normalization.  For 
    355         # solid shapes we can use shell_volume as well since it is equal 
    356         # to form volume. 
    357         combined_scale = scale*volfrac/shell_volume 
    358  
    359         # Combine form factor and structure factor 
    360         #print("beta", beta_mode, F1, F2, S) 
    361         PS = F2 + F1**2*(S-1) if beta_mode else F2*S 
    362         final_result = combined_scale*PS + background 
     521        # shell_volume, so that is needed for number density estimation. 
     522        # For solid shapes we can use shell_volume as well since it is 
     523        # equal to form volume.  If P already has a volfraction parameter, 
     524        # then assume that it is already on absolute scale, and don't 
     525        # include volfrac in the combined_scale. 
     526        combined_scale = scale*(volfrac if not self._volfrac_in_p else 1.0) 
     527        final_result = combined_scale/shell_volume*PS + background 
    363528 
    364529        # Capture intermediate values so user can see them.  These are 
    … …  
    372537        # the results directly rather than through a lazy evaluator. 
    373538        self.results = lambda: _intermediates( 
    374             F1, F2, S, combined_scale, radius_effective, beta_mode) 
     539            F, Fsq, S, combined_scale, radius_effective, beta_mode) 
    375540 
    376541        return final_result 
Note: See TracChangeset for help on using the changeset viewer.