Changeset 9616dfe in sasmodels


Ignore:
Timestamp:
Mar 11, 2018 2:31:50 PM (7 years ago)
Author:
butler
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
331870d
Parents:
dbf1a60 (diff), bf94e6e (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 ticket-896

Files:
7 edited
1 moved

Legend:

Unmodified
Added
Removed

  • example/multiscatfit.py

    rbc248f8 r49d1f8b8  
    3434from sasmodels.data import load_data, set_beam_stop, set_top 
    3535 
    36 from multiscat import MultipleScattering 
     36from sasmodels.multiscat import MultipleScattering 
    3737 
    3838## Load the data 
    3939#data = load_data('DEC07267.DAT') 
    4040#set_beam_stop(data, 0.003, outer=0.025) 
    41 data = load_data('latex_smeared.xml', index=1) 
     41data = load_data('latex_smeared.xml', index=0) 
    4242 
    4343## Define the model 
    … …  
    6666 
    6767# Mulitple scattering probability parameter 
    68 # HACK: the parameter is assigned to model.theta, which is otherwise unused 
    69 # since the dataset is 1D; this won't work for 2D data 
     68# HACK: the probability is stuffed in as an extra parameter to the experiment. 
    7069probability = Parameter(name="probability", value=0.0) 
    7170probability.range(0.0, 0.9) 
    72 model.phi = probability 
    7371 
    74 M = Experiment(data=data, model=model) 
     72M = Experiment(data=data, model=model, extra_pars={'probability': probability}) 
    7573 
    7674# Stack mulitple scattering on top of the existing resolution function. 

  • sasmodels/bumps_model.py

    r2d81cfe r49d1f8b8  
    137137    """ 
    138138    _cache = None # type: Dict[str, np.ndarray] 
    139     def __init__(self, data, model, cutoff=1e-5, name=None): 
     139    def __init__(self, data, model, cutoff=1e-5, name=None, extra_pars=None): 
    140140        # type: (Data, Model, float) -> None 
    141141        # remember inputs so we can inspect from outside 
    … …  
    145145        self._interpret_data(data, model.sasmodel) 
    146146        self._cache = {} 
     147        self.extra_pars = extra_pars 
    147148 
    148149    def update(self): 
    … …  
    166167        Return a dictionary of parameters 
    167168        """ 
    168         return self.model.parameters() 
     169        pars = self.model.parameters() 
     170        if self.extra_pars: 
     171            pars.update(self.extra_pars) 
     172        return pars 
    169173 
    170174    def theory(self): 

  • sasmodels/kerneldll.py

    r2d81cfe rbf94e6e  
    158158        return CC + [source, "-o", output, "-lm"] 
    159159 
    160 # Windows-specific solution 
    161 if os.name == 'nt': 
    162     # Assume the default location of module DLLs is in .sasmodels/compiled_models. 
    163     DLL_PATH = os.path.join(os.path.expanduser("~"), ".sasmodels", "compiled_models") 
    164     if not os.path.exists(DLL_PATH): 
    165         os.makedirs(DLL_PATH) 
    166 else: 
    167     # Set up the default path for compiled modules. 
    168     DLL_PATH = tempfile.gettempdir() 
     160# Assume the default location of module DLLs is in .sasmodels/compiled_models. 
     161DLL_PATH = os.path.join(os.path.expanduser("~"), ".sasmodels", "compiled_models") 
    169162 
    170163ALLOW_SINGLE_PRECISION_DLLS = True 
    … …  
    233226 
    234227    Set *sasmodels.kerneldll.DLL_PATH* to the compiled dll output path. 
    235     The default is the system temporary directory. 
     228    The default is in ~/.sasmodels/compiled_models. 
    236229    """ 
    237230    if dtype == F16: 
    … …  
    250243        need_recompile = dll_time < newest_source 
    251244    if need_recompile: 
     245        # Make sure the DLL path exists 
     246        if not os.path.exists(DLL_PATH): 
     247            os.makedirs(DLL_PATH) 
    252248        basename = splitext(os.path.basename(dll))[0] + "_" 
    253249        system_fd, filename = tempfile.mkstemp(suffix=".c", prefix=basename) 

  • sasmodels/multiscat.py

    rbc248f8 r49d1f8b8  
    144144        frame = np.fft.ifft2(F_convolved) 
    145145        result = scale * _inverse_shift(frame.real, dtype=self.dtype) 
    146         #print("multiscat time", time.time()-t0) 
     146        #print("numpy multiscat time", time.time()-t0) 
    147147        return result 
    148148 
    … …  
    237237        #result = scale * _inverse_shift(frame.real, dtype=self.dtype) 
    238238        result = scale * _inverse_shift(frame.real, dtype=self.dtype) 
    239         #print("multiscat time", time.time()-t0) 
     239        #print("OpenCL multiscat time", time.time()-t0) 
    240240        return result 
    241241 

  • sasmodels/models/core_shell_parallelepiped.c

    re077231 rdbf1a60  
    5959 
    6060    // outer integral (with gauss points), integration limits = 0, 1 
     61    // substitute d_cos_alpha for sin_alpha d_alpha 
    6162    double outer_sum = 0; //initialize integral 
    6263    for( int i=0; i<GAUSS_N; i++) { 
    6364        const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 ); 
    6465        const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha); 
    65  
    66         // inner integral (with gauss points), integration limits = 0, pi/2 
    6766        const double siC = length_c * sas_sinx_x(length_c * cos_alpha * half_q); 
    6867        const double siCt = tC * sas_sinx_x(tC * cos_alpha * half_q); 
     68 
     69        // inner integral (with gauss points), integration limits = 0, 1 
     70        // substitute beta = PI/2 u (so 2/PI * d_(PI/2 * beta) = d_beta) 
    6971        double inner_sum = 0.0; 
    7072        for(int j=0; j<GAUSS_N; j++) { 
    71             const double beta = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
     73            const double u = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
    7274            double sin_beta, cos_beta; 
    73             SINCOS(M_PI_2*beta, sin_beta, cos_beta); 
     75            SINCOS(M_PI_2*u, sin_beta, cos_beta); 
    7476            const double siA = length_a * sas_sinx_x(length_a * mu * sin_beta); 
    7577            const double siB = length_b * sas_sinx_x(length_b * mu * cos_beta); 
    … …  
    9193            inner_sum += GAUSS_W[j] * f * f; 
    9294        } 
     95        // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5 
    9396        inner_sum *= 0.5; 
    9497        // now sum up the outer integral 
    9598        outer_sum += GAUSS_W[i] * inner_sum; 
    9699    } 
     100    // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5 
    97101    outer_sum *= 0.5; 
    98102 

  • sasmodels/models/core_shell_parallelepiped.py

    r97be877 rdbf1a60  
    1111.. math:: 
    1212 
    13     I(q) = \text{scale}\frac{\langle f^2 \rangle}{V} + \text{background} 
     13    I(q) = \frac{\text{scale}}{V} \langle P(q,\alpha,\beta) \rangle  
     14    + \text{background} 
    1415 
    1516where $\langle \ldots \rangle$ is an average over all possible orientations 
    16 of the rectangular solid. 
    17  
    18 The function calculated is the form factor of the rectangular solid below. 
     17of the rectangular solid, and the usual $\Delta \rho^2 \ V^2$ term cannot be 
     18pulled out of the form factor term due to the multiple slds in the model. 
     19 
    1920The core of the solid is defined by the dimensions $A$, $B$, $C$ such that 
    2021$A < B < C$. 
    2122 
    22 .. image:: img/core_shell_parallelepiped_geometry.jpg 
     23.. figure:: img/parallelepiped_geometry.jpg 
     24 
     25   Core of the core shell Parallelepiped with the corresponding definition 
     26   of sides. 
     27 
    2328 
    2429There are rectangular "slabs" of thickness $t_A$ that add to the $A$ dimension 
    2530(on the $BC$ faces). There are similar slabs on the $AC$ $(=t_B)$ and $AB$ 
    26 $(=t_C)$ faces. The projection in the $AB$ plane is then 
    27  
    28 .. image:: img/core_shell_parallelepiped_projection.jpg 
    29  
    30 The volume of the solid is 
     31$(=t_C)$ faces. The projection in the $AB$ plane is 
     32 
     33.. figure:: img/core_shell_parallelepiped_projection.jpg 
     34 
     35   AB cut through the core-shell parllelipiped showing the cross secion of 
     36   four of the six shell slabs. As can be seen This model leaves **"gaps"** 
     37   at the corners of the solid. 
     38 
     39 
     40The total volume of the solid is thus given as 
    3141 
    3242.. math:: 
    3343 
    3444    V = ABC + 2t_ABC + 2t_BAC + 2t_CAB 
    35  
    36 **meaning that there are "gaps" at the corners of the solid.** 
    3745 
    3846The intensity calculated follows the :ref:`parallelepiped` model, with the 
    3947core-shell intensity being calculated as the square of the sum of the 
    40 amplitudes of the core and the slabs on the edges. 
    41  
    42 the scattering amplitude is computed for a particular orientation of the 
    43 core-shell parallelepiped with respect to the scattering vector and then 
    44 averaged over all possible orientations, where $\alpha$ is the angle between 
    45 the $z$ axis and the $C$ axis of the parallelepiped, $\beta$ is 
    46 the angle between projection of the particle in the $xy$ detector plane 
    47 and the $y$ axis. 
    48  
    49 .. math:: 
    50  
    51     F(Q) 
     48amplitudes of the core and the slabs on the edges. The scattering amplitude is 
     49computed for a particular orientation of the core-shell parallelepiped with 
     50respect to the scattering vector and then averaged over all possible 
     51orientations, where $\alpha$ is the angle between the $z$ axis and the $C$ axis 
     52of the parallelepiped, and $\beta$ is the angle between the projection of the 
     53particle in the $xy$ detector plane and the $y$ axis. 
     54 
     55.. math:: 
     56 
     57    P(q)=\frac {\int_{0}^{\pi/2}\int_{0}^{\pi/2}F^2(q,\alpha,\beta) \ sin\alpha 
     58    \ d\alpha \ d\beta} {\int_{0}^{\pi/2} \ sin\alpha \ d\alpha \ d\beta} 
     59 
     60and 
     61 
     62.. math:: 
     63 
     64    F(q,\alpha,\beta) 
    5265    &= (\rho_\text{core}-\rho_\text{solvent}) 
    5366       S(Q_A, A) S(Q_B, B) S(Q_C, C) \\ 
    5467    &+ (\rho_\text{A}-\rho_\text{solvent}) 
    55         \left[S(Q_A, A+2t_A) - S(Q_A, Q)\right] S(Q_B, B) S(Q_C, C) \\ 
     68        \left[S(Q_A, A+2t_A) - S(Q_A, A)\right] S(Q_B, B) S(Q_C, C) \\ 
    5669    &+ (\rho_\text{B}-\rho_\text{solvent}) 
    5770        S(Q_A, A) \left[S(Q_B, B+2t_B) - S(Q_B, B)\right] S(Q_C, C) \\ 
    … …  
    6376.. math:: 
    6477 
    65     S(Q, L) = L \frac{\sin \tfrac{1}{2} Q L}{\tfrac{1}{2} Q L} 
     78    S(Q_X, L) = L \frac{\sin \tfrac{1}{2} Q_X L}{\tfrac{1}{2} Q_X L} 
    6679 
    6780and 
    … …  
    6982.. math:: 
    7083 
    71     Q_A &= \sin\alpha \sin\beta \\ 
    72     Q_B &= \sin\alpha \cos\beta \\ 
    73     Q_C &= \cos\alpha 
     84    Q_A &= q \sin\alpha \sin\beta \\ 
     85    Q_B &= q \sin\alpha \cos\beta \\ 
     86    Q_C &= q \cos\alpha 
    7487 
    7588 
    … …  
    7992is the scattering length of the solvent. 
    8093 
     94.. note::  
     95 
     96   the code actually implements two substitutions: $d(cos\alpha)$ is 
     97   substituted for -$sin\alpha \ d\alpha$ (note that in the 
     98   :ref:`parallelepiped` code this is explicitly implemented with 
     99   $\sigma = cos\alpha$), and $\beta$ is set to $\beta = u \pi/2$ so that 
     100   $du = \pi/2 \ d\beta$.  Thus both integrals go from 0 to 1 rather than 0 
     101   to $\pi/2$. 
     102 
    81103FITTING NOTES 
    82104~~~~~~~~~~~~~ 
    … …  
    95117and length $(C+2t_C)$ values, after appropriately sorting the three dimensions 
    96118to give an oblate or prolate particle, to give an effective radius, 
    97 for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
     119for $S(q)$ when $P(q) * S(q)$ is applied. 
    98120 
    99121For 2d data the orientation of the particle is required, described using 
    … …  
    103125$\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector. 
    104126 
    105 For 2d, constraints must be applied during fitting to ensure that the 
    106 inequality $A < B < C$ is not violated, and hence the correct definition 
    107 of angles is preserved. The calculation will not report an error, 
    108 but the results may be not correct. 
     127.. note:: For 2d, constraints must be applied during fitting to ensure that the 
     128   inequality $A < B < C$ is not violated, and hence the correct definition 
     129   of angles is preserved. The calculation will not report an error, 
     130   but the results may be not correct. 
    109131 
    110132.. figure:: img/parallelepiped_angle_definition.png 

  • sasmodels/models/parallelepiped.c

    r108e70e rdbf1a60  
    3838            inner_total += GAUSS_W[j] * square(si1 * si2); 
    3939        } 
     40        // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5 
    4041        inner_total *= 0.5; 
    4142 
    … …  
    4344        outer_total += GAUSS_W[i] * inner_total * si * si; 
    4445    } 
     46    // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5 
    4547    outer_total *= 0.5; 
    4648 

  • sasmodels/models/parallelepiped.py

    ref07e95 rdbf1a60  
    1010 
    1111 This model calculates the scattering from a rectangular parallelepiped 
    12  (\:numref:`parallelepiped-image`\). 
     12 (:numref:`parallelepiped-image`). 
    1313 If you need to apply polydispersity, see also :ref:`rectangular-prism`. 
    1414 
    … …  
    3939 
    4040    I(q) = \frac{\text{scale}}{V} (\Delta\rho \cdot V)^2 
    41            \left< P(q, \alpha) \right> + \text{background} 
     41           \left< P(q, \alpha, \beta) \right> + \text{background} 
    4242 
    4343where the volume $V = A B C$, the contrast is defined as 
    44 $\Delta\rho = \rho_\text{p} - \rho_\text{solvent}$, 
    45 $P(q, \alpha)$ is the form factor corresponding to a parallelepiped oriented 
    46 at an angle $\alpha$ (angle between the long axis C and $\vec q$), 
    47 and the averaging $\left<\ldots\right>$ is applied over all orientations. 
     44$\Delta\rho = \rho_\text{p} - \rho_\text{solvent}$, $P(q, \alpha, \beta)$ 
     45is the form factor corresponding to a parallelepiped oriented 
     46at an angle $\alpha$ (angle between the long axis C and $\vec q$), and $\beta$ 
     47( the angle between the projection of the particle in the $xy$ detector plane 
     48and the $y$ axis) and the averaging $\left<\ldots\right>$ is applied over all 
     49orientations. 
    4850 
    4951Assuming $a = A/B < 1$, $b = B /B = 1$, and $c = C/B > 1$, the 
    50 form factor is given by (Mittelbach and Porod, 1961) 
     52form factor is given by (Mittelbach and Porod, 1961 [#Mittelbach]_) 
    5153 
    5254.. math:: 
    … …  
    6668    \mu &= qB 
    6769 
    68 The scattering intensity per unit volume is returned in units of |cm^-1|. 
     70where substitution of $\sigma = cos\alpha$ and $\beta = \pi/2 \ u$ have been 
     71applied. 
    6972 
    7073NB: The 2nd virial coefficient of the parallelepiped is calculated based on 
    … …  
    120123.. math:: 
    121124 
    122     P(q_x, q_y) = \left[\frac{\sin(\tfrac{1}{2}qA\cos\alpha)}{(\tfrac{1}{2}qA\cos\alpha)}\right]^2 
    123                   \left[\frac{\sin(\tfrac{1}{2}qB\cos\beta)}{(\tfrac{1}{2}qB\cos\beta)}\right]^2 
    124                   \left[\frac{\sin(\tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1}{2}qC\cos\gamma)}\right]^2 
     125    P(q_x, q_y) = \left[\frac{\sin(\tfrac{1}{2}qA\cos\alpha)}{(\tfrac{1} 
     126                   {2}qA\cos\alpha)}\right]^2 
     127                  \left[\frac{\sin(\tfrac{1}{2}qB\cos\beta)}{(\tfrac{1} 
     128                   {2}qB\cos\beta)}\right]^2 
     129                  \left[\frac{\sin(\tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1} 
     130                   {2}qC\cos\gamma)}\right]^2 
    125131 
    126132with 
    … …  
    160166---------- 
    161167 
    162 P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211 
    163  
    164 R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 
     168.. [#Mittelbach] P Mittelbach and G Porod, *Acta Physica Austriaca*, 
     169   14 (1961) 185-211 
     170.. [#] R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 
    165171 
    166172Authorship and Verification 
Note: See TracChangeset for help on using the changeset viewer.