Changeset a66b004 in sasmodels


Ignore:
Timestamp:
Nov 30, 2017 12:30:08 AM (5 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:
34bbb9c
Parents:
26a6608 (diff), b669b49 (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:
19 added
1 deleted
137 edited

Legend:

Unmodified
Added
Removed
  • LICENSE.txt

    rf68e2a5 r5f8b72b  
    22All rights reserved. 
    33 
    4 Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: 
     4Redistribution and use in source and binary forms, with or without 
     5modification, are permitted provided that the following conditions are met: 
    56 
    6 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 
     71. Redistributions of source code must retain the above copyright notice, 
     8this list of conditions and the following disclaimer. 
    79 
    8 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 
     102. Redistributions in binary form must reproduce the above copyright notice, 
     11this list of conditions and the following disclaimer in the documentation 
     12and/or other materials provided with the distribution. 
    913 
    10 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. 
     143. Neither the name of the copyright holder nor the names of its contributors 
     15may be used to endorse or promote products derived from this software without 
     16specific prior written permission. 
    1117 
    12 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 
     18THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
     19AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
     20IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
     21ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE 
     22LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
     23CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
     24SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 
     25INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
     26CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
     27ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
     28POSSIBILITY OF SUCH DAMAGE. 
  • doc/_extensions/dollarmath.py

    r103ea45 rdd4d95d  
    1212import re 
    1313 
    14 _dollar = re.compile(r"(?:^|(?<=\s|[(]))[$]([^\n]*?)(?<![\\])[$](?:$|(?=\s|[.,;)\\]))") 
     14_dollar = re.compile(r"(?:^|(?<=\s|[-(]))[$]([^\n]*?)(?<![\\])[$](?:$|(?=\s|[-.,;:?\\)]))") 
    1515_notdollar = re.compile(r"\\[$]") 
    1616 
  • doc/developer/calculator.rst

    r870a2f4 r2c108a3  
    77 
    88This document describes the layer between the form factor kernels and the 
    9 model calculator which implements the polydispersity and magnetic SLD 
     9model calculator which implements the dispersity and magnetic SLD 
    1010calculations.  There are three separate implementations of this layer, 
    1111:mod:`kernelcl` for OpenCL, which operates on a single Q value at a time, 
     
    1414 
    1515Each implementation provides three different calls *Iq*, *Iqxy* and *Imagnetic* 
    16 for 1-D, 2-D and 2-D magnetic kernels respectively. 
    17  
    18 The C code is defined in *kernel_iq.c* and *kernel_iq.cl* for DLL and OpenCL 
    19 respectively.  The kernel call looks as follows:: 
     16for 1-D, 2-D and 2-D magnetic kernels respectively. The C code is defined 
     17in *kernel_iq.c*, with the minor differences between OpenCL and DLL handled 
     18by #ifdef statements. 
     19 
     20The kernel call looks as follows:: 
    2021 
    2122  kernel void KERNEL_NAME( 
    2223      int nq,                  // Number of q values in the q vector 
    23       int pd_start,            // Starting position in the polydispersity loop 
    24       int pd_stop,             // Ending position in the polydispersity loop 
    25       ProblemDetails *details, // Polydispersity info 
     24      int pd_start,            // Starting position in the dispersity loop 
     25      int pd_stop,             // Ending position in the dispersity loop 
     26      ProblemDetails *details, // dispersity info 
    2627      double *values,          // Value and weights vector 
    2728      double *q,               // q or (qx,qy) vector 
    2829      double *result,          // returned I(q), with result[nq] = pd_weight 
    29       double cutoff)           // polydispersity weight cutoff 
     30      double cutoff)           // dispersity weight cutoff 
    3031 
    3132The details for OpenCL and the python loop are slightly different, but these 
     
    3435*nq* indicates the number of q values that will be calculated. 
    3536 
    36 The *pd_start* and *pd_stop* parameters set the range of the polydispersity 
    37 loop to compute for the current kernel call.   Give a polydispersity 
     37The *pd_start* and *pd_stop* parameters set the range of the dispersity 
     38loop to compute for the current kernel call.   Give a dispersity 
    3839calculation with 30 weights for length and 30 weights for radius for example, 
    3940there are a total of 900 calls to the form factor required to compute the 
     
    4243the length index to 3 and the radius index to 10 for a position of 3*30+10=100, 
    4344and could then proceed to position 200.  This allows us to interrupt the 
    44 calculation in the middle of a long polydispersity loop without having to 
     45calculation in the middle of a long dispersity loop without having to 
    4546do special tricks with the C code.  More importantly, it stops the OpenCL 
    4647kernel in a reasonable time; because the GPU is used by the operating 
     
    4950 
    5051The *ProblemDetails* structure is a direct map of the 
    51 :class:`details.CallDetails` buffer.  This indicates which parameters are 
    52 polydisperse, and where in the values vector the values and weights can be 
    53 found.  For each polydisperse parameter there is a parameter id, the length 
    54 of the polydispersity loop for that parameter, the offset of the parameter 
     52:class:`details.CallDetails` buffer.  This indicates which parameters have 
     53dispersity, and where in the values vector the values and weights can be 
     54found.  For each parameter with dispersity there is a parameter id, the length 
     55of the dispersity loop for that parameter, the offset of the parameter 
    5556values in the pd value and pd weight vectors and the 'stride' from one index 
    5657to the next, which is used to translate between the position in the 
    57 polydispersity loop and the particular parameter indices.  The *num_eval* 
    58 field is the total size of the polydispersity loop.  *num_weights* is the 
     58dispersity loop and the particular parameter indices.  The *num_eval* 
     59field is the total size of the dispersity loop.  *num_weights* is the 
    5960number of elements in the pd value and pd weight vectors.  *num_active* is 
    60 the number of non-trivial pd loops (polydisperse parameters should be ordered 
    61 by decreasing pd vector length, with a length of 1 meaning no polydispersity). 
     61the number of non-trivial pd loops (parameters with dispersity should be ordered 
     62by decreasing pd vector length, with a length of 1 meaning no dispersity). 
    6263Oriented objects in 2-D need a cos(theta) spherical correction on the angular 
    6364variation in order to preserve the 'surface area' of the weight distribution. 
     
    7273*(Mx, My, Mz)*.  Sample magnetization is translated from *(M, theta, phi)* 
    7374to *(Mx, My, Mz)* before the kernel is called.   After the fixed values comes 
    74 the pd value vector, with the polydispersity values for each parameter 
     75the pd value vector, with the dispersity values for each parameter 
    7576stacked one after the other.  The order isn't important since the location 
    7677for each parameter is stored in the *pd_offset* field of the *ProblemDetails* 
     
    7879values, the pd weight vector is stored, with the same configuration as the 
    7980pd value vector.  Note that the pd vectors can contain values that are not 
    80 in the polydispersity loop; this is used by :class:`mixture.MixtureKernel` 
     81in the dispersity loop; this is used by :class:`mixture.MixtureKernel` 
    8182to make it easier to call the various component kernels. 
    8283 
     
    8788 
    8889The *results* vector contains one slot for each of the *nq* values, plus 
    89 one extra slot at the end for the current polydisperse normalization.  This 
    90 is required when the polydispersity loop is broken across several kernel 
    91 calls. 
     90one extra slot at the end for the weight normalization accumulated across 
     91all points in the dispersity mesh.  This is required when the dispersity 
     92loop is broken across several kernel calls. 
    9293 
    9394*cutoff* is a importance cutoff so that points which contribute negligibly 
     
    9798 
    9899- USE_OPENCL is defined if running in opencl 
    99 - MAX_PD is the maximum depth of the polydispersity loop [model specific] 
     100- MAX_PD is the maximum depth of the dispersity loop [model specific] 
    100101- NUM_PARS is the number of parameter values in the kernel.  This may be 
    101102  more than the number of parameters if some of the parameters are vector 
    102103  values. 
    103104- NUM_VALUES is the number of fixed values, which defines the offset in the 
    104   value list to the polydisperse value and weight vectors. 
     105  value list to the dispersity value and weight vectors. 
    105106- NUM_MAGNETIC is the number of magnetic SLDs 
    106107- MAGNETIC_PARS is a comma separated list of the magnetic SLDs, indicating 
    107108  their locations in the values vector. 
    108 - MAGNETIC_PAR0 to MAGNETIC_PAR2 are the first three magnetic parameter ids 
    109   so we can hard code the setting of magnetic values if there are only a 
    110   few of them. 
    111109- KERNEL_NAME is the name of the function being declared 
    112110- PARAMETER_TABLE is the declaration of the parameters to the kernel: 
     
    152150    Cylinder2D:: 
    153151 
    154         #define CALL_IQ(q, i, var) Iqxy(q[2*i], q[2*i+1], \ 
     152        #define CALL_IQ(q, i, var) Iqxy(qa, qc, \ 
    155153        var.length, \ 
    156154        var.radius, \ 
    157155        var.sld, \ 
    158         var.sld_solvent, \ 
    159         var.theta, \ 
    160         var.phi) 
     156        var.sld_solvent) 
    161157 
    162158- CALL_VOLUME(var) is similar, but for calling the form volume:: 
     
    182178        #define INVALID(var) constrained(var.p1, var.p2, var.p3) 
    183179 
    184 Our design supports a limited number of polydispersity loops, wherein 
    185 we need to cycle through the values of the polydispersity, calculate 
     180Our design supports a limited number of dispersity loops, wherein 
     181we need to cycle through the values of the dispersity, calculate 
    186182the I(q, p) for each combination of parameters, and perform a normalized 
    187183weighted sum across all the weights.  Parameters may be passed to the 
    188 underlying calculation engine as scalars or vectors, but the polydispersity 
     184underlying calculation engine as scalars or vectors, but the dispersity 
    189185calculator treats the parameter set as one long vector. 
    190186 
    191 Let's assume we have 8 parameters in the model, with two polydisperse.  Since 
    192 this is a 1-D model the orientation parameters won't be used:: 
     187Let's assume we have 8 parameters in the model, two of which allow dispersity. 
     188Since this is a 1-D model the orientation parameters won't be used:: 
    193189 
    194190    0: scale        {scl = constant} 
     
    196192    2: radius       {r = vector of 10pts} 
    197193    3: length       {l = vector of 30pts} 
    198     4: sld          {s = constant/(radius**2*length)} 
     194    4: sld          {s1 = constant/(radius**2*length)} 
    199195    5: sld_solvent  {s2 = constant} 
    200196    6: theta        {not used} 
     
    202198 
    203199This generates the following call to the kernel.  Note that parameters 4 and 
    204 5 are treated as polydisperse even though they are not --- this is because 
     2005 are treated as having dispersity even though they don't --- this is because 
    205201it is harmless to do so and simplifies the looping code:: 
    206202 
     
    210206    NUM_MAGNETIC = 2      // two parameters might be magnetic 
    211207    MAGNETIC_PARS = 4, 5  // they are sld and sld_solvent 
    212     MAGNETIC_PAR0 = 4     // sld index 
    213     MAGNETIC_PAR1 = 5     // solvent index 
    214208 
    215209    details { 
     
    218212        pd_offset = {10, 0, 31, 32}   // *length* starts at index 10 in weights 
    219213        pd_stride = {1, 30, 300, 300} // cumulative product of pd length 
    220         num_eval = 300   // 300 values in the polydispersity loop 
     214        num_eval = 300   // 300 values in the dispersity loop 
    221215        num_weights = 42 // 42 values in the pd vector 
    222216        num_active = 2   // only the first two pd are active 
     
    225219 
    226220    values = { scl, bkg,                                  // universal 
    227                r, l, s, s2, theta, phi,                   // kernel pars 
     221               r, l, s1, s2, theta, phi,                  // kernel pars 
    228222               in spin, out spin, spin angle,             // applied magnetism 
    229                mx s, my s, mz s, mx s2, my s2, mz s2,     // magnetic slds 
     223               mx s1, my s1, mz s1, mx s2, my s2, mz s2,  // magnetic slds 
    230224               r0, .., r9, l0, .., l29, s, s2,            // pd values 
    231225               r0, .., r9, l0, .., l29, s, s2}            // pd weights 
     
    235229    result = {r1, ..., r130, pd_norm, x } 
    236230 
    237 The polydisperse parameters are stored in as an array of parameter 
    238 indices, one for each polydisperse parameter, stored in pd_par[n]. 
    239 Non-polydisperse parameters do not appear in this array. Each polydisperse 
     231The dispersity parameters are stored in as an array of parameter 
     232indices, one for each parameter, stored in pd_par[n]. Parameters which do 
     233not support dispersity do not appear in this array. Each dispersity 
    240234parameter has a weight vector whose length is stored in pd_length[n]. 
    241235The weights are stored in a contiguous vector of weights for all 
     
    243237in pd_offset[n].  The values corresponding to the weights are stored 
    244238together in a separate weights[] vector, with offset stored in 
    245 par_offset[pd_par[n]]. Polydisperse parameters should be stored in 
     239par_offset[pd_par[n]]. Dispersity parameters should be stored in 
    246240decreasing order of length for highest efficiency. 
    247241 
    248 We limit the number of polydisperse dimensions to MAX_PD (currently 4), 
    249 though some models may have fewer if they have fewer polydisperse 
     242We limit the number of dispersity dimensions to MAX_PD (currently 4), 
     243though some models may have fewer if they have fewer dispersity 
    250244parameters.  The main reason for the limit is to reduce code size. 
    251 Each additional polydisperse parameter requires a separate polydispersity 
    252 loop.  If more than 4 levels of polydispersity are needed, then kernel_iq.c 
    253 and kernel_iq.cl will need to be extended. 
     245Each additional dispersity parameter requires a separate dispersity 
     246loop.  If more than 4 levels of dispersity are needed, then we need to 
     247switch to a monte carlo importance sampling algorithm with better 
     248performance for high-dimensional integrals. 
    254249 
    255250Constraints between parameters are not supported.  Instead users will 
     
    262257theta since the polar coordinates normalization is tied to this parameter. 
    263258 
    264 If there is no polydispersity we pretend that it is polydisperisty with one 
    265 parameter, pd_start=0 and pd_stop=1.  We may or may not short circuit the 
    266 calculation in this case, depending on how much time it saves. 
     259If there is no dispersity we pretend that we have a disperisty mesh over 
     260a single parameter with a single point in the distribution, giving 
     261pd_start=0 and pd_stop=1. 
    267262 
    268263The problem details structure could be allocated and sent in as an integer 
    269264array using the read-only flag.  This would allow us to copy it once per fit 
    270265along with the weights vector, since features such as the number of 
    271 polydisperity elements per pd parameter or the coordinated won't change 
    272 between function evaluations.  A new parameter vector must be sent for 
    273 each I(q) evaluation.  This is not currently implemented, and would require 
    274 some resturcturing of the :class:`sasview_model.SasviewModel` interface. 
    275  
    276 The results array will be initialized to zero for polydispersity loop 
     266disperity points per pd parameter won't change between function evaluations. 
     267A new parameter vector must be sent for each I(q) evaluation.  This is 
     268not currently implemented, and would require some resturcturing of 
     269the :class:`sasview_model.SasviewModel` interface. 
     270 
     271The results array will be initialized to zero for dispersity loop 
    277272entry zero, and preserved between calls to [start, stop] so that the 
    278273results accumulate by the time the loop has completed.  Background and 
     
    295290 
    296291This will require accumulated error for each I(q) value to be preserved 
    297 between kernel calls to implement this fully.  The kernel_iq.c code, which 
    298 loops over q for each parameter set in the polydispersity loop, will need 
    299 also need the accumalation vector. 
     292between kernel calls to implement this fully.  The *kernel_iq.c* code, which 
     293loops over q for each parameter set in the dispersity loop, will also need 
     294the accumulation vector. 
  • doc/developer/overview.rst

    r870a2f4 r3d40839  
    166166:ref:`Calculator_Interface` 
    167167 
     168.. _orientation_developer: 
     169 
     170Orientation and Numerical Integration 
     171------------------------------------- 
     172 
     173For 2d data from oriented anisotropic particles, the mean particle 
     174orientation is defined by angles $\theta$, $\phi$ and $\Psi$, which are not 
     175in general the same as similarly named angles in many form factors. The 
     176wikipedia page on Euler angles (https://en.wikipedia.org/wiki/Euler_angles) 
     177lists the different conventions available. To quote: "Different authors may 
     178use different sets of rotation axes to define Euler angles, or different 
     179names for the same angles. Therefore, any discussion employing Euler angles 
     180should always be preceded by their definition." 
     181 
     182We are using the $z$-$y$-$z$ convention with extrinsic rotations 
     183$\Psi$-$\theta$-$\phi$ for the particle orientation and $x$-$y$-$z$ 
     184convention with extrinsic rotations $\Psi$-$\theta$-$\phi$ for jitter, with 
     185jitter applied before particle orientation. 
     186 
     187For numerical integration within form factors etc. sasmodels is mostly using 
     188Gaussian quadrature with 20, 76 or 150 points depending on the model. It also 
     189makes use of symmetries such as calculating only over one quadrant rather 
     190than the whole sphere. There is often a U-substitution replacing $\theta$ 
     191with $cos(\theta)$ which changes the limits of integration from 0 to $\pi/2$ 
     192to 0 to 1 and also conveniently absorbs the $sin(\theta)$ scale factor in the 
     193integration. This can cause confusion if checking equations to include in a 
     194paper or thesis! Most models use the same core kernel code expressed in terms 
     195of the rotated view ($q_a$, $q_b$, $q_c$) for both the 1D and the 2D models, 
     196but there are also historical quirks such as the parallelepiped model, which 
     197has a useless transformation representing $j_0(a q_a)$ as $j_0(b q_a a/b)$. 
     198 
     199Useful testing routines include: 
     200 
     201:mod:`asymint` a direct implementation of the surface integral for certain 
     202models to get a more trusted value for the 1D integral using a 
     203reimplementation of the 2D kernel in python and mpmath (which computes math 
     204functions to arbitrary precision). It uses $\theta$ ranging from 0 to $\pi$ 
     205and $\phi$ ranging from 0 to $2\pi$. It perhaps would benefit from including 
     206the U-substitution for $\theta$. 
     207 
     208:mod:`check1d` uses sasmodels 1D integration and compares that with a 
     209rectangle distribution in $\theta$ and $\phi$, with $\theta$ limits set to 
     210$\pm 90/\sqrt(3)$ and $\phi$ limits set to $\pm 180/\sqrt(3)$ [The rectangle 
     211weight function uses the fact that the distribution width column is labelled 
     212sigma to decide that the 1-$\sigma$ width of a rectangular distribution needs to 
     213be multiplied by $\sqrt(3)$ to get the corresponding gaussian equivalent, or 
     214similar reasoning.] This should rotate the sample through the entire 
     215$\theta$-$\phi$ surface according to the pattern that you see in jitter.py when 
     216you modify it to use 'rectangle' rather than 'gaussian' for its distribution 
     217without changing the viewing angle. In order to match the 1-D pattern for 
     218an arbitrary viewing angle on triaxial shapes, we need to integrate 
     219over $\theta$, $\phi$ and $\Psi$. 
     220 
     221When computing the dispersity integral, weights are scaled by 
     222$|\cos(\delta \theta)|$ to account for the points in $\phi$ getting closer 
     223together as $\delta \theta$ increases. 
     224[This will probably change so that instead of adjusting the weights, we will 
     225adjust $\delta\theta$-$\delta\phi$ mesh so that the point density in 
     226$\delta\phi$ is lower at larger $\delta\theta$. The flag USE_SCALED_PHI in 
     227*kernel_iq.c* selects an alternative algorithm.] 
     228 
     229The integrated dispersion is computed at a set of $(qx, qy)$ points $(q 
     230\cos(\alpha), q \sin(\alpha))$ at some angle $\alpha$ (currently angle=0) for 
     231each $q$ used in the 1-D integration. The individual $q$ points should be 
     232equivalent to asymint rect-n when the viewing angle is set to 
     233$(\theta,\phi,\Psi) = (90, 0, 0)$. Such tests can help to validate that 2d 
     234models are consistent with 1d models. 
     235 
     236:mod:`sascomp -sphere=n` uses the same rectangular distribution as check1d to 
     237compute the pattern of the $q_x$-$q_y$ grid. 
     238 
     239The :mod:`sascomp` utility can be used for 2d as well as 1d calculations to 
     240compare results for two sets of parameters or processor types, for example 
     241these two oriented cylinders here should be equivalent. 
     242 
     243:mod:`\./sascomp -2d cylinder theta=0 phi=0,90 theta_pd_type=rectangle phi_pd_type=rectangle phi_pd=10,1 theta_pd=1,10 length=500 radius=10` 
     244 
    168245 
    169246Testing 
  • doc/genapi.py

    r2e66ef5 r706f466  
    5959    #('alignment', 'GPU data alignment [unused]'), 
    6060    ('bumps_model', 'Bumps interface'), 
     61    ('compare', 'Compare models on different compute engines'), 
    6162    ('compare_many', 'Batch compare models on different compute engines'), 
    62     ('compare', 'Compare models on different compute engines'), 
     63    ('conversion_table', 'Model conversion table'), 
    6364    ('convert', 'Sasview to sasmodel converter'), 
    6465    ('core', 'Model access'), 
     
    8283    ('sasview_model', 'Sasview interface'), 
    8384    ('sesans', 'SESANS calculation routines'), 
     85    ('special', 'Special functions library'), 
    8486    ('weights', 'Distribution functions'), 
    8587] 
  • doc/guide/index.rst

    rc0d7ab3 rda5536f  
    1313   resolution.rst 
    1414   magnetism/magnetism.rst 
     15   orientation/orientation.rst 
    1516   sesans/sans_to_sesans.rst 
    1617   sesans/sesans_fitting.rst 
  • doc/guide/magnetism/magnetism.rst

    r1f058ea r4f5afc9  
    55 
    66Models which define a scattering length density parameter can be evaluated 
    7  as magnetic models. In general, the scattering length density (SLD = 
    8  $\beta$) in each region where the SLD is uniform, is a combination of the 
    9  nuclear and magnetic SLDs and, for polarised neutrons, also depends on the 
    10  spin states of the neutrons. 
     7as magnetic models. In general, the scattering length density (SLD = 
     8$\beta$) in each region where the SLD is uniform, is a combination of the 
     9nuclear and magnetic SLDs and, for polarised neutrons, also depends on the 
     10spin states of the neutrons. 
    1111 
    1212For magnetic scattering, only the magnetization component $\mathbf{M_\perp}$ 
    1313perpendicular to the scattering vector $\mathbf{Q}$ contributes to the magnetic 
    1414scattering length. 
    15  
    1615 
    1716.. figure:: 
     
    2827is the Pauli spin. 
    2928 
    30 Assuming that incident neutrons are polarized parallel (+) and anti-parallel (-) 
    31 to the $x'$ axis, the possible spin states after the sample are then 
     29Assuming that incident neutrons are polarized parallel $(+)$ and anti-parallel 
     30$(-)$ to the $x'$ axis, the possible spin states after the sample are then: 
    3231 
    33 No spin-flips (+ +) and (- -) 
     32* Non spin-flip $(+ +)$ and $(- -)$ 
    3433 
    35 Spin-flips    (+ -) and (- +) 
     34* Spin-flip $(+ -)$ and $(- +)$ 
     35 
     36Each measurement is an incoherent mixture of these spin states based on the 
     37fraction of $+$ neutrons before ($u_i$) and after ($u_f$) the sample, 
     38with weighting: 
     39 
     40.. math:: 
     41    -- &= ((1-u_i)(1-u_f))^{1/4} \\ 
     42    -+ &= ((1-u_i)(u_f))^{1/4} \\ 
     43    +- &= ((u_i)(1-u_f))^{1/4} \\ 
     44    ++ &= ((u_i)(u_f))^{1/4} 
     45 
     46Ideally the experiment would measure the pure spin states independently and 
     47perform a simultaneous analysis of the four states, tying all the model 
     48parameters together except $u_i$ and $u_f$. 
    3649 
    3750.. figure:: 
     
    4154$\phi$ and $\theta_{up}$, respectively, then, depending on the spin state of the 
    4255neutrons, the scattering length densities, including the nuclear scattering 
    43 length density ($\beta{_N}$) are 
     56length density $(\beta{_N})$ are 
    4457 
    4558.. math:: 
    4659    \beta_{\pm\pm} =  \beta_N \mp D_M M_{\perp x'} 
    47     \text{ when there are no spin-flips} 
     60    \text{ for non spin-flip states} 
    4861 
    4962and 
     
    5164.. math:: 
    5265    \beta_{\pm\mp} =  -D_M (M_{\perp y'} \pm iM_{\perp z'}) 
    53     \text{ when there are} 
     66    \text{ for spin-flip states} 
    5467 
    5568where 
    5669 
    5770.. math:: 
    58     M_{\perp x'} = M_{0q_x}\cos(\theta_{up})+M_{0q_y}\sin(\theta_{up}) \\ 
    59     M_{\perp y'} = M_{0q_y}\cos(\theta_{up})-M_{0q_x}\sin(\theta_{up}) \\ 
    60     M_{\perp z'} = M_{0z} \\ 
    61     M_{0q_x} = (M_{0x}\cos\phi - M_{0y}\sin\phi)\cos\phi \\ 
    62     M_{0q_y} = (M_{0y}\sin\phi - M_{0x}\cos\phi)\sin\phi 
     71    M_{\perp x'} &= M_{0q_x}\cos(\theta_{up})+M_{0q_y}\sin(\theta_{up}) \\ 
     72    M_{\perp y'} &= M_{0q_y}\cos(\theta_{up})-M_{0q_x}\sin(\theta_{up}) \\ 
     73    M_{\perp z'} &= M_{0z} \\ 
     74    M_{0q_x} &= (M_{0x}\cos\phi - M_{0y}\sin\phi)\cos\phi \\ 
     75    M_{0q_y} &= (M_{0y}\sin\phi - M_{0x}\cos\phi)\sin\phi 
    6376 
    6477Here, $M_{0x}$, $M_{0x}$, $M_{0z}$ are the x, y and z components 
     
    6679 
    6780.. math:: 
    68     M_{0x} = M_0\cos\theta_M\cos\phi_M \\ 
    69     M_{0y} = M_0\sin\theta_M \\ 
    70     M_{0z} = -M_0\cos\theta_M\sin\phi_M 
     81    M_{0x} &= M_0\cos\theta_M\cos\phi_M \\ 
     82    M_{0y} &= M_0\sin\theta_M \\ 
     83    M_{0z} &= -M_0\cos\theta_M\sin\phi_M 
    7184 
    7285and the magnetization angles $\theta_M$ and $\phi_M$ are defined in 
     
    7689 
    7790===========   ================================================================ 
    78  M0_sld        = $D_M M_0$ 
    79  Up_theta      = $\theta_\mathrm{up}$ 
    80  M_theta       = $\theta_M$ 
    81  M_phi         = $\phi_M$ 
    82  Up_frac_i     = (spin up)/(spin up + spin down) neutrons *before* the sample 
    83  Up_frac_f     = (spin up)/(spin up + spin down) neutrons *after* the sample 
     91 M0:sld      $D_M M_0$ 
     92 mtheta:sld   $\theta_M$ 
     93 mphi:sld     $\phi_M$ 
     94 up:angle     $\theta_\mathrm{up}$ 
     95 up:frac_i    $u_i$ = (spin up)/(spin up + spin down) *before* the sample 
     96 up:frac_f    $u_f$ = (spin up)/(spin up + spin down) *after* the sample 
    8497===========   ================================================================ 
    8598 
    8699.. note:: 
    87     The values of the 'Up_frac_i' and 'Up_frac_f' must be in the range 0 to 1. 
     100    The values of the 'up:frac_i' and 'up:frac_f' must be in the range 0 to 1. 
    88101 
    89102*Document History* 
    90103 
    91104| 2015-05-02 Steve King 
    92 | 2017-05-08 Paul Kienzle 
     105| 2017-11-15 Paul Kienzle 
  • doc/guide/pd/polydispersity.rst

    r75e4319 ra66b004  
    66.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
    77 
     8.. _polydispersityhelp: 
     9 
    810Polydispersity Distributions 
    911---------------------------- 
    1012 
    11 With some models in sasmodels we can calculate the average form factor for a 
     13With some models in sasmodels we can calculate the average intensity for a 
    1214population of particles that exhibit size and/or orientational 
    13 polydispersity. The resultant form factor is normalized by the average 
     15polydispersity. The resultant intensity is normalized by the average 
    1416particle volume such that 
    1517 
  • doc/guide/resolution.rst

    r1f058ea r0db85af  
    209209$x'_0 = x_0 \cos(\theta) + y_0 \sin(\theta)$ and 
    210210$y'_0 = -x_0 \sin(\theta) + y_0 \cos(\theta)$. 
    211 Note that the rotation angle is zero for a $x$\ -\ $y$ symmetric 
     211Note that the rotation angle is zero for a $x$-$y$ symmetric 
    212212elliptical Gaussian distribution. The $A$ is a normalization factor. 
    213213 
     
    233233 
    234234Since the weighting factor on each of the bins is known, it is convenient to 
    235 transform $x'$\ -\ $y'$ back to $x$\ -\ $y$ coordinates (by rotating it 
     235transform $x'$-$y'$ back to $x$-$y$ coordinates (by rotating it 
    236236by $-\theta$ around the $z$\ -axis). 
    237237 
     
    254254    y'_0 &= 0 
    255255 
    256 while for a $x$\ -\ $y$ symmetric smear 
     256while for a $x$-$y$ symmetric smear 
    257257 
    258258.. math:: 
  • example/sesans_parameters_sphere.py

    r9217ef8 rfa79f5c  
    4646# DO NOT MODIFY THIS LINE 
    4747problem = sesansfit.sesans_fit(sesans_file, model, initial_vals, custom_params, param_range) 
    48  
  • example/sesans_sphere_2micron.py

    r3330bb4 rfa79f5c  
    11""" 
    2 This is a data file  used to load in sesans data and fit it using the bumps engine 
     2This is a data file used to load in sesans data and fit it using the bumps engine 
    33""" 
    44from bumps.names import * 
     
    3838# Constraints 
    3939# model.param_name = f(other params) 
    40 # EXAMPLE: model.scale = model.radius*model.radius*(1 - phi) - where radius and scale are model functions and phi is 
    41 # a custom parameter 
     40# EXAMPLE: model.scale = model.radius*model.radius*(1 - phi) - where radius 
     41# and scale are model functions and phi is a custom parameter 
    4242model.scale = phi*(1-phi) 
    4343 
     
    4545# DO NOT MODIFY THIS LINE 
    4646problem = sesansfit.sesans_fit(sesans_file, model, initial_vals, custom_params, param_range) 
    47  
  • example/sesansfit.py

    r9217ef8 rfa79f5c  
     1import logging 
     2 
    13from bumps.names import * 
    24from sasmodels import core, bumps_model, sesans 
     
    810    return model 
    911 
    10 def sesans_fit(file, model, initial_vals={}, custom_params={}, param_range=[], acceptance_angle=None): 
     12def sesans_fit(file, model, initial_vals={}, custom_params={}, param_range=[], 
     13               acceptance_angle=None): 
    1114    """ 
    1215 
     
    1922    @return: FitProblem for Bumps usage 
    2023    """ 
     24    logging.basicConfig() 
     25 
    2126    initial_vals['background'] = 0.0 
    2227    try: 
    2328        loader = Loader() 
    24         data = loader.load(file) 
    25         if data is None: raise IOError("Could not load file %r"%(file)) 
     29        data = loader.load(file)[0] 
     30        if data is None: 
     31            raise IOError("Could not load file %r"%(file)) 
    2632 
    27     except: 
     33    except Exception: 
     34        raise 
    2835        # If no loadable data file, generate random data 
    2936        SElength = np.linspace(0, 2400, 61) # [A] 
     
    5057    data.Rmax = 30*radius # [A] 
    5158 
    52     if isinstance(model, basestring): 
     59    if isinstance(model, str): 
    5360        model = get_bumps_model(model) 
    5461 
  • example/spheres2micron.ses

    r6abf703 rfa79f5c  
    1 DataFileTitle   Polystyrene of Markus Strobl,  Full Sine, ++ only  
    2 Sample  Polystyrene 2 um in 53% H2O, 47% D2O  
    3 Settings        D1=D2=20x8 mm,Ds = 16x10 mm (WxH), GF1 =scanning, GF2 = 2.5 A. 2 um polystyrene in 53% H2O, 47% D2O; 8.55% contrast  
    4 Operator        CPD  
    5 Date    ma 7 jul 2014 18:54:43  
    6 ScanType        sine one element scan  
    7 Thickness [mm]  2  
    8 Q_zmax [\AA^-1]         0.05  
    9 Q_ymax [\AA^-1]         0.05  
    10    
    11 spin echo length [nm]    error SEL       wavelength [nm]         error wavelength        polarisation    error pol  
    12 49.778  2.4889  0.211   0.01055 0.99782 0.0044367 
    13 63.041  3.152   0.211   0.01055 1.0026  0.0047862 
    14 76.487  3.8244  0.211   0.01055 0.99601 0.0060598 
    15 89.847  4.4924  0.211   0.01055 0.99175 0.0058257 
    16 103.41  5.1705  0.211   0.01055 0.99543 0.0060966 
    17 116.95  5.8475  0.211   0.01055 0.99512 0.0048106 
    18 130.61  6.5303  0.211   0.01055 0.9975  0.0062594 
    19 144.37  7.2184  0.211   0.01055 0.99473 0.005293 
    20 158.2   7.9102  0.211   0.01055 0.9927  0.0053389 
    21 172.12  8.6062  0.211   0.01055 0.99453 0.0064548 
    22 186.17  9.3087  0.211   0.01055 0.98641 0.0073259 
    23 200.28  10.014  0.211   0.01055 0.98879 0.0078682 
    24 214.46  10.723  0.211   0.01055 0.9895  0.0068058 
    25 228.73  11.436  0.211   0.01055 0.99204 0.0082509 
    26 243.12  12.156  0.211   0.01055 0.99406 0.0051381 
    27 257.55  12.878  0.211   0.01055 0.97965 0.0055447 
    28 272.22  13.611  0.211   0.01055 0.97982 0.0065843 
    29 286.9   14.345  0.211   0.01055 0.97793 0.0071814 
    30 301.73  15.087  0.211   0.01055 0.97835 0.0066197 
    31 316.75  15.838  0.211   0.01055 0.98125 0.0069484 
    32 331.82  16.591  0.211   0.01055 0.97828 0.0068737 
    33 347.16  17.358  0.211   0.01055 0.97682 0.006744 
    34 362.45  18.122  0.211   0.01055 0.98155 0.0059865 
    35 378.09  18.904  0.211   0.01055 0.96446 0.0066788 
    36 393.74  19.687  0.211   0.01055 0.97276 0.0073781 
    37 409.61  20.481  0.211   0.01055 0.97101 0.0071053 
    38 425.55  21.278  0.211   0.01055 0.97501 0.0068244 
    39 441.91  22.096  0.211   0.01055 0.96958 0.0072272 
    40 458.12  22.906  0.211   0.01055 0.95991 0.0064035 
    41 474.77  23.739  0.211   0.01055 0.96219 0.0059227 
    42 491.52  24.576  0.211   0.01055 0.96016 0.0058606 
    43 508.51  25.426  0.211   0.01055 0.9509  0.0049836 
    44 525.68  26.284  0.211   0.01055 0.95558 0.0055218 
    45 543.08  27.154  0.211   0.01055 0.95079 0.0058258 
    46 560.74  28.037  0.211   0.01055 0.9514  0.0043921 
    47 578.69  28.935  0.211   0.01055 0.94562 0.0046962 
    48 596.75  29.838  0.211   0.01055 0.94588 0.0047341 
    49 615.17  30.758  0.211   0.01055 0.93958 0.0045821 
    50 633.77  31.688  0.211   0.01055 0.93731 0.0043108 
    51 652.82  32.641  0.211   0.01055 0.93664 0.0052044 
    52 672.04  33.602  0.211   0.01055 0.93209 0.0035563 
    53 691.57  34.578  0.211   0.01055 0.93385 0.0046002 
    54 711.48  35.574  0.211   0.01055 0.93028 0.0051136 
    55 731.67  36.583  0.211   0.01055 0.93012 0.0035279 
    56 752.17  37.608  0.211   0.01055 0.92594 0.0039458 
    57 773.1   38.655  0.211   0.01055 0.92352 0.0027613 
    58 794.42  39.721  0.211   0.01055 0.91823 0.0034587 
    59 816.14  40.807  0.211   0.01055 0.90862 0.0034262 
    60 838.19  41.91   0.211   0.01055 0.91909 0.003498 
    61 860.69  43.034  0.211   0.01055 0.91013 0.0030537 
    62 883.74  44.187  0.211   0.01055 0.90449 0.0035391 
    63 907.04  45.352  0.211   0.01055 0.90057 0.0029672 
    64 930.95  46.547  0.211   0.01055 0.90316 0.0031245 
    65 955.38  47.769  0.211   0.01055 0.89492 0.0030981 
    66 980.25  49.013  0.211   0.01055 0.89636 0.0036625 
    67 1005.7  50.284  0.211   0.01055 0.88808 0.0022652 
    68 1031.7  51.586  0.211   0.01055 0.88665 0.0044622 
    69 1058.3  52.914  0.211   0.01055 0.88205 0.0030275 
    70 1085.5  54.276  0.211   0.01055 0.88257 0.0028637 
    71 1113.4  55.669  0.211   0.01055 0.87629 0.0030312 
    72 1141.9  57.093  0.211   0.01055 0.87732 0.0026539 
    73 1171    58.55   0.211   0.01055 0.87361 0.0034591 
    74 1201    60.052  0.211   0.01055 0.86595 0.0028129 
    75 1231.6  61.582  0.211   0.01055 0.86804 0.0043311 
    76 1263    63.151  0.211   0.01055 0.86676 0.0028806 
    77 1295.2  64.762  0.211   0.01055 0.85655 0.0028657 
    78 1328.4  66.419  0.211   0.01055 0.85794 0.0034773 
    79 1362.3  68.113  0.211   0.01055 0.85521 0.0032798 
    80 1397.1  69.854  0.211   0.01055 0.84974 0.0044537 
    81 1432.9  71.644  0.211   0.01055 0.84149 0.0032738 
    82 1469.7  73.487  0.211   0.01055 0.84801 0.0037874 
    83 1507.5  75.376  0.211   0.01055 0.84334 0.003206 
    84 1546.4  77.318  0.211   0.01055 0.83749 0.0040707 
    85 1586.4  79.322  0.211   0.01055 0.83921 0.0049938 
    86 1627.6  81.381  0.211   0.01055 0.83434 0.0057337 
    87 1669.9  83.497  0.211   0.01055 0.82878 0.0062417 
    88 1713.6  85.678  0.211   0.01055 0.82309 0.0080179 
    89 1758.5  87.927  0.211   0.01055 0.82433 0.0074191 
    90 1804.8  90.239  0.211   0.01055 0.82023 0.0076431 
    91 1852.5  92.627  0.211   0.01055 0.82756 0.0063064 
    92 1901.7  95.087  0.211   0.01055 0.82584 0.0052944 
    93 1952.4  97.622  0.211   0.01055 0.82799 0.0050662 
    94 2004.8  100.24  0.211   0.01055 0.82345 0.0043127 
    95 2058.8  102.94  0.211   0.01055 0.82296 0.0038 
    96 2114.5  105.72  0.211   0.01055 0.81987 0.0034072 
    97 2172    108.6   0.211   0.01055 0.82021 0.0036752 
    98 2231.4  111.57  0.211   0.01055 0.82765 0.0028851 
    99 2292.8  114.64  0.211   0.01055 0.82664 0.0038942 
    100 2356.2  117.81  0.211   0.01055 0.82702 0.0047371 
    101 2421.8  121.09  0.211   0.01055 0.81593 0.0043772 
    102 2489.6  124.48  0.211   0.01055 0.8251  0.0028026 
    103 2559.5  127.98  0.211   0.01055 0.8292  0.0024574 
    104 2631.9  131.59  0.211   0.01055 0.82626 0.0036198 
    105 2706.8  135.34  0.211   0.01055 0.8208  0.0032314 
    106 2784.3  139.22  0.211   0.01055 0.81959 0.0042731 
    107 2864.5  143.23  0.211   0.01055 0.82653 0.002699 
    108 2947.5  147.38  0.211   0.01055 0.82401 0.0036726 
    109 3033.5  151.67  0.211   0.01055 0.82361 0.0048224 
    110 3122.4  156.12  0.211   0.01055 0.82358 0.0041221 
    111 3214.5  160.73  0.211   0.01055 0.82187 0.0028807 
    112 3310.1  165.5   0.211   0.01055 0.82644 0.003516 
    113 3409    170.45  0.211   0.01055 0.82355 0.0021166 
    114 3511.4  175.57  0.211   0.01055 0.82513 0.0033911 
    115 3617.6  180.88  0.211   0.01055 0.82802 0.0015342 
    116 3727.6  186.38  0.211   0.01055 0.82663 0.0029222 
    117 3841.7  192.08  0.211   0.01055 0.82026 0.0020755 
    118 3960    198     0.211   0.01055 0.83079 0.0026394 
    119 4082.7  204.13  0.211   0.01055 0.82665 0.0027466 
    120 4209.9  210.5   0.211   0.01055 0.82774 0.0025199 
    121 4341.9  217.09  0.211   0.01055 0.83489 0.0030619 
    122 4478.7  223.94  0.211   0.01055 0.81987 0.0020988 
    123 4620.8  231.04  0.211   0.01055 0.8253  0.0023899 
    124 4768.2  238.41  0.211   0.01055 0.82653 0.0022851 
    125 4921.1  246.06  0.211   0.01055 0.82442 0.003383 
    126 5079.8  253.99  0.211   0.01055 0.82827 0.0015979 
    127 5244.7  262.23  0.211   0.01055 0.82494 0.0031129 
    128 5415.7  270.79  0.211   0.01055 0.82183 0.0030149 
    129 5593.3  279.67  0.211   0.01055 0.83217 0.0046976 
    130 5777.8  288.89  0.211   0.01055 0.82227 0.005574 
    131 5969.3  298.46  0.211   0.01055 0.82338 0.0025569 
     1FileFormatVersion       1.0 
     2DataFileTitle           Polystyrene of Markus Strobl,  Full Sine, ++ only 
     3Sample                  Polystyrene 2 um in 53% H2O, 47% D2O 
     4Settings                D1=D2=20x8 mm,Ds = 16x10 mm (WxH), GF1 =scanning, GF2 = 2.5 A. 2 um polystyrene in 53% H2O, 47% D2O; 8.55% contrast 
     5Operator                CPD 
     6Date                    do 10 jul 2014 16:37:30 
     7ScanType                sine one element scan 
     8Thickness               2.00E-01 
     9Thickness_unit          cm 
     10Theta_zmax              0.0168 
     11Theta_zmax_unit         radians 
     12Theta_ymax              0.0168 
     13Theta_ymax_unit         radians 
     14Orientation             Z 
     15SpinEchoLength_unit     A 
     16Depolarisation_unit     A-2 cm-1 
     17Wavelength_unit         A 
     18 
     19BEGIN_DATA 
     20SpinEchoLength Depolarisation Depolarisation_error SpinEchoLength_error Wavelength Wavelength_error Polarisation  Polarisation_error 
     21391.56 0.0041929 0.0036894 19.578 2.11 0.1055 1.0037 0.0032974 
     221564 -0.0046571 0.0038185 78.2 2.11 0.1055 0.99586 0.003386 
     232735.6 -0.017007 0.0038132 136.78 2.11 0.1055 0.98497 0.0033444 
     243907.9 -0.033462 0.0035068 195.39 2.11 0.1055 0.97064 0.0030309 
     255080.2 -0.047483 0.0038208 254.01 2.11 0.1055 0.9586 0.0032613 
     266251.8 -0.070375 0.00376 312.59 2.11 0.1055 0.93926 0.0031446 
     277423.2 -0.092217 0.0037927 371.16 2.11 0.1055 0.92117 0.0031108 
     288595.5 -0.10238 0.004006 429.77 2.11 0.1055 0.91287 0.0032562 
     299767.7 -0.12672 0.0038534 488.39 2.11 0.1055 0.8933 0.0030651 
     3010940 -0.1374 0.004243 546.98 2.11 0.1055 0.88484 0.003343 
     3112112 -0.16072 0.0045837 605.58 2.11 0.1055 0.86666 0.0035372 
     3213284 -0.16623 0.0045613 664.2 2.11 0.1055 0.86242 0.0035027 
     3314456 -0.18468 0.0044918 722.79 2.11 0.1055 0.84837 0.0033931 
     3415628 -0.19143 0.0048967 781.38 2.11 0.1055 0.84328 0.0036768 
     3516800 -0.20029 0.0045421 840.02 2.11 0.1055 0.83666 0.0033837 
     3617971 -0.19798 0.0046642 898.56 2.11 0.1055 0.83838 0.0034819 
     3719143 -0.21442 0.0047052 957.17 2.11 0.1055 0.82619 0.0034614 
     3820316 -0.20885 0.0044931 1015.8 2.11 0.1055 0.8303 0.0033218 
     3921488 -0.21393 0.0049186 1074.4 2.11 0.1055 0.82655 0.00362 
     4022660 -0.20685 0.004423 1133 2.11 0.1055 0.83179 0.0032758 
     4123832 -0.20802 0.0046979 1191.6 2.11 0.1055 0.83092 0.0034758 
     4225003 -0.19848 0.0045953 1250.2 2.11 0.1055 0.838 0.0034289 
     4326175 -0.21117 0.0044567 1308.8 2.11 0.1055 0.82859 0.0032881 
     4427347 -0.21283 0.004137 1367.4 2.11 0.1055 0.82736 0.0030477 
     4528520 -0.2042 0.0044587 1426 2.11 0.1055 0.83375 0.0033101 
     4629692 -0.2112 0.0042852 1484.6 2.11 0.1055 0.82857 0.0031615 
     4730864 -0.20319 0.0043483 1543.2 2.11 0.1055 0.8345 0.003231 
     4832036 -0.20752 0.0044297 1601.8 2.11 0.1055 0.83129 0.0032788 
     4933207 -0.20654 0.0043188 1660.4 2.11 0.1055 0.83201 0.0031995 
     5034380 -0.20126 0.0046375 1719 2.11 0.1055 0.83593 0.0034518 
     5135551 -0.20924 0.0042871 1777.6 2.11 0.1055 0.83001 0.0031684 
     5236724 -0.21323 0.0045471 1836.2 2.11 0.1055 0.82707 0.0033487 
     5337895 -0.21324 0.0045354 1894.7 2.11 0.1055 0.82706 0.00334 
     5439067 -0.19905 0.0044141 1953.4 2.11 0.1055 0.83758 0.003292 
     5540239 -0.1991 0.0047441 2012 2.11 0.1055 0.83754 0.003538 
     5641411 -0.20359 0.0050136 2070.5 2.11 0.1055 0.8342 0.003724 
     5742583 -0.21032 0.0049474 2129.1 2.11 0.1055 0.82922 0.0036529 
     5843755 -0.20689 0.0048203 2187.8 2.11 0.1055 0.83176 0.00357 
     5944927 -0.21075 0.0052337 2246.4 2.11 0.1055 0.8289 0.0038628 
     6046099 -0.19956 0.0047827 2304.9 2.11 0.1055 0.8372 0.0035653 
  • explore/jitter.py

    • Property mode changed from 100644 to 100755
    r85190c2 rff10479  
     1#!/usr/bin/env python 
    12""" 
    23Application to explore the difference between sasview 3.x orientation 
    34dispersity and possible replacement algorithms. 
    45""" 
     6from __future__ import division, print_function 
     7 
     8import sys, os 
     9sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.realpath(__file__)))) 
     10 
     11import argparse 
     12 
    513import mpl_toolkits.mplot3d   # Adds projection='3d' option to subplot 
    614import matplotlib.pyplot as plt 
    715from matplotlib.widgets import Slider, CheckButtons 
    816from matplotlib import cm 
    9  
    1017import numpy as np 
    1118from numpy import pi, cos, sin, sqrt, exp, degrees, radians 
    1219 
    13 def draw_beam(ax): 
     20def draw_beam(ax, view=(0, 0)): 
     21    """ 
     22    Draw the beam going from source at (0, 0, 1) to detector at (0, 0, -1) 
     23    """ 
    1424    #ax.plot([0,0],[0,0],[1,-1]) 
    1525    #ax.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8) 
     
    2232    x = r*np.outer(np.cos(u), np.ones_like(v)) 
    2333    y = r*np.outer(np.sin(u), np.ones_like(v)) 
    24     z = np.outer(np.ones_like(u), v) 
     34    z = 1.3*np.outer(np.ones_like(u), v) 
     35 
     36    theta, phi = view 
     37    shape = x.shape 
     38    points = np.matrix([x.flatten(), y.flatten(), z.flatten()]) 
     39    points = Rz(phi)*Ry(theta)*points 
     40    x, y, z = [v.reshape(shape) for v in points] 
    2541 
    2642    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5) 
    2743 
    28 def draw_shimmy(ax, theta, phi, psi, dtheta, dphi, dpsi): 
    29     size=[0.1, 0.4, 1.0] 
    30     view=[theta, phi, psi] 
    31     shimmy=[0,0,0] 
    32     #draw_shape = draw_parallelepiped 
    33     draw_shape = draw_ellipsoid 
     44def draw_jitter(ax, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0)): 
     45    """ 
     46    Represent jitter as a set of shapes at different orientations. 
     47    """ 
     48    # set max diagonal to 0.95 
     49    scale = 0.95/sqrt(sum(v**2 for v in size)) 
     50    size = tuple(scale*v for v in size) 
     51    draw_shape = draw_parallelepiped 
     52    #draw_shape = draw_ellipsoid 
    3453 
    3554    #np.random.seed(10) 
     
    6483        [ 1,  1,  1], 
    6584    ] 
     85    dtheta, dphi, dpsi = jitter 
    6686    if dtheta == 0: 
    6787        cloud = [v for v in cloud if v[0] == 0] 
     
    7090    if dpsi == 0: 
    7191        cloud = [v for v in cloud if v[2] == 0] 
    72     draw_shape(ax, size, view, shimmy, steps=100, alpha=0.8) 
     92    draw_shape(ax, size, view, [0, 0, 0], steps=100, alpha=0.8) 
     93    scale = {'gaussian':1, 'rectangle':1/sqrt(3), 'uniform':1/3}[dist] 
    7394    for point in cloud: 
    74         shimmy=[dtheta*point[0], dphi*point[1], dpsi*point[2]] 
    75         draw_shape(ax, size, view, shimmy, alpha=0.8) 
     95        delta = [scale*dtheta*point[0], scale*dphi*point[1], scale*dpsi*point[2]] 
     96        draw_shape(ax, size, view, delta, alpha=0.8) 
    7697    for v in 'xyz': 
    7798        a, b, c = size 
     
    80101        getattr(ax, v+'axis').label.set_text(v) 
    81102 
    82 def draw_ellipsoid(ax, size, view, shimmy, steps=25, alpha=1): 
     103def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1): 
     104    """Draw an ellipsoid.""" 
    83105    a,b,c = size 
    84     theta, phi, psi = view 
    85     dtheta, dphi, dpsi = shimmy 
    86  
    87106    u = np.linspace(0, 2 * np.pi, steps) 
    88107    v = np.linspace(0, np.pi, steps) 
     
    90109    y = b*np.outer(np.sin(u), np.sin(v)) 
    91110    z = c*np.outer(np.ones_like(u), np.cos(v)) 
    92  
    93     shape = x.shape 
    94     points = np.matrix([x.flatten(),y.flatten(),z.flatten()]) 
    95     points = Rz(dpsi)*Ry(dtheta)*Rx(dphi)*points 
    96     points = Rz(phi)*Ry(theta)*Rz(psi)*points 
    97     x,y,z = [v.reshape(shape) for v in points] 
     111    x, y, z = transform_xyz(view, jitter, x, y, z) 
    98112 
    99113    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha) 
    100114 
    101 def draw_parallelepiped(ax, size, view, shimmy, alpha=1): 
     115    draw_labels(ax, view, jitter, [ 
     116         ('c+', [ 0, 0, c], [ 1, 0, 0]), 
     117         ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
     118         ('a+', [ a, 0, 0], [ 0, 0, 1]), 
     119         ('a-', [-a, 0, 0], [ 0, 0,-1]), 
     120         ('b+', [ 0, b, 0], [-1, 0, 0]), 
     121         ('b-', [ 0,-b, 0], [-1, 0, 0]), 
     122    ]) 
     123 
     124def draw_parallelepiped(ax, size, view, jitter, steps=None, alpha=1): 
     125    """Draw a parallelepiped.""" 
    102126    a,b,c = size 
    103     theta, phi, psi = view 
    104     dtheta, dphi, dpsi = shimmy 
    105  
    106127    x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1]) 
    107128    y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1]) 
     
    118139    ]) 
    119140 
    120     points = np.matrix([x,y,z]) 
    121     points = Rz(dpsi)*Ry(dtheta)*Rx(dphi)*points 
    122     points = Rz(phi)*Ry(theta)*Rz(psi)*points 
    123  
    124     x,y,z = [np.array(v).flatten() for v in points] 
     141    x, y, z = transform_xyz(view, jitter, x, y, z) 
    125142    ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha) 
    126143 
     144    draw_labels(ax, view, jitter, [ 
     145         ('c+', [ 0, 0, c], [ 1, 0, 0]), 
     146         ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
     147         ('a+', [ a, 0, 0], [ 0, 0, 1]), 
     148         ('a-', [-a, 0, 0], [ 0, 0,-1]), 
     149         ('b+', [ 0, b, 0], [-1, 0, 0]), 
     150         ('b-', [ 0,-b, 0], [-1, 0, 0]), 
     151    ]) 
     152 
    127153def draw_sphere(ax, radius=10., steps=100): 
     154    """Draw a sphere""" 
    128155    u = np.linspace(0, 2 * np.pi, steps) 
    129156    v = np.linspace(0, np.pi, steps) 
     
    134161    ax.plot_surface(x, y, z, rstride=4, cstride=4, color='w') 
    135162 
    136 def draw_mesh_new(ax, theta, dtheta, phi, dphi, flow, radius=10., dist='gauss'): 
    137     theta_center = radians(theta) 
    138     phi_center = radians(phi) 
    139     flow_center = radians(flow) 
    140     dtheta = radians(dtheta) 
    141     dphi = radians(dphi) 
    142  
    143     # 10 point 3-sigma gaussian weights 
    144     t = np.linspace(-3., 3., 11) 
    145     if dist == 'gauss': 
     163PROJECTIONS = [ 
     164    # in order of PROJECTION number; do not change without updating the 
     165    # constants in kernel_iq.c 
     166    'equirectangular', 'sinusoidal', 'guyou', 'azimuthal_equidistance', 
     167] 
     168def draw_mesh(ax, view, jitter, radius=1.2, n=11, dist='gaussian', 
     169              projection='equirectangular'): 
     170    """ 
     171    Draw the dispersion mesh showing the theta-phi orientations at which 
     172    the model will be evaluated. 
     173 
     174    jitter projections 
     175    <https://en.wikipedia.org/wiki/List_of_map_projections> 
     176 
     177    equirectangular (standard latitude-longitude mesh) 
     178        <https://en.wikipedia.org/wiki/Equirectangular_projection> 
     179        Allows free movement in phi (around the equator), but theta is 
     180        limited to +/- 90, and points are cos-weighted. Jitter in phi is 
     181        uniform in weight along a line of latitude.  With small theta and 
     182        phi ranging over +/- 180 this forms a wobbling disk.  With small 
     183        phi and theta ranging over +/- 90 this forms a wedge like a slice 
     184        of an orange. 
     185    azimuthal_equidistance (Postel) 
     186        <https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection> 
     187        Preserves distance from center, and so is an excellent map for 
     188        representing a bivariate gaussian on the surface.  Theta and phi 
     189        operate identically, cutting wegdes from the antipode of the viewing 
     190        angle.  This unfortunately does not allow free movement in either 
     191        theta or phi since the orthogonal wobble decreases to 0 as the body 
     192        rotates through 180 degrees. 
     193    sinusoidal (Sanson-Flamsteed, Mercator equal-area) 
     194        <https://en.wikipedia.org/wiki/Sinusoidal_projection> 
     195        Preserves arc length with latitude, giving bad behaviour at 
     196        theta near +/- 90.  Theta and phi operate somewhat differently, 
     197        so a system with a-b-c dtheta-dphi-dpsi will not give the same 
     198        value as one with b-a-c dphi-dtheta-dpsi, as would be the case 
     199        for azimuthal equidistance.  Free movement using theta or phi 
     200        uniform over +/- 180 will work, but not as well as equirectangular 
     201        phi, with theta being slightly worse.  Computationally it is much 
     202        cheaper for wide theta-phi meshes since it excludes points which 
     203        lie outside the sinusoid near theta +/- 90 rather than packing 
     204        them close together as in equirectangle.  Note that the poles 
     205        will be slightly overweighted for theta > 90 with the circle 
     206        from theta at 90+dt winding backwards around the pole, overlapping 
     207        the circle from theta at 90-dt. 
     208    Guyou (hemisphere-in-a-square) **not weighted** 
     209        <https://en.wikipedia.org/wiki/Guyou_hemisphere-in-a-square_projection> 
     210        With tiling, allows rotation in phi or theta through +/- 180, with 
     211        uniform spacing.  Both theta and phi allow free rotation, with wobble 
     212        in the orthogonal direction reasonably well behaved (though not as 
     213        good as equirectangular phi). The forward/reverse transformations 
     214        relies on elliptic integrals that are somewhat expensive, so the 
     215        behaviour has to be very good to justify the cost and complexity. 
     216        The weighting function for each point has not yet been computed. 
     217        Note: run the module *guyou.py* directly and it will show the forward 
     218        and reverse mappings. 
     219    azimuthal_equal_area  **incomplete** 
     220        <https://en.wikipedia.org/wiki/Lambert_azimuthal_equal-area_projection> 
     221        Preserves the relative density of the surface patches.  Not that 
     222        useful and not completely implemented 
     223    Gauss-Kreuger **not implemented** 
     224        <https://en.wikipedia.org/wiki/Transverse_Mercator_projection#Ellipsoidal_transverse_Mercator> 
     225        Should allow free movement in theta, but phi is distorted. 
     226    """ 
     227    theta, phi, psi = view 
     228    dtheta, dphi, dpsi = jitter 
     229 
     230    t = np.linspace(-1, 1, n) 
     231    weights = np.ones_like(t) 
     232    if dist == 'gaussian': 
     233        t *= 3 
    146234        weights = exp(-0.5*t**2) 
    147     elif dist == 'rect': 
    148         weights = np.ones_like(t) 
     235    elif dist == 'rectangle': 
     236        # Note: uses sasmodels ridiculous definition of rectangle width 
     237        t *= sqrt(3) 
     238    elif dist == 'uniform': 
     239        pass 
    149240    else: 
    150         raise ValueError("expected dist to be 'gauss' or 'rect'") 
    151     theta = dtheta*t 
    152     phi = dphi*t 
    153  
    154     x = radius * np.outer(cos(phi), cos(theta)) 
    155     y = radius * np.outer(sin(phi), cos(theta)) 
    156     z = radius * np.outer(np.ones_like(phi), sin(theta)) 
    157     #w = np.outer(weights, weights*abs(cos(dtheta*t))) 
    158     w = np.outer(weights, weights*abs(cos(theta))) 
    159  
    160     x, y, z, w = [v.flatten() for v in (x,y,z,w)] 
    161     x, y, z = rotate(x, y, z, phi_center, theta_center, flow_center) 
    162  
     241        raise ValueError("expected dist to be gaussian, rectangle or uniform") 
     242 
     243    if projection == 'equirectangular':  #define PROJECTION 1 
     244        def rotate(theta_i, phi_j): 
     245            return Rx(phi_j)*Ry(theta_i) 
     246        def weight(theta_i, phi_j, wi, wj): 
     247            return wi*wj*abs(cos(radians(theta_i))) 
     248    elif projection == 'sinusoidal':  #define PROJECTION 2 
     249        def rotate(theta_i, phi_j): 
     250            latitude = theta_i 
     251            scale = cos(radians(latitude)) 
     252            longitude = phi_j/scale if abs(phi_j) < abs(scale)*180 else 0 
     253            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
     254            return Rx(longitude)*Ry(latitude) 
     255        def weight(theta_i, phi_j, wi, wj): 
     256            latitude = theta_i 
     257            scale = cos(radians(latitude)) 
     258            w = 1 if abs(phi_j) < abs(scale)*180 else 0 
     259            return w*wi*wj 
     260    elif projection == 'guyou':  #define PROJECTION 3  (eventually?) 
     261        def rotate(theta_i, phi_j): 
     262            from guyou import guyou_invert 
     263            #latitude, longitude = guyou_invert([theta_i], [phi_j]) 
     264            longitude, latitude = guyou_invert([phi_j], [theta_i]) 
     265            return Rx(longitude[0])*Ry(latitude[0]) 
     266        def weight(theta_i, phi_j, wi, wj): 
     267            return wi*wj 
     268    elif projection == 'azimuthal_equidistance':  # Note: Rz Ry, not Rx Ry 
     269        def rotate(theta_i, phi_j): 
     270            latitude = sqrt(theta_i**2 + phi_j**2) 
     271            longitude = degrees(np.arctan2(phi_j, theta_i)) 
     272            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
     273            return Rz(longitude)*Ry(latitude) 
     274        def weight(theta_i, phi_j, wi, wj): 
     275            # Weighting for each point comes from the integral: 
     276            #     \int\int I(q, lat, log) sin(lat) dlat dlog 
     277            # We are doing a conformal mapping from disk to sphere, so we need 
     278            # a change of variables g(theta, phi) -> (lat, long): 
     279            #     lat, long = sqrt(theta^2 + phi^2), arctan(phi/theta) 
     280            # giving: 
     281            #     dtheta dphi = det(J) dlat dlong 
     282            # where J is the jacobian from the partials of g. Using 
     283            #     R = sqrt(theta^2 + phi^2), 
     284            # then 
     285            #     J = [[x/R, Y/R], -y/R^2, x/R^2]] 
     286            # and 
     287            #     det(J) = 1/R 
     288            # with the final integral being: 
     289            #    \int\int I(q, theta, phi) sin(R)/R dtheta dphi 
     290            # 
     291            # This does approximately the right thing, decreasing the weight 
     292            # of each point as you go farther out on the disk, but it hasn't 
     293            # yet been checked against the 1D integral results. Prior 
     294            # to declaring this "good enough" and checking that integrals 
     295            # work in practice, we will examine alternative mappings. 
     296            # 
     297            # The issue is that the mapping does not support the case of free 
     298            # rotation about a single axis correctly, with a small deviation 
     299            # in the orthogonal axis independent of the first axis.  Like the 
     300            # usual polar coordiates integration, the integrated sections 
     301            # form wedges, though at least in this case the wedge cuts through 
     302            # the entire sphere, and treats theta and phi identically. 
     303            latitude = sqrt(theta_i**2 + phi_j**2) 
     304            w = sin(radians(latitude))/latitude if latitude != 0 else 1 
     305            return w*wi*wj if latitude < 180 else 0 
     306    elif projection == 'azimuthal_equal_area': 
     307        def rotate(theta_i, phi_j): 
     308            R = min(1, sqrt(theta_i**2 + phi_j**2)/180) 
     309            latitude = 180-degrees(2*np.arccos(R)) 
     310            longitude = degrees(np.arctan2(phi_j, theta_i)) 
     311            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
     312            return Rz(longitude)*Ry(latitude) 
     313        def weight(theta_i, phi_j, wi, wj): 
     314            latitude = sqrt(theta_i**2 + phi_j**2) 
     315            w = sin(radians(latitude))/latitude if latitude != 0 else 1 
     316            return w*wi*wj if latitude < 180 else 0 
     317    else: 
     318        raise ValueError("unknown projection %r"%projection) 
     319 
     320    # mesh in theta, phi formed by rotating z 
     321    z = np.matrix([[0], [0], [radius]]) 
     322    points = np.hstack([rotate(theta_i, phi_j)*z 
     323                        for theta_i in dtheta*t 
     324                        for phi_j in dphi*t]) 
     325    # select just the active points (i.e., those with phi < 180 
     326    w = np.array([weight(theta_i, phi_j, wi, wj) 
     327                  for wi, theta_i in zip(weights, dtheta*t) 
     328                  for wj, phi_j in zip(weights, dphi*t)]) 
     329    #print(max(w), min(w), min(w[w>0])) 
     330    points = points[:, w>0] 
     331    w = w[w>0] 
     332    w /= max(w) 
     333 
     334    if 0: # Kent distribution 
     335        points = np.hstack([Rx(phi_j)*Ry(theta_i)*z for theta_i in 30*t for phi_j in 60*t]) 
     336        xp, yp, zp = [np.array(v).flatten() for v in points] 
     337        kappa = max(1e6, radians(dtheta)/(2*pi)) 
     338        beta = 1/max(1e-6, radians(dphi)/(2*pi))/kappa 
     339        w = exp(kappa*zp) #+ beta*(xp**2 + yp**2) 
     340        print(kappa, dtheta, radians(dtheta), min(w), max(w), sum(w)) 
     341        #w /= abs(cos(radians( 
     342        #w /= sum(w) 
     343 
     344    # rotate relative to beam 
     345    points = orient_relative_to_beam(view, points) 
     346 
     347    x, y, z = [np.array(v).flatten() for v in points] 
     348    #plt.figure(2); plt.clf(); plt.hist(z, bins=np.linspace(-1, 1, 51)) 
    163349    ax.scatter(x, y, z, c=w, marker='o', vmin=0., vmax=1.) 
    164350 
    165 def rotate(x, y, z, phi, theta, psi): 
    166     R = Rz(phi)*Ry(theta)*Rz(psi) 
    167     p = np.vstack([x,y,z]) 
    168     return R*p 
    169  
     351def draw_labels(ax, view, jitter, text): 
     352    """ 
     353    Draw text at a particular location. 
     354    """ 
     355    labels, locations, orientations = zip(*text) 
     356    px, py, pz = zip(*locations) 
     357    dx, dy, dz = zip(*orientations) 
     358 
     359    px, py, pz = transform_xyz(view, jitter, px, py, pz) 
     360    dx, dy, dz = transform_xyz(view, jitter, dx, dy, dz) 
     361 
     362    # TODO: zdir for labels is broken, and labels aren't appearing. 
     363    for label, p, zdir in zip(labels, zip(px, py, pz), zip(dx, dy, dz)): 
     364        zdir = np.asarray(zdir).flatten() 
     365        ax.text(p[0], p[1], p[2], label, zdir=zdir) 
     366 
     367# Definition of rotation matrices comes from wikipedia: 
     368#    https://en.wikipedia.org/wiki/Rotation_matrix#Basic_rotations 
    170369def Rx(angle): 
     370    """Construct a matrix to rotate points about *x* by *angle* degrees.""" 
    171371    a = radians(angle) 
    172     R = [[1., 0., 0.], 
    173          [0.,  cos(a), sin(a)], 
    174          [0., -sin(a), cos(a)]] 
     372    R = [[1, 0, 0], 
     373         [0, +cos(a), -sin(a)], 
     374         [0, +sin(a), +cos(a)]] 
    175375    return np.matrix(R) 
    176376 
    177377def Ry(angle): 
     378    """Construct a matrix to rotate points about *y* by *angle* degrees.""" 
    178379    a = radians(angle) 
    179     R = [[cos(a), 0., -sin(a)], 
    180          [0., 1., 0.], 
    181          [sin(a), 0.,  cos(a)]] 
     380    R = [[+cos(a), 0, +sin(a)], 
     381         [0, 1, 0], 
     382         [-sin(a), 0, +cos(a)]] 
    182383    return np.matrix(R) 
    183384 
    184385def Rz(angle): 
     386    """Construct a matrix to rotate points about *z* by *angle* degrees.""" 
    185387    a = radians(angle) 
    186     R = [[cos(a), -sin(a), 0.], 
    187          [sin(a),  cos(a), 0.], 
    188          [0., 0., 1.]] 
     388    R = [[+cos(a), -sin(a), 0], 
     389         [+sin(a), +cos(a), 0], 
     390         [0, 0, 1]] 
    189391    return np.matrix(R) 
    190392 
    191 def main(): 
     393def transform_xyz(view, jitter, x, y, z): 
     394    """ 
     395    Send a set of (x,y,z) points through the jitter and view transforms. 
     396    """ 
     397    x, y, z = [np.asarray(v) for v in (x, y, z)] 
     398    shape = x.shape 
     399    points = np.matrix([x.flatten(),y.flatten(),z.flatten()]) 
     400    points = apply_jitter(jitter, points) 
     401    points = orient_relative_to_beam(view, points) 
     402    x, y, z = [np.array(v).reshape(shape) for v in points] 
     403    return x, y, z 
     404 
     405def apply_jitter(jitter, points): 
     406    """ 
     407    Apply the jitter transform to a set of points. 
     408 
     409    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple. 
     410    """ 
     411    dtheta, dphi, dpsi = jitter 
     412    points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points 
     413    return points 
     414 
     415def orient_relative_to_beam(view, points): 
     416    """ 
     417    Apply the view transform to a set of points. 
     418 
     419    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple. 
     420    """ 
     421    theta, phi, psi = view 
     422    points = Rz(phi)*Ry(theta)*Rz(psi)*points 
     423    return points 
     424 
     425# translate between number of dimension of dispersity and the number of 
     426# points along each dimension. 
     427PD_N_TABLE = { 
     428    (0, 0, 0): (0, 0, 0),     # 0 
     429    (1, 0, 0): (100, 0, 0),   # 100 
     430    (0, 1, 0): (0, 100, 0), 
     431    (0, 0, 1): (0, 0, 100), 
     432    (1, 1, 0): (30, 30, 0),   # 900 
     433    (1, 0, 1): (30, 0, 30), 
     434    (0, 1, 1): (0, 30, 30), 
     435    (1, 1, 1): (15, 15, 15),  # 3375 
     436} 
     437 
     438def clipped_range(data, portion=1.0, mode='central'): 
     439    """ 
     440    Determine range from data. 
     441 
     442    If *portion* is 1, use full range, otherwise use the center of the range 
     443    or the top of the range, depending on whether *mode* is 'central' or 'top'. 
     444    """ 
     445    if portion == 1.0: 
     446        return data.min(), data.max() 
     447    elif mode == 'central': 
     448        data = np.sort(data.flatten()) 
     449        offset = int(portion*len(data)/2 + 0.5) 
     450        return data[offset], data[-offset] 
     451    elif mode == 'top': 
     452        data = np.sort(data.flatten()) 
     453        offset = int(portion*len(data) + 0.5) 
     454        return data[offset], data[-1] 
     455 
     456def draw_scattering(calculator, ax, view, jitter, dist='gaussian'): 
     457    """ 
     458    Plot the scattering for the particular view. 
     459 
     460    *calculator* is returned from :func:`build_model`.  *ax* are the 3D axes 
     461    on which the data will be plotted.  *view* and *jitter* are the current 
     462    orientation and orientation dispersity.  *dist* is one of the sasmodels 
     463    weight distributions. 
     464    """ 
     465    if dist == 'uniform':  # uniform is not yet in this branch 
     466        dist, scale = 'rectangle', 1/sqrt(3) 
     467    else: 
     468        scale = 1 
     469 
     470    # add the orientation parameters to the model parameters 
     471    theta, phi, psi = view 
     472    theta_pd, phi_pd, psi_pd = [scale*v for v in jitter] 
     473    theta_pd_n, phi_pd_n, psi_pd_n = PD_N_TABLE[(theta_pd>0, phi_pd>0, psi_pd>0)] 
     474    ## increase pd_n for testing jitter integration rather than simple viz 
     475    #theta_pd_n, phi_pd_n, psi_pd_n = [5*v for v in (theta_pd_n, phi_pd_n, psi_pd_n)] 
     476 
     477    pars = dict( 
     478        theta=theta, theta_pd=theta_pd, theta_pd_type=dist, theta_pd_n=theta_pd_n, 
     479        phi=phi, phi_pd=phi_pd, phi_pd_type=dist, phi_pd_n=phi_pd_n, 
     480        psi=psi, psi_pd=psi_pd, psi_pd_type=dist, psi_pd_n=psi_pd_n, 
     481    ) 
     482    pars.update(calculator.pars) 
     483 
     484    # compute the pattern 
     485    qx, qy = calculator._data.x_bins, calculator._data.y_bins 
     486    Iqxy = calculator(**pars).reshape(len(qx), len(qy)) 
     487 
     488    # scale it and draw it 
     489    Iqxy = np.log(Iqxy) 
     490    if calculator.limits: 
     491        # use limits from orientation (0,0,0) 
     492        vmin, vmax = calculator.limits 
     493    else: 
     494        vmax = Iqxy.max() 
     495        vmin = vmax*10**-7 
     496        #vmin, vmax = clipped_range(Iqxy, portion=portion, mode='top') 
     497    #print("range",(vmin,vmax)) 
     498    #qx, qy = np.meshgrid(qx, qy) 
     499    if 0: 
     500        level = np.asarray(255*(Iqxy - vmin)/(vmax - vmin), 'i') 
     501        level[level<0] = 0 
     502        colors = plt.get_cmap()(level) 
     503        ax.plot_surface(qx, qy, -1.1, rstride=1, cstride=1, facecolors=colors) 
     504    elif 1: 
     505        ax.contourf(qx/qx.max(), qy/qy.max(), Iqxy, zdir='z', offset=-1.1, 
     506                    levels=np.linspace(vmin, vmax, 24)) 
     507    else: 
     508        ax.pcolormesh(qx, qy, Iqxy) 
     509 
     510def build_model(model_name, n=150, qmax=0.5, **pars): 
     511    """ 
     512    Build a calculator for the given shape. 
     513 
     514    *model_name* is any sasmodels model.  *n* and *qmax* define an n x n mesh 
     515    on which to evaluate the model.  The remaining parameters are stored in 
     516    the returned calculator as *calculator.pars*.  They are used by 
     517    :func:`draw_scattering` to set the non-orientation parameters in the 
     518    calculation. 
     519 
     520    Returns a *calculator* function which takes a dictionary or parameters and 
     521    produces Iqxy.  The Iqxy value needs to be reshaped to an n x n matrix 
     522    for plotting.  See the :class:`sasmodels.direct_model.DirectModel` class 
     523    for details. 
     524    """ 
     525    from sasmodels.core import load_model_info, build_model 
     526    from sasmodels.data import empty_data2D 
     527    from sasmodels.direct_model import DirectModel 
     528 
     529    model_info = load_model_info(model_name) 
     530    model = build_model(model_info) #, dtype='double!') 
     531    q = np.linspace(-qmax, qmax, n) 
     532    data = empty_data2D(q, q) 
     533    calculator = DirectModel(data, model) 
     534 
     535    # stuff the values for non-orientation parameters into the calculator 
     536    calculator.pars = pars.copy() 
     537    calculator.pars.setdefault('backgound', 1e-3) 
     538 
     539    # fix the data limits so that we can see if the pattern fades 
     540    # under rotation or angular dispersion 
     541    Iqxy = calculator(theta=0, phi=0, psi=0, **calculator.pars) 
     542    Iqxy = np.log(Iqxy) 
     543    vmin, vmax = clipped_range(Iqxy, 0.95, mode='top') 
     544    calculator.limits = vmin, vmax+1 
     545 
     546    return calculator 
     547 
     548def select_calculator(model_name, n=150, size=(10,40,100)): 
     549    """ 
     550    Create a model calculator for the given shape. 
     551 
     552    *model_name* is one of sphere, cylinder, ellipsoid, triaxial_ellipsoid, 
     553    parallelepiped or bcc_paracrystal. *n* is the number of points to use 
     554    in the q range.  *qmax* is chosen based on model parameters for the 
     555    given model to show something intersting. 
     556 
     557    Returns *calculator* and tuple *size* (a,b,c) giving minor and major 
     558    equitorial axes and polar axis respectively.  See :func:`build_model` 
     559    for details on the returned calculator. 
     560    """ 
     561    a, b, c = size 
     562    if model_name == 'sphere': 
     563        calculator = build_model('sphere', n=n, radius=c) 
     564        a = b = c 
     565    elif model_name == 'bcc_paracrystal': 
     566        calculator = build_model('bcc_paracrystal', n=n, dnn=c, 
     567                                  d_factor=0.06, radius=40) 
     568        a = b = c 
     569    elif model_name == 'cylinder': 
     570        calculator = build_model('cylinder', n=n, qmax=0.3, radius=b, length=c) 
     571        a = b 
     572    elif model_name == 'ellipsoid': 
     573        calculator = build_model('ellipsoid', n=n, qmax=1.0, 
     574                                 radius_polar=c, radius_equatorial=b) 
     575        a = b 
     576    elif model_name == 'triaxial_ellipsoid': 
     577        calculator = build_model('triaxial_ellipsoid', n=n, qmax=0.5, 
     578                                 radius_equat_minor=a, 
     579                                 radius_equat_major=b, 
     580                                 radius_polar=c) 
     581    elif model_name == 'parallelepiped': 
     582        calculator = build_model('parallelepiped', n=n, a=a, b=b, c=c) 
     583    else: 
     584        raise ValueError("unknown model %s"%model_name) 
     585 
     586    return calculator, (a, b, c) 
     587 
     588SHAPES = [ 
     589    'parallelepiped', 'triaxial_ellipsoid', 'bcc_paracrystal', 
     590    'cylinder', 'ellipsoid', 
     591    'sphere', 
     592 ] 
     593 
     594DISTRIBUTIONS = [ 
     595    'gaussian', 'rectangle', 'uniform', 
     596] 
     597DIST_LIMITS = { 
     598    'gaussian': 30, 
     599    'rectangle': 90/sqrt(3), 
     600    'uniform': 90, 
     601} 
     602 
     603def run(model_name='parallelepiped', size=(10, 40, 100), 
     604        dist='gaussian', mesh=30, 
     605        projection='equirectangular'): 
     606    """ 
     607    Show an interactive orientation and jitter demo. 
     608 
     609    *model_name* is one of the models available in :func:`select_model`. 
     610    """ 
     611    # projection number according to 1-order position in list, but 
     612    # only 1 and 2 are implemented so far. 
     613    from sasmodels import generate 
     614    generate.PROJECTION = PROJECTIONS.index(projection) + 1 
     615    if generate.PROJECTION > 2: 
     616        print("*** PROJECTION %s not implemented in scattering function ***"%projection) 
     617        generate.PROJECTION = 2 
     618 
     619    # set up calculator 
     620    calculator, size = select_calculator(model_name, n=150, size=size) 
     621 
     622    ## uncomment to set an independent the colour range for every view 
     623    ## If left commented, the colour range is fixed for all views 
     624    calculator.limits = None 
     625 
     626    ## initial view 
     627    #theta, dtheta = 70., 10. 
     628    #phi, dphi = -45., 3. 
     629    #psi, dpsi = -45., 3. 
     630    theta, phi, psi = 0, 0, 0 
     631    dtheta, dphi, dpsi = 0, 0, 0 
     632 
     633    ## create the plot window 
    192634    #plt.hold(True) 
     635    plt.subplots(num=None, figsize=(5.5, 5.5)) 
    193636    plt.set_cmap('gist_earth') 
    194637    plt.clf() 
     638    plt.gcf().canvas.set_window_title(projection) 
    195639    #gs = gridspec.GridSpec(2,1,height_ratios=[4,1]) 
    196640    #ax = plt.subplot(gs[0], projection='3d') 
    197641    ax = plt.axes([0.0, 0.2, 1.0, 0.8], projection='3d') 
    198  
    199     theta, dtheta = 70., 10. 
    200     phi, dphi = -45., 3. 
    201     psi, dpsi = -45., 3. 
    202     theta, phi, psi = 0, 0, 0 
    203     dtheta, dphi, dpsi = 0, 0, 0 
    204     #dist = 'rect' 
    205     dist = 'gauss' 
     642    try:  # CRUFT: not all versions of matplotlib accept 'square' 3d projection 
     643        ax.axis('square') 
     644    except Exception: 
     645        pass 
    206646 
    207647    axcolor = 'lightgoldenrodyellow' 
     648 
     649    ## add control widgets to plot 
    208650    axtheta  = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor) 
    209651    axphi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor) 
     
    212654    sphi = Slider(axphi, 'Phi', -180, 180, valinit=phi) 
    213655    spsi = Slider(axpsi, 'Psi', -180, 180, valinit=psi) 
     656 
    214657    axdtheta  = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor) 
    215658    axdphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor) 
    216659    axdpsi= plt.axes([0.75, 0.05, 0.15, 0.04], axisbg=axcolor) 
    217     sdtheta = Slider(axdtheta, 'dTheta', 0, 30, valinit=dtheta) 
    218     sdphi = Slider(axdphi, 'dPhi', 0, 30, valinit=dphi) 
    219     sdpsi = Slider(axdpsi, 'dPsi', 0, 30, valinit=dphi) 
    220  
     660    # Note: using ridiculous definition of rectangle distribution, whose width 
     661    # in sasmodels is sqrt(3) times the given width.  Divide by sqrt(3) to keep 
     662    # the maximum width to 90. 
     663    dlimit = DIST_LIMITS[dist] 
     664    sdtheta = Slider(axdtheta, 'dTheta', 0, 2*dlimit, valinit=dtheta) 
     665    sdphi = Slider(axdphi, 'dPhi', 0, 2*dlimit, valinit=dphi) 
     666    sdpsi = Slider(axdpsi, 'dPsi', 0, 2*dlimit, valinit=dpsi) 
     667 
     668    ## callback to draw the new view 
    221669    def update(val, axis=None): 
    222         theta, phi, psi = stheta.val, sphi.val, spsi.val 
    223         dtheta, dphi, dpsi = sdtheta.val, sdphi.val, sdpsi.val 
     670        view = stheta.val, sphi.val, spsi.val 
     671        jitter = sdtheta.val, sdphi.val, sdpsi.val 
     672        # set small jitter as 0 if multiple pd dims 
     673        dims = sum(v > 0 for v in jitter) 
     674        limit = [0, 0, 0.5, 5][dims] 
     675        jitter = [0 if v < limit else v for v in jitter] 
    224676        ax.cla() 
    225         draw_beam(ax) 
    226         draw_shimmy(ax, theta, phi, psi, dtheta, dphi, dpsi) 
    227         #if not axis.startswith('d'): 
    228         #    ax.view_init(elev=theta, azim=phi) 
     677        draw_beam(ax, (0, 0)) 
     678        draw_jitter(ax, view, jitter, dist=dist, size=size) 
     679        #draw_jitter(ax, view, (0,0,0)) 
     680        draw_mesh(ax, view, jitter, dist=dist, n=mesh, projection=projection) 
     681        draw_scattering(calculator, ax, view, jitter, dist=dist) 
    229682        plt.gcf().canvas.draw() 
    230683 
     684    ## bind control widgets to view updater 
    231685    stheta.on_changed(lambda v: update(v,'theta')) 
    232686    sphi.on_changed(lambda v: update(v, 'phi')) 
     
    236690    sdpsi.on_changed(lambda v: update(v, 'dpsi')) 
    237691 
     692    ## initialize view 
    238693    update(None, 'phi') 
    239694 
     695    ## go interactive 
    240696    plt.show() 
     697 
     698def main(): 
     699    parser = argparse.ArgumentParser( 
     700        description="Display jitter", 
     701        formatter_class=argparse.ArgumentDefaultsHelpFormatter, 
     702        ) 
     703    parser.add_argument('-p', '--projection', choices=PROJECTIONS, default=PROJECTIONS[0], help='coordinate projection') 
     704    parser.add_argument('-s', '--size', type=str, default='10,40,100', help='a,b,c lengths') 
     705    parser.add_argument('-d', '--distribution', choices=DISTRIBUTIONS, default=DISTRIBUTIONS[0], help='jitter distribution') 
     706    parser.add_argument('-m', '--mesh', type=int, default=30, help='#points in theta-phi mesh') 
     707    parser.add_argument('shape', choices=SHAPES, nargs='?', default=SHAPES[0], help='oriented shape') 
     708    opts = parser.parse_args() 
     709    size = tuple(int(v) for v in opts.size.split(',')) 
     710    run(opts.shape, size=size, 
     711        mesh=opts.mesh, dist=opts.distribution, 
     712        projection=opts.projection) 
    241713 
    242714if __name__ == "__main__": 
  • explore/precision.py

    r237c9cf ra1c5758  
    11#!/usr/bin/env python 
    22r""" 
    3 Show numerical precision of $2 J_1(x)/x$. 
     3Show numerical precision of various expressions. 
     4 
     5Evaluates the same function(s) in single and double precision and compares 
     6the results to 500 digit mpmath evaluation of the same function. 
     7 
     8Note: a quick way to generation C and python code for taylor series 
     9expansions from sympy: 
     10 
     11    import sympy as sp 
     12    x = sp.var("x") 
     13    f = sp.sin(x)/x 
     14    t = sp.series(f, n=12).removeO()  # taylor series with no O(x^n) term 
     15    p = sp.horner(t)   # Horner representation 
     16    p = p.replace(x**2, sp.var("xsq")  # simplify if alternate terms are zero 
     17    p = p.n(15)  # evaluate coefficients to 15 digits (optional) 
     18    c_code = sp.ccode(p, assign_to=sp.var("p"))  # convert to c code 
     19    py_code = c[:-1]  # strip semicolon to convert c to python 
     20 
     21    # mpmath has pade() rational function approximation, which might work 
     22    # better than the taylor series for some functions: 
     23    P, Q = mp.pade(sp.Poly(t.n(15),x).coeffs(), L, M) 
     24    P = sum(a*x**n for n,a in enumerate(reversed(P))) 
     25    Q = sum(a*x**n for n,a in enumerate(reversed(Q))) 
     26    c_code = sp.ccode(sp.horner(P)/sp.horner(Q), assign_to=sp.var("p")) 
     27 
     28    # There are richardson and shanks series accelerators in both sympy 
     29    # and mpmath that may be helpful. 
    430""" 
    531from __future__ import division, print_function 
     
    284310    np_function=scipy.special.erfc, 
    285311    ocl_function=make_ocl("return sas_erfc(q);", "sas_erfc", ["lib/polevl.c", "lib/sas_erf.c"]), 
     312    limits=(-5., 5.), 
     313) 
     314add_function( 
     315    name="expm1(x)", 
     316    mp_function=mp.expm1, 
     317    np_function=np.expm1, 
     318    ocl_function=make_ocl("return expm1(q);", "sas_expm1"), 
    286319    limits=(-5., 5.), 
    287320) 
     
    448481) 
    449482 
     483replacement_expm1 = """\ 
     484      double x = (double)q;  // go back to float for single precision kernels 
     485      // Adapted from the cephes math library. 
     486      // Copyright 1984 - 1992 by Stephen L. Moshier 
     487      if (x != x || x == 0.0) { 
     488         return x; // NaN and +/- 0 
     489      } else if (x < -0.5 || x > 0.5) { 
     490         return exp(x) - 1.0; 
     491      } else { 
     492         const double xsq = x*x; 
     493         const double p = ((( 
     494            +1.2617719307481059087798E-4)*xsq 
     495            +3.0299440770744196129956E-2)*xsq 
     496            +9.9999999999999999991025E-1); 
     497         const double q = (((( 
     498            +3.0019850513866445504159E-6)*xsq 
     499            +2.5244834034968410419224E-3)*xsq 
     500            +2.2726554820815502876593E-1)*xsq 
     501            +2.0000000000000000000897E0); 
     502         double r = x * p; 
     503         r =  r / (q - r); 
     504         return r+r; 
     505       } 
     506""" 
     507add_function( 
     508    name="sas_expm1(x)", 
     509    mp_function=mp.expm1, 
     510    np_function=np.expm1, 
     511    ocl_function=make_ocl(replacement_expm1, "sas_expm1"), 
     512) 
     513 
    450514# Alternate versions of 3 j1(x)/x, for posterity 
    451515def taylor_3j1x_x(x): 
  • extra/pylint.rc

    r823e620 rb669b49  
    2121# List of plugins (as comma separated values of python modules names) to load, 
    2222# usually to register additional checkers. 
    23 load-plugins=pylint_numpy,pylint_pyopencl 
     23load-plugins=pylint_numpy,pylint_pyopencl,pylint_sas 
    2424 
    2525# Use multiple processes to speed up Pylint. 
  • sasmodels/__init__.py

    r997c9ca re65c3ba  
    4747        ] 
    4848    return return_list 
    49  
    50  
  • sasmodels/bumps_model.py

    r74b0495 r2d81cfe  
    2020from .direct_model import DataMixin 
    2121 
     22# pylint: disable=unused-import 
    2223try: 
    2324    from typing import Dict, Union, Tuple, Any 
     
    2829except ImportError: 
    2930    pass 
     31# pylint: enable=unused-import 
    3032 
    3133try: 
     
    3739 
    3840 
    39 def create_parameters(model_info, **kwargs): 
    40     # type: (ModelInfo, **Union[float, str, Parameter]) -> Tuple[Dict[str, Parameter], Dict[str, str]] 
     41def create_parameters(model_info,  # type: ModelInfo 
     42                      **kwargs     # type: Union[float, str, Parameter] 
     43                     ): 
     44    # type: (...) -> Tuple[Dict[str, Parameter], Dict[str, str]] 
    4145    """ 
    4246    Generate Bumps parameters from the model info. 
     
    238242        # pylint: disable=attribute-defined-outside-init 
    239243        self.__dict__ = state 
    240  
  • sasmodels/compare.py

    rced5bd2 r2d81cfe  
    4040from . import core 
    4141from . import kerneldll 
    42 from . import exception 
    4342from .data import plot_theory, empty_data1D, empty_data2D, load_data 
    44 from .direct_model import DirectModel 
    45 from .convert import revert_name, revert_pars, constrain_new_to_old 
     43from .direct_model import DirectModel, get_mesh 
    4644from .generate import FLOAT_RE 
    47  
     45from .weights import plot_weights 
     46 
     47# pylint: disable=unused-import 
    4848try: 
    4949    from typing import Optional, Dict, Any, Callable, Tuple 
    50 except Exception: 
     50except ImportError: 
    5151    pass 
    5252else: 
     
    5454    from .data import Data 
    5555    Calculator = Callable[[float], np.ndarray] 
     56# pylint: enable=unused-import 
    5657 
    5758USAGE = """ 
     
    8384    -pars/-nopars* prints the parameter set or not 
    8485    -default/-demo* use demo vs default parameters 
     86    -sphere[=150] set up spherical integration over theta/phi using n points 
    8587 
    8688    === calculation options === 
    87     -mono*/-poly force monodisperse or allow polydisperse demo parameters 
     89    -mono*/-poly force monodisperse or allow polydisperse random parameters 
    8890    -cutoff=1e-5* cutoff value for including a point in polydispersity 
    8991    -magnetic/-nonmagnetic* suppress magnetism 
     
    9294 
    9395    === precision options === 
    94     -calc=default uses the default calcution precision 
     96    -engine=default uses the default calcution precision 
    9597    -single/-double/-half/-fast sets an OpenCL calculation engine 
    9698    -single!/-double!/-quad! sets an OpenMP calculation engine 
    97     -sasview sets the sasview calculation engine 
    9899 
    99100    === plotting === 
     
    103104    -abs/-rel* plot relative or absolute error 
    104105    -title="note" adds note to the plot title, after the model name 
     106    -weights shows weights plots for the polydisperse parameters 
    105107 
    106108    === output options === 
     
    111113vary from 64-bit to 128-bit, with 80-bit floats being common (1e-19 precision). 
    112114On unix and mac you may need single quotes around the DLL computation 
    113 engines, such as -calc='single!,double!' since !, is treated as a history 
     115engines, such as -engine='single!,double!' since !, is treated as a history 
    114116expansion request in the shell. 
    115117 
     
    123125 
    124126    # compare single and double precision calculation for a barbell 
    125     sascomp barbell -calc=single,double 
     127    sascomp barbell -engine=single,double 
    126128 
    127129    # generate 10 random lorentz models, with seed=27 
     
    132134 
    133135    # model timing test requires multiple evals to perform the estimate 
    134     sascomp pringle -calc=single,double -timing=100,100 -noplot 
     136    sascomp pringle -engine=single,double -timing=100,100 -noplot 
    135137""" 
    136138 
     
    147149kerneldll.ALLOW_SINGLE_PRECISION_DLLS = True 
    148150 
    149 # list of math functions for use in evaluating parameters 
    150 MATH = dict((k,getattr(math, k)) for k in dir(math) if not k.startswith('_')) 
     151def build_math_context(): 
     152    # type: () -> Dict[str, Callable] 
     153    """build dictionary of functions from math module""" 
     154    return dict((k, getattr(math, k)) 
     155                for k in dir(math) if not k.startswith('_')) 
     156 
     157#: list of math functions for use in evaluating parameters 
     158MATH = build_math_context() 
    151159 
    152160# CRUFT python 2.6 
     
    228236        pass 
    229237 
    230     def __exit__(self, exc_type, exc_value, traceback): 
     238    def __exit__(self, exc_type, exc_value, trace): 
    231239        # type: (Any, BaseException, Any) -> None 
    232         # TODO: better typing for __exit__ method 
    233240        np.random.set_state(self._state) 
    234241 
     
    249256    """ 
    250257    Add a beam stop of the given *radius*.  If *outer*, make an annulus. 
    251  
    252     Note: this function does not require sasview 
    253258    """ 
    254259    if hasattr(data, 'qx_data'): 
     
    278283        # orientation in [-180,180], orientation pd in [0,45] 
    279284        if p.endswith('_pd'): 
    280             return 0., 45. 
     285            return 0., 180. 
    281286        else: 
    282287            return -180., 180. 
     
    371376 
    372377def _random_pd(model_info, pars): 
     378    # type: (ModelInfo, Dict[str, float]) -> None 
     379    """ 
     380    Generate a random dispersity distribution for the model. 
     381 
     382    1% no shape dispersity 
     383    85% single shape parameter 
     384    13% two shape parameters 
     385    1% three shape parameters 
     386 
     387    If oriented, then put dispersity in theta, add phi and psi dispersity 
     388    with 10% probability for each. 
     389    """ 
    373390    pd = [p for p in model_info.parameters.kernel_parameters if p.polydisperse] 
    374391    pd_volume = [] 
     
    433450 
    434451 
     452def limit_dimensions(model_info, pars, maxdim): 
     453    # type: (ModelInfo, ParameterSet, float) -> None 
     454    """ 
     455    Limit parameters of units of Ang to maxdim. 
     456    """ 
     457    for p in model_info.parameters.call_parameters: 
     458        value = pars[p.name] 
     459        if p.units == 'Ang' and value > maxdim: 
     460            pars[p.name] = maxdim*10**np.random.uniform(-3, 0) 
     461 
    435462def constrain_pars(model_info, pars): 
    436463    # type: (ModelInfo, ParameterSet) -> None 
     
    477504        if pars['radius'] < pars['thick_string']: 
    478505            pars['radius'], pars['thick_string'] = pars['thick_string'], pars['radius'] 
    479         pass 
    480506 
    481507    elif name == 'rpa': 
     
    495521    Format the parameter list for printing. 
    496522    """ 
     523    is2d = True 
    497524    lines = [] 
    498525    parameters = model_info.parameters 
     
    594621    return pars 
    595622 
    596 def eval_sasview(model_info, data): 
    597     # type: (Modelinfo, Data) -> Calculator 
    598     """ 
    599     Return a model calculator using the pre-4.0 SasView models. 
    600     """ 
    601     # importing sas here so that the error message will be that sas failed to 
    602     # import rather than the more obscure smear_selection not imported error 
    603     import sas 
    604     import sas.models 
    605     from sas.models.qsmearing import smear_selection 
    606     from sas.models.MultiplicationModel import MultiplicationModel 
    607     from sas.models.dispersion_models import models as dispersers 
    608  
    609     def get_model_class(name): 
    610         # type: (str) -> "sas.models.BaseComponent" 
    611         #print("new",sorted(_pars.items())) 
    612         __import__('sas.models.' + name) 
    613         ModelClass = getattr(getattr(sas.models, name, None), name, None) 
    614         if ModelClass is None: 
    615             raise ValueError("could not find model %r in sas.models"%name) 
    616         return ModelClass 
    617  
    618     # WARNING: ugly hack when handling model! 
    619     # Sasview models with multiplicity need to be created with the target 
    620     # multiplicity, so we cannot create the target model ahead of time for 
    621     # for multiplicity models.  Instead we store the model in a list and 
    622     # update the first element of that list with the new multiplicity model 
    623     # every time we evaluate. 
    624  
    625     # grab the sasview model, or create it if it is a product model 
    626     if model_info.composition: 
    627         composition_type, parts = model_info.composition 
    628         if composition_type == 'product': 
    629             P, S = [get_model_class(revert_name(p))() for p in parts] 
    630             model = [MultiplicationModel(P, S)] 
    631         else: 
    632             raise ValueError("sasview mixture models not supported by compare") 
    633     else: 
    634         old_name = revert_name(model_info) 
    635         if old_name is None: 
    636             raise ValueError("model %r does not exist in old sasview" 
    637                             % model_info.id) 
    638         ModelClass = get_model_class(old_name) 
    639         model = [ModelClass()] 
    640     model[0].disperser_handles = {} 
    641  
    642     # build a smearer with which to call the model, if necessary 
    643     smearer = smear_selection(data, model=model) 
    644     if hasattr(data, 'qx_data'): 
    645         q = np.sqrt(data.qx_data**2 + data.qy_data**2) 
    646         index = ((~data.mask) & (~np.isnan(data.data)) 
    647                  & (q >= data.qmin) & (q <= data.qmax)) 
    648         if smearer is not None: 
    649             smearer.model = model  # because smear_selection has a bug 
    650             smearer.accuracy = data.accuracy 
    651             smearer.set_index(index) 
    652             def _call_smearer(): 
    653                 smearer.model = model[0] 
    654                 return smearer.get_value() 
    655             theory = _call_smearer 
    656         else: 
    657             theory = lambda: model[0].evalDistribution([data.qx_data[index], 
    658                                                         data.qy_data[index]]) 
    659     elif smearer is not None: 
    660         theory = lambda: smearer(model[0].evalDistribution(data.x)) 
    661     else: 
    662         theory = lambda: model[0].evalDistribution(data.x) 
    663  
    664     def calculator(**pars): 
    665         # type: (float, ...) -> np.ndarray 
    666         """ 
    667         Sasview calculator for model. 
    668         """ 
    669         oldpars = revert_pars(model_info, pars) 
    670         # For multiplicity models, create a model with the correct multiplicity 
    671         control = oldpars.pop("CONTROL", None) 
    672         if control is not None: 
    673             # sphericalSLD has one fewer multiplicity.  This update should 
    674             # happen in revert_pars, but it hasn't been called yet. 
    675             model[0] = ModelClass(control) 
    676         # paying for parameter conversion each time to keep life simple, if not fast 
    677         for k, v in oldpars.items(): 
    678             if k.endswith('.type'): 
    679                 par = k[:-5] 
    680                 if v == 'gaussian': continue 
    681                 cls = dispersers[v if v != 'rectangle' else 'rectangula'] 
    682                 handle = cls() 
    683                 model[0].disperser_handles[par] = handle 
    684                 try: 
    685                     model[0].set_dispersion(par, handle) 
    686                 except Exception: 
    687                     exception.annotate_exception("while setting %s to %r" 
    688                                                  %(par, v)) 
    689                     raise 
    690  
    691  
    692         #print("sasview pars",oldpars) 
    693         for k, v in oldpars.items(): 
    694             name_attr = k.split('.')  # polydispersity components 
    695             if len(name_attr) == 2: 
    696                 par, disp_par = name_attr 
    697                 model[0].dispersion[par][disp_par] = v 
    698             else: 
    699                 model[0].setParam(k, v) 
    700         return theory() 
    701  
    702     calculator.engine = "sasview" 
    703     return calculator 
    704623 
    705624DTYPE_MAP = { 
     
    795714    than OpenCL. 
    796715    """ 
    797     if dtype == 'sasview': 
    798         return eval_sasview(model_info, data) 
    799     elif dtype is None or not dtype.endswith('!'): 
     716    if dtype is None or not dtype.endswith('!'): 
    800717        return eval_opencl(model_info, data, dtype=dtype, cutoff=cutoff) 
    801718    else: 
     
    815732 
    816733 
    817 def compare(opts, limits=None): 
     734def compare(opts, limits=None, maxdim=np.inf): 
    818735    # type: (Dict[str, Any], Optional[Tuple[float, float]]) -> Tuple[float, float] 
    819736    """ 
     
    826743    as the values are adjusted, making it easier to see the effects of the 
    827744    parameters. 
    828     """ 
    829     limits = np.Inf, -np.Inf 
     745 
     746    *maxdim* is the maximum value for any parameter with units of Angstrom. 
     747    """ 
    830748    for k in range(opts['sets']): 
    831         opts['pars'] = parse_pars(opts) 
     749        if k > 1: 
     750            # print a separate seed for each dataset for better reproducibility 
     751            new_seed = np.random.randint(1000000) 
     752            print("Set %d uses -random=%i"%(k+1, new_seed)) 
     753            np.random.seed(new_seed) 
     754        opts['pars'] = parse_pars(opts, maxdim=maxdim) 
    832755        if opts['pars'] is None: 
    833756            return 
     
    835758        if opts['plot']: 
    836759            limits = plot_models(opts, result, limits=limits, setnum=k) 
     760        if opts['show_weights']: 
     761            base, _ = opts['engines'] 
     762            base_pars, _ = opts['pars'] 
     763            model_info = base._kernel.info 
     764            dim = base._kernel.dim 
     765            plot_weights(model_info, get_mesh(model_info, base_pars, dim=dim)) 
    837766    if opts['plot']: 
    838767        import matplotlib.pyplot as plt 
    839768        plt.show() 
     769    return limits 
    840770 
    841771def run_models(opts, verbose=False): 
    842772    # type: (Dict[str, Any]) -> Dict[str, Any] 
     773    """ 
     774    Process a parameter set, return calculation results and times. 
     775    """ 
    843776 
    844777    base, comp = opts['engines'] 
     
    896829    # work with trimmed data, not the full set 
    897830    sorted_err = np.sort(abs(err.compressed())) 
    898     if len(sorted_err) == 0.: 
     831    if len(sorted_err) == 0: 
    899832        print(label + "  no valid values") 
    900833        return 
     
    912845 
    913846 
    914 def plot_models(opts, result, limits=(np.Inf, -np.Inf), setnum=0): 
     847def plot_models(opts, result, limits=None, setnum=0): 
    915848    # type: (Dict[str, Any], Dict[str, Any], Optional[Tuple[float, float]]) -> Tuple[float, float] 
     849    """ 
     850    Plot the results from :func:`run_model`. 
     851    """ 
     852    import matplotlib.pyplot as plt 
     853 
    916854    base_value, comp_value = result['base_value'], result['comp_value'] 
    917855    base_time, comp_time = result['base_time'], result['comp_time'] 
     
    925863    # Plot if requested 
    926864    view = opts['view'] 
    927     import matplotlib.pyplot as plt 
    928     vmin, vmax = limits 
    929     if have_base: 
    930         vmin = min(vmin, base_value.min()) 
    931         vmax = max(vmax, base_value.max()) 
    932     if have_comp: 
    933         vmin = min(vmin, comp_value.min()) 
    934         vmax = max(vmax, comp_value.max()) 
    935     limits = vmin, vmax 
     865    if limits is None: 
     866        vmin, vmax = np.inf, -np.inf 
     867        if have_base: 
     868            vmin = min(vmin, base_value.min()) 
     869            vmax = max(vmax, base_value.max()) 
     870        if have_comp: 
     871            vmin = min(vmin, comp_value.min()) 
     872            vmax = max(vmax, comp_value.max()) 
     873        limits = vmin, vmax 
    936874 
    937875    if have_base: 
     
    958896                errview = 'linear' 
    959897        if 0:  # 95% cutoff 
    960             sorted = np.sort(err.flatten()) 
    961             cutoff = sorted[int(sorted.size*0.95)] 
     898            sorted_err = np.sort(err.flatten()) 
     899            cutoff = sorted_err[int(sorted_err.size*0.95)] 
    962900            err[err > cutoff] = cutoff 
    963901        #err,errstr = base/comp,"ratio" 
    964         plot_theory(data, None, resid=err, view=view, use_data=use_data) 
    965         plt.yscale(errview) 
     902        plot_theory(data, None, resid=err, view=errview, use_data=use_data) 
     903        plt.xscale('log' if view == 'log' and not opts['is2d'] else 'linear') 
     904        plt.legend(['P%d'%(k+1) for k in range(setnum+1)], loc='best') 
    966905        plt.title("max %s = %.3g"%(errstr, abs(err).max())) 
    967906        #cbar_title = errstr if errview=="linear" else "log "+errstr 
     
    995934OPTIONS = [ 
    996935    # Plotting 
    997     'plot', 'noplot', 
     936    'plot', 'noplot', 'weights', 
    998937    'linear', 'log', 'q4', 
    999938    'rel', 'abs', 
     
    1010949    'demo', 'default',  # TODO: remove demo/default 
    1011950    'nopars', 'pars', 
     951    'sphere', 'sphere=', # integrate over a sphere in 2d with n points 
    1012952 
    1013953    # Calculation options 
     
    1018958 
    1019959    # Precision options 
    1020     'calc=', 
     960    'engine=', 
    1021961    'half', 'fast', 'single', 'double', 'single!', 'double!', 'quad!', 
    1022     'sasview',  # TODO: remove sasview 3.x support 
    1023962 
    1024963    # Output options 
     
    1026965    ] 
    1027966 
    1028 NAME_OPTIONS = set(k for k in OPTIONS if not k.endswith('=')) 
    1029 VALUE_OPTIONS = [k[:-1] for k in OPTIONS if k.endswith('=')] 
     967NAME_OPTIONS = (lambda: set(k for k in OPTIONS if not k.endswith('=')))() 
     968VALUE_OPTIONS = (lambda: [k[:-1] for k in OPTIONS if k.endswith('=')])() 
    1030969 
    1031970 
     
    10751014 
    10761015INTEGER_RE = re.compile("^[+-]?[1-9][0-9]*$") 
    1077 def isnumber(str): 
    1078     match = FLOAT_RE.match(str) 
    1079     isfloat = (match and not str[match.end():]) 
    1080     return isfloat or INTEGER_RE.match(str) 
     1016def isnumber(s): 
     1017    # type: (str) -> bool 
     1018    """Return True if string contains an int or float""" 
     1019    match = FLOAT_RE.match(s) 
     1020    isfloat = (match and not s[match.end():]) 
     1021    return isfloat or INTEGER_RE.match(s) 
    10811022 
    10821023# For distinguishing pairs of models for comparison 
     
    11171058    name = positional_args[-1] 
    11181059 
    1119     # pylint: disable=bad-whitespace 
     1060    # pylint: disable=bad-whitespace,C0321 
    11201061    # Interpret the flags 
    11211062    opts = { 
     
    11451086        'sets'      : 0, 
    11461087        'engine'    : 'default', 
    1147         'evals'     : '1', 
     1088        'count'     : '1', 
     1089        'show_weights' : False, 
     1090        'sphere'    : 0, 
    11481091    } 
    11491092    for arg in flags: 
     
    11681111        elif arg.startswith('-accuracy='): opts['accuracy'] = arg[10:] 
    11691112        elif arg.startswith('-cutoff='):   opts['cutoff'] = arg[8:] 
    1170         elif arg.startswith('-random='):   opts['seed'] = int(arg[8:]) 
    11711113        elif arg.startswith('-title='):    opts['title'] = arg[7:] 
    11721114        elif arg.startswith('-data='):     opts['datafile'] = arg[6:] 
    1173         elif arg.startswith('-calc='):     opts['engine'] = arg[6:] 
    1174         elif arg.startswith('-neval='):    opts['evals'] = arg[7:] 
    1175         elif arg == '-random':  opts['seed'] = np.random.randint(1000000) 
     1115        elif arg.startswith('-engine='):   opts['engine'] = arg[8:] 
     1116        elif arg.startswith('-neval='):    opts['count'] = arg[7:] 
     1117        elif arg.startswith('-random='): 
     1118            opts['seed'] = int(arg[8:]) 
     1119            opts['sets'] = 0 
     1120        elif arg == '-random': 
     1121            opts['seed'] = np.random.randint(1000000) 
     1122            opts['sets'] = 0 
     1123        elif arg.startswith('-sphere'): 
     1124            opts['sphere'] = int(arg[8:]) if len(arg) > 7 else 150 
     1125            opts['is2d'] = True 
    11761126        elif arg == '-preset':  opts['seed'] = -1 
    11771127        elif arg == '-mono':    opts['mono'] = True 
     
    11921142        elif arg == '-double!': opts['engine'] = 'double!' 
    11931143        elif arg == '-quad!':   opts['engine'] = 'quad!' 
    1194         elif arg == '-sasview': opts['engine'] = 'sasview' 
    11951144        elif arg == '-edit':    opts['explore'] = True 
    11961145        elif arg == '-demo':    opts['use_demo'] = True 
    11971146        elif arg == '-default': opts['use_demo'] = False 
     1147        elif arg == '-weights': opts['show_weights'] = True 
    11981148        elif arg == '-html':    opts['html'] = True 
    11991149        elif arg == '-help':    opts['html'] = True 
    1200     # pylint: enable=bad-whitespace 
     1150    # pylint: enable=bad-whitespace,C0321 
    12011151 
    12021152    # Magnetism forces 2D for now 
     
    12321182 
    12331183    if PAR_SPLIT in opts['engine']: 
    1234         engine_types = opts['engine'].split(PAR_SPLIT, 2) 
     1184        opts['engine'] = opts['engine'].split(PAR_SPLIT, 2) 
    12351185        comparison = True 
    12361186    else: 
    1237         engine_types = [opts['engine']]*2 
    1238  
    1239     if PAR_SPLIT in opts['evals']: 
    1240         evals = [int(k) for k in opts['evals'].split(PAR_SPLIT, 2)] 
     1187        opts['engine'] = [opts['engine']]*2 
     1188 
     1189    if PAR_SPLIT in opts['count']: 
     1190        opts['count'] = [int(k) for k in opts['count'].split(PAR_SPLIT, 2)] 
    12411191        comparison = True 
    12421192    else: 
    1243         evals = [int(opts['evals'])]*2 
     1193        opts['count'] = [int(opts['count'])]*2 
    12441194 
    12451195    if PAR_SPLIT in opts['cutoff']: 
    1246         cutoff = [float(k) for k in opts['cutoff'].split(PAR_SPLIT, 2)] 
     1196        opts['cutoff'] = [float(k) for k in opts['cutoff'].split(PAR_SPLIT, 2)] 
    12471197        comparison = True 
    12481198    else: 
    1249         cutoff = [float(opts['cutoff'])]*2 
    1250  
    1251     base = make_engine(model_info[0], data, engine_types[0], cutoff[0]) 
     1199        opts['cutoff'] = [float(opts['cutoff'])]*2 
     1200 
     1201    base = make_engine(model_info[0], data, opts['engine'][0], opts['cutoff'][0]) 
    12521202    if comparison: 
    1253         comp = make_engine(model_info[1], data, engine_types[1], cutoff[1]) 
     1203        comp = make_engine(model_info[1], data, opts['engine'][1], opts['cutoff'][1]) 
    12541204    else: 
    12551205        comp = None 
     
    12601210        'data'      : data, 
    12611211        'name'      : names, 
    1262         'def'       : model_info, 
    1263         'count'     : evals, 
     1212        'info'      : model_info, 
    12641213        'engines'   : [base, comp], 
    12651214        'values'    : values, 
     
    12671216    # pylint: enable=bad-whitespace 
    12681217 
     1218    # Set the integration parameters to the half sphere 
     1219    if opts['sphere'] > 0: 
     1220        set_spherical_integration_parameters(opts, opts['sphere']) 
     1221 
    12691222    return opts 
    12701223 
    1271 def parse_pars(opts): 
    1272     model_info, model_info2 = opts['def'] 
     1224def set_spherical_integration_parameters(opts, steps): 
     1225    # type: (Dict[str, Any], int) -> None 
     1226    """ 
     1227    Set integration parameters for spherical integration over the entire 
     1228    surface in theta-phi coordinates. 
     1229    """ 
     1230    # Set the integration parameters to the half sphere 
     1231    opts['values'].extend([ 
     1232        #'theta=90', 
     1233        'theta_pd=%g'%(90/np.sqrt(3)), 
     1234        'theta_pd_n=%d'%steps, 
     1235        'theta_pd_type=rectangle', 
     1236        #'phi=0', 
     1237        'phi_pd=%g'%(180/np.sqrt(3)), 
     1238        'phi_pd_n=%d'%(2*steps), 
     1239        'phi_pd_type=rectangle', 
     1240        #'background=0', 
     1241    ]) 
     1242    if 'psi' in opts['info'][0].parameters: 
     1243        opts['values'].extend([ 
     1244            #'psi=0', 
     1245            'psi_pd=%g'%(180/np.sqrt(3)), 
     1246            'psi_pd_n=%d'%(2*steps), 
     1247            'psi_pd_type=rectangle', 
     1248        ]) 
     1249 
     1250def parse_pars(opts, maxdim=np.inf): 
     1251    # type: (Dict[str, Any], float) -> Tuple[Dict[str, float], Dict[str, float]] 
     1252    """ 
     1253    Generate a parameter set. 
     1254 
     1255    The default values come from the model, or a randomized model if a seed 
     1256    value is given.  Next, evaluate any parameter expressions, constraining 
     1257    the value of the parameter within and between models.  If *maxdim* is 
     1258    given, limit parameters with units of Angstrom to this value. 
     1259 
     1260    Returns a pair of parameter dictionaries for base and comparison models. 
     1261    """ 
     1262    model_info, model_info2 = opts['info'] 
    12731263 
    12741264    # Get demo parameters from model definition, or use default parameters 
     
    12811271    if opts['seed'] > -1: 
    12821272        pars = randomize_pars(model_info, pars) 
     1273        limit_dimensions(model_info, pars, maxdim) 
    12831274        if model_info != model_info2: 
    12841275            pars2 = randomize_pars(model_info2, pars2) 
     1276            limit_dimensions(model_info, pars2, maxdim) 
    12851277            # Share values for parameters with the same name 
    12861278            for k, v in pars.items(): 
     
    13061298            print("%r invalid; parameters are: %s"%(k, ", ".join(sorted(s)))) 
    13071299            return None 
    1308         v1, v2 = v.split(PAR_SPLIT, 2) if PAR_SPLIT in v else (v,v) 
     1300        v1, v2 = v.split(PAR_SPLIT, 2) if PAR_SPLIT in v else (v, v) 
    13091301        if v1 and k in pars: 
    13101302            presets[k] = float(v1) if isnumber(v1) else v1 
     
    13551347    show html docs for the model 
    13561348    """ 
    1357     import os 
    13581349    from .generate import make_html 
    13591350    from . import rst2html 
    13601351 
    1361     info = opts['def'][0] 
     1352    info = opts['info'][0] 
    13621353    html = make_html(info) 
    13631354    path = os.path.dirname(info.filename) 
    1364     url = "file://"+path.replace("\\","/")[2:]+"/" 
     1355    url = "file://" + path.replace("\\", "/")[2:] + "/" 
    13651356    rst2html.view_html_qtapp(html, url) 
    13661357 
     
    13861377    frame.panel.Layout() 
    13871378    frame.panel.aui.Split(0, wx.TOP) 
    1388     def reset_parameters(event): 
     1379    def _reset_parameters(event): 
    13891380        model.revert_values() 
    13901381        signal.update_parameters(problem) 
    1391     frame.Bind(wx.EVT_TOOL, reset_parameters, frame.ToolBar.GetToolByPos(1)) 
    1392     if is_mac: frame.Show() 
     1382    frame.Bind(wx.EVT_TOOL, _reset_parameters, frame.ToolBar.GetToolByPos(1)) 
     1383    if is_mac: 
     1384        frame.Show() 
    13931385    # If running withing an app, start the main loop 
    13941386    if app: 
     
    14101402        opts['pars'] = list(opts['pars']) 
    14111403        p1, p2 = opts['pars'] 
    1412         m1, m2 = opts['def'] 
     1404        m1, m2 = opts['info'] 
    14131405        self.fix_p2 = m1 != m2 or p1 != p2 
    14141406        model_info = m1 
     
    14291421        self.starting_values = dict((k, v.value) for k, v in pars.items()) 
    14301422        self.pd_types = pd_types 
    1431         self.limits = np.Inf, -np.Inf 
     1423        self.limits = None 
    14321424 
    14331425    def revert_values(self): 
     1426        # type: () -> None 
     1427        """ 
     1428        Restore starting values of the parameters. 
     1429        """ 
    14341430        for k, v in self.starting_values.items(): 
    14351431            self.pars[k].value = v 
    14361432 
    14371433    def model_update(self): 
     1434        # type: () -> None 
     1435        """ 
     1436        Respond to signal that model parameters have been changed. 
     1437        """ 
    14381438        pass 
    14391439 
  • sasmodels/compare_many.py

    rf72d70a r2d81cfe  
    2222from .compare import (randomize_pars, suppress_pd, make_data, 
    2323                      make_engine, get_pars, columnize, 
    24                       constrain_pars, constrain_new_to_old) 
     24                      constrain_pars) 
    2525 
    2626MODELS = core.list_models() 
     
    8080    'double!': 5e-14, 
    8181    'quad!': 5e-18, 
    82     'sasview': 5e-14, 
    8382} 
    8483def compare_instance(name, data, index, N=1, mono=True, cutoff=1e-5, 
    85                      base='sasview', comp='double'): 
     84                     base='single', comp='double'): 
    8685    r""" 
    8786    Compare the model under different calculation engines. 
     
    164163        print("Model %s %d"%(name, k+1), file=sys.stderr) 
    165164        seed = np.random.randint(1e6) 
    166         pars_i = randomize_pars(model_info, pars, seed) 
     165        np.random.seed(seed) 
     166        pars_i = randomize_pars(model_info, pars) 
    167167        constrain_pars(model_info, pars_i) 
    168         if 'sasview' in (base, comp): 
    169             constrain_new_to_old(model_info, pars_i) 
    170168        if mono: 
    171169            pars_i = suppress_pd(pars_i) 
     
    266264        return 
    267265 
    268     data, index = make_data({'qmax':1.0, 'is2d':is2D, 'nq':Nq, 'res':0., 
    269                              'accuracy': 'Low', 'view':'log', 'zero': False}) 
     266    data, index = make_data({ 
     267        'qmin': 0.001, 'qmax': 1.0, 'is2d': is2D, 'nq': Nq, 'res': 0., 
     268        'accuracy': 'Low', 'view':'log', 'zero': False 
     269        }) 
    270270    for model in model_list: 
    271271        compare_instance(model, data, index, N=count, mono=mono, 
  • sasmodels/conversion_table.py

    r505d0ad re65c3ba  
    2929 
    3030CONVERSION_TABLE = { 
    31     (3,1,2) : { 
    32     "adsorbed_layer": [ 
    33         "Core2ndMomentModel", 
    34         { 
    35             "scale": "scale", 
    36             "second_moment": "second_moment", 
    37             "density_shell": "density_poly", 
    38             "sld_solvent": "sld_solv", 
    39             "radius": "radius_core", 
    40             "volfraction": "volf_cores", 
    41             "background": "background", 
    42             "adsorbed_amount": "ads_amount", 
    43             "sld_shell": "sld_poly" 
    44         } 
    45     ], 
    46     "barbell": [ 
    47         "BarBellModel", 
    48         { 
    49             "sld": "sld_barbell", 
    50             "length": "len_bar", 
    51             "radius_bell": "rad_bell", 
    52             "radius": "rad_bar", 
    53             "sld_solvent": "sld_solv" 
    54         } 
    55     ], 
    56     "bcc_paracrystal": [ 
    57         "BCCrystalModel", 
    58         { 
    59             "sld": "sldSph", 
    60             "sld_solvent": "sldSolv" 
    61         } 
    62     ], 
    63     "be_polyelectrolyte": [ 
    64         "BEPolyelectrolyte", 
    65         { 
    66             "ionization_degree": "alpha", 
    67             "polymer_concentration": "c", 
    68             "salt_concentration": "cs", 
    69             "virial_param": "h", 
    70             "background": "background", 
    71             "contrast_factor": "k", 
    72             "bjerrum_length": "lb", 
    73             "monomer_length": "b" 
    74         } 
    75     ], 
    76     "binary_hard_sphere": [ 
    77         "BinaryHSModel", 
    78         { 
    79             "sld_sm": "ss_sld", 
    80             "sld_lg": "ls_sld", 
    81             "volfraction_sm": "vol_frac_ss", 
    82             "radius_lg": "l_radius", 
    83             "radius_sm": "s_radius", 
    84             "volfraction_lg": "vol_frac_ls", 
    85             "sld_solvent": "solvent_sld" 
    86         } 
    87     ], 
    88     "broad_peak": [ 
    89         "BroadPeakModel", 
    90         { 
    91             "peak_pos": "q_peak", 
    92             "scale": None, 
    93             "lorentz_length": "length_l", 
    94             "porod_scale": "scale_p", 
    95             "lorentz_exp": "exponent_l", 
    96             "lorentz_scale": "scale_l", 
    97             "porod_exp": "exponent_p" 
    98         } 
    99     ], 
    100     "capped_cylinder": [ 
    101         "CappedCylinderModel", 
    102         { 
    103             "sld": "sld_capcyl", 
    104             "length": "len_cyl", 
    105             "radius_cap": "rad_cap", 
    106             "radius": "rad_cyl", 
    107             "sld_solvent": "sld_solv" 
    108         } 
    109     ], 
    110     "core_multi_shell": [ 
    111         "CoreMultiShellModel", 
    112         { 
    113             "thickness": "thick_shell", 
    114             "sld": "sld_shell", 
    115             "radius": "rad_core0", 
    116             "sld_core": "sld_core0", 
    117             "sld_solvent": "sld_solv", 
    118             "n": "n_shells", 
    119             "M0:sld_core": "M0_sld_core0", 
    120             "mtheta:sld_core": "M_theta_core0", 
    121             "mphi:sld_core": "M_phi_core0", 
    122             "M0:sld1": "M0_sld_shell1", 
    123             "mtheta:sld1": "M_theta_shell1", 
    124             "mphi:sld1": "M_phi_shell1", 
    125             "M0:sld2": "M0_sld_shell2", 
    126             "mtheta:sld2": "M_theta_shell2", 
    127             "mphi:sld2": "M_phi_shell2", 
    128             "M0:sld3": "M0_sld_shell3", 
    129             "mtheta:sld3": "M_theta_shell3", 
    130             "mphi:sld3": "M_phi_shell3", 
    131             "M0:sld4": "M0_sld_shell4", 
    132             "mtheta:sld4": "M_theta_shell4", 
    133             "mphi:sld4": "M_phi_shell4", 
    134             "M0:sld_solvent": "M0_sld_solv", 
    135             "mtheta:sld_solvent": "M_theta_solv", 
    136             "mphi:sld_solvent": "M_phi_solv", 
    137             "up:frac_i": "Up_frac_i", 
    138             "up:frac_f": "Up_frac_f", 
    139             "up:angle": "Up_theta", 
    140         } 
    141     ], 
    142     "core_shell_bicelle": [ 
    143         "CoreShellBicelleModel", 
    144         { 
    145             "phi": "axis_phi", 
    146             "sld_core": "core_sld", 
    147             "sld_rim": "rim_sld", 
    148             "thick_face": "face_thick", 
    149             "sld_solvent": "solvent_sld", 
    150             "thick_rim": "rim_thick", 
    151             "sld_face": "face_sld", 
    152             "theta": "axis_theta" 
    153         } 
    154     ], 
    155     "core_shell_cylinder": [ 
    156         "CoreShellCylinderModel", 
    157         { 
    158             "theta": "axis_theta", 
    159             "phi": "axis_phi", 
    160             "sld_shell": "shell_sld", 
    161             "sld_solvent": "solvent_sld", 
    162             "sld_core": "core_sld" 
    163         } 
    164     ], 
    165     "core_shell_ellipsoid:1": [ 
    166         "CoreShellEllipsoidModel", 
    167         { 
    168             "sld_core": "sld_core", 
    169             "sld_shell": "sld_shell", 
    170             "sld_solvent": "sld_solvent", 
    171             "radius_equat_core": "equat_core", 
    172             "x_core": "polar_core", 
    173             "thick_shell": "equat_shell", 
    174             "x_polar_shell": "polar_shell", 
    175             "theta": "axis_theta", 
    176             "phi": "axis_phi", 
    177         } 
    178     ], 
    179     "core_shell_ellipsoid": [ 
    180         "CoreShellEllipsoidXTModel", 
    181         { 
    182             "sld_core": "sld_core", 
    183             "sld_shell": "sld_shell", 
    184             "sld_solvent": "sld_solvent", 
    185             "radius_equat_core": "equat_core", 
    186             "thick_shell": "T_shell", 
    187             "x_core": "X_core", 
    188             "x_polar_shell": "XpolarShell", 
    189             "theta": "axis_theta", 
    190             "phi": "axis_phi", 
    191         } 
    192     ], 
    193     "core_shell_parallelepiped": [ 
    194         "CSParallelepipedModel", 
    195         { 
    196             "sld_core": "sld_pcore", 
    197             "sld_a": "sld_rimA", 
    198             "sld_b": "sld_rimB", 
    199             "sld_c": "sld_rimC", 
    200             "sld_solvent": "sld_solv", 
    201             "length_a": "shortA", 
    202             "length_b": "midB", 
    203             "length_c": "longC", 
    204             "thick_rim_a": "rimA", 
    205             "thick_rim_c": "rimC", 
    206             "thick_rim_b": "rimB", 
    207             "theta": "parallel_theta", 
    208             "phi": "parallel_phi", 
    209             "psi": "parallel_psi", 
    210         } 
    211     ], 
    212     "core_shell_sphere": [ 
    213         "CoreShellModel", 
    214         { 
    215             "sld_core": "core_sld", 
    216             "sld_shell": "shell_sld", 
    217             "sld_solvent": "solvent_sld", 
    218             "M0:sld_core": "M0_sld_core", 
    219             "mtheta:sld_core": "M_theta_core", 
    220             "mphi:sld_core": "M_phi_core", 
    221             "M0:sld_shell": "M0_sld_shell", 
    222             "mtheta:sld_shell": "M_theta_shell", 
    223             "mphi:sld_shell": "M_phi_shell", 
    224             "M0:sld_solvent": "M0_sld_solv", 
    225             "mtheta:sld_solvent": "M_theta_solv", 
    226             "mphi:sld_solvent": "M_phi_solv", 
    227             "up:frac_i": "Up_frac_i", 
    228             "up:frac_f": "Up_frac_f", 
    229             "up:angle": "Up_theta" 
    230         } 
    231     ], 
    232     "correlation_length": [ 
    233         "CorrLength", 
    234         { 
    235             "porod_scale": "scale_p", 
    236             "lorentz_scale": "scale_l", 
    237             "porod_exp": "exponent_p", 
    238             "lorentz_exp": "exponent_l", 
    239             "cor_length": "length_l" 
    240         }, 
    241         "CorrLengthModel" 
    242     ], 
    243     "cylinder": [ 
    244         "CylinderModel", 
    245         { 
    246             "sld": "sldCyl", 
    247             "theta": "cyl_theta", 
    248             "phi": "cyl_phi", 
    249             "sld_solvent": "sldSolv", 
    250             "M0:sld": "M0_sld_cyl", 
    251             "mtheta:sld": "M_theta_cyl", 
    252             "mphi:sld": "M_phi_cyl", 
    253             "M0:sld_solvent": "M0_sld_solv", 
    254             "mtheta:sld_solvent": "M_theta_solv", 
    255             "mphi:sld_solvent": "M_phi_solv", 
    256             "up:frac_i": "Up_frac_i", 
    257             "up:frac_f": "Up_frac_f", 
    258             "up:angle": "Up_theta" 
    259         } 
    260     ], 
    261     "dab": [ 
    262         "DABModel", 
    263         { 
    264             "cor_length": "length" 
    265         } 
    266     ], 
    267     "ellipsoid": [ 
    268         "EllipsoidModel", 
    269         { 
    270             "phi": "axis_phi", 
    271             "radius_equatorial": "radius_b", 
    272             "sld": "sldEll", 
    273             "theta": "axis_theta", 
    274             "radius_polar": "radius_a", 
    275             "sld_solvent": "sldSolv" 
    276         } 
    277     ], 
    278     "elliptical_cylinder": [ 
    279         "EllipticalCylinderModel", 
    280         { 
    281             "axis_ratio": "r_ratio", 
    282             "radius_minor": "r_minor", 
    283             "sld": "sldCyl", 
    284             "sld_solvent": "sldSolv", 
    285             "theta": "cyl_theta", 
    286             "phi": "cyl_phi", 
    287             "psi": "cyl_psi", 
    288         } 
    289     ], 
    290     "fcc_paracrystal": [ 
    291         "FCCrystalModel", 
    292         { 
    293             "sld": "sldSph", 
    294             "sld_solvent": "sldSolv" 
    295         } 
    296     ], 
    297     "flexible_cylinder": [ 
    298         "FlexibleCylinderModel", 
    299         { 
    300             "sld": "sldCyl", 
    301             "sld_solvent": "sldSolv" 
    302         } 
    303     ], 
    304     "flexible_cylinder_elliptical": [ 
    305         "FlexCylEllipXModel", 
    306         { 
    307             "sld": "sldCyl", 
    308             "sld_solvent": "sldSolv" 
    309         } 
    310     ], 
    311     "fractal": [ 
    312         "FractalModel", 
    313         { 
    314             "sld_block": "sldBlock", 
    315             "radius": "radius", 
    316             "cor_length": "cor_length", 
    317             "sld_solvent": "sldSolv", 
    318             "fractal_dim": "fractal_dim" 
    319         } 
    320     ], 
    321     "fractal_core_shell": [ 
    322         "FractalCoreShell", 
    323         { 
    324             "sld_core": "core_sld", 
    325             "sld_shell": "shell_sld", 
    326             "sld_solvent": "solvent_sld", 
    327             "radius": "radius", 
    328             "thickness": "thickness", 
    329             "fractal_dim": "frac_dim", 
    330             "cor_length": "cor_length", 
    331             "volfraction": "volfraction", 
    332         }, 
    333         "FractalCoreShellModel" 
    334     ], 
    335     "fuzzy_sphere": [ 
    336         "FuzzySphereModel", 
    337         { 
    338             "sld": "sldSph", 
    339             "fuzziness": "fuzziness", 
    340             "radius": "radius", 
    341             "sld_solvent": "sldSolv" 
    342         } 
    343     ], 
    344     "gauss_lorentz_gel": [ 
    345         "GaussLorentzGel", 
    346         { 
    347             "gauss_scale": "scale_g", 
    348             "cor_length_dynamic": "dyn_colength", 
    349             "cor_length_static": "stat_colength", 
    350             "background": "background", 
    351             "lorentz_scale": "scale_l" 
    352         }, 
    353         "GaussLorentzGelModel" 
    354     ], 
    355     "gaussian_peak": [ 
    356         "Peak Gauss Model", 
    357         { 
    358             "peak_pos": "q0", 
    359             "sigma": "B", 
    360         }, 
    361         "PeakGaussModel", 
    362     ], 
    363     "gel_fit": [ 
    364         "GelFitModel", 
    365         { 
    366             "rg": "radius", 
    367             "lorentz_scale": "lScale", 
    368             "guinier_scale": "gScale", 
    369             "fractal_dim": "FractalExp", 
    370             "cor_length": "zeta", 
    371         } 
    372     ], 
    373     "guinier": [ 
    374         "Guinier", 
    375         { 
    376             "rg": "rg" 
    377         }, 
    378         "GuinierModel", 
    379     ], 
    380     "guinier_porod": [ 
    381         "GuinierPorod", 
    382         { 
    383             "s": "dim", 
    384             "rg": "rg", 
    385             "porod_exp": "m", 
    386             "scale": "scale", 
    387             "background": "background" 
    388         }, 
    389         "GuinierPorodModel", 
    390     ], 
    391     "hardsphere": [ 
    392         "HardsphereStructure", 
    393         { 
    394             "scale": "scale_factor", 
    395             "radius_effective": "effect_radius", 
    396         } 
    397     ], 
    398     "hayter_msa": [ 
    399         "HayterMSAStructure", 
    400         { 
    401             "scale": "scale_factor", 
    402             "radius_effective": "effect_radius", 
    403             "volfraction": "volfraction", 
    404             "charge": "charge", 
    405             "temperature": "temperature", 
    406             "concentration_salt": "saltconc", 
    407             "dielectconst": "dielectconst", 
    408         } 
    409     ], 
    410     "hollow_cylinder": [ 
    411         "HollowCylinderModel", 
    412         { 
    413             "sld": "sldCyl", 
    414             "sld_solvent": "sldSolv", 
    415             "radius": "core_radius", 
    416             "thickness": "radius", 
    417             "length": "length", 
    418             "theta": "axis_theta", 
    419             "phi": "axis_phi", 
    420         } 
    421     ], 
    422     "hollow_rectangular_prism": [ 
    423         "RectangularHollowPrismModel", 
    424         { 
    425             "sld": "sldPipe", 
    426             "sld_solvent": "sldSolv", 
    427             "length_a": "short_side", 
    428             "b2a_ratio": "b2a_ratio", 
    429             "c2a_ratio": "c2a_ratio", 
    430             "thickness": "thickness", 
    431         } 
    432     ], 
    433     "hollow_rectangular_prism_thin_walls": [ 
    434         "RectangularHollowPrismInfThinWallsModel", 
    435         { 
    436             "sld": "sldPipe", 
    437             "sld_solvent": "sldSolv", 
    438             "length_a": "short_side", 
    439             "b2a_ratio": "b2a_ratio", 
    440             "c2a_ratio": "c2a_ratio", 
    441         } 
    442     ], 
    443     "lamellar": [ 
    444         "LamellarModel", 
    445         { 
    446             "sld": "sld_bi", 
    447             "sld_solvent": "sld_sol", 
    448             "thickness": "bi_thick" 
    449         } 
    450     ], 
    451     "lamellar_hg": [ 
    452         "LamellarFFHGModel", 
    453         { 
    454             "sld": "sld_tail", 
    455             "sld_solvent": "sld_solvent", 
    456             "sld_head": "sld_head", 
    457             "length_tail": "t_length", 
    458             "length_head": "h_thickness" 
    459         } 
    460     ], 
    461     "lamellar_hg_stack_caille": [ 
    462         "LamellarPSHGModel", 
    463         { 
    464             "sld": "sld_tail", 
    465             "sld_head": "sld_head", 
    466             "sld_solvent": "sld_solvent", 
    467             "length_tail": "deltaT", 
    468             "length_head": "deltaH", 
    469             "d_spacing": "spacing", 
    470             "Caille_parameter": "caille", 
    471             "Nlayers": "n_plates", 
    472         } 
    473     ], 
    474     "lamellar_stack_caille": [ 
    475         "LamellarPSModel", 
    476         { 
    477             "sld": "sld_bi", 
    478             "sld_solvent": "sld_sol", 
    479             "thickness": "delta", 
    480             "d_spacing": "spacing", 
    481             "Caille_parameter": "caille", 
    482             "Nlayers": "n_plates", 
    483         } 
    484     ], 
    485     "lamellar_stack_paracrystal": [ 
    486         "LamellarPCrystalModel", 
    487         { 
    488             "sld": "sld_layer", 
    489             "sld_solvent": "sld_solvent", 
    490             "thickness": "thickness", 
    491             "d_spacing": "spacing", 
    492             "sigma_d": "pd_spacing", 
    493             "Nlayers": "Nlayers", 
    494         } 
    495     ], 
    496     "line": [ 
    497         "LineModel", 
    498         { 
    499             "slope": "B", 
    500             "scale": None, 
    501             "background": None, 
    502             "intercept": "A" 
    503         } 
    504     ], 
    505     "linear_pearls": [ 
    506         "LinearPearlsModel", 
    507         { 
    508             "sld": "sld_pearl", 
    509             "sld_solvent": "sld_solv", 
    510             "edge_sep": "edge_separation" 
    511         } 
    512     ], 
    513     "lorentz": [ 
    514         "Lorentz", 
    515         { 
    516             "cor_length": "length" 
    517         }, 
    518         "LorentzModel", 
    519     ], 
    520     "mass_fractal": [ 
    521         "MassFractalModel", 
    522         { 
    523             "cutoff_length": "co_length", 
    524             "radius": "radius", 
    525             "fractal_dim_mass": "mass_dim" 
    526         } 
    527     ], 
    528     "mass_surface_fractal": [ 
    529         "MassSurfaceFractal", 
    530         { 
    531             "rg_cluster": "cluster_rg", 
    532             "fractal_dim_mass": "mass_dim", 
    533             "radius": "radius", 
    534             "fractal_dim_surf": "surface_dim", 
    535             "rg_primary": "primary_rg" 
    536         } 
    537     ], 
    538     "mono_gauss_coil": [ 
    539         "Debye", 
    540         { 
    541             "rg": "rg", 
    542             "i_zero": "scale", 
    543             "background": "background", 
    544         }, 
    545         "DebyeModel", 
    546     ], 
    547     "multilayer_vesicle": [ 
    548         "MultiShellModel", 
    549         { 
    550             "radius": "core_radius", 
    551             "sld_solvent": "core_sld", 
    552             "n_shells": "n_pairs", 
    553             "thick_shell": "s_thickness", 
    554             "sld": "shell_sld", 
    555             "thick_solvent": "w_thickness", 
    556         } 
    557     ], 
    558     "onion": [ 
    559         "OnionExpShellModel", 
    560         { 
    561             "n_shells": "n_shells", 
    562             "A": "A_shell", 
    563             "sld_core": "sld_core0", 
    564             "radius_core": "rad_core0", 
    565             "sld_solvent": "sld_solv", 
    566             "thickness": "thick_shell", 
    567             "sld_in": "sld_in_shell", 
    568             "sld_out": "sld_out_shell" 
    569         } 
    570     ], 
    571     "parallelepiped": [ 
    572         "ParallelepipedModel", 
    573         { 
    574             "phi": "parallel_phi", 
    575             "psi": "parallel_psi", 
    576             "sld_solvent": "sldSolv", 
    577             "length_a": "short_a", 
    578             "length_b": "short_b", 
    579             "sld": "sldPipe", 
    580             "theta": "parallel_theta", 
    581             "length_c": "long_c", 
    582             "M0:sld": "M0_sld_pipe", 
    583             "mtheta:sld": "M_theta_pipe", 
    584             "mphi:sld": "M_phi_pipe", 
    585             "M0:sld_solvent": "M0_sld_solv", 
    586             "mtheta:sld_solvent": "M_theta_solv", 
    587             "mphi:sld_solvent": "M_phi_solv", 
    588             "up:frac_i": "Up_frac_i", 
    589             "up:frac_f": "Up_frac_f", 
    590             "up:angle": "Up_theta", 
    591         } 
    592     ], 
    593     "peak_lorentz": [ 
    594         "Peak Lorentz Model", 
    595         { 
    596             "peak_pos": "q0", 
    597             "peak_hwhm": "B" 
    598         }, 
    599         "PeakLorentzModel", 
    600     ], 
    601     "pearl_necklace": [ 
    602         "PearlNecklaceModel", 
    603         { 
    604             "scale": "scale", 
    605             "thick_string": "thick_string", 
    606             "sld_string": "sld_string", 
    607             "sld_solvent": "sld_solv", 
    608             "edge_sep": "edge_separation", 
    609             "num_pearls": "num_pearls", 
    610             "radius": "radius", 
    611             "background": "background", 
    612             "sld": "sld_pearl" 
    613         } 
    614     ], 
    615     "poly_gauss_coil": [ 
    616         "Poly_GaussCoil", 
    617         { 
    618             "rg": "rg", 
    619             "polydispersity": "poly_m", 
    620             "i_zero": "scale", 
    621             "background": "background", 
    622         } 
    623     ], 
    624     "polymer_excl_volume": [ 
    625         "PolymerExclVolume", 
    626         { 
    627             "rg": "rg", 
    628             "scale": "scale", 
    629             "background": "background", 
    630             "porod_exp": "m" 
    631         } 
    632     ], 
    633     "polymer_micelle": [ 
    634         "MicelleSphCoreModel", 
    635         { 
    636             "sld_corona": "rho_corona", 
    637             "sld_solvent": "rho_solv", 
    638             "sld_core": "rho_core", 
    639             "ndensity": "ndensity", 
    640             "v_core": "v_core", 
    641             "v_corona": "v_corona", 
    642             "radius_core": "radius_core", 
    643             "rg": "radius_gyr", 
    644             "d_penetration": "d_penetration", 
    645             "n_aggreg": "n_aggreg", 
    646         } 
    647     ], 
    648     "porod": [ 
    649         "PorodModel", 
    650         { 
    651             "scale": "scale", 
    652             "background": "background" 
    653         } 
    654     ], 
    655     "power_law": [ 
    656         "PowerLawAbsModel", 
    657         { 
    658             "scale": "scale", 
    659             "background": "background", 
    660             "power": "m" 
    661         } 
    662     ], 
    663     "pringle": [ 
    664         "PringlesModel", 
    665         { 
    666             "scale": "scale", 
    667             "sld_solvent": "sld_solvent", 
    668             "thickness": "thickness", 
    669             "beta": "beta", 
    670             "radius": "radius", 
    671             "background": "background", 
    672             "alpha": "alpha", 
    673             "sld": "sld_pringle" 
    674         } 
    675     ], 
    676     "raspberry": [ 
    677         "RaspBerryModel", 
    678         { 
    679             "volfraction_lg": "volf_Lsph", 
    680             "volfraction_sm": "volf_Ssph", 
    681             "radius_sm": "radius_Ssph", 
    682             "radius_lg": "radius_Lsph", 
    683             "sld_lg": "sld_Lsph", 
    684             "sld_sm": "sld_Ssph", 
    685             "sld_solvent": "sld_solv", 
    686             "surface_fraction": "surfrac_Ssph", 
    687             "penetration": "delta_Ssph" 
    688         } 
    689     ], 
    690     "rectangular_prism": [ 
    691         "RectangularPrismModel", 
    692         { 
    693             "sld": "sldPipe", 
    694             "length_a": "short_side", 
    695             "b2a_ratio": "b2a_ratio", 
    696             "c2a_ratio": "c2a_ratio", 
    697             "sld_solvent": "sldSolv" 
    698         } 
    699     ], 
    700     "rpa": [ 
    701         "RPA10Model", 
    702         { 
    703             "K12": "Kab", "K13": "Kac", "K14": "Kad", 
    704             "K23": "Kbc", "K24": "Kbd", "K34": "Kcd", 
    705             "N1": "Na", "N2": "Nb", "N3": "Nc", "N4": "Nd", 
    706             "L1": "La", "L2": "Lb", "L3": "Lc", "L4": "Ld", 
    707             "v1": "va", "v2": "vb", "v3": "vc", "v4": "vd", 
    708             "b1": "ba", "b2": "bb", "b3": "bc", "b4": "bd", 
    709             "Phi1": "Phia", "Phi2": "Phib", "Phi3": "Phic", "Phi4": "Phid", 
    710             "case_num": "lcase_n" 
    711         } 
    712     ], 
    713     "sc_paracrystal": [ 
    714         "SCCrystalModel", 
    715         { 
    716             "sld": "sldSph", 
    717             "sld_solvent": "sldSolv" 
    718         } 
    719     ], 
    720     "sphere": [ 
    721         "SphereModel", 
    722         { 
    723             "sld": "sldSph", 
    724             "radius": "radius", 
    725             "sld_solvent": "sldSolv", 
    726             "M0:sld": "M0_sld_sph", 
    727             "mtheta:sld": "M_theta_sph", 
    728             "mphi:sld": "M_phi_sph", 
    729             "M0:sld_solvent": "M0_sld_solv", 
    730             "mtheta:sld_solvent": "M_theta_solv", 
    731             "mphi:sld_solvent": "M_phi_solv", 
    732             "up:frac_i": "Up_frac_i", 
    733             "up:frac_f": "Up_frac_f", 
    734             "up:angle": "Up_theta" 
    735         } 
    736     ], 
    737     "spherical_sld": [ 
    738         "SphericalSLDModel", 
    739         # Be lazy and use a generator expression to define 
    740         #    sld1: sld_flat0, ... 
    741         #    thickness1: thick_flat0, ... 
    742         #    interface1: thick_inter0, ... 
    743         #    shape1: func_inter0, ... 
    744         #    nu1: nu_inter0, ... 
    745         # but override thickness1 => rad_cor0 and sld1 => sld_core0. 
    746         # Note: explicit key,value pairs given by **{...} override the 
    747         # keys from the gnerator expression ((k,v) for k,v in seq) when 
    748         # used as dict((generator), **{...}) 
    749         dict(((field_new+str(index+1), field_old+str(index)) 
    750               for field_new, field_old in [("sld", "sld_flat"), 
    751                                            ("thickness", "thick_flat"), 
    752                                            ("interface", "thick_inter"), 
    753                                            ("shape", "func_inter"), 
    754                                            ("nu", "nu_inter"),] 
    755               for index in range(11)), 
    756              **{ 
    757                    "n_shells": "n_shells", 
    758                    "n_steps": "npts_inter", 
    759                    "sld_solvent": "sld_solv", 
    760                    "thickness1": "rad_core0", 
    761                    "sld1": "sld_core0", 
    762                }) 
    763     ], 
    764     "squarewell": [ 
    765         "SquareWellStructure", 
    766         { 
    767             "scale": "scale_factor", 
    768             "radius_effective": "effect_radius", 
    769             "wellwidth": "wellwidth", 
    770             "welldepth": "welldepth", 
    771         } 
    772     ], 
    773     "stacked_disks": [ 
    774         "StackedDisksModel", 
    775         { 
    776             "phi": "axis_phi", 
    777             "sld_layer": "layer_sld", 
    778             "sld_core": "core_sld", 
    779             "theta": "axis_theta", 
    780             "sld_solvent": "solvent_sld", 
    781             "n_stacking": "n_stacking", 
    782             "thick_layer": "layer_thick", 
    783             "thick_core": "core_thick", 
    784         } 
    785     ], 
    786     "star_polymer": [ 
    787         "StarPolymer", 
    788         { 
    789             "arms": "arms", 
    790             "rg_squared": "R2" 
    791         } 
    792     ], 
    793     "stickyhardsphere": [ 
    794         "StickyHSStructure", 
    795         { 
    796             "scale": "scale_factor", 
    797             "radius_effective": "effect_radius", 
    798         } 
    799     ], 
    800     "surface_fractal": [ 
    801         "SurfaceFractalModel", 
    802         { 
    803             "cutoff_length": "co_length", 
    804             "radius": "radius", 
    805             "fractal_dim_surf": "surface_dim" 
    806         } 
    807     ], 
    808     "teubner_strey": [ 
    809         "TeubnerStreyModel", 
    810         { 
    811             # Note: parameters are completely rewritten in convert.py 
    812             "volfraction_a": "volfraction_a", 
    813             "sld_a": "sld_a", 
    814             "sld_b": "sld_b", 
    815             "d": "d", 
    816             "xi": "xi", 
    817         } 
    818     ], 
    819     "triaxial_ellipsoid": [ 
    820         "TriaxialEllipsoidModel", 
    821         { 
    822             "phi": "axis_phi", 
    823             "radius_equat_minor": "semi_axisA", 
    824             "radius_polar": "semi_axisC", 
    825             "radius_equat_major": "semi_axisB", 
    826             "sld_solvent": "sldSolv", 
    827             "psi": "axis_psi", 
    828             "sld": "sldEll", 
    829             "theta": "axis_theta" 
    830         } 
    831     ], 
    832     "two_lorentzian": [ 
    833         "TwoLorentzian", 
    834         { 
    835             "lorentz_scale_1": "scale_1", 
    836             "lorentz_scale_2": "scale_2", 
    837             "lorentz_exp_1": "exponent_1", 
    838             "lorentz_exp_2": "exponent_2", 
    839             "lorentz_length_2": "length_2", 
    840             "lorentz_length_1": "length_1", 
    841             "background": "background" 
    842         }, 
    843         "TwoLorentzianModel", 
    844     ], 
    845     "two_power_law": [ 
    846         "TwoPowerLaw", 
    847         { 
    848             "coefficent_1": "coef_A", 
    849             "power_2": "power2", 
    850             "power_1": "power1", 
    851             "background": "background", 
    852             "crossover": "qc" 
    853         }, 
    854         "TwoPowerLawModel", 
    855     ], 
    856     "unified_power_Rg": [ 
    857         "UnifiedPowerRg", 
    858         dict(((field_new+str(index), field_old+str(index)) 
    859               for field_new, field_old in [("rg", "Rg"), 
    860                                            ("power", "power"), 
    861                                            ("G", "G"), 
    862                                            ("B", "B"),] 
    863               for index in range(11)), 
    864              **{ 
    865                    "background": "background", 
    866                    "scale": "scale", 
    867                }), 
    868         "UnifiedPowerRgModel", 
    869     ], 
    870     "vesicle": [ 
    871         "VesicleModel", 
    872         { 
    873             "sld": "shell_sld", 
    874             "sld_solvent": "solv_sld" 
    875         } 
    876     ] 
     31    (3, 1, 2) : { 
     32        "adsorbed_layer": [ 
     33            "Core2ndMomentModel", 
     34            { 
     35                "scale": "scale", 
     36                "second_moment": "second_moment", 
     37                "density_shell": "density_poly", 
     38                "sld_solvent": "sld_solv", 
     39                "radius": "radius_core", 
     40                "volfraction": "volf_cores", 
     41                "background": "background", 
     42                "adsorbed_amount": "ads_amount", 
     43                "sld_shell": "sld_poly" 
     44            } 
     45        ], 
     46        "barbell": [ 
     47            "BarBellModel", 
     48            { 
     49                "sld": "sld_barbell", 
     50                "length": "len_bar", 
     51                "radius_bell": "rad_bell", 
     52                "radius": "rad_bar", 
     53                "sld_solvent": "sld_solv" 
     54            } 
     55        ], 
     56        "bcc_paracrystal": [ 
     57            "BCCrystalModel", 
     58            { 
     59                "sld": "sldSph", 
     60                "sld_solvent": "sldSolv" 
     61            } 
     62        ], 
     63        "be_polyelectrolyte": [ 
     64            "BEPolyelectrolyte", 
     65            { 
     66                "ionization_degree": "alpha", 
     67                "polymer_concentration": "c", 
     68                "salt_concentration": "cs", 
     69                "virial_param": "h", 
     70                "background": "background", 
     71                "contrast_factor": "k", 
     72                "bjerrum_length": "lb", 
     73                "monomer_length": "b" 
     74            } 
     75        ], 
     76        "binary_hard_sphere": [ 
     77            "BinaryHSModel", 
     78            { 
     79                "sld_sm": "ss_sld", 
     80                "sld_lg": "ls_sld", 
     81                "volfraction_sm": "vol_frac_ss", 
     82                "radius_lg": "l_radius", 
     83                "radius_sm": "s_radius", 
     84                "volfraction_lg": "vol_frac_ls", 
     85                "sld_solvent": "solvent_sld" 
     86            } 
     87        ], 
     88        "broad_peak": [ 
     89            "BroadPeakModel", 
     90            { 
     91                "peak_pos": "q_peak", 
     92                "scale": None, 
     93                "lorentz_length": "length_l", 
     94                "porod_scale": "scale_p", 
     95                "lorentz_exp": "exponent_l", 
     96                "lorentz_scale": "scale_l", 
     97                "porod_exp": "exponent_p" 
     98            } 
     99        ], 
     100        "capped_cylinder": [ 
     101            "CappedCylinderModel", 
     102            { 
     103                "sld": "sld_capcyl", 
     104                "length": "len_cyl", 
     105                "radius_cap": "rad_cap", 
     106                "radius": "rad_cyl", 
     107                "sld_solvent": "sld_solv" 
     108            } 
     109        ], 
     110        "core_multi_shell": [ 
     111            "CoreMultiShellModel", 
     112            { 
     113                "thickness": "thick_shell", 
     114                "sld": "sld_shell", 
     115                "radius": "rad_core0", 
     116                "sld_core": "sld_core0", 
     117                "sld_solvent": "sld_solv", 
     118                "n": "n_shells", 
     119                "M0:sld_core": "M0_sld_core0", 
     120                "mtheta:sld_core": "M_theta_core0", 
     121                "mphi:sld_core": "M_phi_core0", 
     122                "M0:sld1": "M0_sld_shell1", 
     123                "mtheta:sld1": "M_theta_shell1", 
     124                "mphi:sld1": "M_phi_shell1", 
     125                "M0:sld2": "M0_sld_shell2", 
     126                "mtheta:sld2": "M_theta_shell2", 
     127                "mphi:sld2": "M_phi_shell2", 
     128                "M0:sld3": "M0_sld_shell3", 
     129                "mtheta:sld3": "M_theta_shell3", 
     130                "mphi:sld3": "M_phi_shell3", 
     131                "M0:sld4": "M0_sld_shell4", 
     132                "mtheta:sld4": "M_theta_shell4", 
     133                "mphi:sld4": "M_phi_shell4", 
     134                "M0:sld_solvent": "M0_sld_solv", 
     135                "mtheta:sld_solvent": "M_theta_solv", 
     136                "mphi:sld_solvent": "M_phi_solv", 
     137                "up:frac_i": "Up_frac_i", 
     138                "up:frac_f": "Up_frac_f", 
     139                "up:angle": "Up_theta", 
     140            } 
     141        ], 
     142        "core_shell_bicelle": [ 
     143            "CoreShellBicelleModel", 
     144            { 
     145                "phi": "axis_phi", 
     146                "sld_core": "core_sld", 
     147                "sld_rim": "rim_sld", 
     148                "thick_face": "face_thick", 
     149                "sld_solvent": "solvent_sld", 
     150                "thick_rim": "rim_thick", 
     151                "sld_face": "face_sld", 
     152                "theta": "axis_theta" 
     153            } 
     154        ], 
     155        "core_shell_cylinder": [ 
     156            "CoreShellCylinderModel", 
     157            { 
     158                "theta": "axis_theta", 
     159                "phi": "axis_phi", 
     160                "sld_shell": "shell_sld", 
     161                "sld_solvent": "solvent_sld", 
     162                "sld_core": "core_sld" 
     163            } 
     164        ], 
     165        "core_shell_ellipsoid:1": [ 
     166            "CoreShellEllipsoidModel", 
     167            { 
     168                "sld_core": "sld_core", 
     169                "sld_shell": "sld_shell", 
     170                "sld_solvent": "sld_solvent", 
     171                "radius_equat_core": "equat_core", 
     172                "x_core": "polar_core", 
     173                "thick_shell": "equat_shell", 
     174                "x_polar_shell": "polar_shell", 
     175                "theta": "axis_theta", 
     176                "phi": "axis_phi", 
     177            } 
     178        ], 
     179        "core_shell_ellipsoid": [ 
     180            "CoreShellEllipsoidXTModel", 
     181            { 
     182                "sld_core": "sld_core", 
     183                "sld_shell": "sld_shell", 
     184                "sld_solvent": "sld_solvent", 
     185                "radius_equat_core": "equat_core", 
     186                "thick_shell": "T_shell", 
     187                "x_core": "X_core", 
     188                "x_polar_shell": "XpolarShell", 
     189                "theta": "axis_theta", 
     190                "phi": "axis_phi", 
     191            } 
     192        ], 
     193        "core_shell_parallelepiped": [ 
     194            "CSParallelepipedModel", 
     195            { 
     196                "sld_core": "sld_pcore", 
     197                "sld_a": "sld_rimA", 
     198                "sld_b": "sld_rimB", 
     199                "sld_c": "sld_rimC", 
     200                "sld_solvent": "sld_solv", 
     201                "length_a": "shortA", 
     202                "length_b": "midB", 
     203                "length_c": "longC", 
     204                "thick_rim_a": "rimA", 
     205                "thick_rim_c": "rimC", 
     206                "thick_rim_b": "rimB", 
     207                "theta": "parallel_theta", 
     208                "phi": "parallel_phi", 
     209                "psi": "parallel_psi", 
     210            } 
     211        ], 
     212        "core_shell_sphere": [ 
     213            "CoreShellModel", 
     214            { 
     215                "sld_core": "core_sld", 
     216                "sld_shell": "shell_sld", 
     217                "sld_solvent": "solvent_sld", 
     218                "M0:sld_core": "M0_sld_core", 
     219                "mtheta:sld_core": "M_theta_core", 
     220                "mphi:sld_core": "M_phi_core", 
     221                "M0:sld_shell": "M0_sld_shell", 
     222                "mtheta:sld_shell": "M_theta_shell", 
     223                "mphi:sld_shell": "M_phi_shell", 
     224                "M0:sld_solvent": "M0_sld_solv", 
     225                "mtheta:sld_solvent": "M_theta_solv", 
     226                "mphi:sld_solvent": "M_phi_solv", 
     227                "up:frac_i": "Up_frac_i", 
     228                "up:frac_f": "Up_frac_f", 
     229                "up:angle": "Up_theta" 
     230            } 
     231        ], 
     232        "correlation_length": [ 
     233            "CorrLength", 
     234            { 
     235                "porod_scale": "scale_p", 
     236                "lorentz_scale": "scale_l", 
     237                "porod_exp": "exponent_p", 
     238                "lorentz_exp": "exponent_l", 
     239                "cor_length": "length_l" 
     240            }, 
     241            "CorrLengthModel" 
     242        ], 
     243        "cylinder": [ 
     244            "CylinderModel", 
     245            { 
     246                "sld": "sldCyl", 
     247                "theta": "cyl_theta", 
     248                "phi": "cyl_phi", 
     249                "sld_solvent": "sldSolv", 
     250                "M0:sld": "M0_sld_cyl", 
     251                "mtheta:sld": "M_theta_cyl", 
     252                "mphi:sld": "M_phi_cyl", 
     253                "M0:sld_solvent": "M0_sld_solv", 
     254                "mtheta:sld_solvent": "M_theta_solv", 
     255                "mphi:sld_solvent": "M_phi_solv", 
     256                "up:frac_i": "Up_frac_i", 
     257                "up:frac_f": "Up_frac_f", 
     258                "up:angle": "Up_theta" 
     259            } 
     260        ], 
     261        "dab": [ 
     262            "DABModel", 
     263            { 
     264                "cor_length": "length" 
     265            } 
     266        ], 
     267        "ellipsoid": [ 
     268            "EllipsoidModel", 
     269            { 
     270                "phi": "axis_phi", 
     271                "radius_equatorial": "radius_b", 
     272                "sld": "sldEll", 
     273                "theta": "axis_theta", 
     274                "radius_polar": "radius_a", 
     275                "sld_solvent": "sldSolv" 
     276            } 
     277        ], 
     278        "elliptical_cylinder": [ 
     279            "EllipticalCylinderModel", 
     280            { 
     281                "axis_ratio": "r_ratio", 
     282                "radius_minor": "r_minor", 
     283                "sld": "sldCyl", 
     284                "sld_solvent": "sldSolv", 
     285                "theta": "cyl_theta", 
     286                "phi": "cyl_phi", 
     287                "psi": "cyl_psi", 
     288            } 
     289        ], 
     290        "fcc_paracrystal": [ 
     291            "FCCrystalModel", 
     292            { 
     293                "sld": "sldSph", 
     294                "sld_solvent": "sldSolv" 
     295            } 
     296        ], 
     297        "flexible_cylinder": [ 
     298            "FlexibleCylinderModel", 
     299            { 
     300                "sld": "sldCyl", 
     301                "sld_solvent": "sldSolv" 
     302            } 
     303        ], 
     304        "flexible_cylinder_elliptical": [ 
     305            "FlexCylEllipXModel", 
     306            { 
     307                "sld": "sldCyl", 
     308                "sld_solvent": "sldSolv" 
     309            } 
     310        ], 
     311        "fractal": [ 
     312            "FractalModel", 
     313            { 
     314                "sld_block": "sldBlock", 
     315                "radius": "radius", 
     316                "cor_length": "cor_length", 
     317                "sld_solvent": "sldSolv", 
     318                "fractal_dim": "fractal_dim" 
     319            } 
     320        ], 
     321        "fractal_core_shell": [ 
     322            "FractalCoreShell", 
     323            { 
     324                "sld_core": "core_sld", 
     325                "sld_shell": "shell_sld", 
     326                "sld_solvent": "solvent_sld", 
     327                "radius": "radius", 
     328                "thickness": "thickness", 
     329                "fractal_dim": "frac_dim", 
     330                "cor_length": "cor_length", 
     331                "volfraction": "volfraction", 
     332            }, 
     333            "FractalCoreShellModel" 
     334        ], 
     335        "fuzzy_sphere": [ 
     336            "FuzzySphereModel", 
     337            { 
     338                "sld": "sldSph", 
     339                "fuzziness": "fuzziness", 
     340                "radius": "radius", 
     341                "sld_solvent": "sldSolv" 
     342            } 
     343        ], 
     344        "gauss_lorentz_gel": [ 
     345            "GaussLorentzGel", 
     346            { 
     347                "gauss_scale": "scale_g", 
     348                "cor_length_dynamic": "dyn_colength", 
     349                "cor_length_static": "stat_colength", 
     350                "background": "background", 
     351                "lorentz_scale": "scale_l" 
     352            }, 
     353            "GaussLorentzGelModel" 
     354        ], 
     355        "gaussian_peak": [ 
     356            "Peak Gauss Model", 
     357            { 
     358                "peak_pos": "q0", 
     359                "sigma": "B", 
     360            }, 
     361            "PeakGaussModel", 
     362        ], 
     363        "gel_fit": [ 
     364            "GelFitModel", 
     365            { 
     366                "rg": "radius", 
     367                "lorentz_scale": "lScale", 
     368                "guinier_scale": "gScale", 
     369                "fractal_dim": "FractalExp", 
     370                "cor_length": "zeta", 
     371            } 
     372        ], 
     373        "guinier": [ 
     374            "Guinier", 
     375            { 
     376                "rg": "rg" 
     377            }, 
     378            "GuinierModel", 
     379        ], 
     380        "guinier_porod": [ 
     381            "GuinierPorod", 
     382            { 
     383                "s": "dim", 
     384                "rg": "rg", 
     385                "porod_exp": "m", 
     386                "scale": "scale", 
     387                "background": "background" 
     388            }, 
     389            "GuinierPorodModel", 
     390        ], 
     391        "hardsphere": [ 
     392            "HardsphereStructure", 
     393            { 
     394                "scale": "scale_factor", 
     395                "radius_effective": "effect_radius", 
     396            } 
     397        ], 
     398        "hayter_msa": [ 
     399            "HayterMSAStructure", 
     400            { 
     401                "scale": "scale_factor", 
     402                "radius_effective": "effect_radius", 
     403                "volfraction": "volfraction", 
     404                "charge": "charge", 
     405                "temperature": "temperature", 
     406                "concentration_salt": "saltconc", 
     407                "dielectconst": "dielectconst", 
     408            } 
     409        ], 
     410        "hollow_cylinder": [ 
     411            "HollowCylinderModel", 
     412            { 
     413                "sld": "sldCyl", 
     414                "sld_solvent": "sldSolv", 
     415                "radius": "core_radius", 
     416                "thickness": "radius", 
     417                "length": "length", 
     418                "theta": "axis_theta", 
     419                "phi": "axis_phi", 
     420            } 
     421        ], 
     422        "hollow_rectangular_prism": [ 
     423            "RectangularHollowPrismModel", 
     424            { 
     425                "sld": "sldPipe", 
     426                "sld_solvent": "sldSolv", 
     427                "length_a": "short_side", 
     428                "b2a_ratio": "b2a_ratio", 
     429                "c2a_ratio": "c2a_ratio", 
     430                "thickness": "thickness", 
     431            } 
     432        ], 
     433        "hollow_rectangular_prism_thin_walls": [ 
     434            "RectangularHollowPrismInfThinWallsModel", 
     435            { 
     436                "sld": "sldPipe", 
     437                "sld_solvent": "sldSolv", 
     438                "length_a": "short_side", 
     439                "b2a_ratio": "b2a_ratio", 
     440                "c2a_ratio": "c2a_ratio", 
     441            } 
     442        ], 
     443        "lamellar": [ 
     444            "LamellarModel", 
     445            { 
     446                "sld": "sld_bi", 
     447                "sld_solvent": "sld_sol", 
     448                "thickness": "bi_thick" 
     449            } 
     450        ], 
     451        "lamellar_hg": [ 
     452            "LamellarFFHGModel", 
     453            { 
     454                "sld": "sld_tail", 
     455                "sld_solvent": "sld_solvent", 
     456                "sld_head": "sld_head", 
     457                "length_tail": "t_length", 
     458                "length_head": "h_thickness" 
     459            } 
     460        ], 
     461        "lamellar_hg_stack_caille": [ 
     462            "LamellarPSHGModel", 
     463            { 
     464                "sld": "sld_tail", 
     465                "sld_head": "sld_head", 
     466                "sld_solvent": "sld_solvent", 
     467                "length_tail": "deltaT", 
     468                "length_head": "deltaH", 
     469                "d_spacing": "spacing", 
     470                "Caille_parameter": "caille", 
     471                "Nlayers": "n_plates", 
     472            } 
     473        ], 
     474        "lamellar_stack_caille": [ 
     475            "LamellarPSModel", 
     476            { 
     477                "sld": "sld_bi", 
     478                "sld_solvent": "sld_sol", 
     479                "thickness": "delta", 
     480                "d_spacing": "spacing", 
     481                "Caille_parameter": "caille", 
     482                "Nlayers": "n_plates", 
     483            } 
     484        ], 
     485        "lamellar_stack_paracrystal": [ 
     486            "LamellarPCrystalModel", 
     487            { 
     488                "sld": "sld_layer", 
     489                "sld_solvent": "sld_solvent", 
     490                "thickness": "thickness", 
     491                "d_spacing": "spacing", 
     492                "sigma_d": "pd_spacing", 
     493                "Nlayers": "Nlayers", 
     494            } 
     495        ], 
     496        "line": [ 
     497            "LineModel", 
     498            { 
     499                "slope": "B", 
     500                "scale": None, 
     501                "background": None, 
     502                "intercept": "A" 
     503            } 
     504        ], 
     505        "linear_pearls": [ 
     506            "LinearPearlsModel", 
     507            { 
     508                "sld": "sld_pearl", 
     509                "sld_solvent": "sld_solv", 
     510                "edge_sep": "edge_separation" 
     511            } 
     512        ], 
     513        "lorentz": [ 
     514            "Lorentz", 
     515            { 
     516                "cor_length": "length" 
     517            }, 
     518            "LorentzModel", 
     519        ], 
     520        "mass_fractal": [ 
     521            "MassFractalModel", 
     522            { 
     523                "cutoff_length": "co_length", 
     524                "radius": "radius", 
     525                "fractal_dim_mass": "mass_dim" 
     526            } 
     527        ], 
     528        "mass_surface_fractal": [ 
     529            "MassSurfaceFractal", 
     530            { 
     531                "rg_cluster": "cluster_rg", 
     532                "fractal_dim_mass": "mass_dim", 
     533                "radius": "radius", 
     534                "fractal_dim_surf": "surface_dim", 
     535                "rg_primary": "primary_rg" 
     536            } 
     537        ], 
     538        "mono_gauss_coil": [ 
     539            "Debye", 
     540            { 
     541                "rg": "rg", 
     542                "i_zero": "scale", 
     543                "background": "background", 
     544            }, 
     545            "DebyeModel", 
     546        ], 
     547        "multilayer_vesicle": [ 
     548            "MultiShellModel", 
     549            { 
     550                "radius": "core_radius", 
     551                "sld_solvent": "core_sld", 
     552                "n_shells": "n_pairs", 
     553                "thick_shell": "s_thickness", 
     554                "sld": "shell_sld", 
     555                "thick_solvent": "w_thickness", 
     556            } 
     557        ], 
     558        "onion": [ 
     559            "OnionExpShellModel", 
     560            { 
     561                "n_shells": "n_shells", 
     562                "A": "A_shell", 
     563                "sld_core": "sld_core0", 
     564                "radius_core": "rad_core0", 
     565                "sld_solvent": "sld_solv", 
     566                "thickness": "thick_shell", 
     567                "sld_in": "sld_in_shell", 
     568                "sld_out": "sld_out_shell" 
     569            } 
     570        ], 
     571        "parallelepiped": [ 
     572            "ParallelepipedModel", 
     573            { 
     574                "phi": "parallel_phi", 
     575                "psi": "parallel_psi", 
     576                "sld_solvent": "sldSolv", 
     577                "length_a": "short_a", 
     578                "length_b": "short_b", 
     579                "sld": "sldPipe", 
     580                "theta": "parallel_theta", 
     581                "length_c": "long_c", 
     582                "M0:sld": "M0_sld_pipe", 
     583                "mtheta:sld": "M_theta_pipe", 
     584                "mphi:sld": "M_phi_pipe", 
     585                "M0:sld_solvent": "M0_sld_solv", 
     586                "mtheta:sld_solvent": "M_theta_solv", 
     587                "mphi:sld_solvent": "M_phi_solv", 
     588                "up:frac_i": "Up_frac_i", 
     589                "up:frac_f": "Up_frac_f", 
     590                "up:angle": "Up_theta", 
     591            } 
     592        ], 
     593        "peak_lorentz": [ 
     594            "Peak Lorentz Model", 
     595            { 
     596                "peak_pos": "q0", 
     597                "peak_hwhm": "B" 
     598            }, 
     599            "PeakLorentzModel", 
     600        ], 
     601        "pearl_necklace": [ 
     602            "PearlNecklaceModel", 
     603            { 
     604                "scale": "scale", 
     605                "thick_string": "thick_string", 
     606                "sld_string": "sld_string", 
     607                "sld_solvent": "sld_solv", 
     608                "edge_sep": "edge_separation", 
     609                "num_pearls": "num_pearls", 
     610                "radius": "radius", 
     611                "background": "background", 
     612                "sld": "sld_pearl" 
     613            } 
     614        ], 
     615        "poly_gauss_coil": [ 
     616            "Poly_GaussCoil", 
     617            { 
     618                "rg": "rg", 
     619                "polydispersity": "poly_m", 
     620                "i_zero": "scale", 
     621                "background": "background", 
     622            } 
     623        ], 
     624        "polymer_excl_volume": [ 
     625            "PolymerExclVolume", 
     626            { 
     627                "rg": "rg", 
     628                "scale": "scale", 
     629                "background": "background", 
     630                "porod_exp": "m" 
     631            } 
     632        ], 
     633        "polymer_micelle": [ 
     634            "MicelleSphCoreModel", 
     635            { 
     636                "sld_corona": "rho_corona", 
     637                "sld_solvent": "rho_solv", 
     638                "sld_core": "rho_core", 
     639                "ndensity": "ndensity", 
     640                "v_core": "v_core", 
     641                "v_corona": "v_corona", 
     642                "radius_core": "radius_core", 
     643                "rg": "radius_gyr", 
     644                "d_penetration": "d_penetration", 
     645                "n_aggreg": "n_aggreg", 
     646            } 
     647        ], 
     648        "porod": [ 
     649            "PorodModel", 
     650            { 
     651                "scale": "scale", 
     652                "background": "background" 
     653            } 
     654        ], 
     655        "power_law": [ 
     656            "PowerLawAbsModel", 
     657            { 
     658                "scale": "scale", 
     659                "background": "background", 
     660                "power": "m" 
     661            } 
     662        ], 
     663        "pringle": [ 
     664            "PringlesModel", 
     665            { 
     666                "scale": "scale", 
     667                "sld_solvent": "sld_solvent", 
     668                "thickness": "thickness", 
     669                "beta": "beta", 
     670                "radius": "radius", 
     671                "background": "background", 
     672                "alpha": "alpha", 
     673                "sld": "sld_pringle" 
     674            } 
     675        ], 
     676        "raspberry": [ 
     677            "RaspBerryModel", 
     678            { 
     679                "volfraction_lg": "volf_Lsph", 
     680                "volfraction_sm": "volf_Ssph", 
     681                "radius_sm": "radius_Ssph", 
     682                "radius_lg": "radius_Lsph", 
     683                "sld_lg": "sld_Lsph", 
     684                "sld_sm": "sld_Ssph", 
     685                "sld_solvent": "sld_solv", 
     686                "surface_fraction": "surfrac_Ssph", 
     687                "penetration": "delta_Ssph" 
     688            } 
     689        ], 
     690        "rectangular_prism": [ 
     691            "RectangularPrismModel", 
     692            { 
     693                "sld": "sldPipe", 
     694                "length_a": "short_side", 
     695                "b2a_ratio": "b2a_ratio", 
     696                "c2a_ratio": "c2a_ratio", 
     697                "sld_solvent": "sldSolv" 
     698            } 
     699        ], 
     700        "rpa": [ 
     701            "RPA10Model", 
     702            { 
     703                "K12": "Kab", "K13": "Kac", "K14": "Kad", 
     704                "K23": "Kbc", "K24": "Kbd", "K34": "Kcd", 
     705                "N1": "Na", "N2": "Nb", "N3": "Nc", "N4": "Nd", 
     706                "L1": "La", "L2": "Lb", "L3": "Lc", "L4": "Ld", 
     707                "v1": "va", "v2": "vb", "v3": "vc", "v4": "vd", 
     708                "b1": "ba", "b2": "bb", "b3": "bc", "b4": "bd", 
     709                "Phi1": "Phia", "Phi2": "Phib", "Phi3": "Phic", "Phi4": "Phid", 
     710                "case_num": "lcase_n" 
     711            } 
     712        ], 
     713        "sc_paracrystal": [ 
     714            "SCCrystalModel", 
     715            { 
     716                "sld": "sldSph", 
     717                "sld_solvent": "sldSolv" 
     718            } 
     719        ], 
     720        "sphere": [ 
     721            "SphereModel", 
     722            { 
     723                "sld": "sldSph", 
     724                "radius": "radius", 
     725                "sld_solvent": "sldSolv", 
     726                "M0:sld": "M0_sld_sph", 
     727                "mtheta:sld": "M_theta_sph", 
     728                "mphi:sld": "M_phi_sph", 
     729                "M0:sld_solvent": "M0_sld_solv", 
     730                "mtheta:sld_solvent": "M_theta_solv", 
     731                "mphi:sld_solvent": "M_phi_solv", 
     732                "up:frac_i": "Up_frac_i", 
     733                "up:frac_f": "Up_frac_f", 
     734                "up:angle": "Up_theta" 
     735            } 
     736        ], 
     737        "spherical_sld": [ 
     738            "SphericalSLDModel", 
     739            # Be lazy and use a generator expression to define 
     740            #    sld1: sld_flat0, ... 
     741            #    thickness1: thick_flat0, ... 
     742            #    interface1: thick_inter0, ... 
     743            #    shape1: func_inter0, ... 
     744            #    nu1: nu_inter0, ... 
     745            # but override thickness1 => rad_cor0 and sld1 => sld_core0. 
     746            # Note: explicit key,value pairs given by **{...} override the 
     747            # keys from the gnerator expression ((k,v) for k,v in seq) when 
     748            # used as dict((generator), **{...}) 
     749            dict(((field_new+str(index+1), field_old+str(index)) 
     750                  for field_new, field_old in [("sld", "sld_flat"), 
     751                                               ("thickness", "thick_flat"), 
     752                                               ("interface", "thick_inter"), 
     753                                               ("shape", "func_inter"), 
     754                                               ("nu", "nu_inter"),] 
     755                  for index in range(11)), 
     756                 **{ 
     757                     "n_shells": "n_shells", 
     758                     "n_steps": "npts_inter", 
     759                     "sld_solvent": "sld_solv", 
     760                     "thickness1": "rad_core0", 
     761                     "sld1": "sld_core0", 
     762                 }) 
     763        ], 
     764        "squarewell": [ 
     765            "SquareWellStructure", 
     766            { 
     767                "scale": "scale_factor", 
     768                "radius_effective": "effect_radius", 
     769                "wellwidth": "wellwidth", 
     770                "welldepth": "welldepth", 
     771            } 
     772        ], 
     773        "stacked_disks": [ 
     774            "StackedDisksModel", 
     775            { 
     776                "phi": "axis_phi", 
     777                "sld_layer": "layer_sld", 
     778                "sld_core": "core_sld", 
     779                "theta": "axis_theta", 
     780                "sld_solvent": "solvent_sld", 
     781                "n_stacking": "n_stacking", 
     782                "thick_layer": "layer_thick", 
     783                "thick_core": "core_thick", 
     784            } 
     785        ], 
     786        "star_polymer": [ 
     787            "StarPolymer", 
     788            { 
     789                "arms": "arms", 
     790                "rg_squared": "R2" 
     791            } 
     792        ], 
     793        "stickyhardsphere": [ 
     794            "StickyHSStructure", 
     795            { 
     796                "scale": "scale_factor", 
     797                "radius_effective": "effect_radius", 
     798            } 
     799        ], 
     800        "surface_fractal": [ 
     801            "SurfaceFractalModel", 
     802            { 
     803                "cutoff_length": "co_length", 
     804                "radius": "radius", 
     805                "fractal_dim_surf": "surface_dim" 
     806            } 
     807        ], 
     808        "teubner_strey": [ 
     809            "TeubnerStreyModel", 
     810            { 
     811                # Note: parameters are completely rewritten in convert.py 
     812                "volfraction_a": "volfraction_a", 
     813                "sld_a": "sld_a", 
     814                "sld_b": "sld_b", 
     815                "d": "d", 
     816                "xi": "xi", 
     817            } 
     818        ], 
     819        "triaxial_ellipsoid": [ 
     820            "TriaxialEllipsoidModel", 
     821            { 
     822                "phi": "axis_phi", 
     823                "radius_equat_minor": "semi_axisA", 
     824                "radius_polar": "semi_axisC", 
     825                "radius_equat_major": "semi_axisB", 
     826                "sld_solvent": "sldSolv", 
     827                "psi": "axis_psi", 
     828                "sld": "sldEll", 
     829                "theta": "axis_theta" 
     830            } 
     831        ], 
     832        "two_lorentzian": [ 
     833            "TwoLorentzian", 
     834            { 
     835                "lorentz_scale_1": "scale_1", 
     836                "lorentz_scale_2": "scale_2", 
     837                "lorentz_exp_1": "exponent_1", 
     838                "lorentz_exp_2": "exponent_2", 
     839                "lorentz_length_2": "length_2", 
     840                "lorentz_length_1": "length_1", 
     841                "background": "background" 
     842            }, 
     843            "TwoLorentzianModel", 
     844        ], 
     845        "two_power_law": [ 
     846            "TwoPowerLaw", 
     847            { 
     848                "coefficent_1": "coef_A", 
     849                "power_2": "power2", 
     850                "power_1": "power1", 
     851                "background": "background", 
     852                "crossover": "qc" 
     853            }, 
     854            "TwoPowerLawModel", 
     855        ], 
     856        "unified_power_Rg": [ 
     857            "UnifiedPowerRg", 
     858            dict(((field_new+str(index), field_old+str(index)) 
     859                  for field_new, field_old in [("rg", "Rg"), 
     860                                               ("power", "power"), 
     861                                               ("G", "G"), 
     862                                               ("B", "B"),] 
     863                  for index in range(11)), 
     864                 **{ 
     865                     "background": "background", 
     866                     "scale": "scale", 
     867                 }), 
     868            "UnifiedPowerRgModel", 
     869        ], 
     870        "vesicle": [ 
     871            "VesicleModel", 
     872            { 
     873                "sld": "shell_sld", 
     874                "sld_solvent": "solv_sld" 
     875            } 
     876        ] 
    877877    } 
    878878} 
  • sasmodels/convert.py

    r07c8d46 r2d81cfe  
    44from __future__ import print_function, division 
    55 
    6 import re 
    76import math 
    87import warnings 
     8 
     9import numpy as np 
    910 
    1011from .conversion_table import CONVERSION_TABLE 
     
    6465    return [pk*scale for pk in par] if isinstance(par, list) else par*scale 
    6566 
    66 def _is_sld(model_info, id): 
     67def _is_sld(model_info, par): 
    6768    """ 
    6869    Return True if parameter is a magnetic magnitude or SLD parameter. 
    6970    """ 
    70     if id.startswith('M0:'): 
     71    if par.startswith('M0:'): 
    7172        return True 
    72     if '_pd' in id or '.' in id: 
     73    if '_pd' in par or '.' in par: 
    7374        return False 
    7475    for p in model_info.parameters.call_parameters: 
    75         if p.id == id: 
     76        if p.id == par: 
    7677            return p.type == 'sld' 
    7778    # check through kernel parameters in case it is a named as a vector 
    7879    for p in model_info.parameters.kernel_parameters: 
    79         if p.id == id: 
     80        if p.id == par: 
    8081            return p.type == 'sld' 
    8182    return False 
     
    8889    *scale=1e-6*.  For forward conversion use *scale=1e6*. 
    8990    """ 
    90     return dict((id, (_rescale(v, scale) if _is_sld(model_info, id) else v)) 
    91                 for id, v in pars.items()) 
    92  
    93  
    94 def _get_translation_table(model_info, version=(3,1,2)): 
     91    return dict((par, (_rescale(v, scale) if _is_sld(model_info, par) else v)) 
     92                for par, v in pars.items()) 
     93 
     94 
     95def _get_translation_table(model_info, version=(3, 1, 2)): 
    9596    conv_param = CONVERSION_TABLE.get(version, {}).get(model_info.id, [None, {}]) 
    9697    translation = conv_param[1].copy() 
     
    130131    newpars = pars.copy() 
    131132    for new, old in mapping.items(): 
    132         if old == new: continue 
    133         if old is None: continue 
    134         for underscore, dot in PD_DOT: 
     133        if old == new: 
     134            continue 
     135        if old is None: 
     136            continue 
     137        for _, dot in PD_DOT: 
    135138            source = old+dot 
    136139            if source in newpars: 
     
    145148    return newpars 
    146149 
    147 def _conversion_target(model_name, version=(3,1,2)): 
     150def _conversion_target(model_name, version=(3, 1, 2)): 
    148151    """ 
    149152    Find the sasmodel name which translates into the sasview name. 
     
    159162    return None 
    160163 
    161 def _hand_convert(name, oldpars, version=(3,1,2)): 
    162     if version == (3,1,2): 
     164def _hand_convert(name, oldpars, version=(3, 1, 2)): 
     165    if version == (3, 1, 2): 
    163166        oldpars = _hand_convert_3_1_2_to_4_1(name, oldpars) 
    164167    return oldpars 
     
    272275        p_scale = oldpars['scale'] 
    273276        p_c1 = oldpars['c1'] 
    274         p_c2= oldpars['c2'] 
     277        p_c2 = oldpars['c2'] 
    275278        i_1 = 0.5*p_c1/p_c2 
    276279        i_2 = math.sqrt(math.fabs(p_scale/p_c2)) 
     
    295298    return oldpars 
    296299 
    297 def convert_model(name, pars, use_underscore=False, model_version=(3,1,2)): 
     300def convert_model(name, pars, use_underscore=False, model_version=(3, 1, 2)): 
    298301    """ 
    299302    Convert model from old style parameter names to new style. 
     
    327330        newpars = _convert_pars(newpars, translation) 
    328331        # TODO: Still not convinced this is the best check 
    329         if not model_info.structure_factor and version == (3,1,2): 
     332        if not model_info.structure_factor and version == (3, 1, 2): 
    330333            newpars = _rescale_sld(model_info, newpars, 1e6) 
    331334        newpars.setdefault('scale', 1.0) 
     
    599602 
    600603    pars = compare.get_pars(model_info, use_demo=False) 
    601     pars = compare.randomize_pars(model_info, pars, seed=seed) 
     604    if seed is not None: 
     605        np.random.seed(seed) 
     606    pars = compare.randomize_pars(model_info, pars) 
    602607    if name == "teubner_strey": 
    603608        # T-S model is underconstrained, so fix the assumptions. 
     
    615620        print("==== %s out ====="%new_name) 
    616621        print(str(compare.parlist(model_info, new_pars, True))) 
    617     assert name==new_name, "%r != %r"%(name, new_name) 
     622    assert name == new_name, "%r != %r"%(name, new_name) 
    618623    for k, v in new_pars.items(): 
    619624        assert k in pars, "%s: %r appeared from conversion"%(name, k) 
    620625        if isinstance(v, float): 
    621             assert abs(v-pars[k])<=abs(1e-12*v), "%s: %r  %s != %s"%(name, k, v, pars[k]) 
     626            assert abs(v-pars[k]) <= abs(1e-12*v), \ 
     627                "%s: %r  %s != %s"%(name, k, v, pars[k]) 
    622628        else: 
    623629            assert v == pars[k], "%s: %r  %s != %s"%(name, k, v, pars[k]) 
  • sasmodels/core.py

    ra85a569 r2d81cfe  
    1010 
    1111import os 
     12from os.path import basename, join as joinpath 
     13from glob import glob 
    1214import re 
    13 from os.path import basename, dirname, join as joinpath 
    14 from glob import glob 
    1515 
    1616import numpy as np # type: ignore 
     
    3535CUSTOM_MODEL_PATH = os.environ.get('SAS_MODELPATH', "") 
    3636if CUSTOM_MODEL_PATH == "": 
    37     path = joinpath(os.path.expanduser("~"), ".sasmodels", "custom_models") 
    38     if not os.path.isdir(path): 
    39         os.makedirs(path) 
    40     CUSTOM_MODEL_PATH = path 
    41  
     37    CUSTOM_MODEL_PATH = joinpath(os.path.expanduser("~"), ".sasmodels", "custom_models") 
     38    if not os.path.isdir(CUSTOM_MODEL_PATH): 
     39        os.makedirs(CUSTOM_MODEL_PATH) 
     40 
     41# pylint: disable=unused-import 
    4242try: 
    4343    from typing import List, Union, Optional, Any 
     
    4646except ImportError: 
    4747    pass 
     48# pylint: enable=unused-import 
    4849 
    4950# TODO: refactor composite model support 
     
    272273    Possible types include 'half', 'single', 'double' and 'quad'.  If the 
    273274    type is 'fast', then this is equivalent to dtype 'single' but using 
    274     fast native functions rather than those with the precision level guaranteed 
    275     by the OpenCL standard. 
     275    fast native functions rather than those with the precision level 
     276    guaranteed by the OpenCL standard.  'default' will choose the appropriate 
     277    default for the model and platform. 
    276278 
    277279    Platform preference can be specfied ("ocl" vs "dll"), with the default 
  • sasmodels/data.py

    rba7302a r2d81cfe  
    3737import numpy as np  # type: ignore 
    3838 
     39# pylint: disable=unused-import 
    3940try: 
    4041    from typing import Union, Dict, List, Optional 
     
    4344else: 
    4445    Data = Union["Data1D", "Data2D", "SesansData"] 
     46# pylint: enable=unused-import 
    4547 
    4648def load_data(filename, index=0): 
     
    6466            data.qmin, data.qmax = data.x.min(), data.x.max() 
    6567            data.mask = (np.isnan(data.y) if data.y is not None 
    66                         else np.zeros_like(data.x, dtype='bool')) 
     68                         else np.zeros_like(data.x, dtype='bool')) 
    6769        elif hasattr(data, 'qx_data'): 
    6870            data.mask = ~data.mask 
     
    137139    *_yaxis*, *_yunit*: label and units for the *y* axis 
    138140    """ 
    139     def __init__(self, x=None, y=None, dx=None, dy=None): 
    140         # type: (Optional[np.ndarray], Optional[np.ndarray], Optional[np.ndarray], Optional[np.ndarray]) -> None 
     141    def __init__(self, 
     142                 x=None,  # type: Optional[np.ndarray] 
     143                 y=None,  # type: Optional[np.ndarray] 
     144                 dx=None, # type: Optional[np.ndarray] 
     145                 dy=None  # type: Optional[np.ndarray] 
     146                ): 
     147        # type: (...) -> None 
    141148        self.x, self.y, self.dx, self.dy = x, y, dx, dy 
    142149        self.dxl = None 
     
    211218    *x_bins*, *y_bins*: grid steps in *x* and *y* directions 
    212219    """ 
    213     def __init__(self, x=None, y=None, z=None, dx=None, dy=None, dz=None): 
    214         # type: (Optional[np.ndarray], Optional[np.ndarray], Optional[np.ndarray], Optional[np.ndarray], Optional[np.ndarray], Optional[np.ndarray]) -> None 
     220    def __init__(self, 
     221                 x=None,   # type: Optional[np.ndarray] 
     222                 y=None,   # type: Optional[np.ndarray] 
     223                 z=None,   # type: Optional[np.ndarray] 
     224                 dx=None,  # type: Optional[np.ndarray] 
     225                 dy=None,  # type: Optional[np.ndarray] 
     226                 dz=None   # type: Optional[np.ndarray] 
     227                ): 
     228        # type: (...) -> None 
    215229        self.qx_data, self.dqx_data = x, dx 
    216230        self.qy_data, self.dqy_data = y, dy 
     
    363377    if hasattr(data, 'isSesans') and data.isSesans: 
    364378        _plot_result_sesans(data, None, None, use_data=True, limits=limits) 
    365     elif hasattr(data, 'qx_data'): 
     379    elif hasattr(data, 'qx_data') and not getattr(data, 'radial', False): 
    366380        _plot_result2D(data, None, None, view, use_data=True, limits=limits) 
    367381    else: 
     
    369383 
    370384 
    371 def plot_theory(data, theory, resid=None, view='log', 
    372                 use_data=True, limits=None, Iq_calc=None): 
    373     # type: (Data, Optional[np.ndarray], Optional[np.ndarray], str, bool, Optional[Tuple[float,float]], Optional[np.ndarray]) -> None 
     385def plot_theory(data,          # type: Data 
     386                theory,        # type: Optional[np.ndarray] 
     387                resid=None,    # type: Optional[np.ndarray] 
     388                view='log',    # type: str 
     389                use_data=True, # type: bool 
     390                limits=None,   # type: Optional[np.ndarray] 
     391                Iq_calc=None   # type: Optional[np.ndarray] 
     392               ): 
     393    # type: (...) -> None 
    374394    """ 
    375395    Plot theory calculation. 
     
    391411    if hasattr(data, 'isSesans') and data.isSesans: 
    392412        _plot_result_sesans(data, theory, resid, use_data=True, limits=limits) 
    393     elif hasattr(data, 'qx_data'): 
     413    elif hasattr(data, 'qx_data') and not getattr(data, 'radial', False): 
    394414        _plot_result2D(data, theory, resid, view, use_data, limits=limits) 
    395415    else: 
     
    417437 
    418438@protect 
    419 def _plot_result1D(data, theory, resid, view, use_data, 
    420                    limits=None, Iq_calc=None): 
    421     # type: (Data1D, Optional[np.ndarray], Optional[np.ndarray], str, bool, Optional[Tuple[float, float]], Optional[np.ndarray]) -> None 
     439def _plot_result1D(data,         # type: Data1D 
     440                   theory,       # type: Optional[np.ndarray] 
     441                   resid,        # type: Optional[np.ndarray] 
     442                   view,         # type: str 
     443                   use_data,     # type: bool 
     444                   limits=None,  # type: Optional[Tuple[float, float]] 
     445                   Iq_calc=None  # type: Optional[np.ndarray] 
     446                  ): 
     447    # type: (...) -> None 
    422448    """ 
    423449    Plot the data and residuals for 1D data. 
     
    425451    import matplotlib.pyplot as plt  # type: ignore 
    426452    from numpy.ma import masked_array, masked  # type: ignore 
     453 
     454    if getattr(data, 'radial', False): 
     455        data.x = data.q_data 
     456        data.y = data.data 
    427457 
    428458    use_data = use_data and data.y is not None 
     
    511541 
    512542@protect 
    513 def _plot_result_sesans(data, theory, resid, use_data, limits=None): 
    514     # type: (SesansData, Optional[np.ndarray], Optional[np.ndarray], bool, Optional[Tuple[float, float]]) -> None 
     543def _plot_result_sesans(data,        # type: SesansData 
     544                        theory,      # type: Optional[np.ndarray] 
     545                        resid,       # type: Optional[np.ndarray] 
     546                        use_data,    # type: bool 
     547                        limits=None  # type: Optional[Tuple[float, float]] 
     548                       ): 
     549    # type: (...) -> None 
    515550    """ 
    516551    Plot SESANS results. 
     
    523558 
    524559    if use_data or use_theory: 
    525         is_tof = (data.lam != data.lam[0]).any() 
     560        is_tof = data.lam is not None and (data.lam != data.lam[0]).any() 
    526561        if num_plots > 1: 
    527562            plt.subplot(1, num_plots, 1) 
     
    556591 
    557592@protect 
    558 def _plot_result2D(data, theory, resid, view, use_data, limits=None): 
    559     # type: (Data2D, Optional[np.ndarray], Optional[np.ndarray], str, bool, Optional[Tuple[float,float]]) -> None 
     593def _plot_result2D(data,         # type: Data2D 
     594                   theory,       # type: Optional[np.ndarray] 
     595                   resid,        # type: Optional[np.ndarray] 
     596                   view,         # type: str 
     597                   use_data,     # type: bool 
     598                   limits=None   # type: Optional[Tuple[float, float]] 
     599                  ): 
     600    # type: (...) -> None 
    560601    """ 
    561602    Plot the data and residuals for 2D data. 
     
    617658 
    618659@protect 
    619 def _plot_2d_signal(data, signal, vmin=None, vmax=None, view='log'): 
    620     # type: (Data2D, np.ndarray, Optional[float], Optional[float], str) -> Tuple[float, float] 
     660def _plot_2d_signal(data,       # type: Data2D 
     661                    signal,     # type: np.ndarray 
     662                    vmin=None,  # type: Optional[float] 
     663                    vmax=None,  # type: Optional[float] 
     664                    view='log'  # type: str 
     665                   ): 
     666    # type: (...) -> Tuple[float, float] 
    621667    """ 
    622668    Plot the target value for the data.  This could be the data itself, 
     
    633679    if view == 'log': 
    634680        valid[valid] = (image[valid] > 0) 
    635         if vmin is None: vmin = image[valid & ~data.mask].min() 
    636         if vmax is None: vmax = image[valid & ~data.mask].max() 
     681        if vmin is None: 
     682            vmin = image[valid & ~data.mask].min() 
     683        if vmax is None: 
     684            vmax = image[valid & ~data.mask].max() 
    637685        image[valid] = np.log10(image[valid]) 
    638686    elif view == 'q4': 
    639687        image[valid] *= (data.qx_data[valid]**2+data.qy_data[valid]**2)**2 
    640         if vmin is None: vmin = image[valid & ~data.mask].min() 
    641         if vmax is None: vmax = image[valid & ~data.mask].max() 
     688        if vmin is None: 
     689            vmin = image[valid & ~data.mask].min() 
     690        if vmax is None: 
     691            vmax = image[valid & ~data.mask].max() 
    642692    else: 
    643         if vmin is None: vmin = image[valid & ~data.mask].min() 
    644         if vmax is None: vmax = image[valid & ~data.mask].max() 
     693        if vmin is None: 
     694            vmin = image[valid & ~data.mask].min() 
     695        if vmax is None: 
     696            vmax = image[valid & ~data.mask].max() 
    645697 
    646698    image[~valid | data.mask] = 0 
     
    651703    ymin, ymax = min(data.qy_data), max(data.qy_data) 
    652704    if view == 'log': 
    653         vmin, vmax = np.log10(vmin), np.log10(vmax) 
     705        vmin_scaled, vmax_scaled = np.log10(vmin), np.log10(vmax) 
     706    else: 
     707        vmin_scaled, vmax_scaled = vmin, vmax 
    654708    plt.imshow(plottable.reshape(len(data.x_bins), len(data.y_bins)), 
    655709               interpolation='nearest', aspect=1, origin='lower', 
    656                extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax) 
     710               extent=[xmin, xmax, ymin, ymax], 
     711               vmin=vmin_scaled, vmax=vmax_scaled) 
    657712    plt.xlabel("$q_x$/A$^{-1}$") 
    658713    plt.ylabel("$q_y$/A$^{-1}$") 
  • sasmodels/details.py

    rccd5f01 r2d81cfe  
    1515 
    1616import numpy as np  # type: ignore 
    17 from numpy import pi, cos, sin 
     17from numpy import cos, sin, radians 
    1818 
    1919try: 
     
    2323    # CRUFT: np.meshgrid requires multiple vectors 
    2424    def meshgrid(*args): 
     25        """See docs from a recent version of numpy""" 
    2526        if len(args) > 1: 
    2627            return np.meshgrid(*args) 
     
    2829            return [np.asarray(v) for v in args] 
    2930 
     31# pylint: disable=unused-import 
    3032try: 
    31     from typing import List 
     33    from typing import List, Tuple, Sequence 
    3234except ImportError: 
    3335    pass 
    3436else: 
    3537    from .modelinfo import ModelInfo 
     38    from .modelinfo import ParameterTable 
     39# pylint: enable=unused-import 
    3640 
    3741 
     
    5357    coordinates, the total circumference decreases as latitude varies from 
    5458    pi r^2 at the equator to 0 at the pole, and the weight associated 
    55     with a range of phi values needs to be scaled by this circumference. 
     59    with a range of latitude values needs to be scaled by this circumference. 
    5660    This scale factor needs to be updated each time the theta value 
    5761    changes.  *theta_par* indicates which of the values in the parameter 
     
    192196    #print("off:",offset) 
    193197 
    194     # Check that we arn't using too many polydispersity loops 
     198    # Check that we aren't using too many polydispersity loops 
    195199    num_active = np.sum(length > 1) 
    196200    max_pd = model_info.parameters.max_pd 
     
    218222 
    219223ZEROS = tuple([0.]*31) 
    220 def make_kernel_args(kernel, pairs): 
    221     # type: (Kernel, Tuple[List[np.ndarray], List[np.ndarray]]) -> Tuple[CallDetails, np.ndarray, bool] 
    222     """ 
    223     Converts (value, weight) pairs into parameters for the kernel call. 
     224def make_kernel_args(kernel, # type: Kernel 
     225                     mesh    # type: Tuple[List[np.ndarray], List[np.ndarray]] 
     226                    ): 
     227    # type: (...) -> Tuple[CallDetails, np.ndarray, bool] 
     228    """ 
     229    Converts (value, dispersity, weight) for each parameter into kernel pars. 
    224230 
    225231    Returns a CallDetails object indicating the polydispersity, a data object 
     
    230236    npars = kernel.info.parameters.npars 
    231237    nvalues = kernel.info.parameters.nvalues 
    232     scalars = [(v[0] if len(v) else np.NaN) for v, w in pairs] 
    233     values, weights = zip(*pairs[2:npars+2]) if npars else ((),()) 
     238    scalars = [value for value, _dispersity, _weight in mesh] 
     239    # skipping scale and background when building values and weights 
     240    _values, dispersity, weights = zip(*mesh[2:npars+2]) if npars else ((), (), ()) 
     241    #weights = correct_theta_weights(kernel.info.parameters, dispersity, weights) 
    234242    length = np.array([len(w) for w in weights]) 
    235243    offset = np.cumsum(np.hstack((0, length))) 
    236244    call_details = make_details(kernel.info, length, offset[:-1], offset[-1]) 
    237245    # Pad value array to a 32 value boundaryd 
    238     data_len = nvalues + 2*sum(len(v) for v in values) 
     246    data_len = nvalues + 2*sum(len(v) for v in dispersity) 
    239247    extra = (32 - data_len%32)%32 
    240     data = np.hstack((scalars,) + values + weights + ZEROS[:extra]) 
     248    data = np.hstack((scalars,) + dispersity + weights + ZEROS[:extra]) 
    241249    data = data.astype(kernel.dtype) 
    242250    is_magnetic = convert_magnetism(kernel.info.parameters, data) 
     
    244252    return call_details, data, is_magnetic 
    245253 
     254def correct_theta_weights(parameters, # type: ParameterTable 
     255                          dispersity, # type: Sequence[np.ndarray] 
     256                          weights     # type: Sequence[np.ndarray] 
     257                         ): 
     258    # type: (...) -> Sequence[np.ndarray] 
     259    """ 
     260    If there is a theta parameter, update the weights of that parameter so that 
     261    the cosine weighting required for polar integration is preserved. 
     262 
     263    Avoid evaluation strictly at the pole, which would otherwise send the 
     264    weight to zero.  This is probably not a problem in practice (if dispersity 
     265    is +/- 90, then you probably should be using a 1-D model of the circular 
     266    average). 
     267 
     268    Note: scale and background parameters are not include in the tuples for 
     269    dispersity and weights, so index is parameters.theta_offset, not 
     270    parameters.theta_offset+2 
     271 
     272    Returns updated weights vectors 
     273    """ 
     274    # TODO: explain in a comment why scale and background are missing 
     275    # Apparently the parameters.theta_offset similarly skips scale and 
     276    # and background, so the indexing works out, but they are still shipped 
     277    # to the kernel, so we need to add two there. 
     278    if parameters.theta_offset >= 0: 
     279        index = parameters.theta_offset 
     280        theta = dispersity[index] 
     281        # TODO: modify the dispersity vector to avoid the theta=-90,90,270,... 
     282        theta_weight = abs(cos(radians(theta))) 
     283        weights = tuple(theta_weight*w if k == index else w 
     284                        for k, w in enumerate(weights)) 
     285    return weights 
     286 
    246287 
    247288def convert_magnetism(parameters, values): 
     289    # type: (ParameterTable, Sequence[np.ndarray]) -> bool 
    248290    """ 
    249291    Convert magnetism values from polar to rectangular coordinates. 
     
    253295    mag = values[parameters.nvalues-3*parameters.nmagnetic:parameters.nvalues] 
    254296    mag = mag.reshape(-1, 3) 
    255     scale = mag[:,0] 
     297    scale = mag[:, 0] 
    256298    if np.any(scale): 
    257         theta, phi = mag[:, 1]*pi/180., mag[:, 2]*pi/180. 
     299        theta, phi = radians(mag[:, 1]), radians(mag[:, 2]) 
    258300        cos_theta = cos(theta) 
    259301        mag[:, 0] = scale*cos_theta*cos(phi)  # mx 
     
    265307 
    266308 
    267 def dispersion_mesh(model_info, pars): 
     309def dispersion_mesh(model_info, mesh): 
    268310    # type: (ModelInfo) -> Tuple[List[np.ndarray], List[np.ndarray]] 
    269311    """ 
    270312    Create a mesh grid of dispersion parameters and weights. 
     313 
     314    *mesh* is a list of (value, dispersity, weights), where the values 
     315    are the individual parameter values, and (dispersity, weights) is 
     316    the distribution of parameter values. 
     317 
     318    Only the volume parameters should be included in this list.  Orientation 
     319    parameters do not affect the calculation of effective radius or volume 
     320    ratio.  This is convenient since it avoids the distinction between 
     321    value and dispersity that is present in orientation parameters but not 
     322    shape parameters. 
    271323 
    272324    Returns [p1,p2,...],w where pj is a vector of values for parameter j 
     
    274326    parameter set in the vector. 
    275327    """ 
    276     value, weight = zip(*pars) 
     328    _, dispersity, weight = zip(*mesh) 
    277329    #weight = [w if len(w)>0 else [1.] for w in weight] 
    278330    weight = np.vstack([v.flatten() for v in meshgrid(*weight)]) 
    279331    weight = np.prod(weight, axis=0) 
    280     value = [v.flatten() for v in meshgrid(*value)] 
     332    dispersity = [v.flatten() for v in meshgrid(*dispersity)] 
    281333    lengths = [par.length for par in model_info.parameters.kernel_parameters 
    282334               if par.type == 'volume'] 
     
    285337        offset = 0 
    286338        for n in lengths: 
    287             pars.append(np.vstack(value[offset:offset+n]) 
    288                         if n > 1 else value[offset]) 
     339            pars.append(np.vstack(dispersity[offset:offset+n]) 
     340                        if n > 1 else dispersity[offset]) 
    289341            offset += n 
    290         value = pars 
    291     return value, weight 
     342        dispersity = pars 
     343    return dispersity, weight 
  • sasmodels/direct_model.py

    rd1ff3a5 r2d81cfe  
    3232from .details import make_kernel_args, dispersion_mesh 
    3333 
     34# pylint: disable=unused-import 
    3435try: 
    3536    from typing import Optional, Dict, Tuple 
     
    4041    from .kernel import Kernel, KernelModel 
    4142    from .modelinfo import Parameter, ParameterSet 
     43# pylint: enable=unused-import 
    4244 
    4345def call_kernel(calculator, pars, cutoff=0., mono=False): 
     
    5557    *mono* is True if polydispersity should be set to none on all parameters. 
    5658    """ 
    57     parameters = calculator.info.parameters 
    58     if mono: 
    59         active = lambda name: False 
    60     elif calculator.dim == '1d': 
    61         active = lambda name: name in parameters.pd_1d 
    62     elif calculator.dim == '2d': 
    63         active = lambda name: name in parameters.pd_2d 
    64     else: 
    65         active = lambda name: True 
    66  
    67     #print("pars",[p.id for p in parameters.call_parameters]) 
    68     vw_pairs = [(get_weights(p, pars) if active(p.name) 
    69                  else ([pars.get(p.name, p.default)], [1.0])) 
    70                 for p in parameters.call_parameters] 
    71  
    72     call_details, values, is_magnetic = make_kernel_args(calculator, vw_pairs) 
     59    mesh = get_mesh(calculator.info, pars, dim=calculator.dim, mono=mono) 
     60    #print("pars", list(zip(*mesh))[0]) 
     61    call_details, values, is_magnetic = make_kernel_args(calculator, mesh) 
    7362    #print("values:", values) 
    7463    return calculator(call_details, values, cutoff, is_magnetic) 
    75  
    7664 
    7765def call_ER(model_info, pars): 
     
    129117    return x, y, model_info.profile_axes 
    130118 
    131  
    132 def get_weights(parameter, values): 
    133     # type: (Parameter, Dict[str, float]) -> Tuple[np.ndarray, np.ndarray] 
     119def get_mesh(model_info, values, dim='1d', mono=False): 
     120    # type: (ModelInfo, Dict[str, float], str, bool) -> List[Tuple[float, np.ndarray, np.ndarry]] 
     121    """ 
     122    Retrieve the dispersity mesh described by the parameter set. 
     123 
     124    Returns a list of *(value, dispersity, weights)* with one tuple for each 
     125    parameter in the model call parameters.  Inactive parameters return the 
     126    default value with a weight of 1.0. 
     127    """ 
     128    parameters = model_info.parameters 
     129    if mono: 
     130        active = lambda name: False 
     131    elif dim == '1d': 
     132        active = lambda name: name in parameters.pd_1d 
     133    elif dim == '2d': 
     134        active = lambda name: name in parameters.pd_2d 
     135    else: 
     136        active = lambda name: True 
     137 
     138    #print("pars",[p.id for p in parameters.call_parameters]) 
     139    mesh = [_get_par_weights(p, values, active(p.name)) 
     140            for p in parameters.call_parameters] 
     141    return mesh 
     142 
     143 
     144def _get_par_weights(parameter, values, active=True): 
     145    # type: (Parameter, Dict[str, float]) -> Tuple[float, np.ndarray, np.ndarray] 
    134146    """ 
    135147    Generate the distribution for parameter *name* given the parameter values 
     
    140152    """ 
    141153    value = float(values.get(parameter.name, parameter.default)) 
    142     relative = parameter.relative_pd 
    143     limits = parameter.limits 
    144     disperser = values.get(parameter.name+'_pd_type', 'gaussian') 
    145154    npts = values.get(parameter.name+'_pd_n', 0) 
    146155    width = values.get(parameter.name+'_pd', 0.0) 
    147     nsigma = values.get(parameter.name+'_pd_nsigma', 3.0) 
    148     if npts == 0 or width == 0: 
    149         return [value], [1.0] 
    150     value, weight = weights.get_weights( 
    151         disperser, npts, width, nsigma, value, limits, relative) 
    152     return value, weight / np.sum(weight) 
    153  
    154  
    155 def _vol_pars(model_info, pars): 
     156    relative = parameter.relative_pd 
     157    if npts == 0 or width == 0.0 or not active: 
     158        # Note: orientation parameters have the viewing angle as the parameter 
     159        # value and the jitter in the distribution, so be sure to set the 
     160        # empty pd for orientation parameters to 0. 
     161        pd = [value if relative or not parameter.polydisperse else 0.0], [1.0] 
     162    else: 
     163        limits = parameter.limits 
     164        disperser = values.get(parameter.name+'_pd_type', 'gaussian') 
     165        nsigma = values.get(parameter.name+'_pd_nsigma', 3.0) 
     166        pd = weights.get_weights(disperser, npts, width, nsigma, 
     167                                 value, limits, relative) 
     168    return value, pd[0], pd[1] 
     169 
     170 
     171def _vol_pars(model_info, values): 
    156172    # type: (ModelInfo, ParameterSet) -> Tuple[np.ndarray, np.ndarray] 
    157     vol_pars = [get_weights(p, pars) 
     173    vol_pars = [_get_par_weights(p, values) 
    158174                for p in model_info.parameters.call_parameters 
    159175                if p.type == 'volume'] 
    160176    #import pylab; pylab.plot(vol_pars[0][0],vol_pars[0][1]); pylab.show() 
    161     value, weight = dispersion_mesh(model_info, vol_pars) 
    162     return value, weight 
     177    dispersity, weight = dispersion_mesh(model_info, vol_pars) 
     178    return dispersity, weight 
     179 
     180 
     181def _make_sesans_transform(data): 
     182    from sas.sascalc.data_util.nxsunit import Converter 
     183 
     184    # Pre-compute the Hankel matrix (H) 
     185    SElength = Converter(data._xunit)(data.x, "A") 
     186 
     187    theta_max = Converter("radians")(data.sample.zacceptance)[0] 
     188    q_max = 2 * np.pi / np.max(data.source.wavelength) * np.sin(theta_max) 
     189    zaccept = Converter("1/A")(q_max, "1/" + data.source.wavelength_unit), 
     190 
     191    Rmax = 10000000 
     192    hankel = sesans.SesansTransform(data.x, SElength, 
     193                                    data.source.wavelength, 
     194                                    zaccept, Rmax) 
     195    return hankel 
    163196 
    164197 
     
    202235 
    203236        if self.data_type == 'sesans': 
    204             q = sesans.make_q(data.sample.zacceptance, data.Rmax) 
     237            res = _make_sesans_transform(data) 
    205238            index = slice(None, None) 
    206             res = None 
    207239            if data.y is not None: 
    208240                Iq, dIq = data.y, data.dy 
     
    210242                Iq, dIq = None, None 
    211243            #self._theory = np.zeros_like(q) 
    212             q_vectors = [q] 
    213             q_mono = sesans.make_all_q(data) 
     244            q_vectors = [res.q_calc] 
    214245        elif self.data_type == 'Iqxy': 
    215246            #if not model.info.parameters.has_2d: 
     
    230261            #self._theory = np.zeros_like(self.Iq) 
    231262            q_vectors = res.q_calc 
    232             q_mono = [] 
    233263        elif self.data_type == 'Iq': 
    234264            index = (data.x >= data.qmin) & (data.x <= data.qmax) 
     
    255285            #self._theory = np.zeros_like(self.Iq) 
    256286            q_vectors = [res.q_calc] 
    257             q_mono = [] 
    258287        elif self.data_type == 'Iq-oriented': 
    259288            index = (data.x >= data.qmin) & (data.x <= data.qmax) 
     
    272301                                      qy_width=data.dxl[index]) 
    273302            q_vectors = res.q_calc 
    274             q_mono = [] 
    275303        else: 
    276304            raise ValueError("Unknown data type") # never gets here 
     
    279307        # so we can save/restore state 
    280308        self._kernel_inputs = q_vectors 
    281         self._kernel_mono_inputs = q_mono 
    282309        self._kernel = None 
    283310        self.Iq, self.dIq, self.index = Iq, dIq, index 
     
    317344        if self._kernel is None: 
    318345            self._kernel = self._model.make_kernel(self._kernel_inputs) 
    319             self._kernel_mono = ( 
    320                 self._model.make_kernel(self._kernel_mono_inputs) 
    321                 if self._kernel_mono_inputs else None) 
    322346 
    323347        Iq_calc = call_kernel(self._kernel, pars, cutoff=cutoff) 
     
    327351        # TODO: refactor so we don't store the result in the model 
    328352        self.Iq_calc = Iq_calc 
    329         if self.data_type == 'sesans': 
    330             Iq_mono = (call_kernel(self._kernel_mono, pars, mono=True) 
    331                        if self._kernel_mono_inputs else None) 
    332             result = sesans.transform(self._data, 
    333                                       self._kernel_inputs[0], Iq_calc, 
    334                                       self._kernel_mono_inputs, Iq_mono) 
    335         else: 
    336             result = self.resolution.apply(Iq_calc) 
    337             if hasattr(self.resolution, 'nx'): 
    338                 self.Iq_calc = ( 
    339                     self.resolution.qx_calc, self.resolution.qy_calc, 
    340                     np.reshape(Iq_calc, (self.resolution.ny, self.resolution.nx)) 
    341                 ) 
     353        result = self.resolution.apply(Iq_calc) 
     354        if hasattr(self.resolution, 'nx'): 
     355            self.Iq_calc = ( 
     356                self.resolution.qx_calc, self.resolution.qy_calc, 
     357                np.reshape(Iq_calc, (self.resolution.ny, self.resolution.nx)) 
     358            ) 
    342359        return result 
    343360 
     
    396413        try: 
    397414            values = [float(v) for v in call.split(',')] 
    398         except Exception: 
     415        except ValueError: 
    399416            values = [] 
    400417        if len(values) == 1: 
  • sasmodels/generate.py

    r6db17bd r2d81cfe  
    163163 
    164164import sys 
    165 from os.path import abspath, dirname, join as joinpath, exists, isdir, getmtime 
     165from os.path import abspath, dirname, join as joinpath, exists, getmtime 
    166166import re 
    167167import string 
    168168from zlib import crc32 
     169from inspect import currentframe, getframeinfo 
    169170 
    170171import numpy as np  # type: ignore 
     
    173174from .custom import load_custom_kernel_module 
    174175 
     176# pylint: disable=unused-import 
    175177try: 
    176178    from typing import Tuple, Sequence, Iterator, Dict 
     
    178180except ImportError: 
    179181    pass 
     182# pylint: enable=unused-import 
     183 
     184# jitter projection to use in the kernel code.  See explore/jitter.py 
     185# for details.  To change it from a program, set generate.PROJECTION. 
     186PROJECTION = 1 
    180187 
    181188def get_data_path(external_dir, target_file): 
     
    649656    line instead. 
    650657    """ 
    651     for path, code in sources: 
     658    for _path, code in sources: 
    652659        if _IQXY_PATTERN.search(code): 
    653660            return True 
    654     else: 
    655         return False 
     661    return False 
    656662 
    657663 
     
    689695    # Load templates and user code 
    690696    kernel_header = load_template('kernel_header.c') 
    691     dll_code = load_template('kernel_iq.c') 
    692     ocl_code = load_template('kernel_iq.cl') 
    693     #ocl_code = load_template('kernel_iq_local.cl') 
     697    kernel_code = load_template('kernel_iq.c') 
    694698    user_code = [(f, open(f).read()) for f in model_sources(model_info)] 
     699 
     700    # What kind of 2D model do we need? 
     701    xy_mode = ('qa' if not _have_Iqxy(user_code) and not isinstance(model_info.Iqxy, str) 
     702               else 'qac' if not partable.is_asymmetric 
     703               else 'qabc') 
    695704 
    696705    # Build initial sources 
     
    717726 
    718727    # Define the parameter table 
    719     # TODO: plug in current line number 
    720     source.append('#line 542 "sasmodels/generate.py"') 
     728    lineno = getframeinfo(currentframe()).lineno + 2 
     729    source.append('#line %d "sasmodels/generate.py"'%lineno) 
     730    #source.append('introduce breakage in generate to test lineno reporting') 
    721731    source.append("#define PARAMETER_TABLE \\") 
    722732    source.append("\\\n".join(p.as_definition() 
     
    734744    source.append(call_volume) 
    735745 
    736     refs = ["_q[_i]"] + _call_pars("_v.", partable.iq_parameters) 
    737     call_iq = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(refs)) 
    738     if _have_Iqxy(user_code) or isinstance(model_info.Iqxy, str): 
    739         # Call 2D model 
    740         refs = ["_q[2*_i]", "_q[2*_i+1]"] + _call_pars("_v.", partable.iqxy_parameters) 
    741         call_iqxy = "#define CALL_IQ(_q,_i,_v) Iqxy(%s)" % (",".join(refs)) 
    742     else: 
    743         # Call 1D model with sqrt(qx^2 + qy^2) 
    744         #warnings.warn("Creating Iqxy = Iq(sqrt(qx^2 + qy^2))") 
    745         # still defined:: refs = ["q[i]"] + _call_pars("v", iq_parameters) 
    746         pars_sqrt = ["sqrt(_q[2*_i]*_q[2*_i]+_q[2*_i+1]*_q[2*_i+1])"] + refs[1:] 
    747         call_iqxy = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(pars_sqrt)) 
     746    model_refs = _call_pars("_v.", partable.iq_parameters) 
     747    pars = ",".join(["_q"] + model_refs) 
     748    call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars 
     749    if xy_mode == 'qabc': 
     750        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs) 
     751        call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqxy(%s)" % pars 
     752        clear_iqxy = "#undef CALL_IQ_ABC" 
     753    elif xy_mode == 'qac': 
     754        pars = ",".join(["_qa", "_qc"] + model_refs) 
     755        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqxy(%s)" % pars