Changeset 26a6608 in sasmodels


Ignore:
Timestamp:
Nov 20, 2017 9:49:03 AM (6 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
a66b004
Parents:
a70959a (diff), fa70e04 (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 'master' into boltzmann

Files:
2 added
6 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/models/core_shell_parallelepiped.c

    r92dfe0c rc69d6d6  
    1 double form_volume(double length_a, double length_b, double length_c,  
     1double form_volume(double length_a, double length_b, double length_c, 
    22                   double thick_rim_a, double thick_rim_b, double thick_rim_c); 
    33double Iq(double q, double core_sld, double arim_sld, double brim_sld, double crim_sld, 
    … …  
    99            double thick_rim_c, double theta, double phi, double psi); 
    1010 
    11 double form_volume(double length_a, double length_b, double length_c,  
     11double form_volume(double length_a, double length_b, double length_c, 
    1212                   double thick_rim_a, double thick_rim_b, double thick_rim_c) 
    1313{ 
    1414    //return length_a * length_b * length_c; 
    15     return length_a * length_b * length_c +  
    16            2.0 * thick_rim_a * length_b * length_c +  
     15    return length_a * length_b * length_c + 
     16           2.0 * thick_rim_a * length_b * length_c + 
    1717           2.0 * thick_rim_b * length_a * length_c + 
    1818           2.0 * thick_rim_c * length_a * length_b; 
    … …  
    3434    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled 
    3535    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 
    36      
     36    //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 
     37 
    3738    const double mu = 0.5 * q * length_b; 
    38      
     39 
    3940    //calculate volume before rescaling (in original code, but not used) 
    40     //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);         
    41     //double vol = length_a * length_b * length_c +  
    42     //       2.0 * thick_rim_a * length_b * length_c +  
     41    //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
     42    //double vol = length_a * length_b * length_c + 
     43    //       2.0 * thick_rim_a * length_b * length_c + 
    4344    //       2.0 * thick_rim_b * length_a * length_c + 
    4445    //       2.0 * thick_rim_c * length_a * length_b; 
    45      
     46 
    4647    // Scale sides by B 
    4748    const double a_scaled = length_a / length_b; 
    4849    const double c_scaled = length_c / length_b; 
    4950 
    50     // ta and tb correspond to the definitions in CSPPKernel, but they don't make sense to me (MG) 
    51     // the a_scaled in the definition of tb was present in CSPPKernel in libCylinder.c, 
    52     // while in cspkernel in csparallelepiped.cpp (used for the 2D), all the definitions 
    53     // for ta, tb, tc use also A + 2*rim_thickness (but not scaled by B!!!) 
    54     double ta = (a_scaled + 2.0*thick_rim_a)/length_b; 
    55     double tb = (a_scaled + 2.0*thick_rim_b)/length_b; 
     51    double ta = a_scaled + 2.0*thick_rim_a/length_b; // incorrect ta = (a_scaled + 2.0*thick_rim_a)/length_b; 
     52    double tb = 1+ 2.0*thick_rim_b/length_b; // incorrect tb = (a_scaled + 2.0*thick_rim_b)/length_b; 
     53    double tc = c_scaled + 2.0*thick_rim_c/length_b; //not present 
    5654 
    5755    double Vin = length_a * length_b * length_c; 
    … …  
    6260    double V1 = (2.0 * thick_rim_a * length_b * length_c);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 
    6361    double V2 = (2.0 * length_a * thick_rim_b * length_c);    // incorrect V2(aa*bb*cc+2*aa*tb*cc) 
     62    double V3 = (2.0 * length_a * length_b * thick_rim_c);    //not present 
     63    double Vot = Vin + V1 + V2 + V3; 
    6464 
    6565    // Scale factors (note that drC is not used later) 
    … …  
    6767    const double drhoA = (arim_sld-solvent_sld); 
    6868    const double drhoB = (brim_sld-solvent_sld); 
    69     //const double drC_Vot = (crim_sld-solvent_sld)*Vot; 
     69    const double drhoC = (crim_sld-solvent_sld);  // incorrect const double drC_Vot = (crim_sld-solvent_sld)*Vot; 
     70 
    7071 
    7172    // Precompute scale factors for combining cross terms from the shape 
    7273    const double scale23 = drhoA*V1; 
    7374    const double scale14 = drhoB*V2; 
    74     const double scale12 = drho0*Vin - scale23 - scale14; 
     75    const double scale24 = drhoC*V3; 
     76    const double scale11 = drho0*Vin; 
     77    const double scale12 = drho0*Vin - scale23 - scale14 - scale24; 
    7578 
    7679    // outer integral (with gauss points), integration limits = 0, 1 
    … …  
    8386        // inner integral (with gauss points), integration limits = 0, 1 
    8487        double inner_total = 0.0; 
     88        double inner_total_crim = 0.0; 
    8589        for(int j=0; j<76; j++) { 
    8690            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
    … …  
    8892            SINCOS(M_PI_2*uu, sin_uu, cos_uu); 
    8993            const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled); 
    90             const double si2 = sas_sinx_x(mu_proj * cos_uu); 
     94            const double si2 = sas_sinx_x(mu_proj * cos_uu ); 
    9195            const double si3 = sas_sinx_x(mu_proj * sin_uu * ta); 
    9296            const double si4 = sas_sinx_x(mu_proj * cos_uu * tb); 
    … …  
    9498            // Expression in libCylinder.c (neither drC nor Vot are used) 
    9599            const double form = scale12*si1*si2 + scale23*si2*si3 + scale14*si1*si4; 
     100            const double form_crim = scale11*si1*si2; 
    96101 
    97             // To note also that in csparallelepiped.cpp, there is a function called 
    98             // cspkernel, which seems to make more sense and has the following comment: 
    99             //   This expression is different from NIST/IGOR package (I strongly believe the IGOR is wrong!!!). 10/15/2010. 
    100             //   tmp =( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3)* 
    101             //   ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3);   //  correct FF : square of sum of phase factors 
    102             // This is the function called by csparallelepiped_analytical_2D_scaled, 
    103             // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c         
    104              
     102 
    105103            //  correct FF : sum of square of phase factors 
    106104            inner_total += Gauss76Wt[j] * form * form; 
     105            inner_total_crim += Gauss76Wt[j] * form_crim * form_crim; 
    107106        } 
    108107        inner_total *= 0.5; 
    109  
     108        inner_total_crim *= 0.5; 
    110109        // now sum up the outer integral 
    111110        const double si = sas_sinx_x(mu * c_scaled * sigma); 
    112         outer_total += Gauss76Wt[i] * inner_total * si * si; 
     111        const double si_crim = sas_sinx_x(mu * tc * sigma); 
     112        outer_total += Gauss76Wt[i] * (inner_total * si * si + inner_total_crim * si_crim * si_crim); 
    113113    } 
    114114    outer_total *= 0.5; 
    … …  
    154154 
    155155    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no 
    156     // the scaling by B. The use of length_a for the 3 of them seems clearly a mistake to me, 
    157     // but for the moment I let it like this until understanding better the code. 
     156    // the scaling by B. 
    158157    double ta = length_a + 2.0*thick_rim_a; 
    159     double tb = length_a + 2.0*thick_rim_b; 
    160     double tc = length_a + 2.0*thick_rim_c; 
     158    double tb = length_b + 2.0*thick_rim_b; 
     159    double tc = length_c + 2.0*thick_rim_c; 
    161160    //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    162161    double siA = sas_sinx_x(0.5*q*length_a*xhat); 
    … …  
    166165    double siBt = sas_sinx_x(0.5*q*tb*yhat); 
    167166    double siCt = sas_sinx_x(0.5*q*tc*zhat); 
    168      
     167 
    169168 
    170169    // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed 
    … …  
    173172               + drA*(siAt-siA)*siB*siC*V1 
    174173               + drB*siA*(siBt-siB)*siC*V2 
    175                + drC*siA*siB*(siCt*siCt-siC)*V3); 
    176     
     174               + drC*siA*siB*(siCt-siC)*V3); 
     175 
    177176    return 1.0e-4 * f * f; 
    178177} 

  • sasmodels/models/core_shell_parallelepiped.py

    r8f04da4 rfa70e04  
    211211 
    212212# rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models 
    213 qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) 
    214 tests = [[{}, 0.2, 0.533149288477], 
    215          [{}, [0.2], [0.533149288477]], 
    216          [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0853299803222], 
    217          [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0853299803222]], 
    218         ] 
    219 del qx, qy  # not necessary to delete, but cleaner 
     213if 0:  # pak: model rewrite; need to update tests 
     214    qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) 
     215    tests = [[{}, 0.2, 0.533149288477], 
     216            [{}, [0.2], [0.533149288477]], 
     217            [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0853299803222], 
     218            [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0853299803222]], 
     219            ] 
     220    del qx, qy  # not necessary to delete, but cleaner 

  • sasmodels/product.py

    r058460c r146793b  
    100100    # Remember the component info blocks so we can build the model 
    101101    model_info.composition = ('product', [p_info, s_info]) 
     102    model_info.control = p_info.control 
     103    model_info.hidden = p_info.hidden 
     104    if getattr(p_info, 'profile', None) is not None: 
     105        profile_pars = set(p.id for p in p_info.parameters.kernel_parameters) 
     106        def profile(**kwargs): 
     107            # extract the profile args 
     108            kwargs = dict((k, v) for k, v in kwargs.items() if k in profile_pars) 
     109            return p_info.profile(**kwargs) 
     110    else: 
     111        profile = None 
     112    model_info.profile = profile 
     113    model_info.profile_axes = p_info.profile_axes 
     114 
    102115    # TODO: delegate random to p_info, s_info 
    103116    #model_info.random = lambda: {} 
    … …  
    129142    def __init__(self, model_info, P, S): 
    130143        # type: (ModelInfo, KernelModel, KernelModel) -> None 
     144        #: Combined info plock for the product model 
    131145        self.info = model_info 
     146        #: Form factor modelling individual particles. 
    132147        self.P = P 
     148        #: Structure factor modelling interaction between particles. 
    133149        self.S = S 
    134         self.dtype = P.dtype 
     150        #: Model precision. This is not really relevant, since it is the 
     151        #: individual P and S models that control the effective dtype, 
     152        #: converting the q-vectors to the correct type when the kernels 
     153        #: for each are created. Ideally this should be set to the more 
     154        #: precise type to avoid loss of precision, but precision in q is 
     155        #: not critical (single is good enough for our purposes), so it just 
     156        #: uses the precision of the form factor. 
     157        self.dtype = P.dtype  # type: np.dtype 
    135158 
    136159    def make_kernel(self, q_vectors): 

  • sasmodels/sasview_model.py

    r9f8ade1 redb0f85  
    205205                                           structure_factor._model_info) 
    206206    ConstructedModel = make_model_from_info(model_info) 
    207     return ConstructedModel() 
     207    return ConstructedModel(form_factor.multiplicity) 
    208208 
    209209 
    … …  
    323323    #: True if model has multiplicity 
    324324    is_multiplicity_model = False 
    325     #: Mulitplicity information 
     325    #: Multiplicity information 
    326326    multiplicity_info = None # type: MultiplicityInfoType 
    327327 
    … …  
    354354        # and lines to plot. 
    355355 
    356         # Get the list of hidden parameters given the mulitplicity 
     356        # Get the list of hidden parameters given the multiplicity 
    357357        # Don't include multiplicity in the list of parameters 
    358358        self.multiplicity = multiplicity 

  • doc/guide/pd/polydispersity.rst

    r1f058ea r75e4319  
    4040calculations are generally more robust with more data points or more angles. 
    4141 
    42 The following five distribution functions are provided: 
     42The following six distribution functions are provided: 
    4343 
    4444*  *Rectangular Distribution* 
     45*  *Uniform Distribution* 
    4546*  *Gaussian Distribution* 
    4647*  *Lognormal Distribution* 
    4748*  *Schulz Distribution* 
    4849*  *Array Distribution* 
     50*  *Boltzmann Distribution* 
    4951 
    5052These are all implemented as *number-average* distributions. 
    … …  
    8385    Rectangular distribution. 
    8486 
     87Uniform Distribution 
     88^^^^^^^^^^^^^^^^^^^^^^^^ 
     89 
     90The Uniform Distribution is defined as 
     91 
     92    .. math:: 
     93 
     94        f(x) = \frac{1}{\text{Norm}} 
     95        \begin{cases} 
     96          1 & \text{for } |x - \bar x| \leq \sigma \\ 
     97          0 & \text{for } |x - \bar x| > \sigma 
     98        \end{cases} 
     99 
     100    where $\bar x$ is the mean of the distribution, $\sigma$ is the half-width, and 
     101    *Norm* is a normalization factor which is determined during the numerical 
     102    calculation. 
     103 
     104    Note that the polydispersity is given by 
     105 
     106    .. math:: \text{PD} = \sigma / \bar x 
     107 
     108    .. figure:: pd_uniform.jpg 
     109 
     110        Uniform distribution. 
     111 
    85112.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
    86113 
    … …  
    181208^^^^^^^^^^^^^^^^^^ 
    182209 
    183 This user-definable distribution should be given as as a simple ASCII text 
     210This user-definable distribution should be given as a simple ASCII text 
    184211file where the array is defined by two columns of numbers: $x$ and $f(x)$. 
    185212The $f(x)$ will be normalized to 1 during the computation. 
    … …  
    200227given for the model will have no affect, and will be ignored when computing 
    201228the average.  This means that any parameter with an array distribution will 
    202 not be fittable. 
     229not be fitable. 
     230 
     231.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
     232 
     233Boltzmann Distribution 
     234^^^^^^^^^^^^^^^^^^^^^^ 
     235 
     236The Boltzmann Distribution is defined as 
     237 
     238.. math:: 
     239 
     240    f(x) = \frac{1}{\text{Norm}} 
     241           \exp\left(-\frac{ | x - \bar x | }{\sigma}\right) 
     242 
     243where $\bar x$ is the mean of the distribution and *Norm* is a normalization 
     244factor which is determined during the numerical calculation. 
     245The width is defined as 
     246 
     247.. math:: \sigma=\frac{k T}{E} 
     248 
     249which is the inverse Boltzmann factor, 
     250where $k$ is the Boltzmann constant, $T$ the temperature in Kelvin and $E$ a 
     251characteristic energy per particle. 
     252 
     253.. figure:: pd_boltzmann.jpg 
     254 
     255    Boltzmann distribution. 
    203256 
    204257.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 

  • sasmodels/weights.py

    r41e7f2e r75e4319  
    5555        """ 
    5656        sigma = self.width * center if relative else self.width 
     57        if not relative: 
     58            # For orientation, the jitter is relative to 0 not the angle 
     59            center = 0 
     60            pass 
    5761        if sigma == 0 or self.npts < 2: 
    5862            if lb <= center <= ub: 
    … …  
    9397        return x, px 
    9498 
     99class UniformDispersion(Dispersion): 
     100    r""" 
     101    Uniform dispersion, with width $\sigma$. 
     102 
     103    .. math:: 
     104 
     105        w = 1 
     106    """ 
     107    type = "uniform" 
     108    default = dict(npts=35, width=0, nsigmas=1) 
     109    def _weights(self, center, sigma, lb, ub): 
     110        x = self._linspace(center, sigma, lb, ub) 
     111        x = x[np.fabs(x-center) <= np.fabs(sigma)] 
     112        return x, np.ones_like(x) 
    95113 
    96114class RectangleDispersion(Dispersion): 
    … …  
    186204        return x, px 
    187205 
     206class BoltzmannDispersion(Dispersion): 
     207    r""" 
     208    Boltzmann dispersion, with $\sigma=k T/E$. 
     209 
     210    .. math:: 
     211 
     212        w = \exp\left( -|x - c|/\sigma\right) 
     213    """ 
     214    type = "boltzmann" 
     215    default = dict(npts=35, width=0, nsigmas=3) 
     216    def _weights(self, center, sigma, lb, ub): 
     217        x = self._linspace(center, sigma, lb, ub) 
     218        px = np.exp(-np.abs(x-center) / np.abs(sigma)) 
     219        return x, px 
    188220 
    189221# dispersion name -> disperser lookup table. 
    … …  
    192224MODELS = OrderedDict((d.type, d) for d in ( 
    193225    RectangleDispersion, 
     226    UniformDispersion, 
    194227    ArrayDispersion, 
    195228    LogNormalDispersion, 
    196229    GaussianDispersion, 
    197230    SchulzDispersion, 
     231    BoltzmannDispersion 
    198232)) 
    199233 
    … …  
    225259    obj = cls(n, width, nsigmas) 
    226260    v, w = obj.get_weights(value, limits[0], limits[1], relative) 
    227     return v, w 
    228  
    229  
    230 def plot_weights(model_info, pairs): 
    231     # type: (ModelInfo, List[Tuple[np.ndarray, np.ndarray]]) -> None 
     261    return v, w/np.sum(w) 
     262 
     263 
     264def plot_weights(model_info, mesh): 
     265    # type: (ModelInfo, List[Tuple[float, np.ndarray, np.ndarray]]) -> None 
    232266    """ 
    233267    Plot the weights returned by :func:`get_weights`. 
    234268 
    235     *model_info* is 
    236     :param model_info: 
    237     :param pairs: 
    238     :return: 
     269    *model_info* defines model parameters, etc. 
     270 
     271    *mesh* is a list of tuples containing (*value*, *dispersity*, *weights*) 
     272    for each parameter, where (*dispersity*, *weights*) pairs are the 
     273    distributions to be plotted. 
    239274    """ 
    240275    import pylab 
    241276 
    242     if any(len(values)>1 for values, weights in pairs): 
     277    if any(len(dispersity)>1 for value, dispersity, weights in mesh): 
    243278        labels = [p.name for p in model_info.parameters.call_parameters] 
    244         pylab.interactive(True) 
     279        #pylab.interactive(True) 
    245280        pylab.figure() 
    246         for (v,w), s in zip(pairs, labels): 
    247             if len(v) > 1: 
    248                 #print("weights for", s, v, w) 
    249                 pylab.plot(v, w, '-o', label=s) 
     281        for (v,x,w), s in zip(mesh, labels): 
     282            if len(x) > 1: 
     283                pylab.plot(x, w, '-o', label=s) 
    250284        pylab.grid(True) 
    251285        pylab.legend() 
Note: See TracChangeset for help on using the changeset viewer.