Changeset 5e2fd41 in sasmodels


Ignore:
Timestamp:
Mar 9, 2017 4:46:17 PM (7 years ago)
Author:
GitHub <noreply@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
4f9e288
Parents:
0011ecc (diff), 9c44b7b (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.
git-author:
Paul Kienzle <pkienzle@…> (03/09/17 16:46:17)
git-committer:
GitHub <noreply@…> (03/09/17 16:46:17)
Message:

Merge pull request #27 from SasView?/ticket-843

ticket-843/870: fix multilayer vesicle polydispersity, docs and parameter names; fixes P*S when parameter names collide. Fixes 843, 870.

Location:
sasmodels
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/conversion_table.py

    r790b16bb r0c2da4b  
    549549            "radius": "core_radius", 
    550550            "sld_solvent": "core_sld", 
    551             "n_pairs": "n_pairs", 
     551            "n_shells": "n_pairs", 
    552552            "thick_shell": "s_thickness", 
    553553            "sld": "shell_sld", 

  • sasmodels/modelinfo.py

    r5124c969 r9c44b7b  
    230230    defined as a sublist with the following elements: 
    231231 
    232     *name* is the name that will be used in the call to the kernel 
    233     function and the name that will be displayed to the user.  Names 
     232    *name* is the name that will be displayed to the user.  Names 
    234233    should be lower case, with words separated by underscore.  If 
    235     acronyms are used, the whole acronym should be upper case. 
     234    acronyms are used, the whole acronym should be upper case. For vector 
     235    parameters, the name will be followed by *[len]* where *len* is an 
     236    integer length of the vector, or the name of the parameter which 
     237    controls the length.  The attribute *id* will be created from name 
     238    without the length. 
    236239 
    237240    *units* should be one of *degrees* for angles, *Ang* for lengths, 
    … …  
    603606        # Using the call_parameters table, we already have expanded forms 
    604607        # for each of the vector parameters; put them in a lookup table 
    605         expanded_pars = dict((p.name, p) for p in self.call_parameters) 
     608        # Note: p.id and p.name are currently identical for the call parameters 
     609        expanded_pars = dict((p.id, p) for p in self.call_parameters) 
    606610 
    607611        def append_group(name): 

  • sasmodels/models/multilayer_vesicle.c

    rc3ccaec rec1d4bc  
    1 static 
    2 double multilayer_vesicle_kernel(double q, 
     1static double 
     2form_volume(double radius, 
     3          double thick_shell, 
     4          double thick_solvent, 
     5          double fp_n_shells) 
     6{ 
     7    int n_shells = (int)(fp_n_shells + 0.5); 
     8    double R_N = radius + n_shells*(thick_shell+thick_solvent) - thick_solvent; 
     9    return M_4PI_3*cube(R_N); 
     10} 
     11 
     12static double 
     13multilayer_vesicle_kernel(double q, 
    314          double volfraction, 
    415          double radius, 
    … …  
    718          double sld_solvent, 
    819          double sld, 
    9           int n_pairs) 
     20          int n_shells) 
    1021{ 
    1122    //calculate with a loop, two shells at a time 
    … …  
    2940 
    3041        //do 2 layers at a time 
    31         ii += 1; 
     42        ii++; 
    3243 
    33     } while(ii <= n_pairs-1);  //change to make 0 < n_pairs < 2 correspond to 
     44    } while(ii <= n_shells-1);  //change to make 0 < n_shells < 2 correspond to 
    3445                               //unilamellar vesicles (C. Glinka, 11/24/03) 
    3546 
    36     fval *= volfraction*1.0e-4*fval/voli; 
    37  
    38     return(fval); 
     47    return 1.0e-4*volfraction*fval*fval;  // Volume normalization happens in caller 
    3948} 
    4049 
    41 static 
    42 double Iq(double q, 
     50static double 
     51Iq(double q, 
    4352          double volfraction, 
    4453          double radius, 
    … …  
    4756          double sld_solvent, 
    4857          double sld, 
    49           double fp_n_pairs) 
     58          double fp_n_shells) 
    5059{ 
    51     int n_pairs = (int)(fp_n_pairs + 0.5); 
     60    int n_shells = (int)(fp_n_shells + 0.5); 
    5261    return multilayer_vesicle_kernel(q, 
    5362           volfraction, 
    … …  
    5766           sld_solvent, 
    5867           sld, 
    59            n_pairs); 
     68           n_shells); 
    6069} 
    6170 

  • sasmodels/models/multilayer_vesicle.py

    rc3ccaec r68f45cb  
    55This model is a trivial extension of the core_shell_sphere function to include 
    66*N* shells where the core is filled with solvent and the shells are interleaved 
    7 with layers of solvent. For *N = 1*, this returns the same as the vesicle model, 
     7with layers of solvent. For $N = 1$, this returns the same as the vesicle model, 
    88except for the normalisation, which here is to outermost volume. 
    99The shell thicknessess and SLD are constant for all shells as expected for 
    … …  
    1919 
    2020.. math:: 
    21     P(q) = \text{scale} \cdot \frac{V_f}{V_t} F^2(q) + \text{background} 
     21    P(q) = \text{scale} \cdot \frac{\phi}{V(R_N)} F^2(q) + \text{background} 
     22 
     23where 
     24 
     25.. math:: 
     26     F(q) = (\rho_\text{shell}-\rho_\text{solv}) \sum_{i=1}^{N} \left[ 
     27     3V(r_i)\frac{\sin(qr_i) - qr_i\cos(qr_i)}{(qr_i)^3} 
     28     - 3V(R_i)\frac{\sin(qR_i) - qR_i\cos(qR_i)}{(qR_i)^3} 
     29     \right] 
    2230 
    2331for 
    2432 
    2533.. math:: 
    26     F(q) = (\rho_\text{shell}-\rho_\text{solv}) \sum_{i=1}^{n_\text{pairs}} 
    27         \left[ 
    28           3V(R_i)\frac{\sin(qR_i)-qR_i\cos(qR_i)}{(qR_i)^3} \\ 
    29           - 3V(R_i+t_s)\frac{\sin(q(R_i+t_s))-q(R_i+t_s)\cos(q(R_i+t_s))}{(q(R_i+t_s))^3} 
    30         \right] 
    3134 
    32 and 
     35     r_i &= r_c + (i-1)(t_s + t_w) && \text{ solvent radius before shell } i \\ 
     36     R_i &= r_i + t_s && \text{ shell radius for shell } i 
    3337 
    34 .. math:: 
    35      R_i = r_c + (i-1)(t_s + t_w) 
     38$\phi$ is the volume fraction of particles, $V(r)$ is the volume of a sphere 
     39of radius $r$, $r_c$ is the radius of the core, $t_s$ is the thickness of 
     40the shell, $t_w$ is the thickness of the solvent layer between the shells, 
     41$\rho_\text{shell}$ is the scattering length density of a shell, and 
     42$\rho_\text{solv}$ is the scattering length density of the solvent. 
    3643 
    37 where $V_f$ is the volume fraction of particles, $V_t$ is the volume of the 
    38 whole particle, $V(r)$ is the volume of a sphere of radius $r$, $r_c$ is the 
    39 radius of the core, $\rho_\text{shell}$ is the scattering length density of a 
    40 shell, $\rho_\text{solv}$ is the scattering length density of the solvent. 
     44The outer-most shell radius $R_N$ is used as the effective radius 
     45for $P(Q)$ when $P(Q) * S(Q)$ is applied. 
    4146 
    42 The outer most radius, $r_o = R_n + t_s$, is used for both the volume fraction 
    43 normalization and for the effective radius for *S(Q)* when $P(Q) * S(Q)$ 
    44 is applied. 
     47For mixed systems in which some vesicles have 1 shell, some have 2, 
     48etc., use polydispersity on $N$ to model the data.  For example, 
     49create a file such as *shell_dist.txt* containing the relative portion 
     50of each vesicle size:: 
     51 
     52    1 20 
     53    2  4 
     54    3  1 
     55 
     56Turn on polydispersity and select an array distribution for the *n_shells* 
     57parameter.  Choose the above *shell_dist.txt* file, and the model will be 
     58computed with 80% 1-shell vesicles, 16% 2-shell vesicles and 4% 
     593-shell vesicles. 
    4560 
    4661The 2D scattering intensity is the same as 1D, regardless of the orientation 
    … …  
    86101    sld_solvent: solvent scattering length density 
    87102    sld: shell scattering length density 
    88     n_pairs:number of "shell plus solvent" layer pairs 
     103    n_shells:number of "shell plus solvent" layer pairs 
    89104    background: incoherent background 
    90105        """ 
    … …  
    95110parameters = [ 
    96111    ["volfraction", "",  0.05, [0.0, 1],  "", "volume fraction of vesicles"], 
    97     ["radius", "Ang", 60.0, [0.0, inf],  "", "radius of solvent filled core"], 
    98     ["thick_shell", "Ang",        10.0, [0.0, inf],  "", "thickness of one shell"], 
    99     ["thick_solvent", "Ang",        10.0, [0.0, inf],  "", "solvent thickness between shells"], 
     112    ["radius", "Ang", 60.0, [0.0, inf],  "volume", "radius of solvent filled core"], 
     113    ["thick_shell", "Ang",        10.0, [0.0, inf],  "volume", "thickness of one shell"], 
     114    ["thick_solvent", "Ang",        10.0, [0.0, inf],  "volume", "solvent thickness between shells"], 
    100115    ["sld_solvent",    "1e-6/Ang^2",  6.4, [-inf, inf], "sld", "solvent scattering length density"], 
    101116    ["sld",   "1e-6/Ang^2",  0.4, [-inf, inf], "sld", "Shell scattering length density"], 
    102     ["n_pairs",     "",            2.0, [1.0, inf],  "", "Number of shell plus solvent layer pairs"], 
     117    ["n_shells",     "",            2.0, [1.0, inf],  "volume", "Number of shell plus solvent layer pairs"], 
    103118    ] 
    104119# pylint: enable=bad-whitespace, line-too-long 
    105120 
     121# TODO: proposed syntax for specifying which parameters can be polydisperse 
     122#polydispersity = ["radius", "thick_shell"] 
     123 
    106124source = ["lib/sas_3j1x_x.c", "multilayer_vesicle.c"] 
    107125 
    108 # TODO: the following line does nothing 
    109 polydispersity = ["radius", "n_pairs"] 
     126def ER(radius, thick_shell, thick_solvent, n_shells): 
     127    n_shells = int(n_shells+0.5) 
     128    return radius + n_shells * (thick_shell + thick_solvent) - thick_solvent 
    110129 
    111130demo = dict(scale=1, background=0, 
    … …  
    116135            sld_solvent=6.4, 
    117136            sld=0.4, 
    118             n_pairs=2.0) 
     137            n_shells=2.0) 
    119138 
    120139tests = [ 
    … …  
    125144      'sld_solvent': 6.4, 
    126145      'sld': 0.4, 
    127       'n_pairs': 2.0, 
     146      'n_shells': 2.0, 
    128147      'scale': 1.0, 
    129148      'background': 0.001, 
    … …  
    136155      'sld_solvent': 6.4, 
    137156      'sld': 0.4, 
    138       'n_pairs': 2.0, 
     157      'n_shells': 2.0, 
    139158      'scale': 1.0, 
    140159      'background': 0.001, 

  • sasmodels/product.py

    r9951a86 rf88e248  
    4545    # structure factor calculator.  Structure factors should not 
    4646    # have any magnetic parameters 
    47     assert(s_info.parameters.kernel_parameters[0].id == ER_ID) 
    48     assert(s_info.parameters.kernel_parameters[1].id == VF_ID) 
    49     assert(s_info.parameters.magnetism_index == []) 
     47    if not s_info.parameters.kernel_parameters[0].id == ER_ID: 
     48        raise TypeError("S needs %s as first parameter"%ER_ID) 
     49    if not s_info.parameters.kernel_parameters[1].id == VF_ID: 
     50        raise TypeError("S needs %s as second parameter"%VF_ID) 
     51    if not s_info.parameters.magnetism_index == []: 
     52        raise TypeError("S should not have SLD parameters") 
    5053    p_id, p_name, p_pars = p_info.id, p_info.name, p_info.parameters 
    5154    s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters 
    52     p_set = set(p.id for p in p_pars.call_parameters) 
    53     s_set = set(p.id for p in s_pars.call_parameters) 
    54  
    55     if p_set & s_set: 
    56         # there is some overlap between the parameter names; tag the 
    57         # overlapping S parameters with name_S. 
    58         # Skip the first parameter of s, which is effective radius 
    59         s_list = [(suffix_parameter(par) if par.id in p_set else par) 
    60                   for par in s_pars.kernel_parameters[1:]] 
    61     else: 
    62         # Skip the first parameter of s, which is effective radius 
    63         s_list = s_pars.kernel_parameters[1:] 
     55 
     56    # Create list of parameters for the combined model.  Skip the first 
     57    # parameter of S, which we verified above is effective radius.  If there 
     58    # are any names in P that overlap with those in S, modify the name in S 
     59    # to distinguish it. 
     60    p_set = set(p.id for p in p_pars.kernel_parameters) 
     61    s_list = [(_tag_parameter(par) if par.id in p_set else par) 
     62              for par in s_pars.kernel_parameters[1:]] 
     63    # Check if still a collision after renaming.  This could happen if for 
     64    # example S has volfrac and P has both volfrac and volfrac_S. 
     65    if any(p.id in p_set for p in s_list): 
     66        raise TypeError("name collision: P has P.name and P.name_S while S has S.name") 
     67 
    6468    translate_name = dict((old.id, new.id) for old, new 
    6569                          in zip(s_pars.kernel_parameters[1:], s_list)) 
    6670    demo = {} 
    67     demo.update((k, v) for k, v in p_info.demo.items() 
    68                 if k not in ("background", "scale")) 
     71    demo.update(p_info.demo.items()) 
    6972    demo.update((translate_name[k], v) for k, v in s_info.demo.items() 
    7073                if k not in ("background", "scale") and not k.startswith(ER_ID)) 
    … …  
    9093    # Remember the component info blocks so we can build the model 
    9194    model_info.composition = ('product', [p_info, s_info]) 
    92     model_info.demo = {} 
     95    model_info.demo = demo 
     96 
     97    ## Show the parameter table with the demo values 
     98    #from .compare import get_pars, parlist 
     99    #print("==== %s ====="%model_info.name) 
     100    #values = get_pars(model_info, use_demo=True) 
     101    #print(parlist(model_info, values, is2d=True)) 
    93102    return model_info 
    94103 
    95 def suffix_parameter(par, suffix): 
     104def _tag_parameter(par): 
     105    """ 
     106    Tag the parameter name with _S to indicate that the parameter comes from 
     107    the structure factor parameter set.  This is only necessary if the 
     108    form factor model includes a parameter of the same name as a parameter 
     109    in the structure factor. 
     110    """ 
    96111    par = copy(par) 
    97     par.name = par.name + " S" 
     112    # Protect against a vector parameter in S by appending the vector length 
     113    # to the renamed parameter.  Note: haven't tested this since no existing 
     114    # structure factor models contain vector parameters. 
     115    vector_length = par.name[len(par.id):] 
    98116    par.id = par.id + "_S" 
     117    par.name = par.id + vector_length 
    99118    return par 
    100119 
Note: See TracChangeset for help on using the changeset viewer.