Changeset fd7291e in sasmodels


Ignore:
Timestamp:
Mar 6, 2019 6:05:52 PM (5 months 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:
c11d09f
Parents:
21c93c3 (diff), 9150036 (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 'beta_approx' into ticket-1022-sum_multiplicity

Files:
50 added
115 edited

Legend:

Unmodified
Added
Removed
  • README.rst

    re30d645 r2a64722  
    1010is available. 
    1111 
    12 Example 
     12Install 
    1313------- 
     14 
     15The easiest way to use sasmodels is from `SasView <http://www.sasview.org/>`_. 
     16 
     17You can also install sasmodels as a standalone package in python. Use 
     18`miniconda <https://docs.conda.io/en/latest/miniconda.html>`_ 
     19or `anaconda <https://www.anaconda.com/>`_ 
     20to create a python environment with the sasmodels dependencies:: 
     21 
     22    $ conda create -n sasmodels -c conda-forge numpy scipy matplotlib pyopencl 
     23 
     24The option ``-n sasmodels`` names the environment sasmodels, and the option 
     25``-c conda-forge`` selects the conda-forge package channel because pyopencl 
     26is not part of the base anaconda distribution. 
     27 
     28Activate the environment and install sasmodels:: 
     29 
     30    $ conda activate sasmodels 
     31    (sasmodels) $ pip install sasmodels 
     32 
     33Install `bumps <https://github.com/bumps/bumps>`_ if you want to use it to fit 
     34your data:: 
     35 
     36    (sasmodels) $ pip install bumps 
     37 
     38Usage 
     39----- 
     40 
     41Check that the works:: 
     42 
     43    (sasmodels) $ python -m sasmodels.compare cylinder 
     44 
     45To show the orientation explorer:: 
     46 
     47    (sasmodels) $ python -m sasmodels.jitter 
     48 
     49Documentation is available online as part of the SasView 
     50`fitting perspective <http://www.sasview.org/docs/index.html>`_ 
     51as well as separate pages for 
     52`individual models <http://www.sasview.org/docs/user/sasgui/perspectives/fitting/models/index.html>`_. 
     53Programming details for sasmodels are available in the 
     54`developer documentation <http://www.sasview.org/docs/dev/dev.html>`_. 
     55 
     56 
     57Fitting Example 
     58--------------- 
    1459 
    1560The example directory contains a radial+tangential data set for an oriented 
    1661rod-like shape. 
    1762 
    18 The data is loaded by sas.dataloader from the sasview package, so sasview 
    19 is needed to run the example. 
     63To load the example data, you will need the SAS data loader from the sasview 
     64package. This is not yet available on PyPI, so you will need a copy of the 
     65SasView source code to run it.  Create a directory somewhere to hold the 
     66sasview and sasmodels source code, which we will refer to as $SOURCE. 
    2067 
    21 To run the example, you need sasview, sasmodels and bumps.  Assuming these 
    22 repositories are installed side by side, change to the sasmodels/example 
    23 directory and enter:: 
     68Use the following to install sasview, and the sasmodels examples:: 
    2469 
    25     PYTHONPATH=..:../../sasview/src ../../bumps/run.py fit.py \ 
    26         cylinder --preview 
     70    (sasmodels) $ cd $SOURCE 
     71    (sasmodels) $ conda install git 
     72    (sasmodels) $ git clone https://github.com/sasview/sasview.git 
     73    (sasmodels) $ git clone https://github.com/sasview/sasmodels.git 
    2774 
    28 See bumps documentation for instructions on running the fit.  With the 
    29 python packages installed, e.g., into a virtual environment, then the 
    30 python path need not be set, and the command would be:: 
     75Set the path to the sasview source on your python path within the sasmodels 
     76environment.  On Windows, this will be:: 
    3177 
    32     bumps fit.py cylinder --preview 
     78    (sasmodels)> set PYTHONPATH="$SOURCE\sasview\src" 
     79    (sasmodels)> cd $SOURCE/sasmodels/example 
     80    (sasmodels)> python -m bumps.cli fit.py cylinder --preview 
     81 
     82On Mac/Linux with the standard shell this will be:: 
     83 
     84    (sasmodels) $ export PYTHONPATH="$SOURCE/sasview/src" 
     85    (sasmodels) $ cd $SOURCE/sasmodels/example 
     86    (sasmodels) $ bumps fit.py cylinder --preview 
    3387 
    3488The fit.py model accepts up to two arguments.  The first argument is the 
     
    3892both radial and tangential simultaneously, use the word "both". 
    3993 
    40 Notes 
    41 ----- 
    42  
    43 cylinder.c + cylinder.py is the cylinder model with renamed variables and 
    44 sld scaled by 1e6 so the numbers are nicer.  The model name is "cylinder" 
    45  
    46 lamellar.py is an example of a single file model with embedded C code. 
     94See `bumps documentation <https://bumps.readthedocs.io/>`_ for detailed 
     95instructions on running the fit. 
    4796 
    4897|TravisStatus|_ 
  • doc/guide/gpu_setup.rst

    r63602b1 r8b31efa  
    9494Device Selection 
    9595================ 
     96**OpenCL drivers** 
     97 
    9698If you have multiple GPU devices you can tell the program which device to use. 
    9799By default, the program looks for one GPU and one CPU device from available 
     
    104106was used to run the model. 
    105107 
    106 **If you don't want to use OpenCL, you can set** *SAS_OPENCL=None* 
    107 **in your environment settings, and it will only use normal programs.** 
    108  
    109 If you want to use one of the other devices, you can run the following 
     108If you want to use a specific driver and devices, you can run the following 
    110109from the python console:: 
    111110 
     
    115114This will provide a menu of different OpenCL drivers available. 
    116115When one is selected, it will say "set PYOPENCL_CTX=..." 
    117 Use that value as the value of *SAS_OPENCL*. 
     116Use that value as the value of *SAS_OPENCL=driver:device*. 
     117 
     118To use the default OpenCL device (rather than CUDA or None), 
     119set *SAS_OPENCL=opencl*. 
     120 
     121In batch queues, you may need to set *XDG_CACHE_HOME=~/.cache*  
     122(Linux only) to a different directory, depending on how the filesystem  
     123is configured.  You should also set *SAS_DLL_PATH* for CPU-only modules. 
     124 
     125    -DSAS_MODELPATH=path sets directory containing custom models 
     126    -DSAS_OPENCL=vendor:device|cuda:device|none sets the target GPU device 
     127    -DXDG_CACHE_HOME=~/.cache sets the pyopencl cache root (linux only) 
     128    -DSAS_COMPILER=tinycc|msvc|mingw|unix sets the DLL compiler 
     129    -DSAS_OPENMP=1 turns on OpenMP for the DLLs 
     130    -DSAS_DLL_PATH=path sets the path to the compiled modules 
     131 
     132 
     133**CUDA drivers** 
     134 
     135If OpenCL drivers are not available on your system, but NVidia CUDA 
     136drivers are available, then set *SAS_OPENCL=cuda* or 
     137*SAS_OPENCL=cuda:n* for a particular device number *n*.  If no device 
     138number is specified, then the CUDA drivers looks for look for 
     139*CUDA_DEVICE=n* or a file ~/.cuda-device containing n for the device number. 
     140 
     141In batch queues, the SLURM command *sbatch --gres=gpu:1 ...* will set 
     142*CUDA_VISIBLE_DEVICES=n*, which ought to set the correct device 
     143number for *SAS_OPENCL=cuda*.  If not, then set 
     144*CUDA_DEVICE=$CUDA_VISIBLE_DEVICES* within the batch script.  You may 
     145need to set the CUDA cache directory to a folder accessible across the 
     146cluster with *PYCUDA_CACHE_DIR* (or *PYCUDA_DISABLE_CACHE* to disable 
     147caching), and you may need to set environment specific compiler flags 
     148with *PYCUDA_DEFAULT_NVCC_FLAGS*.  You should also set *SAS_DLL_PATH*  
     149for CPU-only modules. 
     150 
     151**No GPU support** 
     152 
     153If you don't want to use OpenCL or CUDA, you can set *SAS_OPENCL=None* 
     154in your environment settings, and it will only use normal programs. 
     155 
     156In batch queues, you may need to set *SAS_DLL_PATH* to a directory 
     157accessible on the compute node. 
     158 
    118159 
    119160Device Testing 
     
    154195*Document History* 
    155196 
    156 | 2017-09-27 Paul Kienzle 
     197| 2018-10-15 Paul Kienzle 
  • doc/guide/plugin.rst

    r9d8a027 r9150036  
    302302 
    303303**Note: The order of the parameters in the definition will be the order of the 
    304 parameters in the user interface and the order of the parameters in Iq(), 
    305 Iqac(), Iqabc() and form_volume(). And** *scale* **and** *background* 
    306 **parameters are implicit to all models, so they do not need to be included 
    307 in the parameter table.** 
     304parameters in the user interface and the order of the parameters in Fq(), Iq(), 
     305Iqac(), Iqabc(), form_volume() and shell_volume(). 
     306And** *scale* **and** *background* **parameters are implicit to all models, 
     307so they do not need to be included in the parameter table.** 
    308308 
    309309- **"name"** is the name of the parameter shown on the FitPage. 
     
    374374    scattered intensity. 
    375375 
    376   - "volume" parameters are passed to Iq(), Iqac(), Iqabc() and form_volume(), 
    377     and have polydispersity loops generated automatically. 
     376  - "volume" parameters are passed to Fq(), Iq(), Iqac(), Iqabc(), form_volume() 
     377    and shell_volume(), and have polydispersity loops generated automatically. 
    378378 
    379379  - "orientation" parameters are not passed, but instead are combined with 
     
    503503used. 
    504504 
     505Hollow shapes, where the volume fraction of particle corresponds to the 
     506material in the shell rather than the volume enclosed by the shape, must 
     507also define a *shell_volume(par1, par2, ...)* function.  The parameters 
     508are the same as for *form_volume*.  The *I(q)* calculation should use 
     509*shell_volume* squared as its scale factor for the volume normalization. 
     510The structure factor calculation needs *form_volume* in order to properly 
     511scale the volume fraction parameter, so both functions are required for 
     512hollow shapes. 
     513 
     514Note: Pure python models do not yet support direct computation of the 
     515average of $F(q)$ and $F^2(q)$. 
     516 
    505517Embedded C Models 
    506518................. 
     
    514526This expands into the equivalent C code:: 
    515527 
    516     #include <math.h> 
    517528    double Iq(double q, double par1, double par2, ...); 
    518529    double Iq(double q, double par1, double par2, ...) 
     
    523534*form_volume* defines the volume of the shape. As in python models, it 
    524535includes only the volume parameters. 
     536 
     537*form_volume* defines the volume of the shell for hollow shapes. As in 
     538python models, it includes only the volume parameters. 
    525539 
    526540**source=['fn.c', ...]** includes the listed C source files in the 
     
    559573The INVALID define can go into *Iq*, or *c_code*, or an external C file 
    560574listed in *source*. 
     575 
     576Structure Factors 
     577................. 
     578 
     579Structure factor calculations may need the underlying $<F(q)>$ and $<F^2(q)>$ 
     580rather than $I(q)$.  This is used to compute $\beta = <F(q)>^2/<F^2(q)>$ in 
     581the decoupling approximation to the structure factor. 
     582 
     583Instead of defining the *Iq* function, models can define *Fq* as 
     584something like:: 
     585 
     586    double Fq(double q, double *F1, double *F2, double par1, double par2, ...); 
     587    double Fq(double q, double *F1, double *F2, double par1, double par2, ...) 
     588    { 
     589        // Polar integration loop over all orientations. 
     590        ... 
     591        *F1 = 1e-2 * total_F1 * contrast * volume; 
     592        *F2 = 1e-4 * total_F2 * square(contrast * volume); 
     593        return I(q, par1, par2, ...); 
     594    } 
     595 
     596If the volume fraction scale factor is built into the model (as occurs for 
     597the vesicle model, for example), then scale *F1* by $\surd V_f$ so that 
     598$\beta$ is computed correctly. 
     599 
     600Structure factor calculations are not yet supported for oriented shapes. 
     601 
     602Note: only available as a separate C file listed in *source*, or within 
     603a *c_code* block within the python model definition file. 
    561604 
    562605Oriented Shapes 
     
    862905             - \frac{\sin(x)}{x}\left(\frac{1}{x} - \frac{3!}{x^3} + \frac{5!}{x^5} - \frac{7!}{x^7}\right) 
    863906 
    864         For small arguments , 
     907        For small arguments, 
    865908 
    866909        .. math:: 
     
    10231066          "radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, 
    10241067         0.2, 0.228843], 
    1025         [{"radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, "ER", 120.], 
    1026         [{"radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, "VR", 1.], 
     1068        [{"radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, 
     1069         0.1, None, None, 120., None, 1.],  # q, F, F^2, R_eff, V, form:shell 
     1070        [{"@S": "hardsphere"}, 0.1, None], 
    10271071    ] 
    10281072 
    10291073 
    1030 **tests=[[{parameters}, q, result], ...]** is a list of lists. 
     1074**tests=[[{parameters}, q, Iq], ...]** is a list of lists. 
    10311075Each list is one test and contains, in order: 
    10321076 
     
    10401084- input and output values can themselves be lists if you have several 
    10411085  $q$ values to test for the same model parameters. 
    1042 - for testing *ER* and *VR*, give the inputs as "ER" and "VR" respectively; 
    1043   the output for *VR* should be the sphere/shell ratio, not the individual 
    1044   sphere and shell values. 
     1086- for testing effective radius, volume and form:shell volume ratio, use the 
     1087  extended form of the tests results, with *None, None, R_eff, V, V_r* 
     1088  instead of *Iq*.  This calls the kernel *Fq* function instead of *Iq*. 
     1089- for testing F and F^2 (used for beta approximation) do the same as the 
     1090  effective radius test, but include values for the first two elements, 
     1091  $<F(q)>$ and $<F^2(q)>$. 
     1092- for testing interaction between form factor and structure factor, specify 
     1093  the structure factor name in the parameters as *{"@S": "name", ...}* with 
     1094  the remaining list of parameters defined by the *P@S* product model. 
    10451095 
    10461096.. _Test_Your_New_Model: 
  • doc/guide/scripting.rst

    rbd7630d r23df833  
    188188python kernel.  Once the kernel is in hand, we can then marshal a set of 
    189189parameters into a :class:`sasmodels.details.CallDetails` object and ship it to 
    190 the kernel using the :func:`sansmodels.direct_model.call_kernel` function.  An 
    191 example should help, *example/cylinder_eval.py*:: 
    192  
    193     from numpy import logspace 
     190the kernel using the :func:`sansmodels.direct_model.call_kernel` function.  To 
     191accesses the underlying $<F(q)>$ and $<F^2(q)>$, use 
     192:func:`sasmodels.direct_model.call_Fq` instead. 
     193 
     194The following example should 
     195help, *example/cylinder_eval.py*:: 
     196 
     197    from numpy import logspace, sqrt 
    194198    from matplotlib import pyplot as plt 
    195199    from sasmodels.core import load_model 
    196     from sasmodels.direct_model import call_kernel 
     200    from sasmodels.direct_model import call_kernel, call_Fq 
    197201 
    198202    model = load_model('cylinder') 
    199203    q = logspace(-3, -1, 200) 
    200204    kernel = model.make_kernel([q]) 
    201     Iq = call_kernel(kernel, dict(radius=200.)) 
    202     plt.loglog(q, Iq) 
     205    pars = {'radius': 200, 'radius_pd': 0.1, 'scale': 2} 
     206    Iq = call_kernel(kernel, pars) 
     207    F, Fsq, Reff, V, Vratio = call_Fq(kernel, pars) 
     208 
     209    plt.loglog(q, Iq, label='2 I(q)') 
     210    plt.loglog(q, F**2/V, label='<F(q)>^2/V') 
     211    plt.loglog(q, Fsq/V, label='<F^2(q)>/V') 
     212    plt.xlabel('q (1/A)') 
     213    plt.ylabel('I(q) (1/cm)') 
     214    plt.title('Cylinder with radius 200.') 
     215    plt.legend() 
    203216    plt.show() 
    204217 
    205 On windows, this can be called from the cmd prompt using sasview as:: 
     218.. figure:: direct_call.png 
     219 
     220    Comparison between $I(q)$, $<F(q)>$ and $<F^2(q)>$ for cylinder model. 
     221 
     222This compares $I(q)$ with $<F(q)>$ and $<F^2(q)>$ for a cylinder 
     223with *radius=200 +/- 20* and *scale=2*. Note that *call_Fq* does not 
     224include scale and background, nor does it normalize by the average volume. 
     225The definition of $F = \rho V \hat F$ scaled by the contrast and 
     226volume, compared to the canonical cylinder $\hat F$, with $\hat F(0) = 1$. 
     227Integrating over polydispersity and orientation, the returned values are 
     228$\sum_{r,w\in N(r_o, r_o/10)} \sum_\theta w F(q,r_o,\theta)\sin\theta$ and 
     229$\sum_{r,w\in N(r_o, r_o/10)} \sum_\theta w F^2(q,r_o,\theta)\sin\theta$. 
     230 
     231On windows, this example can be called from the cmd prompt using sasview as 
     232as the python interpreter:: 
    206233 
    207234    SasViewCom example/cylinder_eval.py 
  • example/cylinder_eval.py

    r2e66ef5 r23df833  
    33""" 
    44 
    5 from numpy import logspace 
     5from numpy import logspace, sqrt 
    66from matplotlib import pyplot as plt 
    77from sasmodels.core import load_model 
    8 from sasmodels.direct_model import call_kernel 
     8from sasmodels.direct_model import call_kernel, call_Fq 
    99 
    1010model = load_model('cylinder') 
    1111q = logspace(-3, -1, 200) 
    1212kernel = model.make_kernel([q]) 
    13 Iq = call_kernel(kernel, dict(radius=200.)) 
    14 plt.loglog(q, Iq) 
     13pars = {'radius': 200, 'radius_pd': 0.1, 'scale': 2} 
     14Iq = call_kernel(kernel, pars) 
     15F, Fsq, Reff, V, Vratio = call_Fq(kernel, pars) 
     16plt.loglog(q, Iq, label='2 I(q)') 
     17plt.loglog(q, F**2/V, label='<F(q)>^2/V') 
     18plt.loglog(q, Fsq/V, label='<F^2(q)>/V') 
    1519plt.xlabel('q (1/A)') 
    16 plt.ylabel('I(q)') 
     20plt.ylabel('I(q) (1/cm)') 
    1721plt.title('Cylinder with radius 200.') 
     22plt.legend() 
    1823plt.show() 
  • explore/precision.py

    rfba9ca0 rcd28947  
    207207    return model_info 
    208208 
    209 # Hack to allow second parameter A in two parameter functions 
     209# Hack to allow second parameter A in the gammainc and gammaincc functions. 
     210# Create a 2-D variant of the precision test if we need to handle other two 
     211# parameter functions. 
    210212A = 1 
    211213def parse_extra_pars(): 
     214    """ 
     215    Parse the command line looking for the second parameter "A=..." for the 
     216    gammainc/gammaincc functions. 
     217    """ 
    212218    global A 
    213219 
     
    333339) 
    334340add_function( 
     341    # Note: "a" is given as A=... on the command line via parse_extra_pars 
    335342    name="gammainc(x)", 
    336343    mp_function=lambda x, a=A: mp.gammainc(a, a=0, b=x)/mp.gamma(a), 
     
    339346) 
    340347add_function( 
     348    # Note: "a" is given as A=... on the command line via parse_extra_pars 
    341349    name="gammaincc(x)", 
    342350    mp_function=lambda x, a=A: mp.gammainc(a, a=x, b=mp.inf)/mp.gamma(a), 
     
    403411) 
    404412add_function( 
    405     name="debye", 
     413    name="gauss_coil", 
    406414    mp_function=lambda x: 2*(mp.exp(-x**2) + x**2 - 1)/x**4, 
    407415    np_function=lambda x: 2*(np.expm1(-x**2) + x**2)/x**4, 
    408416    ocl_function=make_ocl(""" 
    409417    const double qsq = q*q; 
    410     if (qsq < 1.0) { // Pade approximation 
     418    // For double: use O(5) Pade with 0.5 cutoff (10 mad + 1 divide) 
     419    // For single: use O(7) Taylor with 0.8 cutoff (7 mad) 
     420    if (qsq < 0.0) { 
    411421        const double x = qsq; 
    412422        if (0) { // 0.36 single 
     
    418428            const double B1=3./8., B2=3./56., B3=1./336.; 
    419429            return (((A3*x + A2)*x + A1)*x + 1.)/(((B3*x + B2)*x + B1)*x + 1.); 
    420         } else if (1) { // 1.0 for single, 0.25 for double 
     430        } else if (0) { // 1.0 for single, 0.25 for double 
    421431            // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}] 
    422432            const double A1=1./15., A2=1./60, A3=0., A4=1./75600.; 
     
    431441                  /(((((B5*x + B4)*x + B3)*x + B2)*x + B1)*x + 1.); 
    432442        } 
    433     } else if (qsq < 1.) { // Taylor series; 0.9 for single, 0.25 for double 
     443    } else if (qsq < 0.8) { 
    434444        const double x = qsq; 
    435445        const double C0 = +1.; 
  • sasmodels/__init__.py

    ra1ec908 r37f38ff  
    1414defining new models. 
    1515""" 
    16 __version__ = "0.98" 
     16__version__ = "0.99" 
    1717 
    1818def data_files(): 
  • sasmodels/compare.py

    r610ef23 rc1799d3  
    4141from . import kerneldll 
    4242from . import kernelcl 
     43from . import kernelcuda 
    4344from .data import plot_theory, empty_data1D, empty_data2D, load_data 
    4445from .direct_model import DirectModel, get_mesh 
     
    115116    === environment variables === 
    116117    -DSAS_MODELPATH=path sets directory containing custom models 
    117     -DSAS_OPENCL=vendor:device|none sets the target OpenCL device 
     118    -DSAS_OPENCL=vendor:device|cuda:device|none sets the target GPU device 
    118119    -DXDG_CACHE_HOME=~/.cache sets the pyopencl cache root (linux only) 
    119120    -DSAS_COMPILER=tinycc|msvc|mingw|unix sets the DLL compiler 
     
    352353 
    353354    # If it is a list of choices, pick one at random with equal probability 
    354     # In practice, the model specific random generator will override. 
    355355    par = model_info.parameters[name] 
    356     if len(par.limits) > 2:  # choice list 
    357         return np.random.randint(len(par.limits)) 
     356    if par.choices:  # choice list 
     357        return np.random.randint(len(par.choices)) 
    358358 
    359359    # If it is a fixed range, pick from it with equal probability. 
     
    725725        set_integration_size(model_info, ngauss) 
    726726 
    727     if dtype != "default" and not dtype.endswith('!') and not kernelcl.use_opencl(): 
     727    if (dtype != "default" and not dtype.endswith('!')  
     728            and not (kernelcl.use_opencl() or kernelcuda.use_cuda())): 
    728729        raise RuntimeError("OpenCL not available " + kernelcl.OPENCL_ERROR) 
    729730 
     
    11511152        'rel_err'   : True, 
    11521153        'explore'   : False, 
    1153         'use_demo'  : True, 
     1154        'use_demo'  : False, 
    11541155        'zero'      : False, 
    11551156        'html'      : False, 
  • sasmodels/core.py

    r2dcd6e7 r07646b6  
    2121from . import mixture 
    2222from . import kernelpy 
     23from . import kernelcuda 
    2324from . import kernelcl 
    2425from . import kerneldll 
     
    205206 
    206207    numpy_dtype, fast, platform = parse_dtype(model_info, dtype, platform) 
    207  
    208208    source = generate.make_source(model_info) 
    209209    if platform == "dll": 
    210210        #print("building dll", numpy_dtype) 
    211211        return kerneldll.load_dll(source['dll'], model_info, numpy_dtype) 
     212    elif platform == "cuda": 
     213        return kernelcuda.GpuModel(source, model_info, numpy_dtype, fast=fast) 
    212214    else: 
    213215        #print("building ocl", numpy_dtype) 
     
    245247    # type: (ModelInfo, str, str) -> (np.dtype, bool, str) 
    246248    """ 
    247     Interpret dtype string, returning np.dtype and fast flag. 
     249    Interpret dtype string, returning np.dtype, fast flag and platform. 
    248250 
    249251    Possible types include 'half', 'single', 'double' and 'quad'.  If the 
     
    253255    default for the model and platform. 
    254256 
    255     Platform preference can be specfied ("ocl" vs "dll"), with the default 
    256     being OpenCL if it is availabe.  If the dtype name ends with '!' then 
    257     platform is forced to be DLL rather than OpenCL. 
     257    Platform preference can be specfied ("ocl", "cuda", "dll"), with the 
     258    default being OpenCL or CUDA if available, otherwise DLL.  If the dtype 
     259    name ends with '!' then platform is forced to be DLL rather than GPU. 
     260    The default platform is set by the environment variable SAS_OPENCL, 
     261    SAS_OPENCL=driver:device for OpenCL, SAS_OPENCL=cuda:device for CUDA 
     262    or SAS_OPENCL=none for DLL. 
    258263 
    259264    This routine ignores the preferences within the model definition.  This 
     
    265270    # Assign default platform, overriding ocl with dll if OpenCL is unavailable 
    266271    # If opencl=False OpenCL is switched off 
    267  
    268272    if platform is None: 
    269273        platform = "ocl" 
    270     if not kernelcl.use_opencl() or not model_info.opencl: 
    271         platform = "dll" 
    272274 
    273275    # Check if type indicates dll regardless of which platform is given 
     
    275277        platform = "dll" 
    276278        dtype = dtype[:-1] 
     279 
     280    # Make sure model allows opencl/gpu 
     281    if not model_info.opencl: 
     282        platform = "dll" 
     283 
     284    # Make sure opencl is available, or fallback to cuda then to dll 
     285    if platform == "ocl" and not kernelcl.use_opencl(): 
     286        platform = "cuda" if kernelcuda.use_cuda() else "dll" 
    277287 
    278288    # Convert special type names "half", "fast", and "quad" 
     
    285295        dtype = "float16" 
    286296 
    287     # Convert dtype string to numpy dtype. 
     297    # Convert dtype string to numpy dtype.  Use single precision for GPU 
     298    # if model allows it, otherwise use double precision. 
    288299    if dtype is None or dtype == "default": 
    289         numpy_dtype = (generate.F32 if platform == "ocl" and model_info.single 
     300        numpy_dtype = (generate.F32 if model_info.single and platform in ("ocl", "cuda") 
    290301                       else generate.F64) 
    291302    else: 
    292303        numpy_dtype = np.dtype(dtype) 
    293304 
    294     # Make sure that the type is supported by opencl, otherwise use dll 
     305    # Make sure that the type is supported by GPU, otherwise use dll 
    295306    if platform == "ocl": 
    296307        env = kernelcl.environment() 
    297         if not env.has_type(numpy_dtype): 
    298             platform = "dll" 
    299             if dtype is None: 
    300                 numpy_dtype = generate.F64 
     308    elif platform == "cuda": 
     309        env = kernelcuda.environment() 
     310    else: 
     311        env = None 
     312    if env is not None and not env.has_type(numpy_dtype): 
     313        platform = "dll" 
     314        if dtype is None: 
     315            numpy_dtype = generate.F64 
    301316 
    302317    return numpy_dtype, fast, platform 
     
    365380    model = load_model("cylinder@hardsphere*sphere") 
    366381    actual = [p.name for p in model.info.parameters.kernel_parameters] 
    367     target = ("sld sld_solvent radius length theta phi volfraction" 
     382    target = ("sld sld_solvent radius length theta phi" 
     383              " radius_effective volfraction " 
     384              " structure_factor_mode radius_effective_mode" 
    368385              " A_sld A_sld_solvent A_radius").split() 
    369386    assert target == actual, "%s != %s"%(target, actual) 
  • sasmodels/data.py

    rbd7630d rc88f983  
    337337    else: 
    338338        dq = resolution * q 
    339  
    340339    data = Data1D(q, Iq, dx=dq, dy=dIq) 
    341340    data.filename = "fake data" 
  • sasmodels/direct_model.py

    r2c4a190 r9150036  
    3232from .details import make_kernel_args, dispersion_mesh 
    3333from .modelinfo import DEFAULT_BACKGROUND 
     34from .product import RADIUS_MODE_ID 
    3435 
    3536# pylint: disable=unused-import 
     
    6465    return calculator(call_details, values, cutoff, is_magnetic) 
    6566 
    66 def call_ER(model_info, pars): 
    67     # type: (ModelInfo, ParameterSet) -> float 
    68     """ 
    69     Call the model ER function using *values*. 
    70  
    71     *model_info* is either *model.info* if you have a loaded model, 
    72     or *kernel.info* if you have a model kernel prepared for evaluation. 
    73     """ 
    74     if model_info.ER is None: 
    75         return 1.0 
    76     elif not model_info.parameters.form_volume_parameters: 
    77         # handle the case where ER is provided but model is not polydisperse 
    78         return model_info.ER() 
    79     else: 
    80         value, weight = _vol_pars(model_info, pars) 
    81         individual_radii = model_info.ER(*value) 
    82         return np.sum(weight*individual_radii) / np.sum(weight) 
    83  
    84  
    85 def call_VR(model_info, pars): 
    86     # type: (ModelInfo, ParameterSet) -> float 
    87     """ 
    88     Call the model VR function using *pars*. 
    89  
    90     *model_info* is either *model.info* if you have a loaded model, 
    91     or *kernel.info* if you have a model kernel prepared for evaluation. 
    92     """ 
    93     if model_info.VR is None: 
    94         return 1.0 
    95     elif not model_info.parameters.form_volume_parameters: 
    96         # handle the case where ER is provided but model is not polydisperse 
    97         return model_info.VR() 
    98     else: 
    99         value, weight = _vol_pars(model_info, pars) 
    100         whole, part = model_info.VR(*value) 
    101         return np.sum(weight*part)/np.sum(weight*whole) 
    102  
    103  
    104 def call_profile(model_info, **pars): 
    105     # type: (ModelInfo, ...) -> Tuple[np.ndarray, np.ndarray, Tuple[str, str]] 
     67def call_Fq(calculator, pars, cutoff=0., mono=False): 
     68    # type: (Kernel, ParameterSet, float, bool) -> np.ndarray 
     69    """ 
     70    Like :func:`call_kernel`, but returning F, F^2, R_eff, V_shell, V_form/V_shell. 
     71 
     72    For solid objects V_shell is equal to V_form and the volume ratio is 1. 
     73 
     74    Use parameter *radius_effective_mode* to select the effective radius 
     75    calculation. 
     76    """ 
     77    R_eff_type = int(pars.pop(RADIUS_MODE_ID, 1.0)) 
     78    mesh = get_mesh(calculator.info, pars, dim=calculator.dim, mono=mono) 
     79    #print("pars", list(zip(*mesh))[0]) 
     80    call_details, values, is_magnetic = make_kernel_args(calculator, mesh) 
     81    #print("values:", values) 
     82    return calculator.Fq(call_details, values, cutoff, is_magnetic, R_eff_type) 
     83 
     84def call_profile(model_info, pars=None): 
     85    # type: (ModelInfo, ParameterSet) -> Tuple[np.ndarray, np.ndarray, Tuple[str, str]] 
    10686    """ 
    10787    Returns the profile *x, y, (xlabel, ylabel)* representing the model. 
    10888    """ 
     89    if pars is None: 
     90        pars = {} 
    10991    args = {} 
    11092    for p in model_info.parameters.kernel_parameters: 
     
    347329 
    348330        # Need to pull background out of resolution for multiple scattering 
    349         background = pars.get('background', DEFAULT_BACKGROUND) 
     331        default_background = self._model.info.parameters.common_parameters[1].default 
     332        background = pars.get('background', default_background) 
    350333        pars = pars.copy() 
    351334        pars['background'] = 0. 
     
    400383        Generate a plottable profile. 
    401384        """ 
    402         return call_profile(self.model.info, **pars) 
     385        return call_profile(self.model.info, pars) 
    403386 
    404387def main(): 
     
    416399    model_name = sys.argv[1] 
    417400    call = sys.argv[2].upper() 
    418     if call != "ER_VR": 
    419         try: 
    420             values = [float(v) for v in call.split(',')] 
    421         except ValueError: 
    422             values = [] 
    423         if len(values) == 1: 
    424             q, = values 
    425             data = empty_data1D([q]) 
    426         elif len(values) == 2: 
    427             qx, qy = values 
    428             data = empty_data2D([qx], [qy]) 
    429         else: 
    430             print("use q or qx,qy or ER or VR") 
    431             sys.exit(1) 
     401    try: 
     402        values = [float(v) for v in call.split(',')] 
     403    except ValueError: 
     404        values = [] 
     405    if len(values) == 1: 
     406        q, = values 
     407        data = empty_data1D([q]) 
     408    elif len(values) == 2: 
     409        qx, qy = values 
     410        data = empty_data2D([qx], [qy]) 
    432411    else: 
    433         data = empty_data1D([0.001])  # Data not used in ER/VR 
     412        print("use q or qx,qy") 
     413        sys.exit(1) 
    434414 
    435415    model_info = load_model_info(model_name) 
     
    439419                for pair in sys.argv[3:] 
    440420                for k, v in [pair.split('=')]) 
    441     if call == "ER_VR": 
    442         ER = call_ER(model_info, pars) 
    443         VR = call_VR(model_info, pars) 
    444         print(ER, VR) 
    445     else: 
    446         Iq = calculator(**pars) 
    447         print(Iq[0]) 
     421    Iq = calculator(**pars) 
     422    print(Iq[0]) 
    448423 
    449424if __name__ == "__main__": 
  • sasmodels/generate.py

    r6e45516 ra8a1d48  
    2424    dimension, or 1.0 if no volume normalization is required. 
    2525 
    26     *ER(p1, p2, ...)* returns the effective radius of the form with 
    27     particular dimensions. 
    28  
    29     *VR(p1, p2, ...)* returns the volume ratio for core-shell style forms. 
     26    *shell_volume(p1, p2, ...)* returns the volume of the shell for forms 
     27    which are hollow. 
     28 
     29    *effective_radius(mode, p1, p2, ...)* returns the effective radius of 
     30    the form with particular dimensions.  Mode determines the type of 
     31    effective radius returned, with mode=1 for equivalent volume. 
    3032 
    3133    #define INVALID(v) (expr)  returns False if v.parameter is invalid 
     
    7274background.  This will work correctly even when polydispersity is off. 
    7375 
    74 *ER* and *VR* are python functions which operate on parameter vectors. 
    75 The constructor code will generate the necessary vectors for computing 
    76 them with the desired polydispersity. 
    7776The kernel module must set variables defining the kernel meta data: 
    7877 
     
    106105    create the OpenCL kernel functions.  The files defining the functions 
    107106    need to be listed before the files which use the functions. 
    108  
    109     *ER* is a python function defining the effective radius.  If it is 
    110     not present, the effective radius is 0. 
    111  
    112     *VR* is a python function defining the volume ratio.  If it is not 
    113     present, the volume ratio is 1. 
    114107 
    115108    *form_volume*, *Iq*, *Iqac*, *Iqabc* are strings containing 
     
    671664 
    672665 
    673 # type in IQXY pattern could be single, float, double, long double, ... 
    674 _IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ab?c|xy))\s*[(]", 
     666_IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ac|abc|xy))\s*[(]", 
    675667                           flags=re.MULTILINE) 
    676668def find_xy_mode(source): 
     
    701693    return 'qa' 
    702694 
     695# Note: not presently used.  Need to know whether Fq is available before 
     696# trying to compile the source.  Leave the code here in case we decide that 
     697# define have_Fq for each form factor is too tedious and error prone. 
     698_FQ_PATTERN = re.compile(r"(^|\s)void\s+Fq[(]", flags=re.MULTILINE) 
     699def contains_Fq(source): 
     700    # type: (List[str]) -> bool 
     701    """ 
     702    Return True if C source defines "void Fq(". 
     703    """ 
     704    for code in source: 
     705        if _FQ_PATTERN.search(code) is not None: 
     706            return True 
     707    return False 
     708 
     709_SHELL_VOLUME_PATTERN = re.compile(r"(^|\s)double\s+shell_volume[(]", flags=re.MULTILINE) 
     710def contains_shell_volume(source): 
     711    # type: (List[str]) -> bool 
     712    """ 
     713    Return True if C source defines "double shell_volume(". 
     714    """ 
     715    for code in source: 
     716        if _SHELL_VOLUME_PATTERN.search(code) is not None: 
     717            return True 
     718    return False 
    703719 
    704720def _add_source(source, code, path, lineno=1): 
     
    730746    # dispersion.  Need to be careful that necessary parameters are available 
    731747    # for computing volume even if we allow non-disperse volume parameters. 
    732  
    733748    partable = model_info.parameters 
    734749 
     
    743758    for path, code in user_code: 
    744759        _add_source(source, code, path) 
    745  
    746760    if model_info.c_code: 
    747761        _add_source(source, model_info.c_code, model_info.filename, 
     
    751765    q, qx, qy, qab, qa, qb, qc \ 
    752766        = [Parameter(name=v) for v in 'q qx qy qab qa qb qc'.split()] 
     767 
    753768    # Generate form_volume function, etc. from body only 
    754769    if isinstance(model_info.form_volume, str): 
    755770        pars = partable.form_volume_parameters 
    756771        source.append(_gen_fn(model_info, 'form_volume', pars)) 
     772    if isinstance(model_info.shell_volume, str): 
     773        pars = partable.form_volume_parameters 
     774        source.append(_gen_fn(model_info, 'shell_volume', pars)) 
    757775    if isinstance(model_info.Iq, str): 
    758776        pars = [q] + partable.iq_parameters 
     
    767785        pars = [qa, qb, qc] + partable.iq_parameters 
    768786        source.append(_gen_fn(model_info, 'Iqabc', pars)) 
     787 
     788    # Check for shell_volume in source 
     789    is_hollow = contains_shell_volume(source) 
    769790 
    770791    # What kind of 2D model do we need?  Is it consistent with the parameters? 
     
    789810    source.append("\\\n".join(p.as_definition() 
    790811                              for p in partable.kernel_parameters)) 
    791  
    792812    # Define the function calls 
     813    call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) 0.0" 
    793814    if partable.form_volume_parameters: 
    794815        refs = _call_pars("_v.", partable.form_volume_parameters) 
    795         call_volume = "#define CALL_VOLUME(_v) form_volume(%s)"%(",".join(refs)) 
     816        if is_hollow: 
     817            call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = form_volume(%s); _shell = shell_volume(%s); } while (0)"%((",".join(refs),)*2) 
     818        else: 
     819            call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = _shell = form_volume(%s); } while (0)"%(",".join(refs)) 
     820        if model_info.effective_radius_type: 
     821            call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) effective_radius(_mode, %s)"%(",".join(refs)) 
    796822    else: 
    797823        # Model doesn't have volume.  We could make the kernel run a little 
    798824        # faster by not using/transferring the volume normalizations, but 
    799825        # the ifdef's reduce readability more than is worthwhile. 
    800         call_volume = "#define CALL_VOLUME(v) 1.0" 
     826        call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = _shell = 1.0; } while (0)" 
    801827    source.append(call_volume) 
    802  
     828    source.append(call_effective_radius) 
    803829    model_refs = _call_pars("_v.", partable.iq_parameters) 
    804     pars = ",".join(["_q"] + model_refs) 
    805     call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars 
     830 
     831    if model_info.have_Fq: 
     832        pars = ",".join(["_q", "&_F1", "&_F2",] + model_refs) 
     833        call_iq = "#define CALL_FQ(_q, _F1, _F2, _v) Fq(%s)" % pars 
     834        clear_iq = "#undef CALL_FQ" 
     835    else: 
     836        pars = ",".join(["_q"] + model_refs) 
     837        call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars 
     838        clear_iq = "#undef CALL_IQ" 
    806839    if xy_mode == 'qabc': 
    807840        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs) 
     
    812845        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqac(%s)" % pars 
    813846        clear_iqxy = "#undef CALL_IQ_AC" 
    814     elif xy_mode == 'qa': 
     847    elif xy_mode == 'qa' and not model_info.have_Fq: 
    815848        pars = ",".join(["_qa"] + model_refs) 
    816849        call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars 
    817850        clear_iqxy = "#undef CALL_IQ_A" 
     851    elif xy_mode == 'qa' and model_info.have_Fq: 
     852        pars = ",".join(["_qa", "&_F1", "&_F2",] + model_refs) 
     853        # Note: uses rare C construction (expr1, expr2) which computes 
     854        # expr1 then expr2 and evaluates to expr2.  This allows us to 
     855        # leave it looking like a function even though it is returning 
     856        # its values by reference. 
     857        call_iqxy = "#define CALL_FQ_A(_qa,_F1,_F2,_v) (Fq(%s),_F2)" % pars 
     858        clear_iqxy = "#undef CALL_FQ_A" 
    818859    elif xy_mode == 'qxy': 
    819860        orientation_refs = _call_pars("_v.", partable.orientation_parameters) 
     
    831872    magpars = [k-2 for k, p in enumerate(partable.call_parameters) 
    832873               if p.type == 'sld'] 
    833  
    834874    # Fill in definitions for numbers of parameters 
    835875    source.append("#define MAX_PD %s"%partable.max_pd) 
     
    839879    source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars)) 
    840880    source.append("#define PROJECTION %d"%PROJECTION) 
    841  
    842881    # TODO: allow mixed python/opencl kernels? 
    843  
    844     ocl = _kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name) 
    845     dll = _kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name) 
     882    ocl = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name) 
     883    dll = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name) 
     884 
    846885    result = { 
    847886        'dll': '\n'.join(source+dll[0]+dll[1]+dll[2]), 
    848887        'opencl': '\n'.join(source+ocl[0]+ocl[1]+ocl[2]), 
    849888    } 
    850  
    851889    return result 
    852890 
    853891 
    854 def _kernels(kernel, call_iq, call_iqxy, clear_iqxy, name): 
     892def _kernels(kernel, call_iq, clear_iq, call_iqxy, clear_iqxy, name): 
    855893    # type: ([str,str], str, str, str) -> List[str] 
    856894    code = kernel[0] 
     
    862900        '#line 1 "%s Iq"' % path, 
    863901        code, 
    864         "#undef CALL_IQ", 
     902        clear_iq, 
    865903        "#undef KERNEL_NAME", 
    866904        ] 
     
    9681006        pars = model_info.parameters.kernel_parameters 
    9691007    else: 
    970         pars = model_info.parameters.COMMON + model_info.parameters.kernel_parameters 
     1008        pars = (model_info.parameters.common_parameters 
     1009                + model_info.parameters.kernel_parameters) 
    9711010    partable = make_partable(pars) 
    9721011    subst = dict(id=model_info.id.replace('_', '-'), 
  • sasmodels/jitter.py

    r1198f90 r7d97437  
    1515    pass 
    1616 
     17import matplotlib as mpl 
    1718import matplotlib.pyplot as plt 
    1819from matplotlib.widgets import Slider 
     
    746747        pass 
    747748 
    748     axcolor = 'lightgoldenrodyellow' 
     749    # CRUFT: use axisbg instead of facecolor for matplotlib<2 
     750    facecolor_prop = 'facecolor' if mpl.__version__ > '2' else 'axisbg' 
     751    props = {facecolor_prop: 'lightgoldenrodyellow'} 
    749752 
    750753    ## add control widgets to plot 
    751     axes_theta = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor) 
    752     axes_phi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor) 
    753     axes_psi = plt.axes([0.1, 0.05, 0.45, 0.04], axisbg=axcolor) 
     754    axes_theta = plt.axes([0.1, 0.15, 0.45, 0.04], **props) 
     755    axes_phi = plt.axes([0.1, 0.1, 0.45, 0.04], **props) 
     756    axes_psi = plt.axes([0.1, 0.05, 0.45, 0.04], **props) 
    754757    stheta = Slider(axes_theta, 'Theta', -90, 90, valinit=theta) 
    755758    sphi = Slider(axes_phi, 'Phi', -180, 180, valinit=phi) 
    756759    spsi = Slider(axes_psi, 'Psi', -180, 180, valinit=psi) 
    757760 
    758     axes_dtheta = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor) 
    759     axes_dphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor) 
    760     axes_dpsi = plt.axes([0.75, 0.05, 0.15, 0.04], axisbg=axcolor) 
     761    axes_dtheta = plt.axes([0.75, 0.15, 0.15, 0.04], **props) 
     762    axes_dphi = plt.axes([0.75, 0.1, 0.15, 0.04], **props) 
     763    axes_dpsi = plt.axes([0.75, 0.05, 0.15, 0.04], **props) 
    761764    # Note: using ridiculous definition of rectangle distribution, whose width 
    762765    # in sasmodels is sqrt(3) times the given width.  Divide by sqrt(3) to keep 
  • sasmodels/kernel.py

    r2d81cfe rcd28947  
    4141    dtype = None  # type: np.dtype 
    4242 
    43     def __call__(self, call_details, values, cutoff, magnetic): 
     43    def Iq(self, call_details, values, cutoff, magnetic): 
    4444        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    45         raise NotImplementedError("need to implement __call__") 
     45        r""" 
     46        Returns I(q) from the polydisperse average scattering. 
     47 
     48        .. math:: 
     49 
     50            I(q) = \text{scale} \cdot P(q) + \text{background} 
     51 
     52        With the correct choice of model and contrast, setting *scale* to 
     53        the volume fraction $V_f$ of particles should match the measured 
     54        absolute scattering.  Some models (e.g., vesicle) have volume fraction 
     55        built into the model, and do not need an additional scale. 
     56        """ 
     57        _, F2, _, shell_volume, _ = self.Fq(call_details, values, cutoff, magnetic, 
     58                              effective_radius_type=0) 
     59        combined_scale = values[0]/shell_volume 
     60        background = values[1] 
     61        return combined_scale*F2 + background 
     62    __call__ = Iq 
     63 
     64    def Fq(self, call_details, values, cutoff, magnetic, effective_radius_type=0): 
     65        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
     66        r""" 
     67        Returns <F(q)>, <F(q)^2>, effective radius, shell volume and 
     68        form:shell volume ratio.  The <F(q)> term may be None if the 
     69        form factor does not support direct computation of $F(q)$ 
     70 
     71        $P(q) = <F^2(q)>/<V>$ is used for structure factor calculations, 
     72 
     73        .. math:: 
     74 
     75            I(q) = \text{scale} \cdot P(q) \cdot S(q) + \text{background} 
     76 
     77        For the beta approximation, this becomes 
     78 
     79        .. math:: 
     80 
     81            I(q) = \text{scale} * P (1 + <F>^2/<F^2> (S - 1)) + \text{background} 
     82                 = \text{scale}/<V> (<F^2> + <F>^2 (S - 1)) + \text{background} 
     83 
     84        $<F(q)>$ and $<F^2(q)>$ are averaged by polydispersity in shape 
     85        and orientation, with each configuration $x_k$ having form factor 
     86        $F(q, x_k)$, weight $w_k$ and volume $V_k$.  The result is: 
     87 
     88        .. math:: 
     89 
     90            P(q) = \frac{\sum w_k F^2(q, x_k) / \sum w_k}{\sum w_k V_k / \sum w_k} 
     91 
     92        The form factor itself is scaled by volume and contrast to compute the 
     93        total scattering.  This is then squared, and the volume weighted 
     94        F^2 is then normalized by volume F.  For a given density, the number 
     95        of scattering centers is assumed to scale linearly with volume.  Later 
     96        scaling the resulting $P(q)$ by the volume fraction of particles 
     97        gives the total scattering on an absolute scale. Most models 
     98        incorporate the volume fraction into the overall scale parameter.  An 
     99        exception is vesicle, which includes the volume fraction parameter in 
     100        the model itself, scaling $F$ by $\surd V_f$ so that the math for 
     101        the beta approximation works out. 
     102 
     103        By scaling $P(q)$ by total weight $\sum w_k$, there is no need to make 
     104        sure that the polydisperisity distributions normalize to one.  In 
     105        particular, any distibution values $x_k$ outside the valid domain 
     106        of $F$ will not be included, and the distribution will be implicitly 
     107        truncated.  This is controlled by the parameter limits defined in the 
     108        model (which truncate the distribution before calling the kernel) as 
     109        well as any region excluded using the *INVALID* macro defined within 
     110        the model itself. 
     111 
     112        The volume used in the polydispersity calculation is the form volume 
     113        for solid objects or the shell volume for hollow objects.  Shell 
     114        volume should be used within $F$ so that the normalizing scale 
     115        represents the volume fraction of the shell rather than the entire 
     116        form.  This corresponds to the volume fraction of shell-forming 
     117        material added to the solvent. 
     118 
     119        The calculation of $S$ requires the effective radius and the 
     120        volume fraction of the particles.  The model can have several 
     121        different ways to compute effective radius, with the 
     122        *effective_radius_type* parameter used to select amongst them.  The 
     123        volume fraction of particles should be determined from the total 
     124        volume fraction of the form, not just the shell volume fraction. 
     125        This makes a difference for hollow shapes, which need to scale 
     126        the volume fraction by the returned volume ratio when computing $S$. 
     127        For solid objects, the shell volume is set to the form volume so 
     128        this scale factor evaluates to one and so can be used for both 
     129        hollow and solid shapes. 
     130        """ 
     131        self._call_kernel(call_details, values, cutoff, magnetic, effective_radius_type) 
     132        #print("returned",self.q_input.q, self.result) 
     133        nout = 2 if self.info.have_Fq and self.dim == '1d' else 1 
     134        total_weight = self.result[nout*self.q_input.nq + 0] 
     135        # Note: total_weight = sum(weight > cutoff), with cutoff >= 0, so it 
     136        # is okay to test directly against zero.  If weight is zero then I(q), 
     137        # etc. must also be zero. 
     138        if total_weight == 0.: 
     139            total_weight = 1. 
     140        # Note: shell_volume == form_volume for solid objects 
     141        form_volume = self.result[nout*self.q_input.nq + 1]/total_weight 
     142        shell_volume = self.result[nout*self.q_input.nq + 2]/total_weight 
     143        effective_radius = self.result[nout*self.q_input.nq + 3]/total_weight 
     144        if shell_volume == 0.: 
     145            shell_volume = 1. 
     146        F1 = self.result[1:nout*self.q_input.nq:nout]/total_weight if nout == 2 else None 
     147        F2 = self.result[0:nout*self.q_input.nq:nout]/total_weight 
     148        return F1, F2, effective_radius, shell_volume, form_volume/shell_volume 
    46149 
    47150    def release(self): 
    48151        # type: () -> None 
    49152        pass 
     153 
     154    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
     155        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
     156        """ 
     157        Call the kernel.  Subclasses defining kernels for particular execution 
     158        engines need to provide an implementation for this. 
     159        """ 
     160        raise NotImplementedError() 
  • sasmodels/kernel_header.c

    r108e70e r07646b6  
    11#ifdef __OPENCL_VERSION__ 
    22# define USE_OPENCL 
     3#elif defined(__CUDACC__) 
     4# define USE_CUDA 
    35#elif defined(_OPENMP) 
    46# define USE_OPENMP 
    57#endif 
     8 
     9// Use SAS_DOUBLE to force the use of double even for float kernels 
     10#define SAS_DOUBLE dou ## ble 
    611 
    712// If opencl is not available, then we are compiling a C function 
    813// Note: if using a C++ compiler, then define kernel as extern "C" 
    914#ifdef USE_OPENCL 
     15 
     16   #define USE_GPU 
     17   #define pglobal global 
     18   #define pconstant constant 
     19 
    1020   typedef int int32_t; 
    11 #  if defined(USE_SINCOS) 
    12 #    define SINCOS(angle,svar,cvar) svar=sincos(angle,&cvar) 
    13 #  else 
    14 #    define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0) 
    15 #  endif 
     21 
     22   #if defined(USE_SINCOS) 
     23   #  define SINCOS(angle,svar,cvar) svar=sincos(angle,&cvar) 
     24   #else 
     25   #  define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0) 
     26   #endif 
    1627   // Intel CPU on Mac gives strange values for erf(); on the verified 
    1728   // platforms (intel, nvidia, amd), the cephes erf() is significantly 
     
    2435   #  define erfcf erfc 
    2536   #endif 
    26 #else // !USE_OPENCL 
    27 // Use SAS_DOUBLE to force the use of double even for float kernels 
    28 #  define SAS_DOUBLE dou ## ble 
    29 #  ifdef __cplusplus 
     37 
     38#elif defined(USE_CUDA) 
     39 
     40   #define USE_GPU 
     41   #define local __shared__ 
     42   #define pglobal 
     43   #define constant __constant__ 
     44   #define pconstant const 
     45   #define kernel extern "C" __global__ 
     46 
     47   // OpenCL powr(a,b) = C99 pow(a,b), b >= 0 
     48   // OpenCL pown(a,b) = C99 pow(a,b), b integer 
     49   #define powr(a,b) pow(a,b) 
     50   #define pown(a,b) pow(a,b) 
     51   //typedef int int32_t; 
     52   #if defined(USE_SINCOS) 
     53   #  define SINCOS(angle,svar,cvar) sincos(angle,&svar,&cvar) 
     54   #else 
     55   #  define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0) 
     56   #endif 
     57 
     58#else // !USE_OPENCL && !USE_CUDA 
     59 
     60   #define local 
     61   #define pglobal 
     62   #define constant const 
     63   #define pconstant const 
     64 
     65   #ifdef __cplusplus 
    3066      #include <cstdio> 
    3167      #include <cmath> 
     
    5187     #endif 
    5288     inline void SINCOS(double angle, double &svar, double &cvar) { svar=sin(angle); cvar=cos(angle); } 
    53 else // !__cplusplus 
     89   #else // !__cplusplus 
    5490     #include <inttypes.h>  // C99 guarantees that int32_t types is here 
    5591     #include <stdio.h> 
     
    68104         #define NEED_EXPM1 
    69105         #define NEED_TGAMMA 
     106         #define NEED_CBRT 
    70107         // expf missing from windows? 
    71108         #define expf exp 
     
    76113     #define kernel 
    77114     #define SINCOS(angle,svar,cvar) do {const double _t_=angle; svar=sin(_t_);cvar=cos(_t_);} while (0) 
    78 #  endif  // !__cplusplus 
    79 #  define global 
    80 #  define local 
    81 #  define constant const 
    82 // OpenCL powr(a,b) = C99 pow(a,b), b >= 0 
    83 // OpenCL pown(a,b) = C99 pow(a,b), b integer 
    84 #  define powr(a,b) pow(a,b) 
    85 #  define pown(a,b) pow(a,b) 
     115   #endif  // !__cplusplus 
     116   // OpenCL powr(a,b) = C99 pow(a,b), b >= 0 
     117   // OpenCL pown(a,b) = C99 pow(a,b), b integer 
     118   #define powr(a,b) pow(a,b) 
     119   #define pown(a,b) pow(a,b) 
     120 
    86121#endif // !USE_OPENCL 
     122 
     123#if defined(NEED_CBRT) 
     124   #define cbrt(_x) pow(_x, 0.33333333333333333333333) 
     125#endif 
    87126 
    88127#if defined(NEED_EXPM1) 
  • sasmodels/kernel_iq.c

    r70530778 r12f4c19  
    2626//  MAGNETIC_PARS : a comma-separated list of indices to the sld 
    2727//      parameters in the parameter table. 
    28 //  CALL_VOLUME(table) : call the form volume function 
     28//  CALL_VOLUME(form, shell, table) : assign form and shell values 
     29//  CALL_EFFECTIVE_RADIUS(type, table) : call the R_eff function 
    2930//  CALL_IQ(q, table) : call the Iq function for 1D calcs. 
    3031//  CALL_IQ_A(q, table) : call the Iq function with |q| for 2D data. 
     32//  CALL_FQ(q, F1, F2, table) : call the Fq function for 1D calcs. 
     33//  CALL_FQ_A(q, F1, F2, table) : call the Iq function with |q| for 2D data. 
    3134//  CALL_IQ_AC(qa, qc, table) : call the Iqxy function for symmetric shapes 
    3235//  CALL_IQ_ABC(qa, qc, table) : call the Iqxy function for asymmetric shapes 
     
    3639//  PROJECTION : equirectangular=1, sinusoidal=2 
    3740//      see explore/jitter.py for definitions. 
     41 
    3842 
    3943#ifndef _PAR_BLOCK_ // protected block so we can include this code twice. 
     
    270274#endif // _QABC_SECTION 
    271275 
    272  
    273276// ==================== KERNEL CODE ======================== 
    274  
    275277kernel 
    276278void KERNEL_NAME( 
    277     int32_t nq,                 // number of q values 
    278     const int32_t pd_start,     // where we are in the dispersity loop 
    279     const int32_t pd_stop,      // where we are stopping in the dispersity loop 
    280     global const ProblemDetails *details, 
    281     global const double *values, 
    282     global const double *q, // nq q values, with padding to boundary 
    283     global double *result,  // nq+1 return values, again with padding 
    284     const double cutoff     // cutoff in the dispersity weight product 
     279    int32_t nq,                   // number of q values 
     280    const int32_t pd_start,       // where we are in the dispersity loop 
     281    const int32_t pd_stop,        // where we are stopping in the dispersity loop 
     282    pglobal const ProblemDetails *details, 
     283    pglobal const double *values, // parameter values and distributions 
     284    pglobal const double *q,      // nq q values, with padding to boundary 
     285    pglobal double *result,       // nq+1 return values, again with padding 
     286    const double cutoff,          // cutoff in the dispersity weight product 
     287    int32_t effective_radius_type // which effective radius to compute 
    285288    ) 
    286289{ 
    287 #ifdef USE_OPENCL 
     290#if defined(USE_GPU) 
    288291  // who we are and what element we are working with 
     292  #if defined(USE_OPENCL) 
    289293  const int q_index = get_global_id(0); 
     294  #else // USE_CUDA 
     295  const int q_index = threadIdx.x + blockIdx.x * blockDim.x; 
     296  #endif 
    290297  if (q_index >= nq) return; 
    291298#else 
     
    338345  // 
    339346  // The code differs slightly between opencl and dll since opencl is only 
    340   // seeing one q value (stored in the variable "this_result") while the dll 
     347  // seeing one q value (stored in the variable "this_F2") while the dll 
    341348  // version must loop over all q. 
    342   #ifdef USE_OPENCL 
    343     double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
    344     double this_result = (pd_start == 0 ? 0.0 : result[q_index]); 
    345   #else // !USE_OPENCL 
    346     double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     349  #if defined(CALL_FQ) 
     350    double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq]); 
     351    double weighted_form = (pd_start == 0 ? 0.0 : result[2*nq+1]); 
     352    double weighted_shell = (pd_start == 0 ? 0.0 : result[2*nq+2]); 
     353    double weighted_radius = (pd_start == 0 ? 0.0 : result[2*nq+3]); 
     354  #else 
     355    double weight_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     356    double weighted_form = (pd_start == 0 ? 0.0 : result[nq+1]); 
     357    double weighted_shell = (pd_start == 0 ? 0.0 : result[nq+2]); 
     358    double weighted_radius = (pd_start == 0 ? 0.0 : result[nq+3]); 
     359  #endif 
     360  #if defined(USE_GPU) 
     361    #if defined(CALL_FQ) 
     362      double this_F2 = (pd_start == 0 ? 0.0 : result[2*q_index+0]); 
     363      double this_F1 = (pd_start == 0 ? 0.0 : result[2*q_index+1]); 
     364    #else 
     365      double this_F2 = (pd_start == 0 ? 0.0 : result[q_index]); 
     366    #endif 
     367  #else // !USE_GPU 
    347368    if (pd_start == 0) { 
    348369      #ifdef USE_OPENMP 
    349370      #pragma omp parallel for 
    350371      #endif 
    351       for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0; 
     372      #if defined(CALL_FQ) 
     373          // 2*nq for F^2,F pairs 
     374          for (int q_index=0; q_index < 2*nq; q_index++) result[q_index] = 0.0; 
     375      #else 
     376          for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0; 
     377      #endif 
    352378    } 
    353379    //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]); 
    354 #endif // !USE_OPENCL 
     380#endif // !USE_GPU 
    355381 
    356382 
     
    375401  const int n4 = pd_length[4]; 
    376402  const int p4 = pd_par[4]; 
    377   global const double *v4 = pd_value + pd_offset[4]; 
    378   global const double *w4 = pd_weight + pd_offset[4]; 
     403  pglobal const double *v4 = pd_value + pd_offset[4]; 
     404  pglobal const double *w4 = pd_weight + pd_offset[4]; 
    379405  int i4 = (pd_start/pd_stride[4])%n4;  // position in level 4 at pd_start 
    380406 
     
    411437          FETCH_Q         // set qx,qy from the q input vector 
    412438          APPLY_ROTATION  // convert qx,qy to qa,qb,qc 
    413           CALL_KERNEL     // scattering = Iqxy(qa, qb, qc, p1, p2, ...) 
     439          CALL_KERNEL     // F2 = Iqxy(qa, qb, qc, p1, p2, ...) 
    414440 
    415441      ++step;  // increment counter representing position in dispersity mesh 
     
    442468// inner loop and defines the macros that use them. 
    443469 
    444 #if defined(CALL_IQ) 
    445   // unoriented 1D 
     470 
     471#if defined(CALL_FQ) // COMPUTE_F1_F2 is true 
     472  // unoriented 1D returning <F> and <F^2> 
     473  // Note that F1 and F2 are returned from CALL_FQ by reference, and the 
     474  // user of the CALL_KERNEL macro below is assuming that F1 and F2 are defined. 
     475  double qk; 
     476  double F1, F2; 
     477  #define FETCH_Q() do { qk = q[q_index]; } while (0) 
     478  #define BUILD_ROTATION() do {} while(0) 
     479  #define APPLY_ROTATION() do {} while(0) 
     480  #define CALL_KERNEL() CALL_FQ(qk,F1,F2,local_values.table) 
     481 
     482#elif defined(CALL_FQ_A) 
     483  // unoriented 2D return <F> and <F^2> 
     484  // Note that the CALL_FQ_A macro is computing _F1_slot and _F2_slot by 
     485  // reference then returning _F2_slot.  We are calling them _F1_slot and 
     486  // _F2_slot here so they don't conflict with _F1 and _F2 in the macro 
     487  // expansion, or with the use of F2 = CALL_KERNEL() when it is used below. 
     488  double qx, qy; 
     489  double _F1_slot, _F2_slot; 
     490  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0) 
     491  #define BUILD_ROTATION() do {} while(0) 
     492  #define APPLY_ROTATION() do {} while(0) 
     493  #define CALL_KERNEL() CALL_FQ_A(sqrt(qx*qx+qy*qy),_F1_slot,_F2_slot,local_values.table) 
     494 
     495#elif defined(CALL_IQ) 
     496  // unoriented 1D return <F^2> 
    446497  double qk; 
    447498  #define FETCH_Q() do { qk = q[q_index]; } while (0) 
    448499  #define BUILD_ROTATION() do {} while(0) 
    449500  #define APPLY_ROTATION() do {} while(0) 
    450   #define CALL_KERNEL() CALL_IQ(qk, local_values.table) 
     501  #define CALL_KERNEL() CALL_IQ(qk,local_values.table) 
    451502 
    452503#elif defined(CALL_IQ_A) 
     
    482533  #define APPLY_ROTATION() qabc_apply(&rotation, qx, qy, &qa, &qb, &qc) 
    483534  #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table) 
     535 
    484536#elif defined(CALL_IQ_XY) 
    485537  // direct call to qx,qy calculator 
     
    495547// logic in the IQ_AC and IQ_ABC branches.  This will become more important 
    496548// if we implement more projections, or more complicated projections. 
    497 #if defined(CALL_IQ) || defined(CALL_IQ_A)  // no orientation 
     549#if defined(CALL_IQ) || defined(CALL_IQ_A) || defined(CALL_FQ) || defined(CALL_FQ_A) 
     550  // no orientation 
    498551  #define APPLY_PROJECTION() const double weight=weight0 
    499552#elif defined(CALL_IQ_XY) // pass orientation to the model 
     
    562615  const int n##_LOOP = details->pd_length[_LOOP]; \ 
    563616  const int p##_LOOP = details->pd_par[_LOOP]; \ 
    564   global const double *v##_LOOP = pd_value + details->pd_offset[_LOOP]; \ 
    565   global const double *w##_LOOP = pd_weight + details->pd_offset[_LOOP]; \ 
     617  pglobal const double *v##_LOOP = pd_value + details->pd_offset[_LOOP]; \ 
     618  pglobal const double *w##_LOOP = pd_weight + details->pd_offset[_LOOP]; \ 
    566619  int i##_LOOP = (pd_start/details->pd_stride[_LOOP])%n##_LOOP; 
    567620 
     
    587640// Pointers to the start of the dispersity and weight vectors, if needed. 
    588641#if MAX_PD>0 
    589   global const double *pd_value = values + NUM_VALUES; 
    590   global const double *pd_weight = pd_value + details->num_weights; 
     642  pglobal const double *pd_value = values + NUM_VALUES; 
     643  pglobal const double *pd_weight = pd_value + details->num_weights; 
    591644#endif 
    592645 
     
    645698    // Note: weight==0 must always be excluded 
    646699    if (weight > cutoff) { 
    647       pd_norm += weight * CALL_VOLUME(local_values.table); 
     700      double form, shell; 
     701      CALL_VOLUME(form, shell, local_values.table); 
     702      weight_norm += weight; 
     703      weighted_form += weight * form; 
     704      weighted_shell += weight * shell; 
     705      if (effective_radius_type != 0) { 
     706        weighted_radius += weight * CALL_EFFECTIVE_RADIUS(effective_radius_type, local_values.table); 
     707      } 
    648708      BUILD_ROTATION(); 
    649709 
    650 #ifndef USE_OPENCL 
     710#if !defined(USE_GPU) 
    651711      // DLL needs to explicitly loop over the q values. 
    652712      #ifdef USE_OPENMP 
     
    654714      #endif 
    655715      for (q_index=0; q_index<nq; q_index++) 
    656 #endif // !USE_OPENCL 
     716#endif // !USE_GPU 
    657717      { 
    658718 
     
    663723        #if defined(MAGNETIC) && NUM_MAGNETIC > 0 
    664724          // Compute the scattering from the magnetic cross sections. 
    665           double scattering = 0.0; 
     725          double F2 = 0.0; 
    666726          const double qsq = qx*qx + qy*qy; 
    667727          if (qsq > 1.e-16) { 
     
    688748//  q_index, qx, qy, xs, sk, local_values.vector[sld_index], px, py, mx, my, mz); 
    689749                } 
    690                 scattering += xs_weight * CALL_KERNEL(); 
     750                F2 += xs_weight * CALL_KERNEL(); 
    691751              } 
    692752            } 
    693753          } 
    694754        #else  // !MAGNETIC 
    695           const double scattering = CALL_KERNEL(); 
     755          #if defined(CALL_FQ) 
     756            CALL_KERNEL(); // sets F1 and F2 by reference 
     757          #else 
     758            const double F2 = CALL_KERNEL(); 
     759          #endif 
    696760        #endif // !MAGNETIC 
    697 //printf("q_index:%d %g %g %g %g\n", q_index, scattering, weight0); 
    698  
    699         #ifdef USE_OPENCL 
    700           this_result += weight * scattering; 
     761//printf("q_index:%d %g %g %g %g\n", q_index, F2, weight0); 
     762 
     763        #if defined(USE_GPU) 
     764          #if defined(CALL_FQ) 
     765            this_F2 += weight * F2; 
     766            this_F1 += weight * F1; 
     767          #else 
     768            this_F2 += weight * F2; 
     769          #endif 
    701770        #else // !USE_OPENCL 
    702           result[q_index] += weight * scattering; 
     771          #if defined(CALL_FQ) 
     772            result[2*q_index+0] += weight * F2; 
     773            result[2*q_index+1] += weight * F1; 
     774          #else 
     775            result[q_index] += weight * F2; 
     776          #endif 
    703777        #endif // !USE_OPENCL 
    704778      } 
    705779    } 
    706780  } 
    707  
    708781// close nested loops 
    709782++step; 
     
    724797#endif 
    725798 
    726 // Remember the current result and the updated norm. 
    727 #ifdef USE_OPENCL 
    728   result[q_index] = this_result; 
    729   if (q_index == 0) result[nq] = pd_norm; 
    730 //if (q_index == 0) printf("res: %g/%g\n", result[0], pd_norm); 
    731 #else // !USE_OPENCL 
    732   result[nq] = pd_norm; 
    733 //printf("res: %g/%g\n", result[0], pd_norm); 
    734 #endif // !USE_OPENCL 
     799// Remember the results and the updated norm. 
     800#if defined(USE_GPU) 
     801  #if defined(CALL_FQ) 
     802  result[2*q_index+0] = this_F2; 
     803  result[2*q_index+1] = this_F1; 
     804  #else 
     805  result[q_index] = this_F2; 
     806  #endif 
     807  if (q_index == 0) 
     808#endif 
     809  { 
     810#if defined(CALL_FQ) 
     811    result[2*nq] = weight_norm; 
     812    result[2*nq+1] = weighted_form; 
     813    result[2*nq+2] = weighted_shell; 
     814    result[2*nq+3] = weighted_radius; 
     815#else 
     816    result[nq] = weight_norm; 
     817    result[nq+1] = weighted_form; 
     818    result[nq+2] = weighted_shell; 
     819    result[nq+3] = weighted_radius; 
     820#endif 
     821  } 
    735822 
    736823// ** clear the macros in preparation for the next kernel ** 
  • sasmodels/kernelcl.py

    rd86f0fc r8064d5e  
    11""" 
    22GPU driver for C kernels 
     3 
     4TODO: docs are out of date 
    35 
    46There should be a single GPU environment running on the system.  This 
     
    5658import time 
    5759 
     60try: 
     61    from time import perf_counter as clock 
     62except ImportError: # CRUFT: python < 3.3 
     63    import sys 
     64    if sys.platform.count("darwin") > 0: 
     65        from time import time as clock 
     66    else: 
     67        from time import clock 
     68 
    5869import numpy as np  # type: ignore 
    5970 
    60  
    61 # Attempt to setup opencl. This may fail if the opencl package is not 
     71# Attempt to setup opencl. This may fail if the pyopencl package is not 
    6272# installed or if it is installed but there are no devices available. 
    6373try: 
     
    7484 
    7585from . import generate 
     86from .generate import F32, F64 
    7687from .kernel import KernelModel, Kernel 
    7788 
     
    131142 
    132143def use_opencl(): 
    133     return HAVE_OPENCL and os.environ.get("SAS_OPENCL", "").lower() != "none" 
     144    sas_opencl = os.environ.get("SAS_OPENCL", "OpenCL").lower() 
     145    return HAVE_OPENCL and sas_opencl != "none" and not sas_opencl.startswith("cuda") 
    134146 
    135147ENV = None 
     
    162174    Return true if device supports the requested precision. 
    163175    """ 
    164     if dtype == generate.F32: 
     176    if dtype == F32: 
    165177        return True 
    166     elif dtype == generate.F64: 
     178    elif dtype == F64: 
    167179        return "cl_khr_fp64" in device.extensions 
    168     elif dtype == generate.F16: 
    169         return "cl_khr_fp16" in device.extensions 
    170180    else: 
     181        # Not supporting F16 type since it isn't accurate enough 
    171182        return False 
    172183 
     
    179190        cl.kernel_work_group_info.PREFERRED_WORK_GROUP_SIZE_MULTIPLE, 
    180191        queue.device) 
    181  
    182 def _stretch_input(vector, dtype, extra=1e-3, boundary=32): 
    183     # type: (np.ndarray, np.dtype, float, int) -> np.ndarray 
    184     """ 
    185     Stretch an input vector to the correct boundary. 
    186  
    187     Performance on the kernels can drop by a factor of two or more if the 
    188     number of values to compute does not fall on a nice power of two 
    189     boundary.   The trailing additional vector elements are given a 
    190     value of *extra*, and so f(*extra*) will be computed for each of 
    191     them.  The returned array will thus be a subset of the computed array. 
    192  
    193     *boundary* should be a power of 2 which is at least 32 for good 
    194     performance on current platforms (as of Jan 2015).  It should 
    195     probably be the max of get_warp(kernel,queue) and 
    196     device.min_data_type_align_size//4. 
    197     """ 
    198     remainder = vector.size % boundary 
    199     if remainder != 0: 
    200         size = vector.size + (boundary - remainder) 
    201         vector = np.hstack((vector, [extra] * (size - vector.size))) 
    202     return np.ascontiguousarray(vector, dtype=dtype) 
    203  
    204192 
    205193def compile_model(context, source, dtype, fast=False): 
     
    239227    """ 
    240228    GPU context, with possibly many devices, and one queue per device. 
     229 
     230    Because the environment can be reset during a live program (e.g., if the 
     231    user changes the active GPU device in the GUI), everything associated 
     232    with the device context must be cached in the environment and recreated 
     233    if the environment changes.  The *cache* attribute is a simple dictionary 
     234    which holds keys and references to objects, such as compiled kernels and 
     235    allocated buffers.  The running program should check in the cache for 
     236    long lived objects and create them if they are not there.  The program 
     237    should not hold onto cached objects, but instead only keep them active 
     238    for the duration of a function call.  When the environment is destroyed 
     239    then the *release* method for each active cache item is called before 
     240    the environment is freed.  This means that each cl buffer should be 
     241    in its own cache entry. 
    241242    """ 
    242243    def __init__(self): 
    243244        # type: () -> None 
    244245        # find gpu context 
    245         #self.context = cl.create_some_context() 
    246  
    247         self.context = None 
    248         if 'SAS_OPENCL' in os.environ: 
    249             #Setting PYOPENCL_CTX as a SAS_OPENCL to create cl context 
    250             os.environ["PYOPENCL_CTX"] = os.environ["SAS_OPENCL"] 
    251         if 'PYOPENCL_CTX' in os.environ: 
    252             self._create_some_context() 
    253  
    254         if not self.context: 
    255             self.context = _get_default_context() 
     246        context_list = _create_some_context() 
     247 
     248        # Find a context for F32 and for F64 (maybe the same one). 
     249        # F16 isn't good enough. 
     250        self.context = {} 
     251        for dtype in (F32, F64): 
     252            for context in context_list: 
     253                if has_type(context.devices[0], dtype): 
     254                    self.context[dtype] = context 
     255                    break 
     256            else: 
     257                self.context[dtype] = None 
     258 
     259        # Build a queue for each context 
     260        self.queue = {} 
     261        context = self.context[F32] 
     262        self.queue[F32] = cl.CommandQueue(context, context.devices[0]) 
     263        if self.context[F64] == self.context[F32]: 
     264            self.queue[F64] = self.queue[F32] 
     265        else: 
     266            context = self.context[F64] 
     267            self.queue[F64] = cl.CommandQueue(context, context.devices[0]) 
    256268 
    257269        # Byte boundary for data alignment 
    258         #self.data_boundary = max(d.min_data_type_align_size 
    259         #                         for d in self.context.devices) 
    260         self.queues = [cl.CommandQueue(context, context.devices[0]) 
    261                        for context in self.context] 
     270        #self.data_boundary = max(context.devices[0].min_data_type_align_size 
     271        #                         for context in self.context.values()) 
     272 
     273        # Cache for compiled programs, and for items in context 
    262274        self.compiled = {} 
     275        self.cache = {} 
    263276 
    264277    def has_type(self, dtype): 
     
    267280        Return True if all devices support a given type. 
    268281        """ 
    269         return any(has_type(d, dtype) 
    270                    for context in self.context 
    271                    for d in context.devices) 
    272  
    273     def get_queue(self, dtype): 
    274         # type: (np.dtype) -> cl.CommandQueue 
    275         """ 
    276         Return a command queue for the kernels of type dtype. 
    277         """ 
    278         for context, queue in zip(self.context, self.queues): 
    279             if all(has_type(d, dtype) for d in context.devices): 
    280                 return queue 
    281  
    282     def get_context(self, dtype): 
    283         # type: (np.dtype) -> cl.Context 
    284         """ 
    285         Return a OpenCL context for the kernels of type dtype. 
    286         """ 
    287         for context in self.context: 
    288             if all(has_type(d, dtype) for d in context.devices): 
    289                 return context 
    290  
    291     def _create_some_context(self): 
    292         # type: () -> cl.Context 
    293         """ 
    294         Protected call to cl.create_some_context without interactivity.  Use 
    295         this if SAS_OPENCL is set in the environment.  Sets the *context* 
    296         attribute. 
    297         """ 
    298         try: 
    299             self.context = [cl.create_some_context(interactive=False)] 
    300         except Exception as exc: 
    301             warnings.warn(str(exc)) 
    302             warnings.warn("pyopencl.create_some_context() failed") 
    303             warnings.warn("the environment variable 'SAS_OPENCL' might not be set correctly") 
     282        return self.context.get(dtype, None) is not None 
    304283 
    305284    def compile_program(self, name, source, dtype, fast, timestamp): 
     
    318297            del self.compiled[key] 
    319298        if key not in self.compiled: 
    320             context = self.get_context(dtype) 
     299            context = self.context[dtype] 
    321300            logging.info("building %s for OpenCL %s", key, 
    322301                         context.devices[0].name.strip()) 
    323             program = compile_model(self.get_context(dtype), 
     302            program = compile_model(self.context[dtype], 
    324303                                    str(source), dtype, fast) 
    325304            self.compiled[key] = (program, timestamp) 
    326305        return program 
     306 
     307    def free_buffer(self, key): 
     308        if key in self.cache: 
     309            self.cache[key].release() 
     310            del self.cache[key] 
     311 
     312    def __del__(self): 
     313        for v in self.cache.values(): 
     314            release = getattr(v, 'release', lambda: None) 
     315            release() 
     316        self.cache = {} 
     317 
     318_CURRENT_ID = 0 
     319def unique_id(): 
     320    global _CURRENT_ID 
     321    _CURRENT_ID += 1 
     322    return _CURRENT_ID 
     323 
     324def _create_some_context(): 
     325    # type: () -> cl.Context 
     326    """ 
     327    Protected call to cl.create_some_context without interactivity. 
     328 
     329    Uses SAS_OPENCL or PYOPENCL_CTX if they are set in the environment, 
     330    otherwise scans for the most appropriate device using 
     331    :func:`_get_default_context`.  Ignore *SAS_OPENCL=OpenCL*, which 
     332    indicates that an OpenCL device should be used without specifying 
     333    which one (and not a CUDA device, or no GPU). 
     334    """ 
     335    # Assume we do not get here if SAS_OPENCL is None or CUDA 
     336    sas_opencl = os.environ.get('SAS_OPENCL', 'opencl') 
     337    if sas_opencl.lower() != 'opencl': 
     338        # Setting PYOPENCL_CTX as a SAS_OPENCL to create cl context 
     339        os.environ["PYOPENCL_CTX"] = sas_opencl 
     340 
     341    if 'PYOPENCL_CTX' in os.environ: 
     342        try: 
     343            return [cl.create_some_context(interactive=False)] 
     344        except Exception as exc: 
     345            warnings.warn(str(exc)) 
     346            warnings.warn("pyopencl.create_some_context() failed") 
     347            warnings.warn("the environment variable 'SAS_OPENCL' or 'PYOPENCL_CTX' might not be set correctly") 
     348 
     349    return _get_default_context() 
    327350 
    328351def _get_default_context(): 
     
    404427        self.dtype = dtype 
    405428        self.fast = fast 
    406         self.program = None # delay program creation 
    407         self._kernels = None 
     429        self.timestamp = generate.ocl_timestamp(self.info) 
     430        self._cache_key = unique_id() 
    408431 
    409432    def __getstate__(self): 
     
    414437        # type: (Tuple[ModelInfo, str, np.dtype, bool]) -> None 
    415438        self.info, self.source, self.dtype, self.fast = state 
    416         self.program = None 
    417439 
    418440    def make_kernel(self, q_vectors): 
    419441        # type: (List[np.ndarray]) -> "GpuKernel" 
    420         if self.program is None: 
    421             compile_program = environment().compile_program 
    422             timestamp = generate.ocl_timestamp(self.info) 
    423             self.program = compile_program( 
     442        return GpuKernel(self, q_vectors) 
     443 
     444    @property 
     445    def Iq(self): 
     446        return self._fetch_kernel('Iq') 
     447 
     448    def fetch_kernel(self, name): 
     449        # type: (str) -> cl.Kernel 
     450        """ 
     451        Fetch the kernel from the environment by name, compiling it if it 
     452        does not already exist. 
     453        """ 
     454        gpu = environment() 
     455        key = self._cache_key 
     456        if key not in gpu.cache: 
     457            program = gpu.compile_program( 
    424458                self.info.name, 
    425459                self.source['opencl'], 
    426460                self.dtype, 
    427461                self.fast, 
    428                 timestamp) 
     462                self.timestamp) 
    429463            variants = ['Iq', 'Iqxy', 'Imagnetic'] 
    430464            names = [generate.kernel_name(self.info, k) for k in variants] 
    431             kernels = [getattr(self.program, k) for k in names] 
    432             self._kernels = dict((k, v) for k, v in zip(variants, kernels)) 
    433         is_2d = len(q_vectors) == 2 
    434         if is_2d: 
    435             kernel = [self._kernels['Iqxy'], self._kernels['Imagnetic']] 
     465            kernels = [getattr(program, k) for k in names] 
     466            data = dict((k, v) for k, v in zip(variants, kernels)) 
     467            # keep a handle to program so GC doesn't collect 
     468            data['program'] = program 
     469            gpu.cache[key] = data 
    436470        else: 
    437             kernel = [self._kernels['Iq']]*2 
    438         return GpuKernel(kernel, self.dtype, self.info, q_vectors) 
    439  
    440     def release(self): 
    441         # type: () -> None 
    442         """ 
    443         Free the resources associated with the model. 
    444         """ 
    445         if self.program is not None: 
    446             self.program = None 
    447  
    448     def __del__(self): 
    449         # type: () -> None 
    450         self.release() 
     471            data = gpu.cache[key] 
     472        return data[name] 
    451473 
    452474# TODO: check that we don't need a destructor for buffers which go out of scope 
     
    473495        # type: (List[np.ndarray], np.dtype) -> None 
    474496        # TODO: do we ever need double precision q? 
    475         env = environment() 
    476497        self.nq = q_vectors[0].size 
    477498        self.dtype = np.dtype(dtype) 
     
    482503        # architectures tested so far. 
    483504        if self.is_2d: 
    484             # Note: 16 rather than 15 because result is 1 longer than input. 
    485             width = ((self.nq+16)//16)*16 
     505            width = ((self.nq+15)//16)*16 
    486506            self.q = np.empty((width, 2), dtype=dtype) 
    487507            self.q[:self.nq, 0] = q_vectors[0] 
    488508            self.q[:self.nq, 1] = q_vectors[1] 
    489509        else: 
    490             # Note: 32 rather than 31 because result is 1 longer than input. 
    491             width = ((self.nq+32)//32)*32 
     510            width = ((self.nq+31)//32)*32 
    492511            self.q = np.empty(width, dtype=dtype) 
    493512            self.q[:self.nq] = q_vectors[0] 
    494513        self.global_size = [self.q.shape[0]] 
    495         context = env.get_context(self.dtype) 
    496         #print("creating inputs of size", self.global_size) 
    497         self.q_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
    498                              hostbuf=self.q) 
     514        self._cache_key = unique_id() 
     515 
     516    @property 
     517    def q_b(self): 
     518        """Lazy creation of q buffer so it can survive context reset""" 
     519        env = environment() 
     520        key = self._cache_key 
     521        if key not in env.cache: 
     522            context = env.context[self.dtype] 
     523            #print("creating inputs of size", self.global_size) 
     524            buffer = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
     525                               hostbuf=self.q) 
     526            env.cache[key] = buffer 
     527        return env.cache[key] 
    499528 
    500529    def release(self): 
    501530        # type: () -> None 
    502531        """ 
    503         Free the memory. 
    504         """ 
    505         if self.q_b is not None: 
    506             self.q_b.release() 
    507             self.q_b = None 
     532        Free the buffer associated with the q value 
     533        """ 
     534        environment().free_buffer(id(self)) 
    508535 
    509536    def __del__(self): 
     
    515542    Callable SAS kernel. 
    516543 
    517     *kernel* is the GpuKernel object to call 
    518  
    519     *model_info* is the module information 
    520  
    521     *q_vectors* is the q vectors at which the kernel should be evaluated 
     544    *model* is the GpuModel object to call 
     545 
     546    The following attributes are defined: 
     547 
     548    *info* is the module information 
    522549 
    523550    *dtype* is the kernel precision 
     551 
     552    *dim* is '1d' or '2d' 
     553 
     554    *result* is a vector to contain the results of the call 
    524555 
    525556    The resulting call method takes the *pars*, a list of values for 
     
    531562    Call :meth:`release` when done with the kernel instance. 
    532563    """ 
    533     def __init__(self, kernel, dtype, model_info, q_vectors): 
     564    def __init__(self, model, q_vectors): 
    534565        # type: (cl.Kernel, np.dtype, ModelInfo, List[np.ndarray]) -> None 
    535         q_input = GpuInput(q_vectors, dtype) 
    536         self.kernel = kernel 
    537         self.info = model_info 
    538         self.dtype = dtype 
    539         self.dim = '2d' if q_input.is_2d else '1d' 
    540         # plus three for the normalization values 
    541         self.result = np.empty(q_input.nq+1, dtype) 
    542  
    543         # Inputs and outputs for each kernel call 
    544         # Note: res may be shorter than res_b if global_size != nq 
     566        dtype = model.dtype 
     567        self.q_input = GpuInput(q_vectors, dtype) 
     568        self._model = model 
     569        # F16 isn't sufficient, so don't support it 
     570        self._as_dtype = np.float64 if dtype == generate.F64 else np.float32 
     571        self._cache_key = unique_id() 
     572 
     573        # attributes accessed from the outside 
     574        self.dim = '2d' if self.q_input.is_2d else '1d' 
     575        self.info = model.info 
     576        self.dtype = model.dtype 
     577 
     578        # holding place for the returned value 
     579        nout = 2 if self.info.have_Fq and self.dim == '1d' else 1 
     580        extra_q = 4  # total weight, form volume, shell volume and R_eff 
     581        self.result = np.empty(self.q_input.nq*nout+extra_q, dtype) 
     582 
     583    @property 
     584    def _result_b(self): 
     585        """Lazy creation of result buffer so it can survive context reset""" 
    545586        env = environment() 
    546         self.queue = env.get_queue(dtype) 
    547  
    548         self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 
    549                                   q_input.global_size[0] * dtype.itemsize) 
    550         self.q_input = q_input # allocated by GpuInput above 
    551  
    552         self._need_release = [self.result_b, self.q_input] 
    553         self.real = (np.float32 if dtype == generate.F32 
    554                      else np.float64 if dtype == generate.F64 
    555                      else np.float16 if dtype == generate.F16 
    556                      else np.float32)  # will never get here, so use np.float32 
    557  
    558     def __call__(self, call_details, values, cutoff, magnetic): 
     587        key = self._cache_key 
     588        if key not in env.cache: 
     589            context = env.context[self.dtype] 
     590            width = ((self.result.size+31)//32)*32 * self.dtype.itemsize 
     591            buffer = cl.Buffer(context, mf.READ_WRITE, width) 
     592            env.cache[key] = buffer 
     593        return env.cache[key] 
     594 
     595    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
    559596        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    560         context = self.queue.context 
    561         # Arrange data transfer to card 
     597        env = environment() 
     598        queue = env.queue[self._model.dtype] 
     599        context = queue.context 
     600 
     601        # Arrange data transfer to/from card 
     602        q_b = self.q_input.q_b 
     603        result_b = self._result_b 
    562604        details_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
    563605                              hostbuf=call_details.buffer) 
     
    565607                             hostbuf=values) 
    566608 
    567         kernel = self.kernel[1 if magnetic else 0] 
    568         args = [ 
     609        name = 'Iq' if self.dim == '1d' else 'Imagnetic' if magnetic else 'Iqxy' 
     610        kernel = self._model.fetch_kernel(name) 
     611        kernel_args = [ 
    569612            np.uint32(self.q_input.nq), None, None, 
    570             details_b, values_b, self.q_input.q_b, self.result_b, 
    571             self.real(cutoff), 
     613            details_b, values_b, q_b, result_b, 
     614            self._as_dtype(cutoff), 
     615            np.uint32(effective_radius_type), 
    572616        ] 
    573617        #print("Calling OpenCL") 
    574618        #call_details.show(values) 
    575         # Call kernel and retrieve results 
     619        #Call kernel and retrieve results 
    576620        wait_for = None 
    577         last_nap = time.clock() 
     621        last_nap = clock() 
    578622        step = 1000000//self.q_input.nq + 1 
    579623        for start in range(0, call_details.num_eval, step): 
    580624            stop = min(start + step, call_details.num_eval) 
    581625            #print("queuing",start,stop) 
    582             args[1:3] = [np.int32(start), np.int32(stop)] 
    583             wait_for = [kernel(self.queue, self.q_input.global_size, None, 
    584                                *args, wait_for=wait_for)] 
     626            kernel_args[1:3] = [np.int32(start), np.int32(stop)] 
     627            wait_for = [kernel(queue, self.q_input.global_size, None, 
     628                               *kernel_args, wait_for=wait_for)] 
    585629            if stop < call_details.num_eval: 
    586630                # Allow other processes to run 
    587631                wait_for[0].wait() 
    588                 current_time = time.clock() 
     632                current_time = clock() 
    589633                if current_time - last_nap > 0.5: 
    590                     time.sleep(0.05) 
     634                    time.sleep(0.001) 
    591635                    last_nap = current_time 
    592         cl.enqueue_copy(self.queue, self.result, self.result_b) 
     636        cl.enqueue_copy(queue, self.result, result_b, wait_for=wait_for) 
    593637        #print("result", self.result) 
    594638 
     
    598642                v.release() 
    599643 
    600         pd_norm = self.result[self.q_input.nq] 
    601         scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
    602         background = values[1] 
    603         #print("scale",scale,values[0],self.result[self.q_input.nq],background) 
    604         return scale*self.result[:self.q_input.nq] + background 
    605         # return self.result[:self.q_input.nq] 
    606  
    607644    def release(self): 
    608645        # type: () -> None 
     
    610647        Release resources associated with the kernel. 
    611648        """ 
    612         for v in self._need_release: 
    613             v.release() 
    614         self._need_release = [] 
     649        environment().free_buffer(id(self)) 
     650        self.q_input.release() 
    615651 
    616652    def __del__(self): 
  • sasmodels/kerneldll.py

    r1a3559f re44432d  
    9999    pass 
    100100# pylint: enable=unused-import 
     101 
     102# Compiler output is a byte stream that needs to be decode in python 3 
     103decode = (lambda s: s) if sys.version_info[0] < 3 else (lambda s: s.decode('utf8')) 
    101104 
    102105if "SAS_DLL_PATH" in os.environ: 
     
    184187        subprocess.check_output(command, shell=shell, stderr=subprocess.STDOUT) 
    185188    except subprocess.CalledProcessError as exc: 
    186         raise RuntimeError("compile failed.\n%s\n%s"%(command_str, exc.output)) 
     189        output = decode(exc.output) 
     190        raise RuntimeError("compile failed.\n%s\n%s"%(command_str, output)) 
    187191    if not os.path.exists(output): 
    188192        raise RuntimeError("compile failed.  File is in %r"%source) 
     
    315319 
    316320        # int, int, int, int*, double*, double*, double*, double*, double 
    317         argtypes = [ct.c_int32]*3 + [ct.c_void_p]*4 + [float_type] 
     321        argtypes = [ct.c_int32]*3 + [ct.c_void_p]*4 + [float_type, ct.c_int32] 
    318322        names = [generate.kernel_name(self.info, variant) 
    319323                 for variant in ("Iq", "Iqxy", "Imagnetic")] 
     
    375379    def __init__(self, kernel, model_info, q_input): 
    376380        # type: (Callable[[], np.ndarray], ModelInfo, PyInput) -> None 
     381        #,model_info,q_input) 
    377382        self.kernel = kernel 
    378383        self.info = model_info 
     
    380385        self.dtype = q_input.dtype 
    381386        self.dim = '2d' if q_input.is_2d else '1d' 
    382         self.result = np.empty(q_input.nq+1, q_input.dtype) 
     387        # leave room for f1/f2 results in case we need to compute beta for 1d models 
     388        nout = 2 if self.info.have_Fq else 1 
     389        # +4 for total weight, shell volume, effective radius, form volume 
     390        self.result = np.empty(q_input.nq*nout + 4, self.dtype) 
    383391        self.real = (np.float32 if self.q_input.dtype == generate.F32 
    384392                     else np.float64 if self.q_input.dtype == generate.F64 
    385393                     else np.float128) 
    386394 
    387     def __call__(self, call_details, values, cutoff, magnetic): 
    388         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    389  
     395    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
     396        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
    390397        kernel = self.kernel[1 if magnetic else 0] 
    391398        args = [ 
     
    394401            None, # pd_stop pd_stride[MAX_PD] 
    395402            call_details.buffer.ctypes.data, # problem 
    396             values.ctypes.data,  #pars 
    397             self.q_input.q.ctypes.data, #q 
     403            values.ctypes.data,  # pars 
     404            self.q_input.q.ctypes.data, # q 
    398405            self.result.ctypes.data,   # results 
    399406            self.real(cutoff), # cutoff 
     407            effective_radius_type, # cutoff 
    400408        ] 
    401409        #print("Calling DLL") 
     
    407415            kernel(*args) # type: ignore 
    408416 
    409         #print("returned",self.q_input.q, self.result) 
    410         pd_norm = self.result[self.q_input.nq] 
    411         scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
    412         background = values[1] 
    413         #print("scale",scale,background) 
    414         return scale*self.result[:self.q_input.nq] + background 
    415  
    416417    def release(self): 
    417418        # type: () -> None 
  • sasmodels/kernelpy.py

    r12eec1e raa8c6e0  
    1212 
    1313import numpy as np  # type: ignore 
     14from numpy import pi 
     15try: 
     16    from numpy import cbrt 
     17except ImportError: 
     18    def cbrt(x): return x ** (1.0/3.0) 
    1419 
    1520from .generate import F64 
     
    158163 
    159164        # Generate a closure which calls the form_volume if it exists. 
    160         form_volume = model_info.form_volume 
    161         self._volume = ((lambda: form_volume(*volume_args)) if form_volume else 
    162                         (lambda: 1.0)) 
    163  
    164     def __call__(self, call_details, values, cutoff, magnetic): 
     165        self._volume_args = volume_args 
     166        volume = model_info.form_volume 
     167        shell = model_info.shell_volume 
     168        radius = model_info.effective_radius 
     169        self._volume = ((lambda: (shell(*volume_args), volume(*volume_args))) if shell and volume 
     170                        else (lambda: [volume(*volume_args)]*2) if volume 
     171                        else (lambda: (1.0, 1.0))) 
     172        self._radius = ((lambda mode: radius(mode, *volume_args)) if radius 
     173                        else (lambda mode: cbrt(0.75/pi*volume(*volume_args))) if volume 
     174                        else (lambda mode: 1.0)) 
     175 
     176 
     177 
     178    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
    165179        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    166180        if magnetic: 
     
    168182        #print("Calling python kernel") 
    169183        #call_details.show(values) 
    170         res = _loops(self._parameter_vector, self._form, self._volume, 
    171                      self.q_input.nq, call_details, values, cutoff) 
    172         return res 
     184        radius = ((lambda: 0.0) if effective_radius_type == 0 
     185                  else (lambda: self._radius(effective_radius_type))) 
     186        self.result = _loops( 
     187            self._parameter_vector, self._form, self._volume, radius, 
     188            self.q_input.nq, call_details, values, cutoff) 
    173189 
    174190    def release(self): 
     
    183199           form,          # type: Callable[[], np.ndarray] 
    184200           form_volume,   # type: Callable[[], float] 
     201           form_radius,   # type: Callable[[], float] 
    185202           nq,            # type: int 
    186203           call_details,  # type: details.CallDetails 
    187204           values,        # type: np.ndarray 
    188            cutoff         # type: float 
     205           cutoff,        # type: float 
    189206          ): 
    190207    # type: (...) -> None 
     
    198215    #                                                              # 
    199216    ################################################################ 
     217 
     218    # WARNING: Trickery ahead 
     219    # The parameters[] vector is embedded in the closures for form(), 
     220    # form_volume() and form_radius().  We set the initial vector from 
     221    # the values for the model parameters. As we loop through the polydispesity 
     222    # mesh, we update the components with the polydispersity values before 
     223    # calling the respective functions. 
    200224    n_pars = len(parameters) 
    201225    parameters[:] = values[2:n_pars+2] 
     226 
    202227    if call_details.num_active == 0: 
    203         pd_norm = float(form_volume()) 
    204         scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
    205         background = values[1] 
    206         return scale*form() + background 
    207  
    208     pd_value = values[2+n_pars:2+n_pars + call_details.num_weights] 
    209     pd_weight = values[2+n_pars + call_details.num_weights:] 
    210  
    211     pd_norm = 0.0 
    212     partial_weight = np.NaN 
    213     weight = np.NaN 
    214  
    215     p0_par = call_details.pd_par[0] 
    216     p0_length = call_details.pd_length[0] 
    217     p0_index = p0_length 
    218     p0_offset = call_details.pd_offset[0] 
    219  
    220     pd_par = call_details.pd_par[:call_details.num_active] 
    221     pd_offset = call_details.pd_offset[:call_details.num_active] 
    222     pd_stride = call_details.pd_stride[:call_details.num_active] 
    223     pd_length = call_details.pd_length[:call_details.num_active] 
    224  
    225     total = np.zeros(nq, 'd') 
    226     for loop_index in range(call_details.num_eval): 
    227         # update polydispersity parameter values 
    228         if p0_index == p0_length: 
    229             pd_index = (loop_index//pd_stride)%pd_length 
    230             parameters[pd_par] = pd_value[pd_offset+pd_index] 
    231             partial_weight = np.prod(pd_weight[pd_offset+pd_index][1:]) 
    232             p0_index = loop_index%p0_length 
    233  
    234         weight = partial_weight * pd_weight[p0_offset + p0_index] 
    235         parameters[p0_par] = pd_value[p0_offset + p0_index] 
    236         p0_index += 1 
    237         if weight > cutoff: 
    238             # Call the scattering function 
    239             # Assume that NaNs are only generated if the parameters are bad; 
    240             # exclude all q for that NaN.  Even better would be to have an 
    241             # INVALID expression like the C models, but that is too expensive. 
    242             Iq = np.asarray(form(), 'd') 
    243             if np.isnan(Iq).any(): 
    244                 continue 
    245  
    246             # update value and norm 
    247             total += weight * Iq 
    248             pd_norm += weight * form_volume() 
    249  
    250     scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
    251     background = values[1] 
    252     return scale*total + background 
     228        total = form() 
     229        weight_norm = 1.0 
     230        weighted_shell, weighted_form = form_volume() 
     231        weighted_radius = form_radius() 
     232 
     233    else: 
     234        pd_value = values[2+n_pars:2+n_pars + call_details.num_weights] 
     235        pd_weight = values[2+n_pars + call_details.num_weights:] 
     236 
     237        weight_norm = 0.0 
     238        weighted_form = 0.0 
     239        weighted_shell = 0.0 
     240        weighted_radius = 0.0 
     241        partial_weight = np.NaN 
     242        weight = np.NaN 
     243 
     244        p0_par = call_details.pd_par[0] 
     245        p0_length = call_details.pd_length[0] 
     246        p0_index = p0_length 
     247        p0_offset = call_details.pd_offset[0] 
     248 
     249        pd_par = call_details.pd_par[:call_details.num_active] 
     250        pd_offset = call_details.pd_offset[:call_details.num_active] 
     251        pd_stride = call_details.pd_stride[:call_details.num_active] 
     252        pd_length = call_details.pd_length[:call_details.num_active] 
     253 
     254        total = np.zeros(nq, 'd') 
     255        for loop_index in range(call_details.num_eval): 
     256            # update polydispersity parameter values 
     257            if p0_index == p0_length: 
     258                pd_index = (loop_index//pd_stride)%pd_length 
     259                parameters[pd_par] = pd_value[pd_offset+pd_index] 
     260                partial_weight = np.prod(pd_weight[pd_offset+pd_index][1:]) 
     261                p0_index = loop_index%p0_length 
     262 
     263            weight = partial_weight * pd_weight[p0_offset + p0_index] 
     264            parameters[p0_par] = pd_value[p0_offset + p0_index] 
     265            p0_index += 1 
     266            if weight > cutoff: 
     267                # Call the scattering function 
     268                # Assume that NaNs are only generated if the parameters are bad; 
     269                # exclude all q for that NaN.  Even better would be to have an 
     270                # INVALID expression like the C models, but that is expensive. 
     271                Iq = np.asarray(form(), 'd') 
     272                if np.isnan(Iq).any(): 
     273                    continue 
     274 
     275                # update value and norm 
     276                total += weight * Iq 
     277                weight_norm += weight 
     278                shell, form = form_volume() 
     279                weighted_form += weight * form 
     280                weighted_shell += weight * shell 
     281                weighted_radius += weight * form_radius() 
     282 
     283    result = np.hstack((total, weight_norm, weighted_form, weighted_shell, weighted_radius)) 
     284    return result 
    253285 
    254286 
  • sasmodels/mixture.py

    rb171acd r39a06c9  
    143143    #model_info.tests = [] 
    144144    #model_info.source = [] 
    145     # Iq, Iqxy, form_volume, ER, VR and sesans 
    146145    # Remember the component info blocks so we can build the model 
    147146    model_info.composition = ('mixture', parts) 
  • sasmodels/model_test.py

    r12eec1e r5024a56  
    55Usage:: 
    66 
    7     python -m sasmodels.model_test [opencl|dll|opencl_and_dll] model1 model2 ... 
     7    python -m sasmodels.model_test [opencl|cuda|dll] model1 model2 ... 
    88 
    99    if model1 is 'all', then all except the remaining models will be tested 
    1010 
    1111Each model is tested using the default parameters at q=0.1, (qx, qy)=(0.1, 0.1), 
    12 and the ER and VR are computed.  The return values at these points are not 
    13 considered.  The test is only to verify that the models run to completion, 
    14 and do not produce inf or NaN. 
     12and Fq is called to make sure R_eff, volume and volume ratio are computed. 
     13The return values at these points are not considered.  The test is only to 
     14verify that the models run to completion, and do not produce inf or NaN. 
    1515 
    1616Tests are defined with the *tests* attribute in the model.py file.  *tests* 
    1717is a list of individual tests to run, where each test consists of the 
    1818parameter values for the test, the q-values and the expected results.  For 
    19 the effective radius test, the q-value should be 'ER'.  For the VR test, 
    20 the q-value should be 'VR'.  For 1-D tests, either specify the q value or 
    21 a list of q-values, and the corresponding I(q) value, or list of I(q) values. 
     19the effective radius test and volume ratio tests, use the extended output 
     20form, which checks each output of kernel.Fq. For 1-D tests, either specify 
     21the q value or a list of q-values, and the corresponding I(q) value, or 
     22list of I(q) values. 
    2223 
    2324That is:: 
     
    3233                        [I(qx1, qy1), I(qx2, qy2), ...]], 
    3334 
    34         [ {parameters}, 'ER', ER(pars) ], 
    35         [ {parameters}, 'VR', VR(pars) ], 
     35        [ {parameters}, q, F(q), F^2(q), R_eff, V, V_r ], 
    3636        ... 
    3737    ] 
     
    6060from . import core 
    6161from .core import list_models, load_model_info, build_model 
    62 from .direct_model import call_kernel, call_ER, call_VR 
     62from .direct_model import call_kernel, call_Fq 
    6363from .exception import annotate_exception 
    6464from .modelinfo import expand_pars 
    6565from .kernelcl import use_opencl 
     66from .kernelcuda import use_cuda 
     67from . import product 
    6668 
    6769# pylint: disable=unused-import 
     
    8082    Construct the pyunit test suite. 
    8183 
    82     *loaders* is the list of kernel drivers to use, which is one of 
    83     *["dll", "opencl"]*, *["dll"]* or *["opencl"]*.  For python models, 
    84     the python driver is always used. 
     84    *loaders* is the list of kernel drivers to use (dll, opencl or cuda). 
     85    For python model the python driver is always used. 
    8586 
    8687    *models* is the list of models to test, or *["all"]* to test all models. 
     
    160161                                    test_method_name, 
    161162                                    platform="ocl", dtype=None, 
     163                                    stash=stash) 
     164            #print("defining", test_name) 
     165            suite.addTest(test) 
     166 
     167        # test using cuda if desired and available 
     168        if 'cuda' in loaders and use_cuda(): 
     169            test_name = "%s-cuda"%model_name 
     170            test_method_name = "test_%s_cuda" % model_info.id 
     171            # Using dtype=None so that the models that are only 
     172            # correct for double precision are not tested using 
     173            # single precision.  The choice is determined by the 
     174            # presence of *single=False* in the model file. 
     175            test = ModelTestCase(test_name, model_info, 
     176                                    test_method_name, 
     177                                    platform="cuda", dtype=None, 
    162178                                    stash=stash) 
    163179            #print("defining", test_name) 
     
    202218                ({}, [0.001, 0.01, 0.1], [None]*3), 
    203219                ({}, [(0.1, 0.1)]*2, [None]*2), 
    204                 # test that ER/VR will run if they exist 
    205                 ({}, 'ER', None), 
    206                 ({}, 'VR', None), 
     220                # test that Fq will run, and return R_eff, V, V_r 
     221                ({}, 0.1, None, None, None, None, None), 
    207222                ] 
    208223            tests = smoke_tests 
     
    210225            if self.info.tests is not None: 
    211226                tests += self.info.tests 
     227            S_tests = [test for test in tests if '@S' in test[0]] 
     228            P_tests = [test for test in tests if '@S' not in test[0]] 
    212229            try: 
    213230                model = build_model(self.info, dtype=self.dtype, 
    214231                                    platform=self.platform) 
    215                 results = [self.run_one(model, test) for test in tests] 
     232                results = [self.run_one(model, test) for test in P_tests] 
     233                for test in S_tests: 
     234                    # pull the S model name out of the test defn 
     235                    pars = test[0].copy() 
     236                    s_name = pars.pop('@S') 
     237                    ps_test = [pars] + list(test[1:]) 
     238                    # build the P@S model 
     239                    s_info = load_model_info(s_name) 
     240                    ps_info = product.make_product_info(self.info, s_info) 
     241                    ps_model = build_model(ps_info, dtype=self.dtype, 
     242                                           platform=self.platform) 
     243                    # run the tests 
     244                    results.append(self.run_one(ps_model, ps_test)) 
     245 
    216246                if self.stash: 
    217247                    for test, target, actual in zip(tests, self.stash[0], results): 
     
    223253 
    224254                # Check for missing tests.  Only do so for the "dll" tests 
    225                 # to reduce noise from both opencl and dll, and because 
     255                # to reduce noise from both opencl and cuda, and because 
    226256                # python kernels use platform="dll". 
    227257                if self.platform == "dll": 
     
    238268        def _find_missing_tests(self): 
    239269            # type: () -> None 
    240             """make sure there are 1D, 2D, ER and VR tests as appropriate""" 
    241             model_has_VR = callable(self.info.VR) 
    242             model_has_ER = callable(self.info.ER) 
     270            """make sure there are 1D and 2D tests as appropriate""" 
    243271            model_has_1D = True 
    244272            model_has_2D = any(p.type == 'orientation' 
     
    248276            single = [test for test in self.info.tests 
    249277                      if not isinstance(test[2], list) and test[2] is not None] 
    250             tests_has_VR = any(test[1] == 'VR' for test in single) 
    251             tests_has_ER = any(test[1] == 'ER' for test in single) 
    252278            tests_has_1D_single = any(isinstance(test[1], float) for test in single) 
    253279            tests_has_2D_single = any(isinstance(test[1], tuple) for test in single) 
     
    262288 
    263289            missing = [] 
    264             if model_has_VR and not tests_has_VR: 
    265                 missing.append("VR") 
    266             if model_has_ER and not tests_has_ER: 
    267                 missing.append("ER") 
    268290            if model_has_1D and not (tests_has_1D_single or tests_has_1D_multiple): 
    269291                missing.append("1D") 
     
    276298            # type: (KernelModel, TestCondition) -> None 
    277299            """Run a single test case.""" 
    278             user_pars, x, y = test 
     300            user_pars, x, y = test[:3] 
    279301            pars = expand_pars(self.info.parameters, user_pars) 
    280302            invalid = invalid_pars(self.info.parameters, pars) 
     
    289311            self.assertEqual(len(y), len(x)) 
    290312 
    291             if x[0] == 'ER': 
    292                 actual = np.array([call_ER(model.info, pars)]) 
    293             elif x[0] == 'VR': 
    294                 actual = np.array([call_VR(model.info, pars)]) 
    295             elif isinstance(x[0], tuple): 
     313            if isinstance(x[0], tuple): 
    296314                qx, qy = zip(*x) 
    297315                q_vectors = [np.array(qx), np.array(qy)] 
    298                 kernel = model.make_kernel(q_vectors) 
    299                 actual = call_kernel(kernel, pars) 
    300316            else: 
    301317                q_vectors = [np.array(x)] 
    302                 kernel = model.make_kernel(q_vectors) 
     318 
     319            kernel = model.make_kernel(q_vectors) 
     320            if len(test) == 3: 
    303321                actual = call_kernel(kernel, pars) 
    304  
    305             self.assertTrue(len(actual) > 0) 
    306             self.assertEqual(len(y), len(actual)) 
    307  
    308             for xi, yi, actual_yi in zip(x, y, actual): 
     322                self._check_vectors(x, y, actual, 'I') 
     323                return actual 
     324            else: 
     325                y1 = y 
     326                y2 = test[3] if not isinstance(test[3], list) else [test[3]] 
     327                F1, F2, R_eff, volume, volume_ratio = call_Fq(kernel, pars) 
     328                if F1 is not None:  # F1 is none for models with Iq instead of Fq 
     329                    self._check_vectors(x, y1, F1, 'F') 
     330                self._check_vectors(x, y2, F2, 'F^2') 
     331                self._check_scalar(test[4], R_eff, 'R_eff') 
     332                self._check_scalar(test[5], volume, 'volume') 
     333                self._check_scalar(test[6], volume_ratio, 'form:shell ratio') 
     334                return F2 
     335 
     336        def _check_scalar(self, target, actual, name): 
     337            if target is None: 
     338                # smoke test --- make sure it runs and produces a value 
     339                self.assertTrue(not np.isnan(actual), 
     340                                'invalid %s: %s' % (name, actual)) 
     341            elif np.isnan(target): 
     342                # make sure nans match 
     343                self.assertTrue(np.isnan(actual), 
     344                                '%s: expected:%s; actual:%s' 
     345                                % (name, target, actual)) 
     346            else: 
     347                # is_near does not work for infinite values, so also test 
     348                # for exact values. 
     349                self.assertTrue(target == actual or is_near(target, actual, 5), 
     350                                '%s: expected:%s; actual:%s' 
     351                                % (name, target, actual)) 
     352 
     353        def _check_vectors(self, x, target, actual, name='I'): 
     354            self.assertTrue(len(actual) > 0, 
     355                            '%s(...) expected return'%name) 
     356            if target is None: 
     357                return 
     358            self.assertEqual(len(target), len(actual), 
     359                             '%s(...) returned wrong length'%name) 
     360            for xi, yi, actual_yi in zip(x, target, actual): 
    309361                if yi is None: 
    310362                    # smoke test --- make sure it runs and produces a value 
    311363                    self.assertTrue(not np.isnan(actual_yi), 
    312                                     'invalid f(%s): %s' % (xi, actual_yi)) 
     364                                    'invalid %s(%s): %s' % (name, xi, actual_yi)) 
    313365                elif np.isnan(yi): 
     366                    # make sure nans match 
    314367                    self.assertTrue(np.isnan(actual_yi), 
    315                                     'f(%s): expected:%s; actual:%s' 
    316                                     % (xi, yi, actual_yi)) 
     368                                    '%s(%s): expected:%s; actual:%s' 
     369                                    % (name, xi, yi, actual_yi)) 
    317370                else: 
    318371                    # is_near does not work for infinite values, so also test 
    319                     # for exact values.  Note that this will not 
     372                    # for exact values. 
    320373                    self.assertTrue(yi == actual_yi or is_near(yi, actual_yi, 5), 
    321                                     'f(%s); expected:%s; actual:%s' 
    322                                     % (xi, yi, actual_yi)) 
    323             return actual 
     374                                    '%s(%s); expected:%s; actual:%s' 
     375                                    % (name, xi, yi, actual_yi)) 
    324376 
    325377    return ModelTestCase 
     
    333385    invalid = [] 
    334386    for par in sorted(pars.keys()): 
     387        # special handling of R_eff mode, which is not a usual parameter 
     388        if par == product.RADIUS_MODE_ID: 
     389            continue 
    335390        parts = par.split('_pd') 
    336391        if len(parts) > 1 and parts[1] not in ("", "_n", "nsigma", "type"): 
     
    383438 
    384439    # Build a test suite containing just the model 
    385     loaders = ['opencl'] if use_opencl() else ['dll'] 
     440    loaders = ['opencl' if use_opencl() else 'cuda' if use_cuda() else 'dll'] 
    386441    suite = unittest.TestSuite() 
    387442    _add_model_to_suite(loaders, suite, model_info) 
     
    444499        loaders = ['opencl'] 
    445500        models = models[1:] 
     501    elif models and models[0] == 'cuda': 
     502        if not use_cuda(): 
     503            print("cuda is not available") 
     504            return 1 
     505        loaders = ['cuda'] 
     506        models = models[1:] 
    446507    elif models and models[0] == 'dll': 
    447508        # TODO: test if compiler is available? 
    448509        loaders = ['dll'] 
    449510        models = models[1:] 
    450     elif models and models[0] == 'opencl_and_dll': 
    451         loaders = ['opencl', 'dll'] if use_opencl() else ['dll'] 
    452         models = models[1:] 
    453511    else: 
    454         loaders = ['opencl', 'dll'] if use_opencl() else ['dll'] 
     512        loaders = ['dll'] 
     513        if use_opencl(): 
     514            loaders.append('opencl') 
     515        if use_cuda(): 
     516            loaders.append('cuda') 
    455517    if not models: 
    456518        print("""\ 
    457519usage: 
    458   python -m sasmodels.model_test [-v] [opencl|dll] model1 model2 ... 
     520  python -m sasmodels.model_test [-v] [opencl|cuda|dll] model1 model2 ... 
    459521 
    460522If -v is included on the command line, then use verbose output. 
    461523 
    462 If neither opencl nor dll is specified, then models will be tested with 
    463 both OpenCL and dll; the compute target is ignored for pure python models. 
     524If no platform is specified, then models will be tested with dll, and 
     525if available, OpenCL and CUDA; the compute target is ignored for pure python models. 
    464526 
    465527If model1 is 'all', then all except the remaining models will be tested. 
     
    481543    Run "nosetests sasmodels" on the command line to invoke it. 
    482544    """ 
    483     loaders = ['opencl', 'dll'] if use_opencl() else ['dll'] 
     545    loaders = ['dll'] 
     546    if use_opencl(): 
     547        loaders.append('opencl') 
     548    if use_cuda(): 
     549        loaders.append('cuda') 
    484550    tests = make_suite(loaders, ['all']) 
    485551    def build_test(test): 
  • sasmodels/modelinfo.py

    r4f4d3e3 rfd7291e  
    164164    parameter.length = length 
    165165    parameter.length_control = control 
    166  
    167166    return parameter 
    168167 
     
    265264    will have a magnitude and a direction, which may be different from 
    266265    other sld parameters. The volume parameters are used for calls 
    267     to form_volume within the kernel (required for volume normalization) 
    268     and for calls to ER and VR for effective radius and volume ratio 
    269     respectively. 
     266    to form_volume within the kernel (required for volume normalization), 
     267    to shell_volume (for hollow shapes), and to effective_radius (for 
     268    structure factor interactions) respectively. 
    270269 
    271270    *description* is a short description of the parameter.  This will 
     
    424423      parameters counted as n individual parameters p1, p2, ... 
    425424 
     425    * *common_parameters* is the list of common parameters, with a unique 
     426      copy for each model so that structure factors can have a default 
     427      background of 0.0. 
     428 
    426429    * *call_parameters* is the complete list of parameters to the kernel, 
    427430      including scale and background, with vector parameters recorded as 
     
    436439    parameters don't use vector notation, and instead use p1, p2, ... 
    437440    """ 
    438     # scale and background are implicit parameters 
    439     COMMON = [Parameter(*p) for p in COMMON_PARAMETERS] 
    440  
    441441    def __init__(self, parameters): 
    442442        # type: (List[Parameter]) -> None 
     443 
     444        # scale and background are implicit parameters 
     445        # Need them to be unique to each model in case they have different 
     446        # properties, such as default=0.0 for structure factor backgrounds. 
     447        self.common_parameters = [Parameter(*p) for p in COMMON_PARAMETERS] 
     448 
    443449        self.kernel_parameters = parameters 
    444450        self._set_vector_lengths() 
    445  
    446451        self.npars = sum(p.length for p in self.kernel_parameters) 
    447452        self.nmagnetic = sum(p.length for p in self.kernel_parameters 
     
    450455        if self.nmagnetic: 
    451456            self.nvalues += 3 + 3*self.nmagnetic 
    452  
    453457        self.call_parameters = self._get_call_parameters() 
    454458        self.defaults = self._get_defaults() 
     
    490494                         if p.polydisperse and p.type not in ('orientation', 'magnetic')) 
    491495        self.pd_2d = set(p.name for p in self.call_parameters if p.polydisperse) 
     496 
     497    def set_zero_background(self): 
     498        """ 
     499        Set the default background to zero for this model.  This is done for 
     500        structure factor models. 
     501        """ 
     502        # type: () -> None 
     503        # Make sure background is the second common parameter. 
     504        assert self.common_parameters[1].id == "background" 
     505        self.common_parameters[1].default = 0.0 
     506        self.defaults = self._get_defaults() 
    492507 
    493508    def check_angles(self): 
     
    589604    def _get_call_parameters(self): 
    590605        # type: () -> List[Parameter] 
    591         full_list = self.COMMON[:] 
     606        full_list = self.common_parameters[:] 
    592607        for p in self.kernel_parameters: 
    593608            if p.length == 1: 
     
    692707 
    693708        # Gather the user parameters in order 
    694         result = control + self.COMMON 
     709        result = control + self.common_parameters 
    695710        for p in self.kernel_parameters: 
    696711            if not is2d and p.type in ('orientation', 'magnetic'): 
     
    741756 
    742757#: Set of variables defined in the model that might contain C code 
    743 C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'c_code'] 
     758C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'shell_volume', 'c_code'] 
    744759 
    745760def _find_source_lines(model_info, kernel_module): 
     
    790805        # Custom sum/multi models 
    791806        return kernel_module.model_info 
     807 
    792808    info = ModelInfo() 
     809 
     810    # Build the parameter table 
    793811    #print("make parameter table", kernel_module.parameters) 
    794812    parameters = make_parameter_table(getattr(kernel_module, 'parameters', [])) 
     813 
     814    # background defaults to zero for structure factor models 
     815    structure_factor = getattr(kernel_module, 'structure_factor', False) 
     816    if structure_factor: 
     817        parameters.set_zero_background() 
     818 
     819    # TODO: remove demo parameters 
     820    # The plots in the docs are generated from the model default values. 
     821    # Sascomp set parameters from the command line, and so doesn't need 
     822    # demo values for testing. 
    795823    demo = expand_pars(parameters, getattr(kernel_module, 'demo', None)) 
     824 
    796825    filename = abspath(kernel_module.__file__).replace('.pyc', '.py') 
    797826    kernel_id = splitext(basename(filename))[0] 
     
    811840    info.category = getattr(kernel_module, 'category', None) 
    812841    info.structure_factor = getattr(kernel_module, 'structure_factor', False) 
     842    # TODO: find Fq by inspection 
     843    info.effective_radius_type = getattr(kernel_module, 'effective_radius_type', None) 
     844    info.have_Fq = getattr(kernel_module, 'have_Fq', False) 
    813845    info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y']) 
    814846    # Note: custom.load_custom_kernel_module assumes the C sources are defined 
     
    816848    info.source = getattr(kernel_module, 'source', []) 
    817849    info.c_code = getattr(kernel_module, 'c_code', None) 
     850    info.effective_radius = getattr(kernel_module, 'effective_radius', None) 
    818851    # TODO: check the structure of the tests 
    819852    info.tests = getattr(kernel_module, 'tests', []) 
    820     info.ER = getattr(kernel_module, 'ER', None) # type: ignore 
    821     info.VR = getattr(kernel_module, 'VR', None) # type: ignore 
    822853    info.form_volume = getattr(kernel_module, 'form_volume', None) # type: ignore 
     854    info.shell_volume = getattr(kernel_module, 'shell_volume', None) # type: ignore 
    823855    info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore 
    824856    info.Iqxy = getattr(kernel_module, 'Iqxy', None) # type: ignore 
     
    844876    info.lineno = {} 
    845877    _find_source_lines(info, kernel_module) 
    846  
    847878    return info 
    848879 
     
    925956    #: provided in the file. 
    926957    structure_factor = None # type: bool 
     958    #: True if the model defines an Fq function with signature 
     959    #: void Fq(double q, double *F1, double *F2, ...) 
     960    have_Fq = False 
     961    #: List of options for computing the effective radius of the shape, 
     962    #: or None if the model is not usable as a form factor model. 
     963    effective_radius_type = None   # type: List[str] 
    927964    #: List of C source files used to define the model.  The source files 
    928965    #: should define the *Iq* function, and possibly *Iqac* or *Iqabc* if the 
    929966    #: model defines orientation parameters. Files containing the most basic 
    930967    #: functions must appear first in the list, followed by the files that 
    931     #: use those functions.  Form factors are indicated by providing 
    932     #: an :attr:`ER` function. 
     968    #: use those functions. 
    933969    source = None           # type: List[str] 
    934     #: The set of tests that must pass.  The format of the tests is described 
    935     #: in :mod:`model_test`. 
    936     tests = None            # type: List[TestCondition] 
    937     #: Returns the effective radius of the model given its volume parameters. 
    938     #: The presence of *ER* indicates that the model is a form factor model 
    939     #: that may be used together with a structure factor to form an implicit 
    940     #: multiplication model. 
    941     #: 
    942     #: The parameters to the *ER* function must be marked with type *volume*. 
    943     #: in the parameter table.  They will appear in the same order as they 
    944     #: do in the table.  The values passed to *ER* will be vectors, with one 
    945     #: value for each polydispersity condition.  For example, if the model 
    946     #: is polydisperse over both length and radius, then both length and 
    947     #: radius will have the same number of values in the vector, with one 
    948     #: value for each *length X radius*.  If only *radius* is polydisperse, 
    949     #: then the value for *length* will be repeated once for each value of 
    950     #: *radius*.  The *ER* function should return one effective radius for 
    951     #: each parameter set.  Multiplicity parameters will be received as 
    952     #: arrays, with one row per polydispersity condition. 
    953     ER = None               # type: Optional[Callable[[np.ndarray], np.ndarray]] 
    954     #: Returns the occupied volume and the total volume for each parameter set. 
    955     #: See :attr:`ER` for details on the parameters. 
    956     VR = None               # type: Optional[Callable[[np.ndarray], Tuple[np.ndarray, np.ndarray]]] 
    957     #: Arbitrary C code containing supporting functions, etc., to be inserted 
    958     #: after everything in source.  This can include Iq and Iqxy functions with 
    959     #: the full function signature, including all parameters. 
    960     c_code = None 
     970    #: inline source code, added after all elements of source 
     971    c_code = None           # type: Optional[str] 
    961972    #: Returns the form volume for python-based models.  Form volume is needed 
    962973    #: for volume normalization in the polydispersity integral.  If no 
     
    966977    #: C code, either defined as a string, or in the sources. 
    967978    form_volume = None      # type: Union[None, str, Callable[[np.ndarray], float]] 
     979    #: Returns the shell volume for python-based models.  Form volume and 
     980    #: shell volume are needed for volume normalization in the polydispersity 
     981    #: integral and structure interactions for hollow shapes.  If no 
     982    #: parameters are *volume* parameters, then shell volume is not needed. 
     983    #: For C-based models, (with :attr:`sources` defined, or with :attr:`Iq` 
     984    #: defined using a string containing C code), shell_volume must also be 
     985    #: C code, either defined as a string, or in the sources. 
     986    shell_volume = None      # type: Union[None, str, Callable[[np.ndarray], float]] 
    968987    #: Returns *I(q, a, b, ...)* for parameters *a*, *b*, etc. defined 
    969988    #: by the parameter table.  *Iq* can be defined as a python function, or 
     
    10001019    #: Line numbers for symbols defining C code 
    10011020    lineno = None           # type: Dict[str, int] 
     1021    #: The set of tests that must pass.  The format of the tests is described 
     1022    #: in :mod:`model_test`. 
     1023    tests = None            # type: List[TestCondition] 
    10021024 
    10031025    def __init__(self): 
  • sasmodels/models/_spherepy.py

    rca4444f r304c775  
    7171    return 1.333333333333333 * pi * radius ** 3 
    7272 
     73def effective_radius(mode, radius): 
     74    return radius 
     75 
    7376def Iq(q, sld, sld_solvent, radius): 
    7477    #print "q",q 
     
    105108sesans.vectorized = True  # sesans accepts an array of z values 
    106109 
    107 def ER(radius): 
    108     return radius 
    109  
    110 # VR defaults to 1.0 
    111  
    112110demo = dict(scale=1, background=0, 
    113111            sld=6, sld_solvent=1, 
  • sasmodels/models/barbell.c

    r108e70e r99658f6  
    6363 
    6464static double 
    65 Iq(double q, double sld, double solvent_sld, 
     65radius_from_excluded_volume(double radius_bell, double radius, double length) 
     66{ 
     67    const double hdist = sqrt(square(radius_bell) - square(radius)); 
     68    const double length_tot = length + 2.0*(hdist+ radius); 
     69    return 0.5*cbrt(0.75*radius_bell*(2.0*radius_bell*length_tot + (radius_bell + length_tot)*(M_PI*radius_bell + length_tot))); 
     70} 
     71 
     72static double 
     73radius_from_volume(double radius_bell, double radius, double length) 
     74{ 
     75    const double vol_barbell = form_volume(radius_bell,radius,length); 
     76    return cbrt(vol_barbell/M_4PI_3); 
     77} 
     78 
     79static double 
     80radius_from_totallength(double radius_bell, double radius, double length) 
     81{ 
     82    const double hdist = sqrt(square(radius_bell) - square(radius)); 
     83    return 0.5*length + hdist + radius_bell; 
     84} 
     85 
     86static double 
     87effective_radius(int mode, double radius_bell, double radius, double length) 
     88{ 
     89    switch (mode) { 
     90    default: 
     91    case 1: // equivalent cylinder excluded volume 
     92        return radius_from_excluded_volume(radius_bell, radius , length); 
     93    case 2: // equivalent volume sphere 
     94        return radius_from_volume(radius_bell, radius , length); 
     95    case 3: // radius 
     96        return radius; 
     97    case 4: // half length 
     98        return 0.5*length; 
     99    case 5: // half total length 
     100        return radius_from_totallength(radius_bell,radius,length); 
     101    } 
     102} 
     103 
     104static void 
     105Fq(double q,double *F1, double *F2, double sld, double solvent_sld, 
    66106    double radius_bell, double radius, double length) 
    67107{ 
     
    72112    const double zm = M_PI_4; 
    73113    const double zb = M_PI_4; 
    74     double total = 0.0; 
     114    double total_F1 = 0.0; 
     115    double total_F2 = 0.0; 
    75116    for (int i = 0; i < GAUSS_N; i++){ 
    76117        const double alpha= GAUSS_Z[i]*zm + zb; 
     
    78119        SINCOS(alpha, sin_alpha, cos_alpha); 
    79120        const double Aq = _fq(q*sin_alpha, q*cos_alpha, h, radius_bell, radius, half_length); 
    80         total += GAUSS_W[i] * Aq * Aq * sin_alpha; 
     121        total_F1 += GAUSS_W[i] * Aq * sin_alpha; 
     122        total_F2 += GAUSS_W[i] * Aq * Aq * sin_alpha; 
    81123    } 
    82124    // translate dx in [-1,1] to dx in [lower,upper] 
    83     const double form = total*zm; 
     125    const double form_avg = total_F1*zm; 
     126    const double form_squared_avg = total_F2*zm; 
    84127 
    85128    //Contrast 
    86129    const double s = (sld - solvent_sld); 
    87     return 1.0e-4 * s * s * form; 
     130    *F1 = 1.0e-2 * s * form_avg; 
     131    *F2 = 1.0e-4 * s * s * form_squared_avg; 
    88132} 
    89  
    90133 
    91134static double 
  • sasmodels/models/barbell.py

    r2d81cfe r99658f6  
    7979.. [#] H Kaya and N R deSouza, *J. Appl. Cryst.*, 37 (2004) 508-509 (addenda 
    8080   and errata) 
     81L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
    8182 
    8283Authorship and Verification 
     
    116117 
    117118source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "barbell.c"] 
     119have_Fq = True  
     120effective_radius_type = [ 
     121    "equivalent cylinder excluded volume","equivalent volume sphere", "radius", "half length", "half total length", 
     122    ] 
    118123 
    119124def random(): 
  • sasmodels/models/binary_hard_sphere.c

    r925ad6e r71b751d  
    11double form_volume(void); 
    22 
    3 double Iq(double q,  
     3double Iq(double q, 
    44    double lg_radius, double sm_radius, 
    55    double lg_vol_frac, double sm_vol_frac, 
    66    double lg_sld, double sm_sld, double solvent_sld 
    77    ); 
    8      
     8 
    99void calculate_psfs(double qval, 
    1010    double r2, double nf2, 
     
    2222    double lg_vol_frac, double sm_vol_frac, 
    2323    double lg_sld, double sm_sld, double solvent_sld) 
    24      
    2524{ 
    2625    double r2,r1,nf2,phi,aa,rho2,rho1,rhos,inten;       //my local names 
     
    2827    double phi1,phi2,phr,a3; 
    2928    double v1,v2,n1,n2,qr1,qr2,b1,b2,sc1,sc2; 
    30      
     29 
    3130    r2 = lg_radius; 
    3231    r1 = sm_radius; 
     
    3635    rho1 = sm_sld; 
    3736    rhos = solvent_sld; 
    38      
    39      
     37 
     38 
    4039    phi = phi1 + phi2; 
    4140    aa = r1/r2; 
     
    4645    // calculate the PSF's here 
    4746    calculate_psfs(q,r2,nf2,aa,phi,&psf11,&psf22,&psf12); 
    48      
     47 
    4948    // /* do form factor calculations  */ 
    50      
     49 
    5150    v1 = M_4PI_3*r1*r1*r1; 
    5251    v2 = M_4PI_3*r2*r2*r2; 
    53      
     52 
    5453    n1 = phi1/v1; 
    5554    n2 = phi2/v2; 
    56      
     55 
    5756    qr1 = r1*q; 
    5857    qr2 = r2*q; 
     
    6867    inten *= 1.0e8; 
    6968    ///*convert rho^2 in 10^-6A to A*/ 
    70     inten *= 1.0e-12;     
     69    inten *= 1.0e-12; 
    7170    return(inten); 
    7271} 
     
    7776    double aa, double phi, 
    7877    double *s11, double *s22, double *s12) 
    79  
    8078{ 
    8179    //  variable qval,r2,nf2,aa,phi,&s11,&s22,&s12 
    82      
     80 
    8381    //   calculate constant terms 
    8482    double s2,v,a3,v1,v2,g11,g12,g22,wmv,wmv3,wmv4; 
     
    8785    double c11,c22,c12,f12,f22,ttt1,ttt2,ttt3,ttt4,yl,y13; 
    8886    double t21,t22,t23,t31,t32,t33,t41,t42,yl3,wma3,y1; 
    89      
     87 
    9088    s2 = 2.0*r2; 
    9189//    s1 = aa*s2;  why is this never used?  check original paper? 
     
    108106    gm1=(v1*a1+a3*v2*a2)*.5; 
    109107    gm12=2.*gm1*(1.-aa)/aa; 
    110     //c   
     108    //c 
    111109    //c   calculate the direct correlation functions and print results 
    112110    //c 
    113111    //  do 20 j=1,npts 
    114      
     112 
    115113    yy=qval*s2; 
    116114    //c   calculate direct correlation functions 
     
    123121    t3=gm1*((4.*ay*ay2-24.*ay)*sin(ay)-(ay2*ay2-12.*ay2+24.)*cos(ay)+24.)/ay3; 
    124122    f11=24.*v1*(t1+t2+t3)/ay3; 
    125      
     123 
    126124    //c ------c22 
    127125    y2=yy*yy; 
     
    131129    tt3=gm1*((4.*y3-24.*yy)*sin(yy)-(y2*y2-12.*y2+24.)*cos(yy)+24.)/ay3; 
    132130    f22=24.*v2*(tt1+tt2+tt3)/y3; 
    133      
     131 
    134132    //c   -----c12 
    135133    yl=.5*yy*(1.-aa); 
     
    151149    ttt4=a1*(t41+t42)/y1; 
    152150    f12=ttt1+24.*v*sqrt(nf2)*sqrt(1.-nf2)*a3*(ttt2+ttt3+ttt4)/(nf2+(1.-nf2)*a3); 
    153      
     151 
    154152    c11=f11; 
    155153    c22=f22; 
    156154    c12=f12; 
    157155    *s11=1./(1.+c11-(c12)*c12/(1.+c22)); 
    158     *s22=1./(1.+c22-(c12)*c12/(1.+c11));  
    159     *s12=-c12/((1.+c11)*(1.+c22)-(c12)*(c12));    
    160      
     156    *s22=1./(1.+c22-(c12)*c12/(1.+c11)); 
     157    *s12=-c12/((1.+c11)*(1.+c22)-(c12)*(c12)); 
     158 
    161159    return; 
    162160} 
    163  
  • sasmodels/models/capped_cylinder.c

    r108e70e r99658f6  
    8585 
    8686static double 
    87 Iq(double q, double sld, double solvent_sld, 
     87radius_from_excluded_volume(double radius, double radius_cap, double length) 
     88{ 
     89    const double hc = radius_cap - sqrt(radius_cap*radius_cap - radius*radius); 
     90    const double length_tot = length + 2.0*hc; 
     91    return 0.5*cbrt(0.75*radius*(2.0*radius*length_tot + (radius + length_tot)*(M_PI*radius + length_tot))); 
     92} 
     93 
     94static double 
     95radius_from_volume(double radius, double radius_cap, double length) 
     96{ 
     97    const double vol_cappedcyl = form_volume(radius,radius_cap,length); 
     98    return cbrt(vol_cappedcyl/M_4PI_3); 
     99} 
     100 
     101static double 
     102radius_from_totallength(double radius, double radius_cap, double length) 
     103{ 
     104    const double hc = radius_cap - sqrt(radius_cap*radius_cap - radius*radius); 
     105    return 0.5*length + hc; 
     106} 
     107 
     108static double 
     109effective_radius(int mode, double radius, double radius_cap, double length) 
     110{ 
     111    switch (mode) { 
     112    default: 
     113    case 1: // equivalent cylinder excluded volume 
     114        return radius_from_excluded_volume(radius, radius_cap, length); 
     115    case 2: // equivalent volume sphere 
     116        return radius_from_volume(radius, radius_cap, length); 
     117    case 3: // radius 
     118        return radius; 
     119    case 4: // half length 
     120        return 0.5*length; 
     121    case 5: // half total length 
     122        return radius_from_totallength(radius, radius_cap,length); 
     123    } 
     124} 
     125 
     126static void 
     127Fq(double q,double *F1, double *F2, double sld, double solvent_sld, 
    88128    double radius, double radius_cap, double length) 
    89129{ 
     
    94134    const double zm = M_PI_4; 
    95135    const double zb = M_PI_4; 
    96     double total = 0.0; 
     136    double total_F1 = 0.0; 
     137    double total_F2 = 0.0; 
    97138    for (int i=0; i<GAUSS_N ;i++) { 
    98139        const double theta = GAUSS_Z[i]*zm + zb; 
     
    103144        const double Aq = _fq(qab, qc, h, radius_cap, radius, half_length); 
    104145        // scale by sin_theta for spherical coord integration 
    105         total += GAUSS_W[i] * Aq * Aq * sin_theta; 
     146        total_F1 += GAUSS_W[i] * Aq * sin_theta; 
     147        total_F2 += GAUSS_W[i] * Aq * Aq * sin_theta; 
    106148    } 
    107149    // translate dx in [-1,1] to dx in [lower,upper] 
    108     const double form = total * zm; 
     150    const double form_avg = total_F1 * zm; 
     151    const double form_squared_avg = total_F2 * zm; 
    109152 
    110153    // Contrast 
    111154    const double s = (sld - solvent_sld); 
    112     return 1.0e-4 * s * s * form; 
     155    *F1 = 1.0e-2 * s * form_avg; 
     156    *F2 = 1.0e-4 * s * s * form_squared_avg; 
    113157} 
    114158 
  • sasmodels/models/capped_cylinder.py

    r2d81cfe r99658f6  
    8282.. [#] H Kaya and N-R deSouza, *J. Appl. Cryst.*, 37 (2004) 508-509 (addenda 
    8383   and errata) 
     84L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
    8485 
    8586Authorship and Verification 
     
    136137 
    137138source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "capped_cylinder.c"] 
     139have_Fq = True 
     140effective_radius_type = [ 
     141    "equivalent cylinder excluded volume", "equivalent volume sphere", "radius", "half length", "half total length", 
     142    ] 
    138143 
    139144def random(): 
  • sasmodels/models/core_multi_shell.c

    rc3ccaec rd42dd4a  
    99 
    1010static double 
    11 form_volume(double core_radius, double fp_n, double thickness[]) 
     11outer_radius(double core_radius, double fp_n, double thickness[]) 
    1212{ 
    1313  double r = core_radius; 
     
    1616    r += thickness[i]; 
    1717  } 
    18   return M_4PI_3 * cube(r); 
     18  return r; 
    1919} 
    2020 
    2121static double 
    22 Iq(double q, double core_sld, double core_radius, 
     22form_volume(double core_radius, double fp_n, double thickness[]) 
     23{ 
     24  return M_4PI_3 * cube(outer_radius(core_radius, fp_n, thickness)); 
     25} 
     26 
     27static double 
     28effective_radius(int mode, double core_radius, double fp_n, double thickness[]) 
     29{ 
     30  switch (mode) { 
     31  default: 
     32  case 1: // outer radius 
     33    return outer_radius(core_radius, fp_n, thickness); 
     34  case 2: // core radius 
     35    return core_radius; 
     36  } 
     37} 
     38 
     39static void 
     40Fq(double q, double *F1, double *F2, double core_sld, double core_radius, 
    2341   double solvent_sld, double fp_n, double sld[], double thickness[]) 
    2442{ 
     
    3452  } 
    3553  f += M_4PI_3 * cube(r) * (solvent_sld - last_sld) * sas_3j1x_x(q*r); 
    36   return f * f * 1.0e-4; 
     54  *F1 = 1e-2 * f; 
     55  *F2 = 1e-4 * f * f; 
    3756} 
  • sasmodels/models/core_multi_shell.py

    r2d81cfe ree60aa7  
    9999 
    100100source = ["lib/sas_3j1x_x.c", "core_multi_shell.c"] 
     101have_Fq = True 
     102effective_radius_type = ["outer radius", "core radius"] 
    101103 
    102104def random(): 
     
    143145    return np.asarray(z), np.asarray(rho) 
    144146 
    145 def ER(radius, n, thickness): 
    146     """Effective radius""" 
    147     n = int(n[0]+0.5)  # n is a control parameter and is not polydisperse 
    148     return np.sum(thickness[:n], axis=0) + radius 
    149  
    150147demo = dict(sld_core=6.4, 
    151148            radius=60, 
  • sasmodels/models/core_shell_bicelle.c

    r108e70e r99658f6  
    3737 
    3838static double 
    39 Iq(double q, 
     39radius_from_excluded_volume(double radius, double thick_rim, double thick_face, double length) 
     40{ 
     41    const double radius_tot = radius + thick_rim; 
     42    const double length_tot = length + 2.0*thick_face; 
     43    return 0.5*cbrt(0.75*radius_tot*(2.0*radius_tot*length_tot + (radius_tot + length_tot)*(M_PI*radius_tot + length_tot))); 
     44} 
     45 
     46static double 
     47radius_from_volume(double radius, double thick_rim, double thick_face, double length) 
     48{ 
     49    const double volume_bicelle = form_volume(radius,thick_rim,thick_face,length); 
     50    return cbrt(volume_bicelle/M_4PI_3); 
     51} 
     52 
     53static double 
     54radius_from_diagonal(double radius, double thick_rim, double thick_face, double length) 
     55{ 
     56    const double radius_tot = radius + thick_rim; 
     57    const double length_tot = length + 2.0*thick_face; 
     58    return sqrt(radius_tot*radius_tot + 0.25*length_tot*length_tot); 
     59} 
     60 
     61static double 
     62effective_radius(int mode, double radius, double thick_rim, double thick_face, double length) 
     63{ 
     64    switch (mode) { 
     65    default: 
     66    case 1: // equivalent cylinder excluded volume 
     67        return radius_from_excluded_volume(radius, thick_rim, thick_face, length); 
     68    case 2: // equivalent sphere 
     69        return radius_from_volume(radius, thick_rim, thick_face, length); 
     70    case 3: // outer rim radius 
     71        return radius + thick_rim; 
     72    case 4: // half outer thickness 
     73        return 0.5*length + thick_face; 
     74    case 5: // half diagonal 
     75        return radius_from_diagonal(radius,thick_rim,thick_face,length); 
     76    } 
     77} 
     78 
     79static void 
     80Fq(double q, 
     81    double *F1, 
     82    double *F2, 
    4083    double radius, 
    4184    double thick_radius, 
     
    5194    const double halflength = 0.5*length; 
    5295 
    53     double total = 0.0; 
     96    double total_F1 = 0.0; 
     97    double total_F2 = 0.0; 
    5498    for(int i=0;i<GAUSS_N;i++) { 
    5599        double theta = (GAUSS_Z[i] + 1.0)*uplim; 
    56100        double sin_theta, cos_theta; // slots to hold sincos function output 
    57101        SINCOS(theta, sin_theta, cos_theta); 
    58         double fq = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face, 
     102        double form = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face, 
    59103                                   halflength, sld_core, sld_face, sld_rim, sld_solvent); 
    60         total += GAUSS_W[i]*fq*fq*sin_theta; 
     104        total_F1 += GAUSS_W[i]*form*sin_theta; 
     105        total_F2 += GAUSS_W[i]*form*form*sin_theta; 
    61106    } 
     107    // Correct for integration range 
     108    total_F1 *= uplim; 
     109    total_F2 *= uplim; 
    62110 
    63     // calculate value of integral to return 
    64     double answer = total*uplim; 
    65     return 1.0e-4*answer; 
     111    *F1 = 1.0e-2*total_F1; 
     112    *F2 = 1.0e-4*total_F2; 
    66113} 
    67114 
  • sasmodels/models/core_shell_bicelle.py

    r2d81cfe r99658f6  
    8989   from Proquest <http://search.proquest.com/docview/304915826?accountid 
    9090   =26379>`_ 
     91    
     92   L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
    9193 
    9294Authorship and Verification 
     
    154156source = ["lib/sas_Si.c", "lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", 
    155157          "core_shell_bicelle.c"] 
     158have_Fq = True 
     159effective_radius_type = [ 
     160    "excluded volume","equivalent volume sphere", "outer rim radius", 
     161    "half outer thickness", "half diagonal", 
     162    ] 
    156163 
    157164def random(): 
  • sasmodels/models/core_shell_bicelle_elliptical.c

    r108e70e r99658f6  
    1111 
    1212static double 
    13 Iq(double q, 
     13radius_from_excluded_volume(double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     14{ 
     15    const double r_equiv     = sqrt((r_minor + thick_rim)*(r_minor*x_core + thick_rim)); 
     16    const double length_tot  = length + 2.0*thick_face; 
     17    return 0.5*cbrt(0.75*r_equiv*(2.0*r_equiv*length_tot + (r_equiv + length_tot)*(M_PI*r_equiv + length_tot))); 
     18} 
     19 
     20static double 
     21radius_from_volume(double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     22{ 
     23    const double volume_bicelle = form_volume(r_minor, x_core, thick_rim,thick_face,length); 
     24    return cbrt(volume_bicelle/M_4PI_3); 
     25} 
     26 
     27static double 
     28radius_from_diagonal(double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     29{ 
     30    const double radius_max = (x_core < 1.0 ? r_minor : x_core*r_minor); 
     31    const double radius_max_tot = radius_max + thick_rim; 
     32    const double length_tot = length + 2.0*thick_face; 
     33    return sqrt(radius_max_tot*radius_max_tot + 0.25*length_tot*length_tot); 
     34} 
     35 
     36static double 
     37effective_radius(int mode, double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     38{ 
     39    switch (mode) { 
     40    default: 
     41    case 1: // equivalent cylinder excluded volume 
     42        return radius_from_excluded_volume(r_minor, x_core, thick_rim, thick_face, length); 
     43    case 2: // equivalent volume sphere 
     44        return radius_from_volume(r_minor, x_core, thick_rim, thick_face, length); 
     45    case 3: // outer rim average radius 
     46        return 0.5*r_minor*(1.0 + x_core) + thick_rim; 
     47    case 4: // outer rim min radius 
     48        return (x_core < 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
     49    case 5: // outer max radius 
     50        return (x_core > 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
     51    case 6: // half outer thickness 
     52        return 0.5*length + thick_face; 
     53    case 7: // half diagonal 
     54        return radius_from_diagonal(r_minor,x_core,thick_rim,thick_face,length); 
     55    } 
     56} 
     57 
     58static void 
     59Fq(double q, 
     60    double *F1, 
     61    double *F2, 
    1462    double r_minor, 
    1563    double x_core, 
     
    3684 
    3785    //initialize integral 
    38     double outer_sum = 0.0; 
     86    double outer_total_F1 = 0.0; 
     87    double outer_total_F2 = 0.0; 
    3988    for(int i=0;i<GAUSS_N;i++) { 
    4089        //setup inner integral over the ellipsoidal cross-section 
    41         //const double va = 0.0; 
    42         //const double vb = 1.0; 
    4390        //const double cos_theta = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 
    4491        const double cos_theta = ( GAUSS_Z[i] + 1.0 )/2.0; 
     
    4895        const double si1 = sas_sinx_x(halfheight*qc); 
    4996        const double si2 = sas_sinx_x((halfheight+thick_face)*qc); 
    50         double inner_sum=0.0; 
     97        double inner_total_F1 = 0; 
     98        double inner_total_F2 = 0; 
    5199        for(int j=0;j<GAUSS_N;j++) { 
    52100            //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe) 
    53             // inner integral limits 
    54             //const double vaj=0.0; 
    55             //const double vbj=M_PI; 
    56             //const double phi = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    57             const double phi = ( GAUSS_Z[j] +1.0)*M_PI_2; 
    58             const double rr = sqrt(r2A - r2B*cos(phi)); 
     101            //const double beta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     102            const double beta = ( GAUSS_Z[j] +1.0)*M_PI_2; 
     103            const double rr = sqrt(r2A - r2B*cos(beta)); 
    59104            const double be1 = sas_2J1x_x(rr*qab); 
    60105            const double be2 = sas_2J1x_x((rr+thick_rim)*qab); 
    61             const double fq = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1; 
     106            const double f = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1; 
    62107 
    63             inner_sum += GAUSS_W[j] * fq * fq; 
     108            inner_total_F1 += GAUSS_W[j] * f; 
     109            inner_total_F2 += GAUSS_W[j] * f * f; 
    64110        } 
    65111        //now calculate outer integral 
    66         outer_sum += GAUSS_W[i] * inner_sum; 
     112        outer_total_F1 += GAUSS_W[i] * inner_total_F1; 
     113        outer_total_F2 += GAUSS_W[i] * inner_total_F2; 
    67114    } 
     115    // now complete change of integration variables (1-0)/(1-(-1))= 0.5 
     116    outer_total_F1 *= 0.25; 
     117    outer_total_F2 *= 0.25; 
    68118 
    69     return outer_sum*2.5e-05; 
     119    //convert from [1e-12 A-1] to [cm-1] 
     120    *F1 = 1e-2*outer_total_F1; 
     121    *F2 = 1e-4*outer_total_F2; 
    70122} 
    71123 
  • sasmodels/models/core_shell_bicelle_elliptical.py

    rfc3ae1b r99658f6  
    100100 
    101101.. [#] 
     102L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
    102103 
    103104Authorship and Verification 
     
    146147source = ["lib/sas_Si.c", "lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", 
    147148          "core_shell_bicelle_elliptical.c"] 
     149have_Fq = True 
     150effective_radius_type = [ 
     151    "equivalent cylinder excluded volume", "equivalent volume sphere", "outer rim average radius", "outer rim min radius", 
     152    "outer max radius", "half outer thickness", "half diagonal", 
     153    ] 
    148154 
    149155def random(): 
     
    179185 
    180186tests = [ 
    181     [{'radius': 30.0, 'x_core': 3.0, 
    182       'thick_rim': 8.0, 'thick_face': 14.0, 'length': 50.0}, 'ER', 1], 
    183     [{'radius': 30.0, 'x_core': 3.0, 
    184       'thick_rim': 8.0, 'thick_face': 14.0, 'length': 50.0}, 'VR', 1], 
     187    #[{'radius': 30.0, 'x_core': 3.0, 
     188    #  'thick_rim': 8.0, 'thick_face': 14.0, 'length': 50.0}, 'ER', 1], 
     189    #[{'radius': 30.0, 'x_core': 3.0, 
     190    #  'thick_rim': 8.0, 'thick_face': 14.0, 'length': 50.0}, 'VR', 1], 
    185191 
    186192    [{'radius': 30.0, 'x_core': 3.0, 
  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

    r108e70e r99658f6  
    1212 
    1313static double 
    14 Iq(double q, 
     14radius_from_excluded_volume(double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     15{ 
     16    const double r_equiv     = sqrt((r_minor + thick_rim)*(r_minor*x_core + thick_rim)); 
     17    const double length_tot  = length + 2.0*thick_face; 
     18    return 0.5*cbrt(0.75*r_equiv*(2.0*r_equiv*length_tot + (r_equiv + length_tot)*(M_PI*r_equiv + length_tot))); 
     19} 
     20 
     21static double 
     22radius_from_volume(double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     23{ 
     24    const double volume_bicelle = form_volume(r_minor, x_core, thick_rim,thick_face,length); 
     25    return cbrt(volume_bicelle/M_4PI_3); 
     26} 
     27 
     28static double 
     29radius_from_diagonal(double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     30{ 
     31    const double radius_max = (x_core < 1.0 ? r_minor : x_core*r_minor); 
     32    const double radius_max_tot = radius_max + thick_rim; 
     33    const double length_tot = length + 2.0*thick_face; 
     34    return sqrt(radius_max_tot*radius_max_tot + 0.25*length_tot*length_tot); 
     35} 
     36 
     37static double 
     38effective_radius(int mode, double r_minor, double x_core, double thick_rim, double thick_face, double length) 
     39{ 
     40    switch (mode) { 
     41    default: 
     42    case 1: // equivalent cylinder excluded volume 
     43        return radius_from_excluded_volume(r_minor, x_core, thick_rim, thick_face, length); 
     44    case 2: // equivalent sphere 
     45        return radius_from_volume(r_minor, x_core, thick_rim, thick_face, length); 
     46    case 3: // outer rim average radius 
     47        return 0.5*r_minor*(1.0 + x_core) + thick_rim; 
     48    case 4: // outer rim min radius 
     49        return (x_core < 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
     50    case 5: // outer max radius 
     51        return (x_core > 1.0 ? x_core*r_minor+thick_rim : r_minor+thick_rim); 
     52    case 6: // half outer thickness 
     53        return 0.5*length + thick_face; 
     54    case 7: // half diagonal (this ignores the missing "corners", so may give unexpected answer if thick_face 
     55            //  or thick_rim is extremely large) 
     56        return radius_from_diagonal(r_minor,x_core,thick_rim,thick_face,length); 
     57    } 
     58} 
     59 
     60static void 
     61Fq(double q, 
     62        double *F1, 
     63        double *F2, 
    1564        double r_minor, 
    1665        double x_core, 
     
    2473        double sigma) 
    2574{ 
    26     double si1,si2,be1,be2; 
    2775     // core_shell_bicelle_elliptical_belt, RKH 5th Oct 2017, core_shell_bicelle_elliptical 
    2876     // tested briefly against limiting cases of cylinder, hollow cylinder & elliptical cylinder models 
     
    3886    const double r2A = 0.5*(square(r_major) + square(r_minor)); 
    3987    const double r2B = 0.5*(square(r_major) - square(r_minor)); 
     88    const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight); 
     89    const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*halfheight; 
     90    const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face); 
    4091    // dr1,2,3 are now for Vcore, Vcore+rim, Vcore+face, 
    41     const double dr1 = (-rhor - rhoh + rhoc + rhosolv) *M_PI*r_minor*r_major* 
    42           2.0*halfheight; 
    43     const double dr2 = (rhor-rhosolv) *M_PI*(r_minor+thick_rim)*( 
    44          r_major+thick_rim)* 2.0*halfheight; 
    45     const double dr3 = (rhoh-rhosolv) *M_PI*r_minor*r_major* 
    46          2.0*(halfheight+thick_face); 
     92    const double dr1 = vol1*(-rhor - rhoh + rhoc + rhosolv); 
     93    const double dr2 = vol2*(rhor-rhosolv); 
     94    const double dr3 = vol3*(rhoh-rhosolv); 
     95 
    4796    //initialize integral 
    48     double outer_sum = 0.0; 
     97    double outer_total_F1 = 0.0; 
     98    double outer_total_F2 = 0.0; 
    4999    for(int i=0;i<GAUSS_N;i++) { 
    50100        //setup inner integral over the ellipsoidal cross-section 
    51101        // since we generate these lots of times, why not store them somewhere? 
    52         //const double cos_alpha = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 
    53         const double cos_alpha = ( GAUSS_Z[i] + 1.0 )/2.0; 
    54         const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 
    55         double inner_sum=0; 
    56         double sinarg1 = q*halfheight*cos_alpha; 
    57         double sinarg2 = q*(halfheight+thick_face)*cos_alpha; 
    58         si1 = sas_sinx_x(sinarg1); 
    59         si2 = sas_sinx_x(sinarg2); 
     102        //const double cos_theta = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 
     103        const double cos_theta = ( GAUSS_Z[i] + 1.0 )/2.0; 
     104        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
     105        const double qab = q*sin_theta; 
     106        const double qc = q*cos_theta; 
     107        const double si1 = sas_sinx_x(halfheight*qc); 
     108        const double si2 = sas_sinx_x((halfheight+thick_face)*qc); 
     109        double inner_total_F1 = 0; 
     110        double inner_total_F2 = 0; 
    60111        for(int j=0;j<GAUSS_N;j++) { 
    61112            //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe) 
     
    63114            const double beta = ( GAUSS_Z[j] +1.0)*M_PI_2; 
    64115            const double rr = sqrt(r2A - r2B*cos(beta)); 
    65             double besarg1 = q*rr*sin_alpha; 
    66             double besarg2 = q*(rr+thick_rim)*sin_alpha; 
    67             be1 = sas_2J1x_x(besarg1); 
    68             be2 = sas_2J1x_x(besarg2); 
    69             inner_sum += GAUSS_W[j] *square(dr1*si1*be1 + 
    70                                               dr2*si1*be2 + 
    71                                               dr3*si2*be1); 
     116            const double be1 = sas_2J1x_x(rr*qab); 
     117            const double be2 = sas_2J1x_x((rr+thick_rim)*qab); 
     118            const double f = dr1*si1*be1 + dr2*si1*be2 + dr3*si2*be1; 
     119 
     120            inner_total_F1 += GAUSS_W[j] * f; 
     121            inner_total_F2 += GAUSS_W[j] * f * f; 
    72122        } 
    73123        //now calculate outer integral 
    74         outer_sum += GAUSS_W[i] * inner_sum; 
     124        outer_total_F1 += GAUSS_W[i] * inner_total_F1; 
     125        outer_total_F2 += GAUSS_W[i] * inner_total_F2; 
    75126    } 
     127    // now complete change of integration variables (1-0)/(1-(-1))= 0.5 
     128    outer_total_F1 *= 0.25; 
     129    outer_total_F2 *= 0.25; 
    76130 
    77     return outer_sum*2.5e-05*exp(-0.5*square(q*sigma)); 
     131    //convert from [1e-12 A-1] to [cm-1] 
     132    *F1 = 1e-2*outer_total_F1*exp(-0.25*square(q*sigma)); 
     133    *F2 = 1e-4*outer_total_F2*exp(-0.5*square(q*sigma)); 
    78134} 
    79135 
     
    111167    const double si1 = sas_sinx_x( halfheight*qc ); 
    112168    const double si2 = sas_sinx_x( (halfheight + thick_face)*qc ); 
    113     const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si1*be2 +  vol3*dr3*si2*be1); 
    114     return 1.0e-4 * Aq*exp(-0.5*(square(qa) + square(qb) + square(qc) )*square(sigma)); 
     169    const double fq = vol1*dr1*si1*be1 + vol2*dr2*si1*be2 +  vol3*dr3*si2*be1; 
     170    const double atten = exp(-0.5*(qa*qa + qb*qb + qc*qc)*(sigma*sigma)); 
     171    return 1.0e-4 * fq*fq * atten; 
    115172} 
  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py

    rfc3ae1b r99658f6  
    112112 
    113113.. [#] 
     114L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
    114115 
    115116Authorship and Verification 
     
    159160source = ["lib/sas_Si.c", "lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", 
    160161          "core_shell_bicelle_elliptical_belt_rough.c"] 
     162have_Fq = True 
     163effective_radius_type = [ 
     164    "equivalent cylinder excluded volume", "equivalent volume sphere", "outer rim average radius", "outer rim min radius", 
     165    "outer max radius", "half outer thickness", "half diagonal", 
     166    ] 
    161167 
    162168demo = dict(scale=1, background=0, 
     
    181187 
    182188tests = [ 
    183     [{'radius': 30.0, 'x_core': 3.0, 'thick_rim':8.0, 'thick_face':14.0, 'length':50.0}, 'ER', 1], 
    184     [{'radius': 30.0, 'x_core': 3.0, 'thick_rim':8.0, 'thick_face':14.0, 'length':50.0}, 'VR', 1], 
     189    #[{'radius': 30.0, 'x_core': 3.0, 'thick_rim':8.0, 'thick_face':14.0, 'length':50.0}, 'ER', 1], 
     190    #[{'radius': 30.0, 'x_core': 3.0, 'thick_rim':8.0, 'thick_face':14.0, 'length':50.0}, 'VR', 1], 
    185191 
    186192    [{'radius': 30.0, 'x_core': 3.0, 'thick_rim':8.0, 'thick_face':14.0, 'length':50.0, 
  • sasmodels/models/core_shell_cylinder.c

    rd86f0fc r99658f6  
    1414 
    1515static double 
    16 Iq(double q, 
     16radius_from_excluded_volume(double radius, double thickness, double length) 
     17{ 
     18    const double radius_tot = radius + thickness; 
     19    const double length_tot = length + 2.0*thickness; 
     20    return 0.5*cbrt(0.75*radius_tot*(2.0*radius_tot*length_tot + (radius_tot + length_tot)*(M_PI*radius_tot + length_tot))); 
     21} 
     22 
     23static double 
     24radius_from_volume(double radius, double thickness, double length) 
     25{ 
     26    const double volume_outer_cyl = form_volume(radius,thickness,length); 
     27    return cbrt(volume_outer_cyl/M_4PI_3); 
     28} 
     29 
     30static double 
     31radius_from_diagonal(double radius, double thickness, double length) 
     32{ 
     33    const double radius_outer = radius + thickness; 
     34    const double length_outer = length + 2.0*thickness; 
     35    return sqrt(radius_outer*radius_outer + 0.25*length_outer*length_outer); 
     36} 
     37 
     38static double 
     39effective_radius(int mode, double radius, double thickness, double length) 
     40{ 
     41    switch (mode) { 
     42    default: 
     43    case 1: //cylinder excluded volume 
     44        return radius_from_excluded_volume(radius, thickness, length); 
     45    case 2: // equivalent volume sphere 
     46        return radius_from_volume(radius, thickness, length); 
     47    case 3: // outer radius 
     48        return radius + thickness; 
     49    case 4: // half outer length 
     50        return 0.5*length + thickness; 
     51    case 5: // half min outer length 
     52        return (radius < 0.5*length ? radius + thickness : 0.5*length + thickness); 
     53    case 6: // half max outer length 
     54        return (radius > 0.5*length ? radius + thickness : 0.5*length + thickness); 
     55    case 7: // half outer diagonal 
     56        return radius_from_diagonal(radius,thickness,length); 
     57    } 
     58} 
     59 
     60static void 
     61Fq(double q, 
     62    double *F1, 
     63    double *F2, 
    1764    double core_sld, 
    1865    double shell_sld, 
     
    2976    const double shell_h = (0.5*length + thickness); 
    3077    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 
    31     double total = 0.0; 
     78    double total_F1 = 0.0; 
     79    double total_F2 = 0.0; 
    3280    for (int i=0; i<GAUSS_N ;i++) { 
    3381        // translate a point in [-1,1] to a point in [0, pi/2] 
     
    4088        const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 
    4189            + _cyl(shell_vd, shell_r*qab, shell_h*qc); 
    42         total += GAUSS_W[i] * fq * fq * sin_theta; 
     90        total_F1 += GAUSS_W[i] * fq * sin_theta; 
     91        total_F2 += GAUSS_W[i] * fq * fq * sin_theta; 
    4392    } 
    4493    // translate dx in [-1,1] to dx in [lower,upper] 
    4594    //const double form = (upper-lower)/2.0*total; 
    46     return 1.0e-4 * total * M_PI_4; 
     95    *F1 = 1.0e-2 * total_F1 * M_PI_4; 
     96    *F2 = 1.0e-4 * total_F2 * M_PI_4; 
    4797} 
    48  
    4998 
    5099static double 
  • sasmodels/models/core_shell_cylinder.py

    re31b19a r99658f6  
    7070   1445-1452 
    7171.. [#kline] S R Kline, *J Appl. Cryst.*, 39 (2006) 895 
     72L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
    7273 
    7374Authorship and Verification 
     
    131132 
    132133source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "core_shell_cylinder.c"] 
    133  
    134 def ER(radius, thickness, length): 
    135     """ 
    136     Returns the effective radius used in the S*P calculation 
    137     """ 
    138     radius = radius + thickness 
    139     length = length + 2 * thickness 
    140     ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius)) 
    141     return 0.5 * (ddd) ** (1. / 3.) 
     134have_Fq = True 
     135effective_radius_type = [ 
     136    "excluded volume", "equivalent volume sphere", "outer radius", "half outer length", 
     137    "half min outer dimension", "half max outer dimension", "half outer diagonal", 
     138    ] 
    142139 
    143140def random(): 
  • sasmodels/models/core_shell_ellipsoid.c

    r108e70e r99658f6  
    3939 
    4040static double 
    41 Iq(double q, 
     41radius_from_volume(double radius_equat_core, double x_core, double thick_shell, double x_polar_shell) 
     42{ 
     43    const double volume_ellipsoid = form_volume(radius_equat_core, x_core, thick_shell, x_polar_shell); 
     44    return cbrt(volume_ellipsoid/M_4PI_3); 
     45} 
     46 
     47static double 
     48radius_from_curvature(double radius_equat_core, double x_core, double thick_shell, double x_polar_shell) 
     49{ 
     50    // Trivial cases 
     51    if (1.0 == x_core && 1.0 == x_polar_shell) return radius_equat_core + thick_shell; 
     52    if ((radius_equat_core + thick_shell)*(radius_equat_core*x_core + thick_shell*x_polar_shell) == 0.)  return 0.; 
     53 
     54    // see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449 
     55    const double radius_equat_tot = radius_equat_core + thick_shell; 
     56    const double radius_polar_tot = radius_equat_core*x_core + thick_shell*x_polar_shell; 
     57    const double ratio = (radius_polar_tot < radius_equat_tot 
     58                          ? radius_polar_tot / radius_equat_tot 
     59                          : radius_equat_tot / radius_polar_tot); 
     60    const double e1 = sqrt(1.0 - ratio*ratio); 
     61    const double b1 = 1.0 + asin(e1) / (e1 * ratio); 
     62    const double bL = (1.0 + e1) / (1.0 - e1); 
     63    const double b2 = 1.0 + 0.5 * ratio * ratio / e1 * log(bL); 
     64    const double delta = 0.75 * b1 * b2; 
     65    const double ddd = 2.0 * (delta + 1.0) * radius_polar_tot * radius_equat_tot * radius_equat_tot; 
     66    return 0.5 * cbrt(ddd); 
     67} 
     68 
     69static double 
     70effective_radius(int mode, double radius_equat_core, double x_core, double thick_shell, double x_polar_shell) 
     71{ 
     72    const double radius_equat_tot = radius_equat_core + thick_shell; 
     73    const double radius_polar_tot = radius_equat_core*x_core + thick_shell*x_polar_shell; 
     74 
     75    switch (mode) { 
     76    default: 
     77    case 1: // average outer curvature 
     78        return radius_from_curvature(radius_equat_core, x_core, thick_shell, x_polar_shell); 
     79    case 2: // equivalent volume sphere 
     80        return radius_from_volume(radius_equat_core, x_core, thick_shell, x_polar_shell); 
     81    case 3: // min outer radius 
     82        return (radius_polar_tot < radius_equat_tot ? radius_polar_tot : radius_equat_tot); 
     83    case 4: // max outer radius 
     84        return (radius_polar_tot > radius_equat_tot ? radius_polar_tot : radius_equat_tot); 
     85    } 
     86} 
     87 
     88static void 
     89Fq(double q, 
     90    double *F1, 
     91    double *F2, 
    4292    double radius_equat_core, 
    4393    double x_core, 
     
    58108    const double m = 0.5; 
    59109    const double b = 0.5; 
    60     double total = 0.0;     //initialize intergral 
     110    double total_F1 = 0.0;     //initialize intergral 
     111    double total_F2 = 0.0;     //initialize intergral 
    61112    for(int i=0;i<GAUSS_N;i++) { 
    62113        const double cos_theta = GAUSS_Z[i]*m + b; 
     
    66117            equat_shell, polar_shell, 
    67118            sld_core_shell, sld_shell_solvent); 
    68         total += GAUSS_W[i] * fq * fq; 
     119        total_F1 += GAUSS_W[i] * fq; 
     120        total_F2 += GAUSS_W[i] * fq * fq; 
    69121    } 
    70     total *= m; 
     122    total_F1 *= m; 
     123    total_F2 *= m; 
    71124 
    72125    // convert to [cm-1] 
    73     return 1.0e-4 * total; 
     126    *F1 = 1.0e-2 * total_F1; 
     127    *F2 = 1.0e-4 * total_F2; 
    74128} 
     129 
    75130 
    76131static double 
  • sasmodels/models/core_shell_ellipsoid.py

    r2d81cfe r99658f6  
    145145 
    146146source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "core_shell_ellipsoid.c"] 
    147  
    148 def ER(radius_equat_core, x_core, thick_shell, x_polar_shell): 
    149     """ 
    150         Returns the effective radius used in the S*P calculation 
    151     """ 
    152     from .ellipsoid import ER as ellipsoid_ER 
    153     polar_outer = radius_equat_core*x_core + thick_shell*x_polar_shell 
    154     equat_outer = radius_equat_core + thick_shell 
    155     return ellipsoid_ER(polar_outer, equat_outer) 
     147have_Fq = True 
     148effective_radius_type = [ 
     149    "average outer curvature", "equivalent volume sphere",      
     150    "min outer radius", "max outer radius", 
     151    ] 
    156152 
    157153def random(): 
  • sasmodels/models/core_shell_parallelepiped.c

    rdbf1a60 r99658f6  
    2828 
    2929static double 
    30 Iq(double q, 
     30radius_from_excluded_volume(double length_a, double length_b, double length_c, 
     31                   double thick_rim_a, double thick_rim_b, double thick_rim_c) 
     32{ 
     33    double r_equiv, length; 
     34    double lengths[3] = {length_a+thick_rim_a, length_b+thick_rim_b, length_c+thick_rim_c}; 
     35    double lengthmax = fmax(lengths[0],fmax(lengths[1],lengths[2])); 
     36    double length_1 = lengthmax; 
     37    double length_2 = lengthmax; 
     38    double length_3 = lengthmax; 
     39 
     40    for(int ilen=0; ilen<3; ilen++) { 
     41        if (lengths[ilen] < length_1) { 
     42            length_2 = length_1; 
     43            length_1 = lengths[ilen]; 
     44            } else { 
     45                if (lengths[ilen] < length_2) { 
     46                        length_2 = lengths[ilen]; 
     47                } 
     48            } 
     49    } 
     50    if(length_2-length_1 > length_3-length_2) { 
     51        r_equiv = sqrt(length_2*length_3/M_PI); 
     52        length  = length_1; 
     53    } else  { 
     54        r_equiv = sqrt(length_1*length_2/M_PI); 
     55        length  = length_3; 
     56    } 
     57 
     58    return 0.5*cbrt(0.75*r_equiv*(2.0*r_equiv*length + (r_equiv + length)*(M_PI*r_equiv + length))); 
     59} 
     60 
     61static double 
     62radius_from_volume(double length_a, double length_b, double length_c, 
     63                   double thick_rim_a, double thick_rim_b, double thick_rim_c) 
     64{ 
     65    const double volume = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
     66    return cbrt(volume/M_4PI_3); 
     67} 
     68 
     69static double 
     70radius_from_crosssection(double length_a, double length_b, double thick_rim_a, double thick_rim_b) 
     71{ 
     72    const double area_xsec_paral = length_a*length_b + 2.0*thick_rim_a*length_b + 2.0*thick_rim_b*length_a; 
     73    return sqrt(area_xsec_paral/M_PI); 
     74} 
     75 
     76static double 
     77effective_radius(int mode, double length_a, double length_b, double length_c, 
     78                 double thick_rim_a, double thick_rim_b, double thick_rim_c) 
     79{ 
     80    switch (mode) { 
     81    default: 
     82    case 1: // equivalent cylinder excluded volume 
     83        return radius_from_excluded_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
     84    case 2: // equivalent volume sphere 
     85        return radius_from_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
     86    case 3: // half outer length a 
     87        return 0.5 * length_a + thick_rim_a; 
     88    case 4: // half outer length b 
     89        return 0.5 * length_b + thick_rim_b; 
     90    case 5: // half outer length c 
     91        return 0.5 * length_c + thick_rim_c; 
     92    case 6: // equivalent circular cross-section 
     93        return radius_from_crosssection(length_a, length_b, thick_rim_a, thick_rim_b); 
     94    case 7: // half outer ab diagonal 
     95        return 0.5*sqrt(square(length_a+ 2.0*thick_rim_a) + square(length_b+ 2.0*thick_rim_b)); 
     96    case 8: // half outer diagonal 
     97        return 0.5*sqrt(square(length_a+ 2.0*thick_rim_a) + square(length_b+ 2.0*thick_rim_b) + square(length_c+ 2.0*thick_rim_c)); 
     98    } 
     99} 
     100 
     101static void 
     102Fq(double q, 
     103    double *F1, 
     104    double *F2, 
    31105    double core_sld, 
    32106    double arim_sld, 
     
    60134    // outer integral (with gauss points), integration limits = 0, 1 
    61135    // substitute d_cos_alpha for sin_alpha d_alpha 
    62     double outer_sum = 0; //initialize integral 
     136    double outer_sum_F1 = 0; //initialize integral 
     137    double outer_sum_F2 = 0; //initialize integral 
    63138    for( int i=0; i<GAUSS_N; i++) { 
    64139        const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 ); 
     
    69144        // inner integral (with gauss points), integration limits = 0, 1 
    70145        // substitute beta = PI/2 u (so 2/PI * d_(PI/2 * beta) = d_beta) 
    71         double inner_sum = 0.0; 
     146        double inner_sum_F1 = 0.0; 
     147        double inner_sum_F2 = 0.0; 
    72148        for(int j=0; j<GAUSS_N; j++) { 
    73149            const double u = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
     
    91167#endif 
    92168 
    93             inner_sum += GAUSS_W[j] * f * f; 
     169            inner_sum_F1 += GAUSS_W[j] * f; 
     170            inner_sum_F2 += GAUSS_W[j] * f * f; 
    94171        } 
    95172        // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5 
    96         inner_sum *= 0.5; 
    97         // now sum up the outer integral 
    98         outer_sum += GAUSS_W[i] * inner_sum; 
     173        // and sum up the outer integral 
     174        outer_sum_F1 += GAUSS_W[i] * inner_sum_F1 * 0.5; 
     175        outer_sum_F2 += GAUSS_W[i] * inner_sum_F2 * 0.5; 
    99176    } 
    100177    // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5 
    101     outer_sum *= 0.5; 
     178    outer_sum_F1 *= 0.5; 
     179    outer_sum_F2 *= 0.5; 
    102180 
    103181    //convert from [1e-12 A-1] to [cm-1] 
    104     return 1.0e-4 * outer_sum; 
     182    *F1 = 1.0e-2 * outer_sum_F1; 
     183    *F2 = 1.0e-4 * outer_sum_F2; 
    105184} 
    106185 
  • sasmodels/models/core_shell_parallelepiped.py

    rf89ec96 r99658f6  
    9898is the scattering length of the solvent. 
    9999 
    100 .. note::  
     100.. note:: 
    101101 
    102102   the code actually implements two substitutions: $d(cos\alpha)$ is 
     
    173173   from Proquest <http://search.proquest.com/docview/304915826?accountid 
    174174   =26379>`_ 
     175L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
    175176 
    176177Authorship and Verification 
     
    226227 
    227228source = ["lib/gauss76.c", "core_shell_parallelepiped.c"] 
    228  
    229  
    230 def ER(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c): 
    231     """ 
    232         Return equivalent radius (ER) 
    233     """ 
    234     from .parallelepiped import ER as ER_p 
    235  
    236     a = length_a + 2*thick_rim_a 
    237     b = length_b + 2*thick_rim_b 
    238     c = length_c + 2*thick_rim_c 
    239     return ER_p(a, b, c) 
    240  
    241 # VR defaults to 1.0 
     229have_Fq = True 
     230effective_radius_type = [ 
     231    "equivalent cylinder excluded volume",  
     232    "equivalent volume sphere", 
     233    "half outer length_a", "half outer length_b", "half outer length_c", 
     234    "equivalent circular cross-section", 
     235    "half outer ab diagonal", "half outer diagonal", 
     236    ] 
    242237 
    243238def random(): 
  • sasmodels/models/core_shell_sphere.c

    r3a48772 rd42dd4a  
    1 double form_volume(double radius, double thickness); 
    2 double Iq(double q, double radius, double thickness, double core_sld, double shell_sld, double solvent_sld); 
     1static double 
     2form_volume(double radius, double thickness) 
     3{ 
     4    return M_4PI_3 * cube(radius + thickness); 
     5} 
    36 
    4 double Iq(double q, double radius, double thickness, double core_sld, double shell_sld, double solvent_sld) { 
     7static double 
     8effective_radius(int mode, double radius, double thickness) 
     9{ 
     10    switch (mode) { 
     11    default: 
     12    case 1: // outer radius 
     13        return radius + thickness; 
     14    case 2: // core radius 
     15        return radius; 
     16    } 
     17} 
    518 
    6  
    7     double intensity = core_shell_kernel(q, 
     19static void 
     20Fq(double q, double *F1, double *F2, double radius, 
     21   double thickness, double core_sld, double shell_sld, double solvent_sld) { 
     22    double form = core_shell_fq(q, 
    823                              radius, 
    924                              thickness, 
     
    1126                              shell_sld, 
    1227                              solvent_sld); 
    13     return intensity; 
     28    *F1 = 1.0e-2*form; 
     29    *F2 = 1.0e-4*form*form; 
    1430} 
    15  
    16 double form_volume(double radius, double thickness) 
    17 { 
    18     return M_4PI_3 * cube(radius + thickness); 
    19 } 
  • sasmodels/models/core_shell_sphere.py

    rda1c8d1 r304c775  
    5858title = "Form factor for a monodisperse spherical particle with particle with a core-shell structure." 
    5959description = """ 
    60     F^2(q) = 3/V_s [V_c (sld_core-sld_shell) (sin(q*radius)-q*radius*cos(q*radius))/(q*radius)^3 
    61                    + V_s (sld_shell-sld_solvent) (sin(q*r_s)-q*r_s*cos(q*r_s))/(q*r_s)^3] 
     60    F(q) = [V_c (sld_core-sld_shell) 3 (sin(q*radius)-q*radius*cos(q*radius))/(q*radius)^3 
     61            + V_s (sld_shell-sld_solvent) 3 (sin(q*r_s)-q*r_s*cos(q*r_s))/(q*r_s)^3] 
    6262 
    6363            V_s: Volume of the sphere shell 
     
    7777 
    7878source = ["lib/sas_3j1x_x.c", "lib/core_shell.c", "core_shell_sphere.c"] 
     79have_Fq = True 
     80effective_radius_type = ["outer radius", "core radius"] 
    7981 
    8082demo = dict(scale=1, background=0, radius=60, thickness=10, 
    8183            sld_core=1.0, sld_shell=2.0, sld_solvent=0.0) 
    82  
    83 def ER(radius, thickness): 
    84     """ 
    85         Equivalent radius 
    86         @param radius: core radius 
    87         @param thickness: shell thickness 
    88     """ 
    89     return radius + thickness 
    9084 
    9185def random(): 
     
    10296 
    10397tests = [ 
    104     [{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0], 
     98    [{'radius': 20.0, 'thickness': 10.0}, 0.1, None, None, 30.0, 4.*pi/3*30**3, 1.0], 
    10599 
    106100    # The SasView test result was 0.00169, with a background of 0.001 
  • sasmodels/models/cylinder.c

    r108e70e r99658f6  
    88 
    99static double 
    10 fq(double qab, double qc, double radius, double length) 
     10_fq(double qab, double qc, double radius, double length) 
    1111{ 
    1212    return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length); 
     
    1414 
    1515static double 
    16 orient_avg_1D(double q, double radius, double length) 
     16radius_from_excluded_volume(double radius, double length) 
     17{ 
     18    return 0.5*cbrt(0.75*radius*(2.0*radius*length + (radius + length)*(M_PI*radius + length))); 
     19} 
     20 
     21static double 
     22radius_from_volume(double radius, double length) 
     23{ 
     24    return cbrt(M_PI*radius*radius*length/M_4PI_3); 
     25} 
     26 
     27static double 
     28radius_from_diagonal(double radius, double length) 
     29{ 
     30    return sqrt(radius*radius + 0.25*length*length); 
     31} 
     32 
     33static double 
     34effective_radius(int mode, double radius, double length) 
     35{ 
     36    switch (mode) { 
     37    default: 
     38    case 1: 
     39        return radius_from_excluded_volume(radius, length); 
     40    case 2: 
     41        return radius_from_volume(radius, length); 
     42    case 3: 
     43        return radius; 
     44    case 4: 
     45        return 0.5*length; 
     46    case 5: 
     47        return (radius < 0.5*length ? radius : 0.5*length); 
     48    case 6: 
     49        return (radius > 0.5*length ? radius : 0.5*length); 
     50    case 7: 
     51        return radius_from_diagonal(radius,length); 
     52    } 
     53} 
     54 
     55static void 
     56Fq(double q, 
     57    double *F1, 
     58    double *F2, 
     59    double sld, 
     60    double solvent_sld, 
     61    double radius, 
     62    double length) 
    1763{ 
    1864    // translate a point in [-1,1] to a point in [0, pi/2] 
     
    2066    const double zb = M_PI_4; 
    2167 
    22     double total = 0.0; 
     68    double total_F1 = 0.0; 
     69    double total_F2 = 0.0; 
    2370    for (int i=0; i<GAUSS_N ;i++) { 
    2471        const double theta = GAUSS_Z[i]*zm + zb; 
     
    2673        // theta (theta,phi) the projection of the cylinder on the detector plane 
    2774        SINCOS(theta , sin_theta, cos_theta); 
    28         const double form = fq(q*sin_theta, q*cos_theta, radius, length); 
    29         total += GAUSS_W[i] * form * form * sin_theta; 
     75        const double form = _fq(q*sin_theta, q*cos_theta, radius, length); 
     76        total_F1 += GAUSS_W[i] * form * sin_theta; 
     77        total_F2 += GAUSS_W[i] * form * form * sin_theta; 
    3078    } 
    3179    // translate dx in [-1,1] to dx in [lower,upper] 
    32     return total*zm; 
     80    total_F1 *= zm; 
     81    total_F2 *= zm; 
     82    const double s = (sld - solvent_sld) * form_volume(radius, length); 
     83    *F1 = 1e-2 * s * total_F1; 
     84    *F2 = 1e-4 * s * s * total_F2; 
    3385} 
    3486 
    35 static double 
    36 Iq(double q, 
    37     double sld, 
    38     double solvent_sld, 
    39     double radius, 
    40     double length) 
    41 { 
    42     const double s = (sld - solvent_sld) * form_volume(radius, length); 
    43     return 1.0e-4 * s * s * orient_avg_1D(q, radius, length); 
    44 } 
     87 
    4588 
    4689static double 
     
    5194    double length) 
    5295{ 
     96    const double form = _fq(qab, qc, radius, length); 
    5397    const double s = (sld-solvent_sld) * form_volume(radius, length); 
    54     const double form = fq(qab, qc, radius, length); 
    5598    return 1.0e-4 * square(s * form); 
    5699} 
     100 
  • sasmodels/models/cylinder.py

    r2d81cfe r99658f6  
    9898J. S. Pedersen, Adv. Colloid Interface Sci. 70, 171-210 (1997). 
    9999G. Fournet, Bull. Soc. Fr. Mineral. Cristallogr. 74, 39-113 (1951). 
     100L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
    100101""" 
    101102 
     
    138139 
    139140source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "cylinder.c"] 
    140  
    141 def ER(radius, length): 
    142     """ 
    143         Return equivalent radius (ER) 
    144     """ 
    145     ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius)) 
    146     return 0.5 * (ddd) ** (1. / 3.) 
     141have_Fq = True 
     142effective_radius_type = [ 
     143    "excluded volume", "equivalent volume sphere", "radius", 
     144    "half length", "half min dimension", "half max dimension", "half diagonal", 
     145    ] 
    147146 
    148147def random(): 
     
    169168            phi_pd=10, phi_pd_n=5) 
    170169 
     170# pylint: disable=bad-whitespace, line-too-long 
    171171qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) 
    172172# After redefinition of angles, find new tests values.  Was 10 10 in old coords 
     
    182182] 
    183183del qx, qy  # not necessary to delete, but cleaner 
     184 
     185# Default radius and length 
     186radius, length = parameters[2][2], parameters[3][2] 
     187tests.extend([ 
     188    ({'radius_effective_mode': 0}, 0.1, None, None, 0., pi*radius**2*length, 1.0),    
     189    ({'radius_effective_mode': 1}, 0.1, None, None, 0.5*(0.75*radius*(2.0*radius*length + (radius + length)*(pi*radius + length)))**(1./3.), None, None),     
     190    ({'radius_effective_mode': 2}, 0.1, None, None, (0.75*radius**2*length)**(1./3.), None, None), 
     191    ({'radius_effective_mode': 3}, 0.1, None, None, radius, None, None), 
     192    ({'radius_effective_mode': 4}, 0.1, None, None, length/2., None, None), 
     193    ({'radius_effective_mode': 5}, 0.1, None, None, min(radius, length/2.), None, None), 
     194    ({'radius_effective_mode': 6}, 0.1, None, None, max(radius, length/2.), None, None), 
     195    ({'radius_effective_mode': 7}, 0.1, None, None, np.sqrt(4*radius**2 + length**2)/2., None, None), 
     196]) 
     197del radius, length 
     198# pylint: enable=bad-whitespace, line-too-long 
     199 
    184200# ADDED by:  RKH  ON: 18Mar2016 renamed sld's etc 
  • sasmodels/models/dab.py

    r2d81cfe r304c775  
    7474    return pars 
    7575 
    76 # ER defaults to 1.0 
    77  
    78 # VR defaults to 1.0 
    79  
    8076demo = dict(scale=1, background=0, cor_length=50) 
  • sasmodels/models/ellipsoid.c

    r108e70e r99658f6  
    55} 
    66 
    7 static  double 
    8 Iq(double q, 
     7static double 
     8radius_from_volume(double radius_polar, double radius_equatorial) 
     9{ 
     10    return cbrt(radius_polar*radius_equatorial*radius_equatorial); 
     11} 
     12 
     13static double 
     14radius_from_curvature(double radius_polar, double radius_equatorial) 
     15{ 
     16    // Trivial cases 
     17    if (radius_polar == radius_equatorial) return radius_polar; 
     18    if (radius_polar * radius_equatorial == 0.)  return 0.; 
     19 
     20    // see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449 
     21    const double ratio = (radius_polar < radius_equatorial 
     22                          ? radius_polar / radius_equatorial 
     23                          : radius_equatorial / radius_polar); 
     24    const double e1 = sqrt(1.0 - ratio*ratio); 
     25    const double b1 = 1.0 + asin(e1) / (e1 * ratio); 
     26    const double bL = (1.0 + e1) / (1.0 - e1); 
     27    const double b2 = 1.0 + 0.5 * ratio * ratio / e1 * log(bL); 
     28    const double delta = 0.75 * b1 * b2; 
     29    const double ddd = 2.0 * (delta + 1.0) * radius_polar * radius_equatorial * radius_equatorial; 
     30    return 0.5 * cbrt(ddd); 
     31} 
     32 
     33static double 
     34effective_radius(int mode, double radius_polar, double radius_equatorial) 
     35{ 
     36    switch (mode) { 
     37    default: 
     38    case 1: // average curvature 
     39        return radius_from_curvature(radius_polar, radius_equatorial); 
     40    case 2: // equivalent volume sphere 
     41        return radius_from_volume(radius_polar, radius_equatorial); 
     42    case 3: // min radius 
     43        return (radius_polar < radius_equatorial ? radius_polar : radius_equatorial); 
     44    case 4: // max radius 
     45        return (radius_polar > radius_equatorial ? radius_polar : radius_equatorial); 
     46    } 
     47} 
     48 
     49 
     50static void 
     51Fq(double q, 
     52    double *F1, 
     53    double *F2, 
    954    double sld, 
    1055    double sld_solvent, 
     
    2065    //     i(h) = int_0^1 Phi^2(h a sqrt(1 + u^2(v^2-1)) du 
    2166    const double v_square_minus_one = square(radius_polar/radius_equatorial) - 1.0; 
    22  
    2367    // translate a point in [-1,1] to a point in [0, 1] 
    2468    // const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper+lower)/2; 
    2569    const double zm = 0.5; 
    2670    const double zb = 0.5; 
    27     double total = 0.0; 
     71    double total_F2 = 0.0; 
     72    double total_F1 = 0.0; 
    2873    for (int i=0;i<GAUSS_N;i++) { 
    2974        const double u = GAUSS_Z[i]*zm + zb; 
    3075        const double r = radius_equatorial*sqrt(1.0 + u*u*v_square_minus_one); 
    3176        const double f = sas_3j1x_x(q*r); 
    32         total += GAUSS_W[i] * f * f; 
     77        total_F2 += GAUSS_W[i] * f * f; 
     78        total_F1 += GAUSS_W[i] * f; 
    3379    } 
    3480    // translate dx in [-1,1] to dx in [lower,upper] 
    35     const double form = total*zm; 
     81    total_F1 *= zm; 
     82    total_F2 *= zm; 
    3683    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    37     return 1.0e-4 * s * s * form; 
     84    *F1 = 1e-2 * s * total_F1; 
     85    *F2 = 1e-4 * s * s * total_F2; 
    3886} 
    3987 
  • sasmodels/models/ellipsoid.py

    r0168844 r99658f6  
    166166             ] 
    167167 
     168 
    168169source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "ellipsoid.c"] 
    169  
    170 def ER(radius_polar, radius_equatorial): 
    171     # see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449 
    172     ee = np.empty_like(radius_polar) 
    173     idx = radius_polar > radius_equatorial 
    174     ee[idx] = (radius_polar[idx] ** 2 - radius_equatorial[idx] ** 2) / radius_polar[idx] ** 2 
    175     idx = radius_polar < radius_equatorial 
    176     ee[idx] = (radius_equatorial[idx] ** 2 - radius_polar[idx] ** 2) / radius_equatorial[idx] ** 2 
    177     valid = (radius_polar * radius_equatorial != 0) & (radius_polar != radius_equatorial) 
    178     bd = 1.0 - ee[valid] 
    179     e1 = np.sqrt(ee[valid]) 
    180     b1 = 1.0 + np.arcsin(e1) / (e1 * np.sqrt(bd)) 
    181     bL = (1.0 + e1) / (1.0 - e1) 
    182     b2 = 1.0 + bd / 2 / e1 * np.log(bL) 
    183     delta = 0.75 * b1 * b2 
    184     ddd = 2.0 * (delta + 1.0) * (radius_polar * radius_equatorial**2)[valid] 
    185  
    186     r = np.zeros_like(radius_polar) 
    187     r[valid] = 0.5 * cbrt(ddd) 
    188     idx = radius_polar == radius_equatorial 
    189     r[idx] = radius_polar[idx] 
    190     return r 
     170have_Fq = True 
     171effective_radius_type = [ 
     172    "average curvature", "equivalent volume sphere", "min radius", "max radius", 
     173    ] 
    191174 
    192175def random(): 
  • sasmodels/models/elliptical_cylinder.c

    r108e70e r99658f6  
    66 
    77static double 
    8 Iq(double q, double radius_minor, double r_ratio, double length, 
     8radius_from_excluded_volume(double radius_minor, double r_ratio, double length) 
     9{ 
     10    const double r_equiv = sqrt(radius_minor*radius_minor*r_ratio); 
     11    return 0.5*cbrt(0.75*r_equiv*(2.0*r_equiv*length + (r_equiv + length)*(M_PI*r_equiv + length))); 
     12} 
     13 
     14static double 
     15radius_from_volume(double radius_minor, double r_ratio, double length) 
     16{ 
     17    const double volume_ellcyl = form_volume(radius_minor,r_ratio,length); 
     18    return cbrt(volume_ellcyl/M_4PI_3); 
     19} 
     20 
     21static double 
     22radius_from_min_dimension(double radius_minor, double r_ratio, double hlength) 
     23{ 
     24    const double rad_min = (r_ratio > 1.0 ? radius_minor : r_ratio*radius_minor); 
     25    return (rad_min < hlength ? rad_min : hlength); 
     26} 
     27 
     28static double 
     29radius_from_max_dimension(double radius_minor, double r_ratio, double hlength) 
     30{ 
     31    const double rad_max = (r_ratio < 1.0 ? radius_minor : r_ratio*radius_minor); 
     32    return (rad_max > hlength ? rad_max : hlength); 
     33} 
     34 
     35static double 
     36radius_from_diagonal(double radius_minor, double r_ratio, double length) 
     37{ 
     38    const double radius_max = (r_ratio > 1.0 ? radius_minor*r_ratio : radius_minor); 
     39    return sqrt(radius_max*radius_max + 0.25*length*length); 
     40} 
     41 
     42static double 
     43effective_radius(int mode, double radius_minor, double r_ratio, double length) 
     44{ 
     45    switch (mode) { 
     46    default: 
     47    case 1: // equivalent cylinder excluded volume 
     48        return radius_from_excluded_volume(radius_minor, r_ratio, length); 
     49    case 2: // equivalent volume sphere 
     50        return radius_from_volume(radius_minor, r_ratio, length); 
     51    case 3: // average radius 
     52        return 0.5*radius_minor*(1.0 + r_ratio); 
     53    case 4: // min radius 
     54        return (r_ratio > 1.0 ? radius_minor : r_ratio*radius_minor); 
     55    case 5: // max radius 
     56        return (r_ratio < 1.0 ? radius_minor : r_ratio*radius_minor); 
     57    case 6: // equivalent circular cross-section 
     58        return sqrt(radius_minor*radius_minor*r_ratio); 
     59    case 7: // half length 
     60        return 0.5*length; 
     61    case 8: // half min dimension 
     62        return radius_from_min_dimension(radius_minor,r_ratio,0.5*length); 
     63    case 9: // half max dimension 
     64        return radius_from_max_dimension(radius_minor,r_ratio,0.5*length); 
     65    case 10: // half diagonal 
     66        return radius_from_diagonal(radius_minor,r_ratio,length); 
     67    } 
     68} 
     69 
     70static void 
     71Fq(double q, double *F1, double *F2, double radius_minor, double r_ratio, double length, 
    972   double sld, double solvent_sld) 
    1073{ 
     
    2184 
    2285    //initialize integral 
    23     double outer_sum = 0.0; 
     86    double outer_sum_F1 = 0.0; 
     87    double outer_sum_F2 = 0.0; 
    2488    for(int i=0;i<GAUSS_N;i++) { 
    2589        //setup inner integral over the ellipsoidal cross-section 
     
    2791        const double sin_val = sqrt(1.0 - cos_val*cos_val); 
    2892        //const double arg = radius_minor*sin_val; 
    29         double inner_sum=0; 
     93        double inner_sum_F1 = 0.0; 
     94        double inner_sum_F2 = 0.0; 
    3095        for(int j=0;j<GAUSS_N;j++) { 
    3196            const double theta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    3297            const double r = sin_val*sqrt(rA - rB*cos(theta)); 
    3398            const double be = sas_2J1x_x(q*r); 
    34             inner_sum += GAUSS_W[j] * be * be; 
     99            inner_sum_F1 += GAUSS_W[j] * be; 
     100            inner_sum_F2 += GAUSS_W[j] * be * be; 
    35101        } 
    36102        //now calculate the value of the inner integral 
    37         inner_sum *= 0.5*(vbj-vaj); 
     103        inner_sum_F1 *= 0.5*(vbj-vaj); 
     104        inner_sum_F2 *= 0.5*(vbj-vaj); 
    38105 
    39106        //now calculate outer integral 
    40107        const double si = sas_sinx_x(q*0.5*length*cos_val); 
    41         outer_sum += GAUSS_W[i] * inner_sum * si * si; 
     108        outer_sum_F1 += GAUSS_W[i] * inner_sum_F1 * si; 
     109        outer_sum_F2 += GAUSS_W[i] * inner_sum_F2 * si * si; 
    42110    } 
    43     outer_sum *= 0.5*(vb-va); 
    44  
    45     //divide integral by Pi 
    46     const double form = outer_sum/M_PI; 
     111    // correct limits and divide integral by pi 
     112    outer_sum_F1 *= 0.5*(vb-va)/M_PI; 
     113    outer_sum_F2 *= 0.5*(vb-va)/M_PI; 
    47114 
    48115    // scale by contrast and volume, and convert to to 1/cm units 
    49     const double vol = form_volume(radius_minor, r_ratio, length); 
    50     const double delrho = sld - solvent_sld; 
    51     return 1.0e-4*square(delrho*vol)*form; 
     116    const double volume = form_volume(radius_minor, r_ratio, length); 
     117    const double contrast = sld - solvent_sld; 
     118    const double s = contrast*volume; 
     119    *F1 = 1.0e-2*s*outer_sum_F1; 
     120    *F2 = 1.0e-4*s*s*outer_sum_F2; 
    52121} 
    53122 
     
    63132    const double be = sas_2J1x_x(qr); 
    64133    const double si = sas_sinx_x(qc*0.5*length); 
    65     const double Aq = be * si; 
    66     const double delrho = sld - solvent_sld; 
    67     const double vol = form_volume(radius_minor, r_ratio, length); 
    68     return 1.0e-4 * square(delrho * vol * Aq); 
     134    const double fq = be * si; 
     135    const double contrast = sld - solvent_sld; 
     136    const double volume = form_volume(radius_minor, r_ratio, length); 
     137    return 1.0e-4 * square(contrast * volume * fq); 
    69138} 
  • sasmodels/models/elliptical_cylinder.py

    ra261a83 r99658f6  
    8888L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and 
    8989Neutron Scattering*, Plenum, New York, (1987) [see table 3.4] 
     90L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
    9091 
    9192Authorship and Verification 
     
    122123 
    123124source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "elliptical_cylinder.c"] 
     125have_Fq = True 
     126effective_radius_type = [ 
     127    "equivalent cylinder excluded volume", "equivalent volume sphere", "average radius", "min radius", "max radius", 
     128    "equivalent circular cross-section", 
     129    "half length", "half min dimension", "half max dimension", "half diagonal", 
     130    ] 
    124131 
    125132demo = dict(scale=1, background=0, radius_minor=100, axis_ratio=1.5, length=400.0, 
    126133            sld=4.0, sld_solvent=1.0, theta=10.0, phi=20, psi=30, 
    127134            theta_pd=10, phi_pd=2, psi_pd=3) 
    128  
    129 def ER(radius_minor, axis_ratio, length): 
    130     """ 
    131         Equivalent radius 
    132         @param radius_minor: Ellipse minor radius 
    133         @param axis_ratio: Ratio of major radius over minor radius 
    134         @param length: Length of the cylinder 
    135     """ 
    136     radius = sqrt(radius_minor * radius_minor * axis_ratio) 
    137     ddd = 0.75 * radius * (2 * radius * length 
    138                            + (length + radius) * (length + pi * radius)) 
    139     return 0.5 * (ddd) ** (1. / 3.) 
    140135 
    141136def random(): 
     
    161156 
    162157tests = [ 
    163     [{'radius_minor': 20.0, 'axis_ratio': 1.5, 'length':400.0}, 'ER', 79.89245454155024], 
    164     [{'radius_minor': 20.0, 'axis_ratio': 1.2, 'length':300.0}, 'VR', 1], 
     158#    [{'radius_minor': 20.0, 'axis_ratio': 1.5, 'length':400.0}, 'ER', 79.89245454155024], 
     159#    [{'radius_minor': 20.0, 'axis_ratio': 1.2, 'length':300.0}, 'VR', 1], 
    165160 
    166161    # The SasView test result was 0.00169, with a background of 0.001 
  • sasmodels/models/fcc_paracrystal.c

    r108e70e r71b751d  
    7979} 
    8080 
    81  
    8281static double Iqabc(double qa, double qb, double qc, 
    8382    double dnn, double d_factor, double radius, 
  • sasmodels/models/flexible_cylinder_elliptical.c

    r74768cb r71b751d  
    7272    return result; 
    7373} 
    74  
  • sasmodels/models/fractal.c

    r4788822 r71b751d  
    1919 
    2020    //calculate P(q) for the spherical subunits 
    21     const double pq = square(form_volume(radius) * (sld_block-sld_solvent) 
    22                       *sas_3j1x_x(q*radius)); 
     21    const double fq = form_volume(radius) * (sld_block-sld_solvent) 
     22                      *sas_3j1x_x(q*radius); 
    2323 
    2424    // scale to units cm-1 sr-1 (assuming data on absolute scale) 
     
    2727    //    combined: 1e-4 * I(q) 
    2828 
    29     return 1.e-4 * volfraction * sq * pq; 
     29    return 1.e-4 * volfraction * sq * fq * fq; 
    3030} 
    31  
  • sasmodels/models/fractal_core_shell.c

    rbdd08df rbad3093  
    1616   double cor_length) 
    1717{ 
    18     //The radius for the building block of the core shell particle that is  
     18    //The radius for the building block of the core shell particle that is 
    1919    //needed by the Teixeira fractal S(q) is the radius of the whole particle. 
    2020    const double cs_radius = radius + thickness; 
    2121    const double sq = fractal_sq(q, cs_radius, fractal_dim, cor_length); 
    22     const double pq = core_shell_kernel(q, radius, thickness, 
    23                                         core_sld, shell_sld, solvent_sld); 
     22    const double fq = core_shell_fq(q, radius, thickness, 
     23                                    core_sld, shell_sld, solvent_sld); 
    2424 
    25     // Note: core_shell_kernel already performs the 1e-4 unit conversion 
    26     return volfraction * sq * pq; 
     25    return 1.0e-4 * volfraction * sq * fq * fq; 
    2726} 
    28  
  • sasmodels/models/fractal_core_shell.py

    reb3eb38 r2cc8aa2  
    134134    return radius + thickness 
    135135 
    136 tests = [[{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0], 
    137  
     136#tests = [[{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0], 
     137tests = [ 
    138138#         # At some point the SasView 3.x test result was deemed incorrect. The 
    139139          #following tests were verified against NIST IGOR macros ver 7.850. 
  • sasmodels/models/fuzzy_sphere.py

    r2d81cfe ree60aa7  
    7878              ["sld_solvent", "1e-6/Ang^2",  3, [-inf, inf], "sld",    "Solvent scattering length density"], 
    7979              ["radius",      "Ang",        60, [0, inf],    "volume", "Sphere radius"], 
    80               ["fuzziness",   "Ang",        10, [0, inf],    "",       "std deviation of Gaussian convolution for interface (must be << radius)"], 
     80              ["fuzziness",   "Ang",        10, [0, inf],    "volume",       "std deviation of Gaussian convolution for interface (must be << radius)"], 
    8181             ] 
    8282# pylint: enable=bad-whitespace,line-too-long 
    8383 
    84 source = ["lib/sas_3j1x_x.c"] 
    85  
    86 # No volume normalization despite having a volume parameter 
    87 # This should perhaps be volume normalized? 
    88 form_volume = """ 
    89     return M_4PI_3*cube(radius); 
    90     """ 
    91  
    92 Iq = """ 
    93     const double qr = q*radius; 
    94     const double bes = sas_3j1x_x(qr); 
    95     const double qf = q*fuzziness; 
    96     const double fq = bes * (sld - sld_solvent) * form_volume(radius) * exp(-0.5*qf*qf); 
    97     return 1.0e-4*fq*fq; 
    98     """ 
    99  
    100 def ER(radius): 
    101     """ 
    102     Return radius 
    103     """ 
    104     return radius 
    105  
    106 # VR defaults to 1.0 
     84source = ["lib/sas_3j1x_x.c","fuzzy_sphere.c"] 
     85have_Fq = True 
     86effective_radius_type = ["radius", "radius + fuzziness"] 
    10787 
    10888def random(): 
  • sasmodels/models/gaussian_peak.py

    r2d81cfe r304c775  
    5050    """ 
    5151 
    52 # VR defaults to 1.0 
    53  
    5452def random(): 
    5553    peak_pos = 10**np.random.uniform(-3, -1) 
  • sasmodels/models/hardsphere.py

    r2d81cfe rc1799d3  
    6262 
    6363#             ["name", "units", default, [lower, upper], "type","description"], 
    64 parameters = [["radius_effective", "Ang", 50.0, [0, inf], "volume", 
     64parameters = [["radius_effective", "Ang", 50.0, [0, inf], "", 
    6565               "effective radius of hard sphere"], 
    6666              ["volfraction", "", 0.2, [0, 0.74], "", 
    6767               "volume fraction of hard spheres"], 
    6868             ] 
    69  
    70 # No volume normalization despite having a volume parameter 
    71 # This should perhaps be volume normalized? 
    72 form_volume = """ 
    73     return 1.0; 
    74     """ 
    7569 
    7670Iq = r""" 
     
    168162    return pars 
    169163 
    170 # ER defaults to 0.0 
    171 # VR defaults to 1.0 
    172  
    173 demo = dict(radius_effective=200, volfraction=0.2, 
    174             radius_effective_pd=0.1, radius_effective_pd_n=40) 
    175164# Q=0.001 is in the Taylor series, low Q part, so add Q=0.1, 
    176165# assuming double precision sasview is correct 
  • sasmodels/models/hayter_msa.py

    r2d81cfe r304c775  
    8989    return 1.0; 
    9090    """ 
    91 # ER defaults to 0.0 
    92 # VR defaults to 1.0 
    9391 
    9492def random(): 
  • sasmodels/models/hollow_cylinder.c

    r108e70e r99658f6  
    11//#define INVALID(v) (v.radius_core >= v.radius) 
    22 
    3 // From Igor library 
    43static double 
    5 _hollow_cylinder_scaling(double integrand, double delrho, double volume) 
     4shell_volume(double radius, double thickness, double length) 
    65{ 
    7     return 1.0e-4 * square(volume * delrho) * integrand; 
     6    return M_PI*length*(square(radius+thickness) - radius*radius); 
     7} 
     8 
     9static double 
     10form_volume(double radius, double thickness, double length) 
     11{ 
     12    return M_PI*length*square(radius+thickness); 
     13} 
     14 
     15static double 
     16radius_from_excluded_volume(double radius, double thickness, double length) 
     17{ 
     18    const double radius_tot = radius + thickness; 
     19    return 0.5*cbrt(0.75*radius_tot*(2.0*radius_tot*length + (radius_tot + length)*(M_PI*radius_tot + length))); 
     20} 
     21 
     22static double 
     23radius_from_volume(double radius, double thickness, double length) 
     24{ 
     25    const double volume_outer_cyl = M_PI*square(radius + thickness)*length; 
     26    return cbrt(volume_outer_cyl/M_4PI_3); 
     27} 
     28 
     29static double 
     30radius_from_diagonal(double radius, double thickness, double length) 
     31{ 
     32    return sqrt(square(radius + thickness) + 0.25*square(length)); 
     33} 
     34 
     35static double 
     36effective_radius(int mode, double radius, double thickness, double length) 
     37{ 
     38    switch (mode) { 
     39    default: 
     40    case 1: // excluded volume 
     41        return radius_from_excluded_volume(radius, thickness, length); 
     42    case 2: // equivalent volume sphere 
     43        return radius_from_volume(radius, thickness, length); 
     44    case 3: // outer radius 
     45        return radius + thickness; 
     46    case 4: // half length 
     47        return 0.5*length; 
     48    case 5: // half outer min dimension 
     49        return (radius + thickness < 0.5*length ? radius + thickness : 0.5*length); 
     50    case 6: // half outer max dimension 
     51        return (radius + thickness > 0.5*length ? radius + thickness : 0.5*length); 
     52    case 7: // half outer diagonal 
     53        return radius_from_diagonal(radius,thickness,length); 
     54    } 
    855} 
    956 
     
    2269} 
    2370 
    24 static double 
    25 form_volume(double radius, double thickness, double length) 
    26 { 
    27     double v_shell = M_PI*length*(square(radius+thickness) - radius*radius); 
    28     return v_shell; 
    29 } 
    30  
    31  
    32 static double 
    33 Iq(double q, double radius, double thickness, double length, 
     71static void 
     72Fq(double q, double *F1, double *F2, double radius, double thickness, double length, 
    3473    double sld, double solvent_sld) 
    3574{ 
     
    3776    const double upper = 1.0;        //limits of numerical integral 
    3877 
    39     double summ = 0.0;            //initialize intergral 
     78    double total_F1 = 0.0;            //initialize intergral 
     79    double total_F2 = 0.0; 
    4080    for (int i=0;i<GAUSS_N;i++)