Changeset 29c455a in sasmodels


Ignore:
Timestamp:
Mar 19, 2019 10:49:29 AM (2 months ago)
Author:
GitHub <noreply@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
0d362b7
Parents:
37f38ff (diff), edc783e (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
git-author:
Paul Butler <butlerpd@…> (03/19/19 10:49:29)
git-committer:
GitHub <noreply@…> (03/19/19 10:49:29)
Message:

Merge pull request #90 from SasView?/beta_approx

Beta approx

closes #1201
closes #1202
closes #822
closes #1076

Files:
50 added
115 edited

Legend:

Unmodified
Added
Removed
  • .travis.yml

    r5c36bf1 rf3767c2  
    99  - os: linux 
    1010    env: 
    11     - PY=3.6 
     11    - PY=3.7 
    1212  - os: osx 
    1313    language: generic 
     
    1717    language: generic 
    1818    env: 
    19     - PY=3.5 
     19    - PY=3.7 
    2020branches: 
    2121  only: 
  • 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

    r94bfa42 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 
     
    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/asymint.py

    rc462169 r66dbbfb  
    2727import os, sys 
    2828sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.realpath(__file__)))) 
     29import warnings 
    2930 
    3031import numpy as np 
     
    3637 
    3738import sasmodels.special as sp 
    38  
    39 # Need to parse shape early since it determines the kernel function 
    40 # that will be used for the various integrators 
    41 shape = 'parallelepiped' if len(sys.argv) < 2 else sys.argv[1] 
    42 Qstr = '0.005' if len(sys.argv) < 3 else sys.argv[2] 
    4339 
    4440DTYPE = 'd' 
     
    201197    return norm, Fq 
    202198 
    203 if shape == 'sphere': 
    204     RADIUS = 50  # integer for the sake of mpf 
    205     NORM, KERNEL = make_sphere(radius=RADIUS) 
    206     NORM_MP, KERNEL_MP = make_sphere(radius=RADIUS, env=MPenv) 
    207 elif shape == 'cylinder': 
    208     #RADIUS, LENGTH = 10, 100000 
    209     RADIUS, LENGTH = 10, 300  # integer for the sake of mpf 
    210     NORM, KERNEL = make_cylinder(radius=RADIUS, length=LENGTH) 
    211     NORM_MP, KERNEL_MP = make_cylinder(radius=RADIUS, length=LENGTH, env=MPenv) 
    212 elif shape == 'triaxial_ellipsoid': 
    213     #A, B, C = 4450, 14000, 47 
    214     A, B, C = 445, 140, 47  # integer for the sake of mpf 
    215     NORM, KERNEL = make_triaxial_ellipsoid(A, B, C) 
    216     NORM_MP, KERNEL_MP = make_triaxial_ellipsoid(A, B, C, env=MPenv) 
    217 elif shape == 'parallelepiped': 
    218     #A, B, C = 4450, 14000, 47 
    219     A, B, C = 445, 140, 47  # integer for the sake of mpf 
    220     NORM, KERNEL = make_parallelepiped(A, B, C) 
    221     NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv) 
    222 elif shape == 'core_shell_parallelepiped': 
    223     #A, B, C = 4450, 14000, 47 
    224     #A, B, C = 445, 140, 47  # integer for the sake of mpf 
    225     A, B, C = 114, 1380, 6800 
    226     DA, DB, DC = 21, 58, 2300 
    227     SLDA, SLDB, SLDC = "5", "-0.3", "11.5" 
    228     ## default parameters from sasmodels 
    229     #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 400,75,35,10,10,10,2,4,2 
    230     ## swap A-B-C to C-B-A 
    231     #A, B, C, DA, DB, DC, SLDA, SLDB, SLDC = C, B, A, DC, DB, DA, SLDC, SLDB, SLDA 
    232     #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 10,20,30,100,200,300,1,2,3 
    233     #SLD_SOLVENT,CONTRAST = 0, 4 
    234     if 1: # C shortest 
    235         B, C = C, B 
    236         DB, DC = DC, DB 
    237         SLDB, SLDC = SLDC, SLDB 
    238     elif 0: # C longest 
    239         A, C = C, A 
    240         DA, DC = DC, DA 
    241         SLDA, SLDC = SLDC, SLDA 
    242     #NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 
    243     NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 
    244     NORM_MP, KERNEL_MP = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC, env=MPenv) 
    245 elif shape == 'paracrystal': 
    246     LATTICE = 'bcc' 
    247     #LATTICE = 'fcc' 
    248     #LATTICE = 'sc' 
    249     DNN, D_FACTOR = 220, '0.06'  # mpmath needs to initialize floats from string 
    250     RADIUS = 40  # integer for the sake of mpf 
    251     NORM, KERNEL = make_paracrystal( 
    252         radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE) 
    253     NORM_MP, KERNEL_MP = make_paracrystal( 
    254         radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE, env=MPenv) 
    255 else: 
    256     raise ValueError("Unknown shape %r"%shape) 
     199NORM = 1.0  # type: float 
     200KERNEL = None  # type: CALLABLE[[ndarray, ndarray, ndarray], ndarray] 
     201NORM_MP = 1  # type: mpf 
     202KERNEL = None  # type: CALLABLE[[mpf, mpf, mpf], mpf] 
     203 
     204SHAPES = [ 
     205    'sphere', 
     206    'cylinder', 
     207    'triaxial_ellipsoid', 
     208    'parallelepiped', 
     209    'core_shell_parallelepiped', 
     210    'fcc_paracrystal', 
     211    'bcc_paracrystal', 
     212    'sc_paracrystal', 
     213] 
     214def build_shape(shape, **pars): 
     215    global NORM, KERNEL 
     216    global NORM_MP, KERNEL_MP 
     217 
     218    # Note: using integer or string defaults for the sake of mpf 
     219    if shape == 'sphere': 
     220        RADIUS = pars.get('radius', 50) 
     221        NORM, KERNEL = make_sphere(radius=RADIUS) 
     222        NORM_MP, KERNEL_MP = make_sphere(radius=RADIUS, env=MPenv) 
     223    elif shape == 'cylinder': 
     224        #RADIUS, LENGTH = 10, 100000 
     225        RADIUS = pars.get('radius', 10) 
     226        LENGTH = pars.get('radius', 300) 
     227        NORM, KERNEL = make_cylinder(radius=RADIUS, length=LENGTH) 
     228        NORM_MP, KERNEL_MP = make_cylinder(radius=RADIUS, length=LENGTH, env=MPenv) 
     229    elif shape == 'triaxial_ellipsoid': 
     230        #A, B, C = 4450, 14000, 47 
     231        A = pars.get('a', 445) 
     232        B = pars.get('b', 140) 
     233        C = pars.get('c', 47) 
     234        NORM, KERNEL = make_triaxial_ellipsoid(A, B, C) 
     235        NORM_MP, KERNEL_MP = make_triaxial_ellipsoid(A, B, C, env=MPenv) 
     236    elif shape == 'parallelepiped': 
     237        #A, B, C = 4450, 14000, 47 
     238        A = pars.get('a', 445) 
     239        B = pars.get('b', 140) 
     240        C = pars.get('c', 47) 
     241        NORM, KERNEL = make_parallelepiped(A, B, C) 
     242        NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv) 
     243    elif shape == 'core_shell_parallelepiped': 
     244        #A, B, C = 4450, 14000, 47 
     245        #A, B, C = 445, 140, 47  # integer for the sake of mpf 
     246        A = pars.get('a', 114) 
     247        B = pars.get('b', 1380) 
     248        C = pars.get('c', 6800) 
     249        DA = pars.get('da', 21) 
     250        DB = pars.get('db', 58) 
     251        DC = pars.get('dc', 2300) 
     252        SLDA = pars.get('slda', "5") 
     253        SLDB = pars.get('sldb', "-0.3") 
     254        SLDC = pars.get('sldc', "11.5") 
     255        ## default parameters from sasmodels 
     256        #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 400,75,35,10,10,10,2,4,2 
     257        ## swap A-B-C to C-B-A 
     258        #A, B, C, DA, DB, DC, SLDA, SLDB, SLDC = C, B, A, DC, DB, DA, SLDC, SLDB, SLDA 
     259        #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 10,20,30,100,200,300,1,2,3 
     260        #SLD_SOLVENT,CONTRAST = 0, 4 
     261        if 1: # C shortest 
     262            B, C = C, B 
     263            DB, DC = DC, DB 
     264            SLDB, SLDC = SLDC, SLDB 
     265        elif 0: # C longest 
     266            A, C = C, A 
     267            DA, DC = DC, DA 
     268            SLDA, SLDC = SLDC, SLDA 
     269        #NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 
     270        NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 
     271        NORM_MP, KERNEL_MP = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC, env=MPenv) 
     272    elif shape.endswith('paracrystal'): 
     273        LATTICE, _ = shape.split('_') 
     274        DNN = pars.get('dnn', 220) 
     275        D_FACTOR = pars.get('d_factor', '0.06') 
     276        RADIUS = pars.get('radius', 40) 
     277        NORM, KERNEL = make_paracrystal( 
     278            radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE) 
     279        NORM_MP, KERNEL_MP = make_paracrystal( 
     280            radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE, env=MPenv) 
     281    else: 
     282        raise ValueError("Unknown shape %r"%shape) 
    257283 
    258284# Note: hardcoded in mp_quad 
     
    273299 
    274300# 2D integration functions 
    275 def mp_quad_2d(q, shape): 
     301def mp_quad_2d(q): 
    276302    evals = [0] 
    277303    def integrand(theta, phi): 
     
    393419    print("simpson-%d"%n, n**2, simps(simps(Zq, dx=dx), dx=dy)*SCALE/(4*pi)) 
    394420    print("romb-%d"%n, n**2, romb(romb(Zq, dx=dx), dx=dy)*SCALE/(4*pi)) 
     421 
     422def quadpy_method(q, rule): 
     423    """ 
     424    Use *rule*="name:index" where name and index are chosen from below. 
     425 
     426    Available rule names and the corresponding indices:: 
     427 
     428        AlbrechtCollatz: [1-5] 
     429        BazantOh: 9, 11, 13 
     430        HeoXu: 13, 15, 17, 19-[1-2], 21-[1-6], 23-[1-3], 25-[1-2], 27-[1-3], 
     431            29, 31, 33, 35, 37, 39-[1-2] 
     432        FliegeMaier: 4, 9, 16, 25 
     433        Lebedev: 3[a-c], 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 35, 
     434            41, 47, 53, 59, 65, 71, 77 83, 89, 95, 101, 107, 113, 119, 125, 131 
     435        McLaren: [1-10] 
     436        Stroud: U3 3-1, U3 5-[1-5], U3 7-[1-2], U3 8-1, U3 9-[1-3], 
     437            U3 11-[1-3], U3 14-1 
     438    """ 
     439    try: 
     440        import quadpy 
     441    except ImportError: 
     442        warnings.warn("use 'pip install quadpy' to enable quadpy.sphere tests") 
     443        return 
     444 
     445    from quadpy.sphere import (AlbrechtCollatz, BazantOh, HeoXu, 
     446        FliegeMaier, Lebedev, McLaren, Stroud, integrate_spherical) 
     447    RULES = { 
     448        'AlbrechtCollatz': AlbrechtCollatz, 
     449        'BazantOh': BazantOh, 
     450        'HeoXu': HeoXu, 
     451        'FliegeMaier': FliegeMaier, 
     452        'Lebedev': Lebedev, 
     453        'McLaren': McLaren, 
     454        'Stroud': Stroud, 
     455    } 
     456    int_index = 'AlbrechtCollatz', 'McLaren' 
     457 
     458    rule_name, rule_index = rule.split(':') 
     459    index = int(rule_index) if rule_name in int_index else rule_index 
     460    rule_obj = RULES[rule_name](index) 
     461    fn = lambda azimuthal, polar: kernel_2d(q=q, theta=polar, phi=azimuthal) 
     462    Iq = integrate_spherical(fn, rule=rule_obj)/(4*pi) 
     463    print("%s degree=%d points=%s => %.15g" 
     464          % (rule, rule_obj.degree, len(rule_obj.points), Iq)) 
    395465 
    396466def plot_2d(q, n=300): 
     
    414484    pylab.show() 
    415485 
    416 def main(Qstr): 
    417     Q = float(Qstr) 
    418     if shape == 'sphere': 
     486def main(): 
     487    import argparse 
     488 
     489    parser = argparse.ArgumentParser( 
     490        description="asymmetric integration explorer", 
     491        formatter_class=argparse.ArgumentDefaultsHelpFormatter, 
     492        ) 
     493    parser.add_argument('-s', '--shape', choices=SHAPES, 
     494                        default='parallelepiped', 
     495                        help='oriented shape') 
     496    parser.add_argument('-q', '--q_value', type=str, default='0.005', 
     497                        help='Q value to evaluate') 
     498    parser.add_argument('pars', type=str, nargs='*', default=[], 
     499                        help='p=val for p in shape parameters') 
     500    opts = parser.parse_args() 
     501    pars = {k: v for par in opts.pars for k, v in [par.split('=')]} 
     502    build_shape(opts.shape, **pars) 
     503 
     504    Q = float(opts.q_value) 
     505    if opts.shape == 'sphere': 
    419506        print("exact", NORM*sp.sas_3j1x_x(Q*RADIUS)**2) 
    420     print("gauss-20", *gauss_quad_2d(Q, n=20)) 
    421     print("gauss-76", *gauss_quad_2d(Q, n=76)) 
    422     print("gauss-150", *gauss_quad_2d(Q, n=150)) 
    423     print("gauss-500", *gauss_quad_2d(Q, n=500)) 
    424     print("gauss-1025", *gauss_quad_2d(Q, n=1025)) 
    425     print("gauss-2049", *gauss_quad_2d(Q, n=2049)) 
    426     print("gauss-20 usub", *gauss_quad_usub(Q, n=20)) 
    427     print("gauss-76 usub", *gauss_quad_usub(Q, n=76)) 
    428     print("gauss-150 usub", *gauss_quad_usub(Q, n=150)) 
     507 
     508    # Methods from quadpy, if quadpy is available 
     509    #  AlbrechtCollatz: [1-5] 
     510    #  BazantOh: 9, 11, 13 
     511    #  HeoXu: 13, 15, 17, 19-[1-2], 21-[1-6], 23-[1-3], 25-[1-2], 27-[1-3], 
     512    #     29, 31, 33, 35, 37, 39-[1-2] 
     513    #  FliegeMaier: 4, 9, 16, 25 
     514    #  Lebedev: 3[a-c], 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 35, 
     515    #     41, 47, 53, 59, 65, 71, 77 83, 89, 95, 101, 107, 113, 119, 125, 131 
     516    #  McLaren: [1-10] 
     517    #  Stroud: U3 3-1, U3 5-[1-5], U3 7-[1-2], U3 8-1, U3 9-[1-3], 
     518    #     U3 11-[1-3], U3 14-1 
     519    quadpy_method(Q, "AlbrechtCollatz:5") 
     520    quadpy_method(Q, "HeoXu:39-2") 
     521    quadpy_method(Q, "FliegeMaier:25") 
     522    quadpy_method(Q, "Lebedev:19") 
     523    quadpy_method(Q, "Lebedev:131") 
     524    quadpy_method(Q, "McLaren:10") 
     525    quadpy_method(Q, "Stroud:U3 14-1") 
     526 
     527    print("gauss-20 points=%d => %.15g" % gauss_quad_2d(Q, n=20)) 
     528    print("gauss-76 points=%d => %.15g" % gauss_quad_2d(Q, n=76)) 
     529    print("gauss-150 points=%d => %.15g" % gauss_quad_2d(Q, n=150)) 
     530    print("gauss-500 points=%d => %.15g" % gauss_quad_2d(Q, n=500)) 
     531    print("gauss-1025 points=%d => %.15g" % gauss_quad_2d(Q, n=1025)) 
     532    print("gauss-2049 points=%d => %.15g" % gauss_quad_2d(Q, n=2049)) 
     533    print("gauss-20 usub points=%d => %.15g" % gauss_quad_usub(Q, n=20)) 
     534    print("gauss-76 usub points=%d => %.15g" % gauss_quad_usub(Q, n=76)) 
     535    print("gauss-150 usub points=%d => %.15g" % gauss_quad_usub(Q, n=150)) 
     536 
    429537    #gridded_2d(Q, n=2**8+1) 
    430538    gridded_2d(Q, n=2**10+1) 
    431539    #gridded_2d(Q, n=2**12+1) 
    432540    #gridded_2d(Q, n=2**15+1) 
    433     if shape not in ('paracrystal', 'core_shell_parallelepiped'): 
    434         # adaptive forms on models for which the calculations are fast enough 
     541    # adaptive forms on models for which the calculations are fast enough 
     542    SLOW_SHAPES = { 
     543        'fcc_paracrystal', 'bcc_paracrystal', 'sc_paracrystal', 
     544        'core_shell_parallelepiped', 
     545    } 
     546    if opts.shape not in SLOW_SHAPES: 
    435547        print("dblquad", *scipy_dblquad_2d(Q)) 
    436548        print("semi-romberg-100", *semi_romberg_2d(Q, n=100)) 
    437549        print("romberg", *scipy_romberg_2d(Q)) 
    438550        with mp.workprec(100): 
    439             print("mpmath", *mp_quad_2d(mp.mpf(Qstr), shape)) 
     551            print("mpmath", *mp_quad_2d(mp.mpf(opts.q_value))) 
    440552    plot_2d(Q, n=200) 
    441553 
    442554if __name__ == "__main__": 
    443     main(Qstr) 
     555    main() 
  • 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/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/convert.py

    r610ef23 r21c93c3  
    105105                    translation[newid+str(k)] = oldid+str(k) 
    106106    # Remove control parameter from the result 
    107     if model_info.control: 
    108         translation[model_info.control] = "CONTROL" 
     107    control_pars = [p.id for p in model_info.parameters.kernel_parameters 
     108                    if p.is_control] 
     109    if control_pars: 
     110        control_id = control_pars[0] 
     111        translation[control_id] = "CONTROL" 
    109112    return translation 
    110113 
  • sasmodels/core.py

    r2dcd6e7 rd92182f  
    2121from . import mixture 
    2222from . import kernelpy 
     23from . import kernelcuda 
    2324from . import kernelcl 
    2425from . import kerneldll 
     
    6364        * all: all models 
    6465        * py: python models only 
    65         * c: compiled models only 
    66         * single: models which support single precision 
    67         * double: models which require double precision 
    68         * opencl: controls if OpenCL is supperessed 
    69         * 1d: models which are 1D only, or 2D using abs(q) 
    70         * 2d: models which can be 2D 
    71         * magnetic: models with an sld 
    72         * nommagnetic: models without an sld 
     66        * c: c models only 
     67        * single: c models which support single precision 
     68        * double: c models which require double precision 
     69        * opencl: c models which run in opencl 
     70        * dll: c models which do not run in opencl 
     71        * 1d: models without orientation 
     72        * 2d: models with orientation 
     73        * magnetic: models supporting magnetic sld 
     74        * nommagnetic: models without magnetic parameter 
    7375 
    7476    For multiple conditions, combine with plus.  For example, *c+single+2d* 
     
    9496    info = load_model_info(name) 
    9597    pars = info.parameters.kernel_parameters 
    96     if kind == "py" and callable(info.Iq): 
    97         return True 
    98     elif kind == "c" and not callable(info.Iq): 
    99         return True 
    100     elif kind == "double" and not info.single: 
    101         return True 
    102     elif kind == "single" and info.single: 
    103         return True 
    104     elif kind == "opencl" and info.opencl: 
    105         return True 
    106     elif kind == "2d" and any(p.type == 'orientation' for p in pars): 
    107         return True 
    108     elif kind == "1d" and all(p.type != 'orientation' for p in pars): 
    109         return True 
    110     elif kind == "magnetic" and any(p.type == 'sld' for p in pars): 
    111         return True 
    112     elif kind == "nonmagnetic" and any(p.type != 'sld' for p in pars): 
    113         return True 
     98    # TODO: may be adding Fq to the list at some point 
     99    is_pure_py = callable(info.Iq) 
     100    if kind == "py": 
     101        return is_pure_py 
     102    elif kind == "c": 
     103        return not is_pure_py 
     104    elif kind == "double": 
     105        return not info.single and not is_pure_py 
     106    elif kind == "single": 
     107        return info.single and not is_pure_py 
     108    elif kind == "opencl": 
     109        return info.opencl 
     110    elif kind == "dll": 
     111        return not info.opencl and not is_pure_py 
     112    elif kind == "2d": 
     113        return any(p.type == 'orientation' for p in pars) 
     114    elif kind == "1d": 
     115        return all(p.type != 'orientation' for p in pars) 
     116    elif kind == "magnetic": 
     117        return any(p.type == 'sld' for p in pars) 
     118    elif kind == "nonmagnetic": 
     119        return not any(p.type == 'sld' for p in pars) 
    114120    return False 
    115121 
     
    205211 
    206212    numpy_dtype, fast, platform = parse_dtype(model_info, dtype, platform) 
    207  
    208213    source = generate.make_source(model_info) 
    209214    if platform == "dll": 
    210215        #print("building dll", numpy_dtype) 
    211216        return kerneldll.load_dll(source['dll'], model_info, numpy_dtype) 
     217    elif platform == "cuda": 
     218        return kernelcuda.GpuModel(source, model_info, numpy_dtype, fast=fast) 
    212219    else: 
    213220        #print("building ocl", numpy_dtype) 
     
    245252    # type: (ModelInfo, str, str) -> (np.dtype, bool, str) 
    246253    """ 
    247     Interpret dtype string, returning np.dtype and fast flag. 
     254    Interpret dtype string, returning np.dtype, fast flag and platform. 
    248255 
    249256    Possible types include 'half', 'single', 'double' and 'quad'.  If the 
     
    253260    default for the model and platform. 
    254261 
    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. 
     262    Platform preference can be specfied ("ocl", "cuda", "dll"), with the 
     263    default being OpenCL or CUDA if available, otherwise DLL.  If the dtype 
     264    name ends with '!' then platform is forced to be DLL rather than GPU. 
     265    The default platform is set by the environment variable SAS_OPENCL, 
     266    SAS_OPENCL=driver:device for OpenCL, SAS_OPENCL=cuda:device for CUDA 
     267    or SAS_OPENCL=none for DLL. 
    258268 
    259269    This routine ignores the preferences within the model definition.  This 
     
    265275    # Assign default platform, overriding ocl with dll if OpenCL is unavailable 
    266276    # If opencl=False OpenCL is switched off 
    267  
    268277    if platform is None: 
    269278        platform = "ocl" 
    270     if not kernelcl.use_opencl() or not model_info.opencl: 
    271         platform = "dll" 
    272279 
    273280    # Check if type indicates dll regardless of which platform is given 
     
    275282        platform = "dll" 
    276283        dtype = dtype[:-1] 
     284 
     285    # Make sure model allows opencl/gpu 
     286    if not model_info.opencl: 
     287        platform = "dll" 
     288 
     289    # Make sure opencl is available, or fallback to cuda then to dll 
     290    if platform == "ocl" and not kernelcl.use_opencl(): 
     291        platform = "cuda" if kernelcuda.use_cuda() else "dll" 
    277292 
    278293    # Convert special type names "half", "fast", and "quad" 
     
    285300        dtype = "float16" 
    286301 
    287     # Convert dtype string to numpy dtype. 
     302    # Convert dtype string to numpy dtype.  Use single precision for GPU 
     303    # if model allows it, otherwise use double precision. 
    288304    if dtype is None or dtype == "default": 
    289         numpy_dtype = (generate.F32 if platform == "ocl" and model_info.single 
     305        numpy_dtype = (generate.F32 if model_info.single and platform in ("ocl", "cuda") 
    290306                       else generate.F64) 
    291307    else: 
    292308        numpy_dtype = np.dtype(dtype) 
    293309 
    294     # Make sure that the type is supported by opencl, otherwise use dll 
     310    # Make sure that the type is supported by GPU, otherwise use dll 
    295311    if platform == "ocl": 
    296312        env = kernelcl.environment() 
    297         if not env.has_type(numpy_dtype): 
    298             platform = "dll" 
    299             if dtype is None: 
    300                 numpy_dtype = generate.F64 
     313    elif platform == "cuda": 
     314        env = kernelcuda.environment() 
     315    else: 
     316        env = None 
     317    if env is not None and not env.has_type(numpy_dtype): 
     318        platform = "dll" 
     319        if dtype is None: 
     320            numpy_dtype = generate.F64 
    301321 
    302322    return numpy_dtype, fast, platform 
    303  
    304 def list_models_main(): 
    305     # type: () -> None 
    306     """ 
    307     Run list_models as a main program.  See :func:`list_models` for the 
    308     kinds of models that can be requested on the command line. 
    309     """ 
    310     import sys 
    311     kind = sys.argv[1] if len(sys.argv) > 1 else "all" 
    312     print("\n".join(list_models(kind))) 
    313323 
    314324def test_composite_order(): 
     
    365375    model = load_model("cylinder@hardsphere*sphere") 
    366376    actual = [p.name for p in model.info.parameters.kernel_parameters] 
    367     target = ("sld sld_solvent radius length theta phi volfraction" 
     377    target = ("sld sld_solvent radius length theta phi" 
     378              " radius_effective volfraction " 
     379              " structure_factor_mode radius_effective_mode" 
    368380              " A_sld A_sld_solvent A_radius").split() 
    369381    assert target == actual, "%s != %s"%(target, actual) 
    370382 
     383def list_models_main(): 
     384    # type: () -> None 
     385    """ 
     386    Run list_models as a main program.  See :func:`list_models` for the 
     387    kinds of models that can be requested on the command line. 
     388    """ 
     389    import sys 
     390    kind = sys.argv[1] if len(sys.argv) > 1 else "all" 
     391    try: 
     392        models = list_models(kind) 
     393    except Exception as exc: 
     394        print(list_models.__doc__) 
     395        return 1 
     396 
     397    print("\n".join(list_models(kind))) 
     398 
    371399if __name__ == "__main__": 
    372400    list_models_main() 
  • 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 rcff2939  
    11#!/usr/bin/env python 
     2# -*- coding: utf-8 -*- 
    23""" 
    34Jitter Explorer 
     
    56 
    67Application to explore orientation angle and angular dispersity. 
     8 
     9From the command line:: 
     10 
     11    # Show docs 
     12    python -m sasmodels.jitter --help 
     13 
     14    # Guyou projection jitter, uniform over 20 degree theta and 10 in phi 
     15    python -m sasmodels.jitter --projection=guyou --dist=uniform --jitter=20,10,0 
     16 
     17From a jupyter cell:: 
     18 
     19    import ipyvolume as ipv 
     20    from sasmodels import jitter 
     21    import importlib; importlib.reload(jitter) 
     22    jitter.set_plotter("ipv") 
     23 
     24    size = (10, 40, 100) 
     25    view = (20, 0, 0) 
     26 
     27    #size = (15, 15, 100) 
     28    #view = (60, 60, 0) 
     29 
     30    dview = (0, 0, 0) 
     31    #dview = (5, 5, 0) 
     32    #dview = (15, 180, 0) 
     33    #dview = (180, 15, 0) 
     34 
     35    projection = 'equirectangular' 
     36    #projection = 'azimuthal_equidistance' 
     37    #projection = 'guyou' 
     38    #projection = 'sinusoidal' 
     39    #projection = 'azimuthal_equal_area' 
     40 
     41    dist = 'uniform' 
     42    #dist = 'gaussian' 
     43 
     44    jitter.run(size=size, view=view, jitter=dview, dist=dist, projection=projection) 
     45    #filename = projection+('_theta' if dview[0] == 180 else '_phi' if dview[1] == 180 else '') 
     46    #ipv.savefig(filename+'.png') 
    747""" 
    848from __future__ import division, print_function 
     
    1050import argparse 
    1151 
    12 try: # CRUFT: travis-ci does not support mpl_toolkits.mplot3d 
    13     import mpl_toolkits.mplot3d  # Adds projection='3d' option to subplot 
    14 except ImportError: 
    15     pass 
    16  
    17 import matplotlib.pyplot as plt 
    18 from matplotlib.widgets import Slider 
    1952import numpy as np 
    2053from numpy import pi, cos, sin, sqrt, exp, degrees, radians 
    2154 
    22 def draw_beam(axes, view=(0, 0)): 
     55def draw_beam(axes, view=(0, 0), alpha=0.5, steps=25): 
    2356    """ 
    2457    Draw the beam going from source at (0, 0, 1) to detector at (0, 0, -1) 
    2558    """ 
    2659    #axes.plot([0,0],[0,0],[1,-1]) 
    27     #axes.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8) 
    28  
    29     steps = 25 
     60    #axes.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=alpha) 
     61 
    3062    u = np.linspace(0, 2 * np.pi, steps) 
    31     v = np.linspace(-1, 1, steps) 
     63    v = np.linspace(-1, 1, 2) 
    3264 
    3365    r = 0.02 
     
    4173    points = Rz(phi)*Ry(theta)*points 
    4274    x, y, z = [v.reshape(shape) for v in points] 
    43  
    44     axes.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5) 
     75    axes.plot_surface(x, y, z, color='yellow', alpha=alpha) 
     76 
     77    # TODO: draw endcaps on beam 
     78    ## Drawing tiny balls on the end will work 
     79    #draw_sphere(axes, radius=0.02, center=(0, 0, 1.3), color='yellow', alpha=alpha) 
     80    #draw_sphere(axes, radius=0.02, center=(0, 0, -1.3), color='yellow', alpha=alpha) 
     81    ## The following does not work 
     82    #triangles = [(0, i+1, i+2) for i in range(steps-2)] 
     83    #x_cap, y_cap = x[:, 0], y[:, 0] 
     84    #for z_cap in z[:, 0], z[:, -1]: 
     85    #    axes.plot_trisurf(x_cap, y_cap, z_cap, triangles, 
     86    #                      color='yellow', alpha=alpha) 
     87 
    4588 
    4689def draw_ellipsoid(axes, size, view, jitter, steps=25, alpha=1): 
     
    5497    x, y, z = transform_xyz(view, jitter, x, y, z) 
    5598 
    56     axes.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha) 
     99    axes.plot_surface(x, y, z, color='w', alpha=alpha) 
    57100 
    58101    draw_labels(axes, view, jitter, [ 
     
    123166    return atoms 
    124167 
    125 def draw_parallelepiped(axes, size, view, jitter, steps=None, alpha=1): 
     168def draw_box(axes, size, view): 
     169    a, b, c = size 
     170    x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1]) 
     171    y = b*np.array([+1, +1, -1, -1, +1, +1, -1, -1]) 
     172    z = c*np.array([+1, +1, +1, +1, -1, -1, -1, -1]) 
     173    x, y, z = transform_xyz(view, None, x, y, z) 
     174    def draw(i, j): 
     175        axes.plot([x[i],x[j]], [y[i], y[j]], [z[i], z[j]], color='black') 
     176    draw(0, 1) 
     177    draw(0, 2) 
     178    draw(0, 3) 
     179    draw(7, 4) 
     180    draw(7, 5) 
     181    draw(7, 6) 
     182 
     183def draw_parallelepiped(axes, size, view, jitter, steps=None, 
     184                        color=(0.6, 1.0, 0.6), alpha=1): 
    126185    """Draw a parallelepiped.""" 
    127186    a, b, c = size 
     
    141200 
    142201    x, y, z = transform_xyz(view, jitter, x, y, z) 
    143     axes.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha) 
    144  
    145     # Draw pink face on box. 
     202    axes.plot_trisurf(x, y, triangles=tri, Z=z, color=color, alpha=alpha, 
     203                      linewidth=0) 
     204 
     205    # Colour the c+ face of the box. 
    146206    # Since I can't control face color, instead draw a thin box situated just 
    147207    # in front of the "c+" face.  Use the c face so that rotations about psi 
    148208    # rotate that face. 
    149     if 1: 
     209    if 0: 
     210        color = (1, 0.6, 0.6)  # pink 
    150211        x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1]) 
    151212        y = b*np.array([+1, +1, -1, -1, +1, +1, -1, -1]) 
    152213        z = c*np.array([+1, +1, +1, +1, -1, -1, -1, -1]) 
    153214        x, y, z = transform_xyz(view, jitter, x, y, abs(z)+0.001) 
    154         axes.plot_trisurf(x, y, triangles=tri, Z=z, color=[1, 0.6, 0.6], alpha=alpha) 
     215        axes.plot_trisurf(x, y, triangles=tri, Z=z, color=color, alpha=alpha) 
    155216 
    156217    draw_labels(axes, view, jitter, [ 
     
    163224    ]) 
    164225 
    165 def draw_sphere(axes, radius=10., steps=100): 
     226def draw_sphere(axes, radius=0.5, steps=25, center=(0,0,0), color='w', alpha=1.): 
    166227    """Draw a sphere""" 
    167228    u = np.linspace(0, 2 * np.pi, steps) 
    168229    v = np.linspace(0, np.pi, steps) 
    169230 
    170     x = radius * np.outer(np.cos(u), np.sin(v)) 
    171     y = radius * np.outer(np.sin(u), np.sin(v)) 
    172     z = radius * np.outer(np.ones(np.size(u)), np.cos(v)) 
    173     axes.plot_surface(x, y, z, rstride=4, cstride=4, color='w') 
    174  
    175 def draw_jitter(axes, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0), 
    176                 draw_shape=draw_parallelepiped): 
     231    x = radius * np.outer(np.cos(u), np.sin(v)) + center[0] 
     232    y = radius * np.outer(np.sin(u), np.sin(v)) + center[1] 
     233    z = radius * np.outer(np.ones(np.size(u)), np.cos(v)) + center[2] 
     234    axes.plot_surface(x, y, z, color=color, alpha=alpha) 
     235    #axes.plot_wireframe(x, y, z) 
     236 
     237def draw_axes(axes, origin=(-1, -1, -1), length=(2, 2, 2)): 
     238    x, y, z = origin 
     239    dx, dy, dz = length 
     240    axes.plot([x, x+dx], [y, y], [z, z], color='black') 
     241    axes.plot([x, x], [y, y+dy], [z, z], color='black') 
     242    axes.plot([x, x], [y, y], [z, z+dz], color='black') 
     243 
     244def draw_person_on_sphere(axes, view, height=0.5, radius=0.5): 
     245    limb_offset = height * 0.05 
     246    head_radius = height * 0.10 
     247    head_height = height - head_radius 
     248    neck_length = head_radius * 0.50 
     249    shoulder_height = height - 2*head_radius - neck_length 
     250    torso_length = shoulder_height * 0.55 
     251    torso_radius = torso_length * 0.30 
     252    leg_length = shoulder_height - torso_length 
     253    arm_length = torso_length * 0.90 
     254 
     255    def _draw_part(x, z): 
     256        y = np.zeros_like(x) 
     257        xp, yp, zp = transform_xyz(view, None, x, y, z + radius) 
     258        axes.plot(xp, yp, zp, color='k') 
     259 
     260    # circle for head 
     261    u = np.linspace(0, 2 * np.pi, 40) 
     262    x = head_radius * np.cos(u) 
     263    z = head_radius * np.sin(u) + head_height 
     264    _draw_part(x, z) 
     265 
     266    # rectangle for body 
     267    x = np.array([-torso_radius, torso_radius, torso_radius, -torso_radius, -torso_radius]) 
     268    z = np.array([0., 0, torso_length, torso_length, 0]) + leg_length 
     269    _draw_part(x, z) 
     270 
     271    # arms 
     272    x = np.array([-torso_radius - limb_offset, -torso_radius - limb_offset, -torso_radius]) 
     273    z = np.array([shoulder_height - arm_length, shoulder_height, shoulder_height]) 
     274    _draw_part(x, z) 
     275    _draw_part(-x, z) 
     276 
     277    # legs 
     278    x = np.array([-torso_radius + limb_offset, -torso_radius + limb_offset]) 
     279    z = np.array([0, leg_length]) 
     280    _draw_part(x, z) 
     281    _draw_part(-x, z) 
     282 
     283    limits = [-radius-height, radius+height] 
     284    axes.set_xlim(limits) 
     285    axes.set_ylim(limits) 
     286    axes.set_zlim(limits) 
     287    axes.set_axis_off() 
     288 
     289def draw_jitter(axes, view, jitter, dist='gaussian', 
     290                size=(0.1, 0.4, 1.0), 
     291                draw_shape=draw_parallelepiped, 
     292                projection='equirectangular', 
     293                alpha=0.8, 
     294                views=None): 
    177295    """ 
    178296    Represent jitter as a set of shapes at different orientations. 
    179297    """ 
     298    project, project_weight = get_projection(projection) 
     299 
    180300    # set max diagonal to 0.95 
    181301    scale = 0.95/sqrt(sum(v**2 for v in size)) 
    182302    size = tuple(scale*v for v in size) 
    183303 
    184     #np.random.seed(10) 
    185     #cloud = np.random.randn(10,3) 
    186     cloud = [ 
    187         [-1, -1, -1], 
    188         [-1, -1, +0], 
    189         [-1, -1, +1], 
    190         [-1, +0, -1], 
    191         [-1, +0, +0], 
    192         [-1, +0, +1], 
    193         [-1, +1, -1], 
    194         [-1, +1, +0], 
    195         [-1, +1, +1], 
    196         [+0, -1, -1], 
    197         [+0, -1, +0], 
    198         [+0, -1, +1], 
    199         [+0, +0, -1], 
    200         [+0, +0, +0], 
    201         [+0, +0, +1], 
    202         [+0, +1, -1], 
    203         [+0, +1, +0], 
    204         [+0, +1, +1], 
    205         [+1, -1, -1], 
    206         [+1, -1, +0], 
    207         [+1, -1, +1], 
    208         [+1, +0, -1], 
    209         [+1, +0, +0], 
    210         [+1, +0, +1], 
    211         [+1, +1, -1], 
    212         [+1, +1, +0], 
    213         [+1, +1, +1], 
    214     ] 
    215304    dtheta, dphi, dpsi = jitter 
    216     if dtheta == 0: 
    217         cloud = [v for v in cloud if v[0] == 0] 
    218     if dphi == 0: 
    219         cloud = [v for v in cloud if v[1] == 0] 
    220     if dpsi == 0: 
    221         cloud = [v for v in cloud if v[2] == 0] 
    222     draw_shape(axes, size, view, [0, 0, 0], steps=100, alpha=0.8) 
    223     scale = {'gaussian':1, 'rectangle':1/sqrt(3), 'uniform':1/3}[dist] 
    224     for point in cloud: 
    225         delta = [scale*dtheta*point[0], scale*dphi*point[1], scale*dpsi*point[2]] 
    226         draw_shape(axes, size, view, delta, alpha=0.8) 
     305    base = {'gaussian':3, 'rectangle':sqrt(3), 'uniform':1}[dist] 
     306    def steps(delta): 
     307        if views is None: 
     308            n = max(3, min(25, 2*int(base*delta/5))) 
     309        else: 
     310            n = views 
     311        return base*delta*np.linspace(-1, 1, n) if delta > 0 else [0.] 
     312    for theta in steps(dtheta): 
     313        for phi in steps(dphi): 
     314            for psi in steps(dpsi): 
     315                w = project_weight(theta, phi, 1.0, 1.0) 
     316                if w > 0: 
     317                    dview = project(theta, phi, psi) 
     318                    draw_shape(axes, size, view, dview, alpha=alpha) 
    227319    for v in 'xyz': 
    228320        a, b, c = size 
    229321        lim = np.sqrt(a**2 + b**2 + c**2) 
    230322        getattr(axes, 'set_'+v+'lim')([-lim, lim]) 
    231         getattr(axes, v+'axis').label.set_text(v) 
     323        #getattr(axes, v+'axis').label.set_text(v) 
    232324 
    233325PROJECTIONS = [ 
     
    237329    'azimuthal_equal_area', 
    238330] 
    239 def draw_mesh(axes, view, jitter, radius=1.2, n=11, dist='gaussian', 
    240               projection='equirectangular'): 
    241     """ 
    242     Draw the dispersion mesh showing the theta-phi orientations at which 
    243     the model will be evaluated. 
    244  
     331def get_projection(projection): 
     332 
     333    """ 
    245334    jitter projections 
    246335    <https://en.wikipedia.org/wiki/List_of_map_projections> 
     
    298387    # TODO: try Kent distribution instead of a gaussian warped by projection 
    299388 
    300     dist_x = np.linspace(-1, 1, n) 
    301     weights = np.ones_like(dist_x) 
    302     if dist == 'gaussian': 
    303         dist_x *= 3 
    304         weights = exp(-0.5*dist_x**2) 
    305     elif dist == 'rectangle': 
    306         # Note: uses sasmodels ridiculous definition of rectangle width 
    307         dist_x *= sqrt(3) 
    308     elif dist == 'uniform': 
    309         pass 
    310     else: 
    311         raise ValueError("expected dist to be gaussian, rectangle or uniform") 
    312  
    313389    if projection == 'equirectangular':  #define PROJECTION 1 
    314         def _rotate(theta_i, phi_j): 
    315             return Rx(phi_j)*Ry(theta_i) 
     390        def _project(theta_i, phi_j, psi): 
     391            latitude, longitude = theta_i, phi_j 
     392            return latitude, longitude, psi 
     393            #return Rx(phi_j)*Ry(theta_i) 
    316394        def _weight(theta_i, phi_j, w_i, w_j): 
    317395            return w_i*w_j*abs(cos(radians(theta_i))) 
    318396    elif projection == 'sinusoidal':  #define PROJECTION 2 
    319         def _rotate(theta_i, phi_j): 
     397        def _project(theta_i, phi_j, psi): 
    320398            latitude = theta_i 
    321399            scale = cos(radians(latitude)) 
    322400            longitude = phi_j/scale if abs(phi_j) < abs(scale)*180 else 0 
    323401            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
    324             return Rx(longitude)*Ry(latitude) 
    325         def _weight(theta_i, phi_j, w_i, w_j): 
     402            return latitude, longitude, psi 
     403            #return Rx(longitude)*Ry(latitude) 
     404        def _project(theta_i, phi_j, w_i, w_j): 
    326405            latitude = theta_i 
    327406            scale = cos(radians(latitude)) 
     
    329408            return active*w_i*w_j 
    330409    elif projection == 'guyou':  #define PROJECTION 3  (eventually?) 
    331         def _rotate(theta_i, phi_j): 
     410        def _project(theta_i, phi_j, psi): 
    332411            from .guyou import guyou_invert 
    333412            #latitude, longitude = guyou_invert([theta_i], [phi_j]) 
    334413            longitude, latitude = guyou_invert([phi_j], [theta_i]) 
    335             return Rx(longitude[0])*Ry(latitude[0]) 
     414            return latitude, longitude, psi 
     415            #return Rx(longitude[0])*Ry(latitude[0]) 
    336416        def _weight(theta_i, phi_j, w_i, w_j): 
    337417            return w_i*w_j 
    338     elif projection == 'azimuthal_equidistance':  # Note: Rz Ry, not Rx Ry 
    339         def _rotate(theta_i, phi_j): 
     418    elif projection == 'azimuthal_equidistance': 
     419        # Note that calculates angles for Rz Ry rather than Rx Ry 
     420        def _project(theta_i, phi_j, psi): 
    340421            latitude = sqrt(theta_i**2 + phi_j**2) 
    341422            longitude = degrees(np.arctan2(phi_j, theta_i)) 
    342423            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
    343             return Rz(longitude)*Ry(latitude) 
     424            return latitude, longitude, psi-longitude, 'zyz' 
     425            #R = Rz(longitude)*Ry(latitude)*Rz(psi) 
     426            #return R_to_xyz(R) 
     427            #return Rz(longitude)*Ry(latitude) 
    344428        def _weight(theta_i, phi_j, w_i, w_j): 
    345429            # Weighting for each point comes from the integral: 
     
    375459            return weight*w_i*w_j if latitude < 180 else 0 
    376460    elif projection == 'azimuthal_equal_area': 
    377         def _rotate(theta_i, phi_j): 
     461        # Note that calculates angles for Rz Ry rather than Rx Ry 
     462        def _project(theta_i, phi_j, psi): 
    378463            radius = min(1, sqrt(theta_i**2 + phi_j**2)/180) 
    379464            latitude = 180-degrees(2*np.arccos(radius)) 
    380465            longitude = degrees(np.arctan2(phi_j, theta_i)) 
    381466            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
    382             return Rz(longitude)*Ry(latitude) 
     467            return latitude, longitude, psi, "zyz" 
     468            #R = Rz(longitude)*Ry(latitude)*Rz(psi) 
     469            #return R_to_xyz(R) 
     470            #return Rz(longitude)*Ry(latitude) 
    383471        def _weight(theta_i, phi_j, w_i, w_j): 
    384472            latitude = sqrt(theta_i**2 + phi_j**2) 
     
    388476        raise ValueError("unknown projection %r"%projection) 
    389477 
     478    return _project, _weight 
     479 
     480def R_to_xyz(R): 
     481    """ 
     482    Return phi, theta, psi Tait-Bryan angles corresponding to the given rotation matrix. 
     483 
     484    Extracting Euler Angles from a Rotation Matrix 
     485    Mike Day, Insomniac Games 
     486    https://d3cw3dd2w32x2b.cloudfront.net/wp-content/uploads/2012/07/euler-angles1.pdf 
     487    Based on: Shoemake’s "Euler Angle Conversion", Graphics Gems IV, pp.  222-229 
     488    """ 
     489    phi = np.arctan2(R[1, 2], R[2, 2]) 
     490    theta = np.arctan2(-R[0, 2], np.sqrt(R[0, 0]**2 + R[0, 1]**2)) 
     491    psi = np.arctan2(R[0, 1], R[0, 0]) 
     492    return np.degrees(phi), np.degrees(theta), np.degrees(psi) 
     493 
     494def draw_mesh(axes, view, jitter, radius=1.2, n=11, dist='gaussian', 
     495              projection='equirectangular'): 
     496    """ 
     497    Draw the dispersion mesh showing the theta-phi orientations at which 
     498    the model will be evaluated. 
     499    """ 
     500 
     501    _project, _weight = get_projection(projection) 
     502    def _rotate(theta, phi, z): 
     503        dview = _project(theta, phi, 0.) 
     504        if len(dview) == 4: # hack for zyz coords 
     505            return Rz(dview[1])*Ry(dview[0])*z 
     506        else: 
     507            return Rx(dview[1])*Ry(dview[0])*z 
     508 
     509 
     510    dist_x = np.linspace(-1, 1, n) 
     511    weights = np.ones_like(dist_x) 
     512    if dist == 'gaussian': 
     513        dist_x *= 3 
     514        weights = exp(-0.5*dist_x**2) 
     515    elif dist == 'rectangle': 
     516        # Note: uses sasmodels ridiculous definition of rectangle width 
     517        dist_x *= sqrt(3) 
     518    elif dist == 'uniform': 
     519        pass 
     520    else: 
     521        raise ValueError("expected dist to be gaussian, rectangle or uniform") 
     522 
    390523    # mesh in theta, phi formed by rotating z 
    391524    dtheta, dphi, dpsi = jitter 
    392525    z = np.matrix([[0], [0], [radius]]) 
    393     points = np.hstack([_rotate(theta_i, phi_j)*z 
     526    points = np.hstack([_rotate(theta_i, phi_j, z) 
    394527                        for theta_i in dtheta*dist_x 
    395528                        for phi_j in dphi*dist_x]) 
     
    469602    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple. 
    470603    """ 
    471     dtheta, dphi, dpsi = jitter 
    472     points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points 
     604    if jitter is None: 
     605        return points 
     606    # Hack to deal with the fact that azimuthal_equidistance uses euler angles 
     607    if len(jitter) == 4: 
     608        dtheta, dphi, dpsi, _ = jitter 
     609        points = Rz(dphi)*Ry(dtheta)*Rz(dpsi)*points 
     610    else: 
     611        dtheta, dphi, dpsi = jitter 
     612        points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points 
    473613    return points 
    474614 
     
    480620    """ 
    481621    theta, phi, psi = view 
    482     points = Rz(phi)*Ry(theta)*Rz(psi)*points 
     622    points = Rz(phi)*Ry(theta)*Rz(psi)*points # viewing angle 
     623    #points = Rz(psi)*Ry(pi/2-theta)*Rz(phi)*points # 1-D integration angles 
     624    #points = Rx(phi)*Ry(theta)*Rz(psi)*points  # angular dispersion angle 
    483625    return points 
     626 
     627def orient_relative_to_beam_quaternion(view, points): 
     628    """ 
     629    Apply the view transform to a set of points. 
     630 
     631    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple. 
     632 
     633    This variant uses quaternions rather than rotation matrices for the 
     634    computation.  It works but it is not used because it doesn't solve 
     635    any problems.  The challenge of mapping theta/phi/psi to SO(3) does 
     636    not disappear by calculating the transform differently. 
     637    """ 
     638    theta, phi, psi = view 
     639    x, y, z = [1, 0, 0], [0, 1, 0], [0, 0, 1] 
     640    q = Quaternion(1, [0, 0, 0]) 
     641    ## Compose a rotation about the three axes by rotating 
     642    ## the unit vectors before applying the rotation. 
     643    #q = Quaternion.from_angle_axis(theta, q.rot(x)) * q 
     644    #q = Quaternion.from_angle_axis(phi, q.rot(y)) * q 
     645    #q = Quaternion.from_angle_axis(psi, q.rot(z)) * q 
     646    ## The above turns out to be equivalent to reversing 
     647    ## the order of application, so ignore it and use below. 
     648    q = q * Quaternion.from_angle_axis(theta, x) 
     649    q = q * Quaternion.from_angle_axis(phi, y) 
     650    q = q * Quaternion.from_angle_axis(psi, z) 
     651    ## Reverse the order by post-multiply rather than pre-multiply 
     652    #q = Quaternion.from_angle_axis(theta, x) * q 
     653    #q = Quaternion.from_angle_axis(phi, y) * q 
     654    #q = Quaternion.from_angle_axis(psi, z) * q 
     655    #print("axes psi", q.rot(np.matrix([x, y, z]).T)) 
     656    return q.rot(points) 
     657#orient_relative_to_beam = orient_relative_to_beam_quaternion 
     658 
     659# Simple stand-alone quaternion class 
     660import numpy as np 
     661from copy import copy 
     662class Quaternion(object): 
     663    def __init__(self, w, r): 
     664         self.w = w 
     665         self.r = np.asarray(r, dtype='d') 
     666    @staticmethod 
     667    def from_angle_axis(theta, r): 
     668         theta = np.radians(theta)/2 
     669         r = np.asarray(r) 
     670         w = np.cos(theta) 
     671         r = np.sin(theta)*r/np.dot(r,r) 
     672         return Quaternion(w, r) 
     673    def __mul__(self, other): 
     674        if isinstance(other, Quaternion): 
     675            w = self.w*other.w - np.dot(self.r, other.r) 
     676            r = self.w*other.r + other.w*self.r + np.cross(self.r, other.r) 
     677            return Quaternion(w, r) 
     678    def rot(self, v): 
     679        v = np.asarray(v).T 
     680        use_transpose = (v.shape[-1] != 3) 
     681        if use_transpose: v = v.T 
     682        v = v + np.cross(2*self.r, np.cross(self.r, v) + self.w*v) 
     683        #v = v + 2*self.w*np.cross(self.r, v) + np.cross(2*self.r, np.cross(self.r, v)) 
     684        if use_transpose: v = v.T 
     685        return v.T 
     686    def conj(self): 
     687        return Quaternion(self.w, -self.r) 
     688    def inv(self): 
     689        return self.conj()/self.norm()**2 
     690    def norm(self): 
     691        return np.sqrt(self.w**2 + np.sum(self.r**2)) 
     692    def __str__(self): 
     693        return "%g%+gi%+gj%+gk"%(self.w, self.r[0], self.r[1], self.r[2]) 
     694def test_qrot(): 
     695    # Define rotation of 60 degrees around an axis in y-z that is 60 degrees from y. 
     696    # The rotation axis is determined by rotating the point [0, 1, 0] about x. 
     697    ax = Quaternion.from_angle_axis(60, [1, 0, 0]).rot([0, 1, 0]) 
     698    q = Quaternion.from_angle_axis(60, ax) 
     699    # Set the point to be rotated, and its expected rotated position. 
     700    p = [1, -1, 2] 
     701    target = [(10+4*np.sqrt(3))/8, (1+2*np.sqrt(3))/8, (14-3*np.sqrt(3))/8] 
     702    #print(q, q.rot(p) - target) 
     703    assert max(abs(q.rot(p) - target)) < 1e-14 
     704#test_qrot() 
     705#import sys; sys.exit() 
    484706 
    485707# translate between number of dimension of dispersity and the number of 
     
    555777        vmin = vmax*10**-7 
    556778        #vmin, vmax = clipped_range(Iqxy, portion=portion, mode='top') 
     779    #vmin, vmax = Iqxy.min(), Iqxy.max() 
    557780    #print("range",(vmin,vmax)) 
    558781    #qx, qy = np.meshgrid(qx, qy) 
    559782    if 0: 
     783        from matplotlib import cm 
    560784        level = np.asarray(255*(Iqxy - vmin)/(vmax - vmin), 'i') 
    561785        level[level < 0] = 0 
    562786        colors = plt.get_cmap()(level) 
    563         axes.plot_surface(qx, qy, -1.1, rstride=1, cstride=1, facecolors=colors) 
     787        #colors = cm.coolwarm(level) 
     788        #colors = cm.gist_yarg(level) 
     789        #colors = cm.Wistia(level) 
     790        colors[level<=0, 3] = 0.  # set floor to transparent 
     791        x, y = np.meshgrid(qx/qx.max(), qy/qy.max()) 
     792        axes.plot_surface(x, y, -1.1*np.ones_like(x), facecolors=colors) 
    564793    elif 1: 
    565794        axes.contourf(qx/qx.max(), qy/qy.max(), Iqxy, zdir='z', offset=-1.1, 
     
    691920} 
    692921 
     922 
    693923def run(model_name='parallelepiped', size=(10, 40, 100), 
     924        view=(0, 0, 0), jitter=(0, 0, 0), 
    694925        dist='gaussian', mesh=30, 
    695926        projection='equirectangular'): 
     
    701932 
    702933    *size* gives the dimensions (a, b, c) of the shape. 
     934 
     935    *view* gives the initial view (theta, phi, psi) of the shape. 
     936 
     937    *view* gives the initial jitter (dtheta, dphi, dpsi) of the shape. 
    703938 
    704939    *dist* is the type of dispersition: gaussian, rectangle, or uniform. 
     
    720955    calculator, size = select_calculator(model_name, n=150, size=size) 
    721956    draw_shape = DRAW_SHAPES.get(model_name, draw_parallelepiped) 
     957    #draw_shape = draw_fcc 
    722958 
    723959    ## uncomment to set an independent the colour range for every view 
     
    725961    calculator.limits = None 
    726962 
    727     ## initial view 
    728     #theta, dtheta = 70., 10. 
    729     #phi, dphi = -45., 3. 
    730     #psi, dpsi = -45., 3. 
    731     theta, phi, psi = 0, 0, 0 
    732     dtheta, dphi, dpsi = 0, 0, 0 
     963    PLOT_ENGINE(calculator, draw_shape, size, view, jitter, dist, mesh, projection) 
     964 
     965def mpl_plot(calculator, draw_shape, size, view, jitter, dist, mesh, projection): 
     966    # Note: travis-ci does not support mpl_toolkits.mplot3d, but this shouldn't be 
     967    # an issue since we are lazy-loading the package on a path that isn't tested. 
     968    import mpl_toolkits.mplot3d  # Adds projection='3d' option to subplot 
     969    import matplotlib as mpl 
     970    import matplotlib.pyplot as plt 
     971    from matplotlib.widgets import Slider 
    733972 
    734973    ## create the plot window 
     
    746985        pass 
    747986 
    748     axcolor = 'lightgoldenrodyellow' 
     987    # CRUFT: use axisbg instead of facecolor for matplotlib<2 
     988    facecolor_prop = 'facecolor' if mpl.__version__ > '2' else 'axisbg' 
     989    props = {facecolor_prop: 'lightgoldenrodyellow'} 
    749990 
    750991    ## 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     stheta = Slider(axes_theta, 'Theta', -90, 90, valinit=theta) 
    755     sphi = Slider(axes_phi, 'Phi', -180, 180, valinit=phi) 
    756     spsi = Slider(axes_psi, 'Psi', -180, 180, valinit=psi) 
    757  
    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) 
     992    axes_theta = plt.axes([0.05, 0.15, 0.50, 0.04], **props) 
     993    axes_phi = plt.axes([0.05, 0.10, 0.50, 0.04], **props) 
     994    axes_psi = plt.axes([0.05, 0.05, 0.50, 0.04], **props) 
     995    stheta = Slider(axes_theta, u'Ξ', -90, 90, valinit=0) 
     996    sphi = Slider(axes_phi, u'φ', -180, 180, valinit=0) 
     997    spsi = Slider(axes_psi, u'ψ', -180, 180, valinit=0) 
     998 
     999    axes_dtheta = plt.axes([0.70, 0.15, 0.20, 0.04], **props) 
     1000    axes_dphi = plt.axes([0.70, 0.1, 0.20, 0.04], **props) 
     1001    axes_dpsi = plt.axes([0.70, 0.05, 0.20, 0.04], **props) 
     1002 
    7611003    # Note: using ridiculous definition of rectangle distribution, whose width 
    7621004    # in sasmodels is sqrt(3) times the given width.  Divide by sqrt(3) to keep 
    7631005    # the maximum width to 90. 
    7641006    dlimit = DIST_LIMITS[dist] 
    765     sdtheta = Slider(axes_dtheta, 'dTheta', 0, 2*dlimit, valinit=dtheta) 
    766     sdphi = Slider(axes_dphi, 'dPhi', 0, 2*dlimit, valinit=dphi) 
    767     sdpsi = Slider(axes_dpsi, 'dPsi', 0, 2*dlimit, valinit=dpsi) 
    768  
     1007    sdtheta = Slider(axes_dtheta, u'Δξ', 0, 2*dlimit, valinit=0) 
     1008    sdphi = Slider(axes_dphi, u'Δφ', 0, 2*dlimit, valinit=0) 
     1009    sdpsi = Slider(axes_dpsi, u'Δψ', 0, 2*dlimit, valinit=0) 
     1010 
     1011    ## initial view and jitter 
     1012    theta, phi, psi = view 
     1013    stheta.set_val(theta) 
     1014    sphi.set_val(phi) 
     1015    spsi.set_val(psi) 
     1016    dtheta, dphi, dpsi = jitter 
     1017    sdtheta.set_val(dtheta) 
     1018    sdphi.set_val(dphi) 
     1019    sdpsi.set_val(dpsi) 
    7691020 
    7701021    ## callback to draw the new view 
     
    7741025        # set small jitter as 0 if multiple pd dims 
    7751026        dims = sum(v > 0 for v in jitter) 
    776         limit = [0, 0.5, 5][dims] 
     1027        limit = [0, 0.5, 5, 5][dims] 
    7771028        jitter = [0 if v < limit else v for v in jitter] 
    7781029        axes.cla() 
    779         draw_beam(axes, (0, 0)) 
    780         draw_jitter(axes, view, jitter, dist=dist, size=size, draw_shape=draw_shape) 
    781         #draw_jitter(axes, view, (0,0,0)) 
     1030 
     1031        ## Visualize as person on globe 
     1032        #draw_sphere(axes) 
     1033        #draw_person_on_sphere(axes, view) 
     1034 
     1035        ## Move beam instead of shape 
     1036        #draw_beam(axes, -view[:2]) 
     1037        #draw_jitter(axes, (0,0,0), (0,0,0), views=3) 
     1038 
     1039        ## Move shape and draw scattering 
     1040        draw_beam(axes, (0, 0), alpha=1.) 
     1041        draw_jitter(axes, view, jitter, dist=dist, size=size, alpha=1., 
     1042                    draw_shape=draw_shape, projection=projection, views=3) 
    7821043        draw_mesh(axes, view, jitter, dist=dist, n=mesh, projection=projection) 
    7831044        draw_scattering(calculator, axes, view, jitter, dist=dist) 
     1045 
    7841046        plt.gcf().canvas.draw() 
    7851047 
     
    7971059    ## go interactive 
    7981060    plt.show() 
     1061 
     1062 
     1063def map_colors(z, kw): 
     1064    from matplotlib import cm 
     1065 
     1066    cmap = kw.pop('cmap', cm.coolwarm) 
     1067    alpha = kw.pop('alpha', None) 
     1068    vmin = kw.pop('vmin', z.min()) 
     1069    vmax = kw.pop('vmax', z.max()) 
     1070    c = kw.pop('c', None) 
     1071    color = kw.pop('color', c) 
     1072    if color is None: 
     1073        znorm = ((z - vmin) / (vmax - vmin)).clip(0, 1) 
     1074        color = cmap(znorm) 
     1075    elif isinstance(color, np.ndarray) and color.shape == z.shape: 
     1076        color = cmap(color) 
     1077    if alpha is None: 
     1078        if isinstance(color, np.ndarray): 
     1079            color = color[..., :3] 
     1080    else: 
     1081        color[..., 3] = alpha 
     1082    kw['color'] = color 
     1083 
     1084def make_vec(*args): 
     1085    #return [np.asarray(v, 'd').flatten() for v in args] 
     1086    return [np.asarray(v, 'd') for v in args] 
     1087 
     1088def make_image(z, kw): 
     1089    import PIL.Image 
     1090    from matplotlib import cm 
     1091 
     1092    cmap = kw.pop('cmap', cm.coolwarm) 
     1093 
     1094    znorm = (z-z.min())/z.ptp() 
     1095    c = cmap(znorm) 
     1096    c = c[..., :3] 
     1097    rgb = np.asarray(c*255, 'u1') 
     1098    image = PIL.Image.fromarray(rgb, mode='RGB') 
     1099    return image 
     1100 
     1101 
     1102_IPV_MARKERS = { 
     1103    'o': 'sphere', 
     1104} 
     1105_IPV_COLORS = { 
     1106    'w': 'white', 
     1107    'k': 'black', 
     1108    'c': 'cyan', 
     1109    'm': 'magenta', 
     1110    'y': 'yellow', 
     1111    'r': 'red', 
     1112    'g': 'green', 
     1113    'b': 'blue', 
     1114} 
     1115def ipv_fix_color(kw): 
     1116    alpha = kw.pop('alpha', None) 
     1117    color = kw.get('color', None) 
     1118    if isinstance(color, str): 
     1119        color = _IPV_COLORS.get(color, color) 
     1120        kw['color'] = color 
     1121    if alpha is not None: 
     1122        color = kw['color'] 
     1123        #TODO: convert color to [r, g, b, a] if not already 
     1124        if isinstance(color, (tuple, list)): 
     1125            if len(color) == 3: 
     1126                color = (color[0], color[1], color[2], alpha) 
     1127            else: 
     1128                color = (color[0], color[1], color[2], alpha*color[3]) 
     1129            color = np.array(color) 
     1130        if isinstance(color, np.ndarray) and color.shape[-1] == 4: 
     1131            color[..., 3] = alpha 
     1132            kw['color'] = color 
     1133 
     1134def ipv_set_transparency(kw, obj): 
     1135    color = kw.get('color', None) 
     1136    if (isinstance(color, np.ndarray) 
     1137            and color.shape[-1] == 4 
     1138            and (color[..., 3] != 1.0).any()): 
     1139        obj.material.transparent = True 
     1140        obj.material.side = "FrontSide" 
     1141 
     1142def ipv_axes(): 
     1143    import ipyvolume as ipv 
     1144 
     1145    class Plotter: 
     1146        # transparency can be achieved by setting the following: 
     1147        #    mesh.color = [r, g, b, alpha] 
     1148        #    mesh.material.transparent = True 
     1149        #    mesh.material.side = "FrontSide" 
     1150        # smooth(ish) rotation can be achieved by setting: 
     1151        #    slide.continuous_update = True 
     1152        #    figure.animation = 0. 
     1153        #    mesh.material.x = x 
     1154        #    mesh.material.y = y 
     1155        #    mesh.material.z = z 
     1156        # maybe need to synchronize update of x/y/z to avoid shimmy when moving 
     1157        def plot(self, x, y, z, **kw): 
     1158            ipv_fix_color(kw) 
     1159            x, y, z = make_vec(x, y, z) 
     1160            ipv.plot(x, y, z, **kw) 
     1161        def plot_surface(self, x, y, z, **kw): 
     1162            facecolors = kw.pop('facecolors', None) 
     1163            if facecolors is not None: 
     1164                kw['color'] = facecolors 
     1165            ipv_fix_color(kw) 
     1166            x, y, z = make_vec(x, y, z) 
     1167            h = ipv.plot_surface(x, y, z, **kw) 
     1168            ipv_set_transparency(kw, h) 
     1169            #h.material.side = "DoubleSide" 
     1170            return h 
     1171        def plot_trisurf(self, x, y, triangles=None, Z=None, **kw): 
     1172            kw.pop('linewidth', None) 
     1173            ipv_fix_color(kw) 
     1174            x, y, z = make_vec(x, y, Z) 
     1175            if triangles is not None: 
     1176                triangles = np.asarray(triangles) 
     1177            h = ipv.plot_trisurf(x, y, z, triangles=triangles, **kw) 
     1178            ipv_set_transparency(kw, h) 
     1179            return h 
     1180        def scatter(self, x, y, z, **kw): 
     1181            x, y, z = make_vec(x, y, z) 
     1182            map_colors(z, kw) 
     1183            marker = kw.get('marker', None) 
     1184            kw['marker'] = _IPV_MARKERS.get(marker, marker) 
     1185            h = ipv.scatter(x, y, z, **kw) 
     1186            ipv_set_transparency(kw, h) 
     1187            return h 
     1188        def contourf(self, x, y, v, zdir='z', offset=0, levels=None, **kw): 
     1189            # Don't use contour for now (although we might want to later) 
     1190            self.pcolor(x, y, v, zdir='z', offset=offset, **kw) 
     1191        def pcolor(self, x, y, v, zdir='z', offset=0, **kw): 
     1192            x, y, v = make_vec(x, y, v) 
     1193            image = make_image(v, kw) 
     1194            xmin, xmax = x.min(), x.max() 
     1195            ymin, ymax = y.min(), y.max() 
     1196            x = np.array([[xmin, xmax], [xmin, xmax]]) 
     1197            y = np.array([[ymin, ymin], [ymax, ymax]]) 
     1198            z = x*0 + offset 
     1199            u = np.array([[0., 1], [0, 1]]) 
     1200            v = np.array([[0., 0], [1, 1]]) 
     1201            h = ipv.plot_mesh(x, y, z, u=u, v=v, texture=image, wireframe=False) 
     1202            ipv_set_transparency(kw, h) 
     1203            h.material.side = "DoubleSide" 
     1204            return h 
     1205        def text(self, *args, **kw): 
     1206            pass 
     1207        def set_xlim(self, limits): 
     1208            ipv.xlim(*limits) 
     1209        def set_ylim(self, limits): 
     1210            ipv.ylim(*limits) 
     1211        def set_zlim(self, limits): 
     1212            ipv.zlim(*limits) 
     1213        def set_axes_on(self): 
     1214            ipv.style.axis_on() 
     1215        def set_axis_off(self): 
     1216            ipv.style.axes_off() 
     1217    return Plotter() 
     1218 
     1219def ipv_plot(calculator, draw_shape, size, view, jitter, dist, mesh, projection): 
     1220    import ipywidgets as widgets 
     1221    import ipyvolume as ipv 
     1222 
     1223    axes = ipv_axes() 
     1224 
     1225    def draw(view, jitter): 
     1226        camera = ipv.gcf().camera 
     1227        #print(ipv.gcf().__dict__.keys()) 
     1228        #print(dir(ipv.gcf())) 
     1229        ipv.figure(animation=0.)  # no animation when updating object mesh 
     1230 
     1231        # set small jitter as 0 if multiple pd dims 
     1232        dims = sum(v > 0 for v in jitter) 
     1233        limit = [0, 0.5, 5, 5][dims] 
     1234        jitter = [0 if v < limit else v for v in jitter] 
     1235 
     1236        ## Visualize as person on globe 
     1237        #draw_beam(axes, (0, 0)) 
     1238        #draw_sphere(axes) 
     1239        #draw_person_on_sphere(axes, view) 
     1240 
     1241        ## Move beam instead of shape 
     1242        #draw_beam(axes, view=(-view[0], -view[1])) 
     1243        #draw_jitter(axes, view=(0,0,0), jitter=(0,0,0)) 
     1244 
     1245        ## Move shape and draw scattering 
     1246        draw_beam(axes, (0, 0), steps=25) 
     1247        draw_jitter(axes, view, jitter, dist=dist, size=size, alpha=1.0, 
     1248                    draw_shape=draw_shape, projection=projection) 
     1249        draw_mesh(axes, view, jitter, dist=dist, n=mesh, radius=0.95, 
     1250                  projection=projection) 
     1251        draw_scattering(calculator, axes, view, jitter, dist=dist) 
     1252 
     1253        draw_axes(axes, origin=(-1, -1, -1.1)) 
     1254        ipv.style.box_off() 
     1255        ipv.style.axes_off() 
     1256        ipv.xyzlabel(" ", " ", " ") 
     1257 
     1258        ipv.gcf().camera = camera 
     1259        ipv.show() 
     1260 
     1261 
     1262    trange, prange = (-180., 180., 1.), (-180., 180., 1.) 
     1263    dtrange, dprange = (0., 180., 1.), (0., 180., 1.) 
     1264 
     1265    ## Super simple interfaca, but uses non-ascii variable namese 
     1266    # Ξ φ ψ Δξ Δφ Δψ 
     1267    #def update(**kw): 
     1268    #    view = kw['Ξ'], kw['φ'], kw['ψ'] 
     1269    #    jitter = kw['Δξ'], kw['Δφ'], kw['Δψ'] 
     1270    #    draw(view, jitter) 
     1271    #widgets.interact(update, Ξ=trange, φ=prange, ψ=prange, Δξ=dtrange, Δφ=dprange, Δψ=dprange) 
     1272 
     1273    def update(theta, phi, psi, dtheta, dphi, dpsi): 
     1274        draw(view=(theta, phi, psi), jitter=(dtheta, dphi, dpsi)) 
     1275 
     1276    def slider(name, slice, init=0.): 
     1277        return widgets.FloatSlider( 
     1278            value=init, 
     1279            min=slice[0], 
     1280            max=slice[1], 
     1281            step=slice[2], 
     1282            description=name, 
     1283            disabled=False, 
     1284            #continuous_update=True, 
     1285            continuous_update=False, 
     1286            orientation='horizontal', 
     1287            readout=True, 
     1288            readout_format='.1f', 
     1289            ) 
     1290    theta = slider(u'Ξ', trange, view[0]) 
     1291    phi = slider(u'φ', prange, view[1]) 
     1292    psi = slider(u'ψ', prange, view[2]) 
     1293    dtheta = slider(u'Δξ', dtrange, jitter[0]) 
     1294    dphi = slider(u'Δφ', dprange, jitter[1]) 
     1295    dpsi = slider(u'Δψ', dprange, jitter[2]) 
     1296    fields = { 
     1297        'theta': theta, 'phi': phi, 'psi': psi, 
     1298        'dtheta': dtheta, 'dphi': dphi, 'dpsi': dpsi, 
     1299    } 
     1300    ui = widgets.HBox([ 
     1301        widgets.VBox([theta, phi, psi]), 
     1302        widgets.VBox([dtheta, dphi, dpsi]) 
     1303    ]) 
     1304 
     1305    out = widgets.interactive_output(update, fields) 
     1306    display(ui, out) 
     1307 
     1308 
     1309_ENGINES = { 
     1310    "matplotlib": mpl_plot, 
     1311    "mpl": mpl_plot, 
     1312    #"plotly": plotly_plot, 
     1313    "ipvolume": ipv_plot, 
     1314    "ipv": ipv_plot, 
     1315} 
     1316PLOT_ENGINE = _ENGINES["matplotlib"] 
     1317def set_plotter(name): 
     1318    global PLOT_ENGINE 
     1319    PLOT_ENGINE = _ENGINES[name] 
    7991320 
    8001321def main(): 
     
    8081329    parser.add_argument('-s', '--size', type=str, default='10,40,100', 
    8091330                        help='a,b,c lengths') 
     1331    parser.add_argument('-v', '--view', type=str, default='0,0,0', 
     1332                        help='initial view angles') 
     1333    parser.add_argument('-j', '--jitter', type=str, default='0,0,0', 
     1334                        help='initial angular dispersion') 
    8101335    parser.add_argument('-d', '--distribution', choices=DISTRIBUTIONS, 
    8111336                        default=DISTRIBUTIONS[0], 
     
    8161341                        help='oriented shape') 
    8171342    opts = parser.parse_args() 
    818     size = tuple(int(v) for v in opts.size.split(',')) 
    819     run(opts.shape, size=size, 
     1343    size = tuple(float(v) for v in opts.size.split(',')) 
     1344    view = tuple(float(v) for v in opts.view.split(',')) 
     1345    jitter = tuple(float(v) for v in opts.jitter.split(',')) 
     1346    run(opts.shape, size=size, view=view, jitter=jitter, 
    8201347        mesh=opts.mesh, dist=opts.distribution, 
    8211348        projection=opts.projection) 
  • sasmodels/kernel.py

    r2d81cfe r36a2418  
    2323# pylint: enable=unused-import 
    2424 
     25 
    2526class KernelModel(object): 
    2627    info = None  # type: ModelInfo 
     
    3435        pass 
    3536 
     37 
    3638class Kernel(object): 
    3739    #: kernel dimension, either "1d" or "2d" 
     
    4143    dtype = None  # type: np.dtype 
    4244 
    43     def __call__(self, call_details, values, cutoff, magnetic): 
     45    def Iq(self, call_details, values, cutoff, magnetic): 
    4446        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    45         raise NotImplementedError("need to implement __call__") 
     47        r""" 
     48        Returns I(q) from the polydisperse average scattering. 
     49 
     50        .. math:: 
     51 
     52            I(q) = \text{scale} \cdot P(q) + \text{background} 
     53 
     54        With the correct choice of model and contrast, setting *scale* to 
     55        the volume fraction $V_f$ of particles should match the measured 
     56        absolute scattering.  Some models (e.g., vesicle) have volume fraction 
     57        built into the model, and do not need an additional scale. 
     58        """ 
     59        _, F2, _, shell_volume, _ = self.Fq(call_details, values, cutoff, magnetic, 
     60                              effective_radius_type=0) 
     61        combined_scale = values[0]/shell_volume 
     62        background = values[1] 
     63        return combined_scale*F2 + background 
     64    __call__ = Iq 
     65 
     66    def Fq(self, call_details, values, cutoff, magnetic, effective_radius_type=0): 
     67        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
     68        r""" 
     69        Returns <F(q)>, <F(q)^2>, effective radius, shell volume and 
     70        form:shell volume ratio.  The <F(q)> term may be None if the 
     71        form factor does not support direct computation of $F(q)$ 
     72 
     73        $P(q) = <F^2(q)>/<V>$ is used for structure factor calculations, 
     74 
     75        .. math:: 
     76 
     77            I(q) = \text{scale} \cdot P(q) \cdot S(q) + \text{background} 
     78 
     79        For the beta approximation, this becomes 
     80 
     81        .. math:: 
     82 
     83            I(q) = \text{scale} * P (1 + <F>^2/<F^2> (S - 1)) + \text{background} 
     84                 = \text{scale}/<V> (<F^2> + <F>^2 (S - 1)) + \text{background} 
     85 
     86        $<F(q)>$ and $<F^2(q)>$ are averaged by polydispersity in shape 
     87        and orientation, with each configuration $x_k$ having form factor 
     88        $F(q, x_k)$, weight $w_k$ and volume $V_k$.  The result is: 
     89 
     90        .. math:: 
     91 
     92            P(q) = \frac{\sum w_k F^2(q, x_k) / \sum w_k}{\sum w_k V_k / \sum w_k} 
     93 
     94        The form factor itself is scaled by volume and contrast to compute the 
     95        total scattering.  This is then squared, and the volume weighted 
     96        F^2 is then normalized by volume F.  For a given density, the number 
     97        of scattering centers is assumed to scale linearly with volume.  Later 
     98        scaling the resulting $P(q)$ by the volume fraction of particles 
     99        gives the total scattering on an absolute scale. Most models 
     100        incorporate the volume fraction into the overall scale parameter.  An 
     101        exception is vesicle, which includes the volume fraction parameter in 
     102        the model itself, scaling $F$ by $\surd V_f$ so that the math for 
     103        the beta approximation works out. 
     104 
     105        By scaling $P(q)$ by total weight $\sum w_k$, there is no need to make 
     106        sure that the polydisperisity distributions normalize to one.  In 
     107        particular, any distibution values $x_k$ outside the valid domain 
     108        of $F$ will not be included, and the distribution will be implicitly 
     109        truncated.  This is controlled by the parameter limits defined in the 
     110        model (which truncate the distribution before calling the kernel) as 
     111        well as any region excluded using the *INVALID* macro defined within 
     112        the model itself. 
     113 
     114        The volume used in the polydispersity calculation is the form volume 
     115        for solid objects or the shell volume for hollow objects.  Shell 
     116        volume should be used within $F$ so that the normalizing scale 
     117        represents the volume fraction of the shell rather than the entire 
     118        form.  This corresponds to the volume fraction of shell-forming 
     119        material added to the solvent. 
     120 
     121        The calculation of $S$ requires the effective radius and the 
     122        volume fraction of the particles.  The model can have several 
     123        different ways to compute effective radius, with the 
     124        *effective_radius_type* parameter used to select amongst them.  The 
     125        volume fraction of particles should be determined from the total 
     126        volume fraction of the form, not just the shell volume fraction. 
     127        This makes a difference for hollow shapes, which need to scale 
     128        the volume fraction by the returned volume ratio when computing $S$. 
     129        For solid objects, the shell volume is set to the form volume so 
     130        this scale factor evaluates to one and so can be used for both 
     131        hollow and solid shapes. 
     132        """ 
     133        self._call_kernel(call_details, values, cutoff, magnetic, effective_radius_type) 
     134        #print("returned",self.q_input.q, self.result) 
     135        nout = 2 if self.info.have_Fq and self.dim == '1d' else 1 
     136        total_weight = self.result[nout*self.q_input.nq + 0] 
     137        # Note: total_weight = sum(weight > cutoff), with cutoff >= 0, so it 
     138        # is okay to test directly against zero.  If weight is zero then I(q), 
     139        # etc. must also be zero. 
     140        if total_weight == 0.: 
     141            total_weight = 1. 
     142        # Note: shell_volume == form_volume for solid objects 
     143        form_volume = self.result[nout*self.q_input.nq + 1]/total_weight 
     144        shell_volume = self.result[nout*self.q_input.nq + 2]/total_weight 
     145        effective_radius = self.result[nout*self.q_input.nq + 3]/total_weight 
     146        if shell_volume == 0.: 
     147            shell_volume = 1. 
     148        F1 = self.result[1:nout*self.q_input.nq:nout]/total_weight if nout == 2 else None 
     149        F2 = self.result[0:nout*self.q_input.nq:nout]/total_weight 
     150        return F1, F2, effective_radius, shell_volume, form_volume/shell_volume 
    46151 
    47152    def release(self): 
    48153        # type: () -> None 
    49154        pass 
     155 
     156    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
     157        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
     158        """ 
     159        Call the kernel.  Subclasses defining kernels for particular execution 
     160        engines need to provide an implementation for this. 
     161        """ 
     162        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 r0d26e91  
    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: 
     
    6575    from pyopencl import mem_flags as mf 
    6676    from pyopencl.characterize import get_fast_inaccurate_build_options 
    67     # Ask OpenCL for the default context so that we know that one exists 
     77    # Ask OpenCL for the default context so that we know that one exists. 
    6878    cl.create_some_context(interactive=False) 
    6979    HAVE_OPENCL = True 
     
    7484 
    7585from . import generate 
     86from .generate import F32, F64 
    7687from .kernel import KernelModel, Kernel 
    7788 
     
    8596# pylint: enable=unused-import 
    8697 
    87 # CRUFT: pyopencl < 2017.1  (as of June 2016 needs quotes around include path) 
     98 
     99# CRUFT: pyopencl < 2017.1 (as of June 2016 needs quotes around include path). 
    88100def quote_path(v): 
    89101    """ 
     
    96108    return '"'+v+'"' if v and ' ' in v and not v[0] in "\"'-" else v 
    97109 
     110 
    98111def fix_pyopencl_include(): 
    99112    """ 
     
    102115    import pyopencl as cl 
    103116    if hasattr(cl, '_DEFAULT_INCLUDE_OPTIONS'): 
    104         cl._DEFAULT_INCLUDE_OPTIONS = [quote_path(v) for v in cl._DEFAULT_INCLUDE_OPTIONS] 
     117        cl._DEFAULT_INCLUDE_OPTIONS = [ 
     118            quote_path(v) for v in cl._DEFAULT_INCLUDE_OPTIONS 
     119            ] 
     120 
    105121 
    106122if HAVE_OPENCL: 
     
    115131MAX_LOOPS = 2048 
    116132 
    117  
    118133# Pragmas for enable OpenCL features.  Be sure to protect them so that they 
    119134# still compile even if OpenCL is not present. 
     
    130145""" 
    131146 
     147 
    132148def use_opencl(): 
    133     return HAVE_OPENCL and os.environ.get("SAS_OPENCL", "").lower() != "none" 
     149    sas_opencl = os.environ.get("SAS_OPENCL", "OpenCL").lower() 
     150    return HAVE_OPENCL and sas_opencl != "none" and not sas_opencl.startswith("cuda") 
     151 
    134152 
    135153ENV = None 
     
    140158    global ENV 
    141159    ENV = GpuEnvironment() if use_opencl() else None 
     160 
    142161 
    143162def environment(): 
     
    157176    return ENV 
    158177 
     178 
    159179def has_type(device, dtype): 
    160180    # type: (cl.Device, np.dtype) -> bool 
     
    162182    Return true if device supports the requested precision. 
    163183    """ 
    164     if dtype == generate.F32: 
     184    if dtype == F32: 
    165185        return True 
    166     elif dtype == generate.F64: 
     186    elif dtype == F64: 
    167187        return "cl_khr_fp64" in device.extensions 
    168     elif dtype == generate.F16: 
    169         return "cl_khr_fp16" in device.extensions 
    170188    else: 
     189        # Not supporting F16 type since it isn't accurate enough. 
    171190        return False 
     191 
    172192 
    173193def get_warp(kernel, queue): 
     
    179199        cl.kernel_work_group_info.PREFERRED_WORK_GROUP_SIZE_MULTIPLE, 
    180200        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) 
    203201 
    204202 
     
    223221        source_list.insert(0, _F64_PRAGMA) 
    224222 
    225     # Note: USE_SINCOS makes the intel cpu slower under opencl 
     223    # Note: USE_SINCOS makes the Intel CPU slower under OpenCL. 
    226224    if context.devices[0].type == cl.device_type.GPU: 
    227225        source_list.insert(0, "#define USE_SINCOS\n") 
     
    230228    source = "\n".join(source_list) 
    231229    program = cl.Program(context, source).build(options=options) 
     230 
    232231    #print("done with "+program) 
    233232    return program 
    234233 
    235234 
    236 # for now, this returns one device in the context 
    237 # TODO: create a context that contains all devices on all platforms 
     235# For now, this returns one device in the context. 
     236# TODO: Create a context that contains all devices on all platforms. 
    238237class GpuEnvironment(object): 
    239238    """ 
    240     GPU context, with possibly many devices, and one queue per device. 
     239    GPU context for OpenCL, with possibly many devices and one queue per device. 
    241240    """ 
    242241    def __init__(self): 
    243242        # type: () -> None 
    244         # 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() 
    256  
    257         # 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] 
     243        # Find gpu context. 
     244        context_list = _create_some_context() 
     245 
     246        # Find a context for F32 and for F64 (maybe the same one). 
     247        # F16 isn't good enough. 
     248        self.context = {} 
     249        for dtype in (F32, F64): 
     250            for context in context_list: 
     251                if has_type(context.devices[0], dtype): 
     252                    self.context[dtype] = context 
     253                    break 
     254            else: 
     255                self.context[dtype] = None 
     256 
     257        # Build a queue for each context. 
     258        self.queue = {} 
     259        context = self.context[F32] 
     260        self.queue[F32] = cl.CommandQueue(context, context.devices[0]) 
     261        if self.context[F64] == self.context[F32]: 
     262            self.queue[F64] = self.queue[F32] 
     263        else: 
     264            context = self.context[F64] 
     265            self.queue[F64] = cl.CommandQueue(context, context.devices[0]) 
     266 
     267        ## Byte boundary for data alignment. 
     268        #self.data_boundary = max(context.devices[0].min_data_type_align_size 
     269        #                         for context in self.context.values()) 
     270 
     271        # Cache for compiled programs, and for items in context. 
    262272        self.compiled = {} 
    263273 
     
    267277        Return True if all devices support a given type. 
    268278        """ 
    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") 
     279        return self.context.get(dtype, None) is not None 
    304280 
    305281    def compile_program(self, name, source, dtype, fast, timestamp): 
     
    309285        """ 
    310286        # Note: PyOpenCL caches based on md5 hash of source, options and device 
    311         # so we don't really need to cache things for ourselves.  I'll do so 
    312         # anyway just to save some data munging time. 
     287        # but I'll do so as well just to save some data munging time. 
    313288        tag = generate.tag_source(source) 
    314289        key = "%s-%s-%s%s"%(name, dtype, tag, ("-fast" if fast else "")) 
    315         # Check timestamp on program 
     290        # Check timestamp on program. 
    316291        program, program_timestamp = self.compiled.get(key, (None, np.inf)) 
    317292        if program_timestamp < timestamp: 
    318293            del self.compiled[key] 
    319294        if key not in self.compiled: 
    320             context = self.get_context(dtype) 
     295            context = self.context[dtype] 
    321296            logging.info("building %s for OpenCL %s", key, 
    322297                         context.devices[0].name.strip()) 
    323             program = compile_model(self.get_context(dtype), 
     298            program = compile_model(self.context[dtype], 
    324299                                    str(source), dtype, fast) 
    325300            self.compiled[key] = (program, timestamp) 
    326301        return program 
     302 
     303 
     304def _create_some_context(): 
     305    # type: () -> cl.Context 
     306    """ 
     307    Protected call to cl.create_some_context without interactivity. 
     308 
     309    Uses SAS_OPENCL or PYOPENCL_CTX if they are set in the environment, 
     310    otherwise scans for the most appropriate device using 
     311    :func:`_get_default_context`.  Ignore *SAS_OPENCL=OpenCL*, which 
     312    indicates that an OpenCL device should be used without specifying 
     313    which one (and not a CUDA device, or no GPU). 
     314    """ 
     315    # Assume we do not get here if SAS_OPENCL is None or CUDA. 
     316    sas_opencl = os.environ.get('SAS_OPENCL', 'opencl') 
     317    if sas_opencl.lower() != 'opencl': 
     318        # Setting PYOPENCL_CTX as a SAS_OPENCL to create cl context. 
     319        os.environ["PYOPENCL_CTX"] = sas_opencl 
     320 
     321    if 'PYOPENCL_CTX' in os.environ: 
     322        try: 
     323            return [cl.create_some_context(interactive=False)] 
     324        except Exception as exc: 
     325            warnings.warn(str(exc)) 
     326            warnings.warn("pyopencl.create_some_context() failed.  The " 
     327                "environment variable 'SAS_OPENCL' or 'PYOPENCL_CTX' might " 
     328                "not be set correctly") 
     329 
     330    return _get_default_context() 
     331 
    327332 
    328333def _get_default_context(): 
     
    337342    # is running may increase throughput. 
    338343    # 
    339     # Macbook pro, base install: 
     344    # MacBook Pro, base install: 
    340345    #     {'Apple': [Intel CPU, NVIDIA GPU]} 
    341     # Macbook pro, base install: 
     346    # MacBook Pro, base install: 
    342347    #     {'Apple': [Intel CPU, Intel GPU]} 
    343     # 2 x nvidia 295 with Intel and NVIDIA opencl drivers installed 
     348    # 2 x NVIDIA 295 with Intel and NVIDIA opencl drivers install: 
    344349    #     {'Intel': [CPU], 'NVIDIA': [GPU, GPU, GPU, GPU]} 
    345350    gpu, cpu = None, None 
     
    364369            else: 
    365370                # System has cl.device_type.ACCELERATOR or cl.device_type.CUSTOM 
    366                 # Intel Phi for example registers as an accelerator 
     371                # Intel Phi for example registers as an accelerator. 
    367372                # Since the user installed a custom device on their system 
    368373                # and went through the pain of sorting out OpenCL drivers for 
     
    371376                gpu = device 
    372377 
    373     # order the devices by gpu then by cpu; when searching for an available 
     378    # Order the devices by gpu then by cpu; when searching for an available 
    374379    # device by data type they will be checked in this order, which means 
    375380    # that if the gpu supports double then the cpu will never be used (though 
     
    398403    that the compiler is allowed to take shortcuts. 
    399404    """ 
     405    info = None  # type: ModelInfo 
     406    source = ""  # type: str 
     407    dtype = None  # type: np.dtype 
     408    fast = False  # type: bool 
     409    _program = None  # type: cl.Program 
     410    _kernels = None  # type: Dict[str, cl.Kernel] 
     411 
    400412    def __init__(self, source, model_info, dtype=generate.F32, fast=False): 
    401413        # type: (Dict[str,str], ModelInfo, np.dtype, bool) -> None 
     
    404416        self.dtype = dtype 
    405417        self.fast = fast 
    406         self.program = None # delay program creation 
    407         self._kernels = None 
    408418 
    409419    def __getstate__(self): 
     
    414424        # type: (Tuple[ModelInfo, str, np.dtype, bool]) -> None 
    415425        self.info, self.source, self.dtype, self.fast = state 
    416         self.program = None 
     426        self._program = self._kernels = None 
    417427 
    418428    def make_kernel(self, q_vectors): 
    419429        # 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( 
    424                 self.info.name, 
    425                 self.source['opencl'], 
    426                 self.dtype, 
    427                 self.fast, 
    428                 timestamp) 
    429             variants = ['Iq', 'Iqxy', 'Imagnetic'] 
    430             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']] 
    436         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() 
    451  
    452 # TODO: check that we don't need a destructor for buffers which go out of scope 
     430        return GpuKernel(self, q_vectors) 
     431 
     432    def get_function(self, name): 
     433        # type: (str) -> cl.Kernel 
     434        """ 
     435        Fetch the kernel from the environment by name, compiling it if it 
     436        does not already exist. 
     437        """ 
     438        if self._program is None: 
     439            self._prepare_program() 
     440        return self._kernels[name] 
     441 
     442    def _prepare_program(self): 
     443        # type: (str) -> None 
     444        env = environment() 
     445        timestamp = generate.ocl_timestamp(self.info) 
     446        program = env.compile_program( 
     447            self.info.name, 
     448            self.source['opencl'], 
     449            self.dtype, 
     450            self.fast, 
     451            timestamp) 
     452        variants = ['Iq', 'Iqxy', 'Imagnetic'] 
     453        names = [generate.kernel_name(self.info, k) for k in variants] 
     454        functions = [getattr(program, k) for k in names] 
     455        self._kernels = {k: v for k, v in zip(variants, functions)} 
     456        # Keep a handle to program so GC doesn't collect. 
     457        self._program = program 
     458 
     459 
     460# TODO: Check that we don't need a destructor for buffers which go out of scope. 
    453461class GpuInput(object): 
    454462    """ 
     
    472480    def __init__(self, q_vectors, dtype=generate.F32): 
    473481        # type: (List[np.ndarray], np.dtype) -> None 
    474         # TODO: do we ever need double precision q? 
    475         env = environment() 
     482        # TODO: Do we ever need double precision q? 
    476483        self.nq = q_vectors[0].size 
    477484        self.dtype = np.dtype(dtype) 
    478485        self.is_2d = (len(q_vectors) == 2) 
    479         # TODO: stretch input based on get_warp() 
    480         # not doing it now since warp depends on kernel, which is not known 
     486        # TODO: Stretch input based on get_warp(). 
     487        # Not doing it now since warp depends on kernel, which is not known 
    481488        # at this point, so instead using 32, which is good on the set of 
    482489        # architectures tested so far. 
    483490        if self.is_2d: 
    484             # Note: 16 rather than 15 because result is 1 longer than input. 
    485             width = ((self.nq+16)//16)*16 
     491            width = ((self.nq+15)//16)*16 
    486492            self.q = np.empty((width, 2), dtype=dtype) 
    487493            self.q[:self.nq, 0] = q_vectors[0] 
    488494            self.q[:self.nq, 1] = q_vectors[1] 
    489495        else: 
    490             # Note: 32 rather than 31 because result is 1 longer than input. 
    491             width = ((self.nq+32)//32)*32 
     496            width = ((self.nq+31)//32)*32 
    492497            self.q = np.empty(width, dtype=dtype) 
    493498            self.q[:self.nq] = q_vectors[0] 
    494499        self.global_size = [self.q.shape[0]] 
    495         context = env.get_context(self.dtype) 
    496500        #print("creating inputs of size", self.global_size) 
     501 
     502        # Transfer input value to GPU. 
     503        env = environment() 
     504        context = env.context[self.dtype] 
    497505        self.q_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
    498506                             hostbuf=self.q) 
     
    501509        # type: () -> None 
    502510        """ 
    503         Free the memory. 
     511        Free the buffer associated with the q value. 
    504512        """ 
    505513        if self.q_b is not None: 
     
    511519        self.release() 
    512520 
     521 
    513522class GpuKernel(Kernel): 
    514523    """ 
    515524    Callable SAS kernel. 
    516525 
    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 
    522  
    523     *dtype* is the kernel precision 
    524  
    525     The resulting call method takes the *pars*, a list of values for 
    526     the fixed parameters to the kernel, and *pd_pars*, a list of (value,weight) 
    527     vectors for the polydisperse parameters.  *cutoff* determines the 
    528     integration limits: any points with combined weight less than *cutoff* 
    529     will not be calculated. 
     526    *model* is the GpuModel object to call 
     527 
     528    The kernel is derived from :class:`Kernel`, providing the 
     529    :meth:`call_kernel` method to evaluate the kernel for a given set of 
     530    parameters.  Because of the need to move the q values to the GPU before 
     531    evaluation, the kernel is instantiated for a particular set of q vectors, 
     532    and can be called many times without transfering q each time. 
    530533 
    531534    Call :meth:`release` when done with the kernel instance. 
    532535    """ 
    533     def __init__(self, kernel, dtype, model_info, q_vectors): 
    534         # 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 
     536    #: SAS model information structure. 
     537    info = None  # type: ModelInfo 
     538    #: Kernel precision. 
     539    dtype = None  # type: np.dtype 
     540    #: Kernel dimensions (1d or 2d). 
     541    dim = ""  # type: str 
     542    #: Calculation results, updated after each call to :meth:`_call_kernel`. 
     543    result = None  # type: np.ndarray 
     544 
     545    def __init__(self, model, q_vectors): 
     546        # type: (GpuModel, List[np.ndarray]) -> None 
     547        dtype = model.dtype 
     548        self.q_input = GpuInput(q_vectors, dtype) 
     549        self._model = model 
     550 
     551        # Attributes accessed from the outside. 
     552        self.dim = '2d' if self.q_input.is_2d else '1d' 
     553        self.info = model.info 
    538554        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 
     555 
     556        # Converter to translate input to target type. 
     557        self._as_dtype = np.float64 if dtype == generate.F64 else np.float32 
     558 
     559        # Holding place for the returned value. 
     560        nout = 2 if self.info.have_Fq and self.dim == '1d' else 1 
     561        extra_q = 4  # Total weight, form volume, shell volume and R_eff. 
     562        self.result = np.empty(self.q_input.nq*nout + extra_q, dtype) 
     563 
     564        # Allocate result value on GPU. 
    545565        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): 
    559         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    560         context = self.queue.context 
    561         # Arrange data transfer to card 
     566        context = env.context[self.dtype] 
     567        width = ((self.result.size+31)//32)*32 * self.dtype.itemsize 
     568        self._result_b = cl.Buffer(context, mf.READ_WRITE, width) 
     569 
     570    def _call_kernel(self, call_details, values, cutoff, magnetic, 
     571                     effective_radius_type): 
     572        # type: (CallDetails, np.ndarray, float, bool, int) -> np.ndarray 
     573        env = environment() 
     574        queue = env.queue[self._model.dtype] 
     575        context = queue.context 
     576 
     577        # Arrange data transfer to card. 
    562578        details_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
    563579                              hostbuf=call_details.buffer) 
     
    565581                             hostbuf=values) 
    566582 
    567         kernel = self.kernel[1 if magnetic else 0] 
    568         args = [ 
    569             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), 
     583        # Setup kernel function and arguments. 
     584        name = 'Iq' if self.dim == '1d' else 'Imagnetic' if magnetic else 'Iqxy' 
     585        kernel = self._model.get_function(name) 
     586        kernel_args = [ 
     587            np.uint32(self.q_input.nq),  # Number of inputs. 
     588            None,  # Placeholder for pd_start. 
     589            None,  # Placeholder for pd_stop. 
     590            details_b,  # Problem definition. 
     591            values_b,  # Parameter values. 
     592            self.q_input.q_b,  # Q values. 
     593            self._result_b,   # Result storage. 
     594            self._as_dtype(cutoff),  # Probability cutoff. 
     595            np.uint32(effective_radius_type),  # R_eff mode. 
    572596        ] 
     597 
     598        # Call kernel and retrieve results. 
    573599        #print("Calling OpenCL") 
    574600        #call_details.show(values) 
    575         # Call kernel and retrieve results 
    576601        wait_for = None 
    577         last_nap = time.clock() 
     602        last_nap = clock() 
    578603        step = 1000000//self.q_input.nq + 1 
    579604        for start in range(0, call_details.num_eval, step): 
    580605            stop = min(start + step, call_details.num_eval) 
    581606            #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)] 
     607            kernel_args[1:3] = [np.int32(start), np.int32(stop)] 
     608            wait_for = [kernel(queue, self.q_input.global_size, None, 
     609                               *kernel_args, wait_for=wait_for)] 
    585610            if stop < call_details.num_eval: 
    586                 # Allow other processes to run 
     611                # Allow other processes to run. 
    587612                wait_for[0].wait() 
    588                 current_time = time.clock() 
     613                current_time = clock() 
    589614                if current_time - last_nap > 0.5: 
    590                     time.sleep(0.05) 
     615                    time.sleep(0.001) 
    591616                    last_nap = current_time 
    592         cl.enqueue_copy(self.queue, self.result, self.result_b) 
     617        cl.enqueue_copy(queue, self.result, self._result_b, wait_for=wait_for) 
    593618        #print("result", self.result) 
    594619 
    595         # Free buffers 
    596         for v in (details_b, values_b): 
    597             if v is not None: 
    598                 v.release() 
    599  
    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] 
     620        # Free buffers. 
     621        details_b.release() 
     622        values_b.release() 
    606623 
    607624    def release(self): 
     
    610627        Release resources associated with the kernel. 
    611628        """ 
    612         for v in self._need_release: 
    613             v.release() 
    614         self._need_release = [] 
     629        self.q_input.release() 
     630        if self._result_b is not None: 
     631            self._result_b.release() 
     632            self._result_b = None 
    615633 
    616634    def __del__(self): 
  • sasmodels/kerneldll.py

    r1a3559f r3199b17  
    100100# pylint: enable=unused-import 
    101101 
     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')) 
     104 
    102105if "SAS_DLL_PATH" in os.environ: 
    103106    SAS_DLL_PATH = os.environ["SAS_DLL_PATH"] 
     
    112115        COMPILER = "tinycc" 
    113116    elif "VCINSTALLDIR" in os.environ: 
    114         # If vcvarsall.bat has been called, then VCINSTALLDIR is in the environment 
    115         # and we can use the MSVC compiler.  Otherwise, if tinycc is available 
    116         # the use it.  Otherwise, hope that mingw is available. 
     117        # If vcvarsall.bat has been called, then VCINSTALLDIR is in the 
     118        # environment and we can use the MSVC compiler.  Otherwise, if 
     119        # tinycc is available then use it.  Otherwise, hope that mingw 
     120        # is available. 
    117121        COMPILER = "msvc" 
    118122    else: 
     
    121125    COMPILER = "unix" 
    122126 
    123 ARCH = "" if ct.sizeof(ct.c_void_p) > 4 else "x86"  # 4 byte pointers on x86 
     127ARCH = "" if ct.sizeof(ct.c_void_p) > 4 else "x86"  # 4 byte pointers on x86. 
    124128if COMPILER == "unix": 
    125     # Generic unix compile 
    126     # On mac users will need the X code command line tools installed 
     129    # Generic unix compile. 
     130    # On Mac users will need the X code command line tools installed. 
    127131    #COMPILE = "gcc-mp-4.7 -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm -lgomp" 
    128132    CC = "cc -shared -fPIC -std=c99 -O2 -Wall".split() 
    129     # add openmp support if not running on a mac 
     133    # Add OpenMP support if not running on a Mac. 
    130134    if sys.platform != "darwin": 
    131         # OpenMP seems to be broken on gcc 5.4.0 (ubuntu 16.04.9) 
     135        # OpenMP seems to be broken on gcc 5.4.0 (ubuntu 16.04.9). 
    132136        # Shut it off for all unix until we can investigate. 
    133137        #CC.append("-fopenmp") 
     
    141145    # vcomp90.dll on the path.  One may be found here: 
    142146    #       C:/Windows/winsxs/x86_microsoft.vc90.openmp*/vcomp90.dll 
    143     # Copy this to the python directory and uncomment the OpenMP COMPILE 
    144     # TODO: remove intermediate OBJ file created in the directory 
    145     # TODO: maybe don't use randomized name for the c file 
    146     # TODO: maybe ask distutils to find MSVC 
     147    # Copy this to the python directory and uncomment the OpenMP COMPILE. 
     148    # TODO: Remove intermediate OBJ file created in the directory. 
     149    # TODO: Maybe don't use randomized name for the c file. 
     150    # TODO: Maybe ask distutils to find MSVC. 
    147151    CC = "cl /nologo /Ox /MD /W3 /GS- /DNDEBUG".split() 
    148152    if "SAS_OPENMP" in os.environ: 
     
    169173ALLOW_SINGLE_PRECISION_DLLS = True 
    170174 
     175 
    171176def compile(source, output): 
    172177    # type: (str, str) -> None 
     
    180185    logging.info(command_str) 
    181186    try: 
    182         # need shell=True on windows to keep console box from popping up 
     187        # Need shell=True on windows to keep console box from popping up. 
    183188        shell = (os.name == 'nt') 
    184189        subprocess.check_output(command, shell=shell, stderr=subprocess.STDOUT) 
    185190    except subprocess.CalledProcessError as exc: 
    186         raise RuntimeError("compile failed.\n%s\n%s"%(command_str, exc.output)) 
     191        output = decode(exc.output) 
     192        raise RuntimeError("compile failed.\n%s\n%s"%(command_str, output)) 
    187193    if not os.path.exists(output): 
    188194        raise RuntimeError("compile failed.  File is in %r"%source) 
     195 
    189196 
    190197def dll_name(model_info, dtype): 
     
    198205    basename += ARCH + ".so" 
    199206 
    200     # Hack to find precompiled dlls 
     207    # Hack to find precompiled dlls. 
    201208    path = joinpath(generate.DATA_PATH, '..', 'compiled_models', basename) 
    202209    if os.path.exists(path): 
     
    238245        raise ValueError("16 bit floats not supported") 
    239246    if dtype == F32 and not ALLOW_SINGLE_PRECISION_DLLS: 
    240         dtype = F64  # Force 64-bit dll 
    241     # Note: dtype may be F128 for long double precision 
     247        dtype = F64  # Force 64-bit dll. 
     248    # Note: dtype may be F128 for long double precision. 
    242249 
    243250    dll = dll_path(model_info, dtype) 
     
    250257        need_recompile = dll_time < newest_source 
    251258    if need_recompile: 
    252         # Make sure the DLL path exists 
     259        # Make sure the DLL path exists. 
    253260        if not os.path.exists(SAS_DLL_PATH): 
    254261            os.makedirs(SAS_DLL_PATH) 
     
    259266            file_handle.write(source) 
    260267        compile(source=filename, output=dll) 
    261         # comment the following to keep the generated c file 
    262         # Note: if there is a syntax error then compile raises an error 
     268        # Comment the following to keep the generated C file. 
     269        # Note: If there is a syntax error then compile raises an error 
    263270        # and the source file will not be deleted. 
    264271        os.unlink(filename) 
     
    299306        self.dllpath = dllpath 
    300307        self._dll = None  # type: ct.CDLL 
    301         self._kernels = None # type: List[Callable, Callable] 
     308        self._kernels = None  # type: List[Callable, Callable] 
    302309        self.dtype = np.dtype(dtype) 
    303310 
     
    315322 
    316323        # int, int, int, int*, double*, double*, double*, double*, double 
    317         argtypes = [ct.c_int32]*3 + [ct.c_void_p]*4 + [float_type] 
     324        argtypes = [ct.c_int32]*3 + [ct.c_void_p]*4 + [float_type, ct.c_int32] 
    318325        names = [generate.kernel_name(self.info, variant) 
    319326                 for variant in ("Iq", "Iqxy", "Imagnetic")] 
     
    334341        # type: (List[np.ndarray]) -> DllKernel 
    335342        q_input = PyInput(q_vectors, self.dtype) 
    336         # Note: pickle not supported for DllKernel 
     343        # Note: DLL is lazy loaded. 
    337344        if self._dll is None: 
    338345            self._load_dll() 
     
    354361        self._dll = None 
    355362 
     363 
    356364class DllKernel(Kernel): 
    357365    """ 
     
    375383    def __init__(self, kernel, model_info, q_input): 
    376384        # type: (Callable[[], np.ndarray], ModelInfo, PyInput) -> None 
     385        dtype = q_input.dtype 
     386        self.q_input = q_input 
    377387        self.kernel = kernel 
     388 
     389        # Attributes accessed from the outside. 
     390        self.dim = '2d' if q_input.is_2d else '1d' 
    378391        self.info = model_info 
    379         self.q_input = q_input 
    380         self.dtype = q_input.dtype 
    381         self.dim = '2d' if q_input.is_2d else '1d' 
    382         self.result = np.empty(q_input.nq+1, q_input.dtype) 
    383         self.real = (np.float32 if self.q_input.dtype == generate.F32 
    384                      else np.float64 if self.q_input.dtype == generate.F64 
    385                      else np.float128) 
    386  
    387     def __call__(self, call_details, values, cutoff, magnetic): 
    388         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    389  
     392        self.dtype = dtype 
     393 
     394        # Converter to translate input to target type. 
     395        self._as_dtype = (np.float32 if dtype == generate.F32 
     396                          else np.float64 if dtype == generate.F64 
     397                          else np.float128) 
     398 
     399        # Holding place for the returned value. 
     400        nout = 2 if self.info.have_Fq else 1 
     401        extra_q = 4  # Total weight, form volume, shell volume and R_eff. 
     402        self.result = np.empty(self.q_input.nq*nout + extra_q, dtype) 
     403 
     404    def _call_kernel(self, call_details, values, cutoff, magnetic, 
     405                     effective_radius_type): 
     406        # type: (CallDetails, np.ndarray, float, bool, int) -> np.ndarray 
     407 
     408        # Setup kernel function and arguments. 
    390409        kernel = self.kernel[1 if magnetic else 0] 
    391         args = [ 
    392             self.q_input.nq, # nq 
    393             None, # pd_start 
    394             None, # pd_stop pd_stride[MAX_PD] 
    395             call_details.buffer.ctypes.data, # problem 
    396             values.ctypes.data,  #pars 
    397             self.q_input.q.ctypes.data, #q 
    398             self.result.ctypes.data,   # results 
    399             self.real(cutoff), # cutoff 
     410        kernel_args = [ 
     411            self.q_input.nq,  # Number of inputs. 
     412            None,  # Placeholder for pd_start. 
     413            None,  # Placeholder for pd_stop. 
     414            call_details.buffer.ctypes.data,  # Problem definition. 
     415            values.ctypes.data,  # Parameter values. 
     416            self.q_input.q.ctypes.data,  # Q values. 
     417            self.result.ctypes.data,   # Result storage. 
     418            self._as_dtype(cutoff),  # Probability cutoff. 
     419            effective_radius_type,  # R_eff mode. 
    400420        ] 
     421 
     422        # Call kernel and retrieve results. 
    401423        #print("Calling DLL") 
    402424        #call_details.show(values) 
    403425        step = 100 
     426        # TODO: Do we need the explicit sleep like the OpenCL and CUDA loops? 
    404427        for start in range(0, call_details.num_eval, step): 
    405428            stop = min(start + step, call_details.num_eval) 
    406             args[1:3] = [start, stop] 
    407             kernel(*args) # type: ignore 
    408  
    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 
     429            kernel_args[1:3] = [start, stop] 
     430            kernel(*kernel_args) # type: ignore 
    415431 
    416432    def release(self): 
    417433        # type: () -> None 
    418434        """ 
    419         Release any resources associated with the kernel. 
     435        Release resources associated with the kernel. 
    420436        """ 
    421         self.q_input.release() 
     437        # TODO: OpenCL/CUDA allocate q_input in __init__ and free it in release. 
     438        # Should we be doing the same for DLL? 
     439        #self.q_input.release() 
     440        pass 
     441 
     442    def __del__(self): 
     443        # type: () -> None 
     444        self.release() 
  • sasmodels/kernelpy.py

    r12eec1e r3199b17  
    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 
     
    2833logger = logging.getLogger(__name__) 
    2934 
     35 
    3036class PyModel(KernelModel): 
    3137    """ 
     
    3339    """ 
    3440    def __init__(self, model_info): 
    35         # Make sure Iq is available and vectorized 
     41        # Make sure Iq is available and vectorized. 
    3642        _create_default_functions(model_info) 
    3743        self.info = model_info 
     
    4854        """ 
    4955        pass 
     56 
    5057 
    5158class PyInput(object): 
     
    8693        self.q = None 
    8794 
     95 
    8896class PyKernel(Kernel): 
    8997    """ 
     
    126134        parameter_vector = np.empty(len(partable.call_parameters)-2, 'd') 
    127135 
    128         # Create views into the array to hold the arguments 
     136        # Create views into the array to hold the arguments. 
    129137        offset = 0 
    130138        kernel_args, volume_args = [], [] 
     
    158166 
    159167        # 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): 
     168        self._volume_args = volume_args 
     169        volume = model_info.form_volume 
     170        shell = model_info.shell_volume 
     171        radius = model_info.effective_radius 
     172        self._volume = ((lambda: (shell(*volume_args), volume(*volume_args))) if shell and volume 
     173                        else (lambda: [volume(*volume_args)]*2) if volume 
     174                        else (lambda: (1.0, 1.0))) 
     175        self._radius = ((lambda mode: radius(mode, *volume_args)) if radius 
     176                        else (lambda mode: cbrt(0.75/pi*volume(*volume_args))) if volume 
     177                        else (lambda mode: 1.0)) 
     178 
     179    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
    165180        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    166181        if magnetic: 
     
    168183        #print("Calling python kernel") 
    169184        #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 
     185        radius = ((lambda: 0.0) if effective_radius_type == 0 
     186                  else (lambda: self._radius(effective_radius_type))) 
     187        self.result = _loops( 
     188            self._parameter_vector, self._form, self._volume, radius, 
     189            self.q_input.nq, call_details, values, cutoff) 
    173190 
    174191    def release(self): 
     
    179196        self.q_input.release() 
    180197        self.q_input = None 
     198 
    181199 
    182200def _loops(parameters,    # type: np.ndarray 
    183201           form,          # type: Callable[[], np.ndarray] 
    184202           form_volume,   # type: Callable[[], float] 
     203           form_radius,   # type: Callable[[], float] 
    185204           nq,            # type: int 
    186205           call_details,  # type: details.CallDetails 
    187206           values,        # type: np.ndarray 
    188            cutoff         # type: float 
     207           cutoff,        # type: float 
    189208          ): 
    190209    # type: (...) -> None 
     
    198217    #                                                              # 
    199218    ################################################################ 
     219 
     220    # WARNING: Trickery ahead 
     221    # The parameters[] vector is embedded in the closures for form(), 
     222    # form_volume() and form_radius().  We set the initial vector from 
     223    # the values for the model parameters. As we loop through the polydispesity 
     224    # mesh, we update the components with the polydispersity values before 
     225    # calling the respective functions. 
    200226    n_pars = len(parameters) 
    201227    parameters[:] = values[2:n_pars+2] 
     228 
    202229    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 
     230        total = form() 
     231        weight_norm = 1.0 
     232        weighted_shell, weighted_form = form_volume() 
     233        weighted_radius = form_radius() 
     234 
     235    else: 
     236        pd_value = values[2+n_pars:2+n_pars + call_details.num_weights] 
     237        pd_weight = values[2+n_pars + call_details.num_weights:] 
     238 
     239        weight_norm = 0.0 
     240        weighted_form = 0.0 
     241        weighted_shell = 0.0 
     242        weighted_radius = 0.0 
     243        partial_weight = np.NaN 
     244        weight = np.NaN 
     245 
     246        p0_par = call_details.pd_par[0] 
     247        p0_length = call_details.pd_length[0] 
     248        p0_index = p0_length 
     249        p0_offset = call_details.pd_offset[0] 
     250 
     251        pd_par = call_details.pd_par[:call_details.num_active] 
     252        pd_offset = call_details.pd_offset[:call_details.num_active] 
     253        pd_stride = call_details.pd_stride[:call_details.num_active] 
     254        pd_length = call_details.pd_length[:call_details.num_active] 
     255 
     256        total = np.zeros(nq, 'd') 
     257        for loop_index in range(call_details.num_eval): 
     258            # Update polydispersity parameter values. 
     259            if p0_index == p0_length: 
     260                pd_index = (loop_index//pd_stride)%pd_length 
     261                parameters[pd_par] = pd_value[pd_offset+pd_index] 
     262                partial_weight = np.prod(pd_weight[pd_offset+pd_index][1:]) 
     263                p0_index = loop_index%p0_length 
     264 
     265            weight = partial_weight * pd_weight[p0_offset + p0_index] 
     266            parameters[p0_par] = pd_value[p0_offset + p0_index] 
     267            p0_index += 1 
     268            if weight > cutoff: 
     269                # Call the scattering function. 
     270                # Assume that NaNs are only generated if the parameters are bad; 
     271                # exclude all q for that NaN.  Even better would be to have an 
     272                # INVALID expression like the C models, but that is expensive. 
     273                Iq = np.asarray(form(), 'd') 
     274                if np.isnan(Iq).any(): 
     275                    continue 
     276 
     277                # Update value and norm. 
     278                total += weight * Iq 
     279                weight_norm += weight 
     280                shell, form = form_volume() 
     281                weighted_form += weight * form 
     282                weighted_shell += weight * shell 
     283                weighted_radius += weight * form_radius() 
     284 
     285    result = np.hstack((total, weight_norm, weighted_form, weighted_shell, weighted_radius)) 
     286    return result 
    253287 
    254288 
     
    261295    any functions that are not already marked as vectorized. 
    262296    """ 
    263     # Note: must call create_vector_Iq before create_vector_Iqxy 
     297    # Note: Must call create_vector_Iq before create_vector_Iqxy. 
    264298    _create_vector_Iq(model_info) 
    265299    _create_vector_Iqxy(model_info) 
  • sasmodels/mixture.py

    recf895e 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 rd92182f  
    55Usage:: 
    66 
    7     python -m sasmodels.model_test [opencl|dll|opencl_and_dll] model1 model2 ... 
    8  
    9     if model1 is 'all', then all except the remaining models will be tested 
     7    python -m sasmodels.model_test [opencl|cuda|dll|all] model1 model2 ... 
     8 
     9If model1 is 'all', then all except the remaining models will be tested. 
     10Subgroups are also possible, such as 'py', 'single' or '1d'.  See 
     11:func:`core.list_models` for details. 
    1012 
    1113Each 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. 
     14and Fq is called to make sure R_eff, volume and volume ratio are computed. 
     15The return values at these points are not considered.  The test is only to 
     16verify that the models run to completion, and do not produce inf or NaN. 
    1517 
    1618Tests are defined with the *tests* attribute in the model.py file.  *tests* 
    1719is a list of individual tests to run, where each test consists of the 
    1820parameter 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. 
     21the effective radius test and volume ratio tests, use the extended output 
     22form, which checks each output of kernel.Fq. For 1-D tests, either specify 
     23the q value or a list of q-values, and the corresponding I(q) value, or 
     24list of I(q) values. 
    2225 
    2326That is:: 
     
    3235                        [I(qx1, qy1), I(qx2, qy2), ...]], 
    3336 
    34         [ {parameters}, 'ER', ER(pars) ], 
    35         [ {parameters}, 'VR', VR(pars) ], 
     37        [ {parameters}, q, F(q), F^2(q), R_eff, V, V_r ], 
    3638        ... 
    3739    ] 
     
    4547from __future__ import print_function 
    4648 
     49import argparse 
    4750import sys 
    4851import unittest 
     
    6063from . import core 
    6164from .core import list_models, load_model_info, build_model 
    62 from .direct_model import call_kernel, call_ER, call_VR 
     65from .direct_model import call_kernel, call_Fq 
    6366from .exception import annotate_exception 
    6467from .modelinfo import expand_pars 
    6568from .kernelcl import use_opencl 
     69from .kernelcuda import use_cuda 
     70from . import product 
    6671 
    6772# pylint: disable=unused-import 
     
    8085    Construct the pyunit test suite. 
    8186 
    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. 
     87    *loaders* is the list of kernel drivers to use (dll, opencl or cuda). 
     88    For python model the python driver is always used. 
    8589 
    8690    *models* is the list of models to test, or *["all"]* to test all models. 
     
    8892    suite = unittest.TestSuite() 
    8993 
    90     if models[0] in core.KINDS: 
     94    try: 
     95        # See if the first model parses as a model group 
     96        group = list_models(models[0]) 
    9197        skip = models[1:] 
    92         models = list_models(models[0]) 
    93     else: 
     98        models = group 
     99    except Exception: 
    94100        skip = [] 
    95101    for model_name in models: 
     
    160166                                    test_method_name, 
    161167                                    platform="ocl", dtype=None, 
     168                                    stash=stash) 
     169            #print("defining", test_name) 
     170            suite.addTest(test) 
     171 
     172        # test using cuda if desired and available 
     173        if 'cuda' in loaders and use_cuda(): 
     174            test_name = "%s-cuda" % model_info.id 
     175            test_method_name = "test_%s_cuda" % model_info.id 
     176            # Using dtype=None so that the models that are only 
     177            # correct for double precision are not tested using 
     178            # single precision.  The choice is determined by the 
     179            # presence of *single=False* in the model file. 
     180            test = ModelTestCase(test_name, model_info, 
     181                                    test_method_name, 
     182                                    platform="cuda", dtype=None, 
    162183                                    stash=stash) 
    163184            #print("defining", test_name) 
     
    202223                ({}, [0.001, 0.01, 0.1], [None]*3), 
    203224                ({}, [(0.1, 0.1)]*2, [None]*2), 
    204                 # test that ER/VR will run if they exist 
    205                 ({}, 'ER', None), 
    206                 ({}, 'VR', None), 
     225                # test that Fq will run, and return R_eff, V, V_r 
     226                ({}, 0.1, None, None, None, None, None), 
    207227                ] 
    208228            tests = smoke_tests 
     
    210230            if self.info.tests is not None: 
    211231                tests += self.info.tests 
     232            S_tests = [test for test in tests if '@S' in test[0]] 
     233            P_tests = [test for test in tests if '@S' not in test[0]] 
    212234            try: 
    213235                model = build_model(self.info, dtype=self.dtype, 
    214236                                    platform=self.platform) 
    215                 results = [self.run_one(model, test) for test in tests] 
     237                results = [self.run_one(model, test) for test in P_tests] 
     238                for test in S_tests: 
     239                    # pull the S model name out of the test defn 
     240                    pars = test[0].copy() 
     241                    s_name = pars.pop('@S') 
     242                    ps_test = [pars] + list(test[1:]) 
     243                    # build the P@S model 
     244                    s_info = load_model_info(s_name) 
     245                    ps_info = product.make_product_info(self.info, s_info) 
     246                    ps_model = build_model(ps_info, dtype=self.dtype, 
     247                                           platform=self.platform) 
     248                    # run the tests 
     249                    results.append(self.run_one(ps_model, ps_test)) 
     250 
    216251                if self.stash: 
    217252                    for test, target, actual in zip(tests, self.stash[0], results): 
     
    223258 
    224259                # Check for missing tests.  Only do so for the "dll" tests 
    225                 # to reduce noise from both opencl and dll, and because 
     260                # to reduce noise from both opencl and cuda, and because 
    226261                # python kernels use platform="dll". 
    227262                if self.platform == "dll": 
     
    238273        def _find_missing_tests(self): 
    239274            # 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) 
     275            """make sure there are 1D and 2D tests as appropriate""" 
    243276            model_has_1D = True 
    244277            model_has_2D = any(p.type == 'orientation' 
     
    248281            single = [test for test in self.info.tests 
    249282                      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) 
    252283            tests_has_1D_single = any(isinstance(test[1], float) for test in single) 
    253284            tests_has_2D_single = any(isinstance(test[1], tuple) for test in single) 
     
    262293 
    263294            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") 
    268295            if model_has_1D and not (tests_has_1D_single or tests_has_1D_multiple): 
    269296                missing.append("1D") 
     
    276303            # type: (KernelModel, TestCondition) -> None 
    277304            """Run a single test case.""" 
    278             user_pars, x, y = test 
     305            user_pars, x, y = test[:3] 
    279306            pars = expand_pars(self.info.parameters, user_pars) 
    280307            invalid = invalid_pars(self.info.parameters, pars) 
     
    289316            self.assertEqual(len(y), len(x)) 
    290317 
    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): 
     318            if isinstance(x[0], tuple): 
    296319                qx, qy = zip(*x) 
    297320                q_vectors = [np.array(qx), np.array(qy)] 
    298                 kernel = model.make_kernel(q_vectors) 
    299                 actual = call_kernel(kernel, pars) 
    300321            else: 
    301322                q_vectors = [np.array(x)] 
    302                 kernel = model.make_kernel(q_vectors) 
     323 
     324            kernel = model.make_kernel(q_vectors) 
     325            if len(test) == 3: 
    303326                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): 
     327                self._check_vectors(x, y, actual, 'I') 
     328                return actual 
     329            else: 
     330                y1 = y 
     331                y2 = test[3] if not isinstance(test[3], list) else [test[3]] 
     332                F1, F2, R_eff, volume, volume_ratio = call_Fq(kernel, pars) 
     333                if F1 is not None:  # F1 is none for models with Iq instead of Fq 
     334                    self._check_vectors(x, y1, F1, 'F') 
     335                self._check_vectors(x, y2, F2, 'F^2') 
     336                self._check_scalar(test[4], R_eff, 'R_eff') 
     337                self._check_scalar(test[5], volume, 'volume') 
     338                self._check_scalar(test[6], volume_ratio, 'form:shell ratio') 
     339                return F2 
     340 
     341        def _check_scalar(self, target, actual, name): 
     342            if target is None: 
     343                # smoke test --- make sure it runs and produces a value 
     344                self.assertTrue(not np.isnan(actual), 
     345                                'invalid %s: %s' % (name, actual)) 
     346            elif np.isnan(target): 
     347                # make sure nans match 
     348                self.assertTrue(np.isnan(actual), 
     349                                '%s: expected:%s; actual:%s' 
     350                                % (name, target, actual)) 
     351            else: 
     352                # is_near does not work for infinite values, so also test 
     353                # for exact values. 
     354                self.assertTrue(target == actual or is_near(target, actual, 5), 
     355                                '%s: expected:%s; actual:%s' 
     356                                % (name, target, actual)) 
     357 
     358        def _check_vectors(self, x, target, actual, name='I'): 
     359            self.assertTrue(len(actual) > 0, 
     360                            '%s(...) expected return'%name) 
     361            if target is None: 
     362                return 
     363            self.assertEqual(len(target), len(actual), 
     364                             '%s(...) returned wrong length'%name) 
     365            for xi, yi, actual_yi in zip(x, target, actual): 
    309366                if yi is None: 
    310367                    # smoke test --- make sure it runs and produces a value 
    311368                    self.assertTrue(not np.isnan(actual_yi), 
    312                                     'invalid f(%s): %s' % (xi, actual_yi)) 
     369                                    'invalid %s(%s): %s' % (name, xi, actual_yi)) 
    313370                elif np.isnan(yi): 
     371                    # make sure nans match 
    314372                    self.assertTrue(np.isnan(actual_yi), 
    315                                     'f(%s): expected:%s; actual:%s' 
    316                                     % (xi, yi, actual_yi)) 
     373                                    '%s(%s): expected:%s; actual:%s' 
     374                                    % (name, xi, yi, actual_yi)) 
    317375                else: 
    318376                    # is_near does not work for infinite values, so also test 
    319                     # for exact values.  Note that this will not 
     377                    # for exact values. 
    320378                    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 
     379                                    '%s(%s); expected:%s; actual:%s' 
     380                                    % (name, xi, yi, actual_yi)) 
    324381 
    325382    return ModelTestCase 
     
    333390    invalid = [] 
    334391    for par in sorted(pars.keys()): 
     392        # special handling of R_eff mode, which is not a usual parameter 
     393        if par == product.RADIUS_MODE_ID: 
     394            continue 
    335395        parts = par.split('_pd') 
    336396        if len(parts) > 1 and parts[1] not in ("", "_n", "nsigma", "type"): 
     
    383443 
    384444    # Build a test suite containing just the model 
    385     loaders = ['opencl'] if use_opencl() else ['dll'] 
     445    loaders = ['opencl' if use_opencl() else 'cuda' if use_cuda() else 'dll'] 
    386446    suite = unittest.TestSuite() 
    387447    _add_model_to_suite(loaders, suite, model_info) 
     
    419479 
    420480 
    421 def main(*models): 
    422     # type: (*str) -> int 
    423     """ 
    424     Run tests given is models. 
    425  
    426     Returns 0 if success or 1 if any tests fail. 
    427     """ 
    428     try: 
    429         from xmlrunner import XMLTestRunner as TestRunner 
    430         test_args = {'output': 'logs'} 
    431     except ImportError: 
    432         from unittest import TextTestRunner as TestRunner 
    433         test_args = {} 
    434  
    435     if models and models[0] == '-v': 
    436         verbosity = 2 
    437         models = models[1:] 
    438     else: 
    439         verbosity = 1 
    440     if models and models[0] == 'opencl': 
    441         if not use_opencl(): 
    442             print("opencl is not available") 
    443             return 1 
    444         loaders = ['opencl'] 
    445         models = models[1:] 
    446     elif models and models[0] == 'dll': 
    447         # TODO: test if compiler is available? 
    448         loaders = ['dll'] 
    449         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:] 
    453     else: 
    454         loaders = ['opencl', 'dll'] if use_opencl() else ['dll'] 
    455     if not models: 
    456         print("""\ 
    457 usage: 
    458   python -m sasmodels.model_test [-v] [opencl|dll] model1 model2 ... 
    459  
    460 If -v is included on the command line, then use verbose output. 
    461  
    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. 
    464  
    465 If model1 is 'all', then all except the remaining models will be tested. 
    466  
    467 """) 
    468  
    469         return 1 
    470  
    471     runner = TestRunner(verbosity=verbosity, **test_args) 
    472     result = runner.run(make_suite(loaders, models)) 
    473     return 1 if result.failures or result.errors else 0 
    474  
    475  
    476481def model_tests(): 
    477482    # type: () -> Iterator[Callable[[], None]] 
     
    481486    Run "nosetests sasmodels" on the command line to invoke it. 
    482487    """ 
    483     loaders = ['opencl', 'dll'] if use_opencl() else ['dll'] 
     488    loaders = ['dll'] 
     489    if use_opencl(): 
     490        loaders.append('opencl') 
     491    if use_cuda(): 
     492        loaders.append('cuda') 
    484493    tests = make_suite(loaders, ['all']) 
    485494    def build_test(test): 
     
    504513 
    505514 
     515def main(): 
     516    # type: (*str) -> int 
     517    """ 
     518    Run tests given is models. 
     519 
     520    Returns 0 if success or 1 if any tests fail. 
     521    """ 
     522    try: 
     523        from xmlrunner import XMLTestRunner as TestRunner 
     524        test_args = {'output': 'logs'} 
     525    except ImportError: 
     526        from unittest import TextTestRunner as TestRunner 
     527        test_args = {} 
     528 
     529    parser = argparse.ArgumentParser(description="Test SasModels Models") 
     530    parser.add_argument("-v", "--verbose", action="store_const", 
     531                        default=1, const=2, help="Use verbose output") 
     532    parser.add_argument("-e", "--engine", default="all", 
     533                        help="Engines on which to run the test.  " 
     534                        "Valid values are opencl, cuda, dll, and all. " 
     535                        "Defaults to all if no value is given") 
     536    parser.add_argument("models", nargs="*", 
     537                        help="The names of the models to be tested.  " 
     538                        "If the first model is 'all', then all but the listed " 
     539                        "models will be tested.  See core.list_models() for " 
     540                        "names of other groups, such as 'py' or 'single'.") 
     541    args, models = parser.parse_known_args() 
     542 
     543    if args.engine == "opencl": 
     544        if not use_opencl(): 
     545            print("opencl is not available") 
     546            return 1 
     547        loaders = ['opencl'] 
     548    elif args.engine == "dll": 
     549        loaders = ["dll"] 
     550    elif args.engine == "cuda": 
     551        if not use_cuda(): 
     552            print("cuda is not available") 
     553            return 1 
     554        loaders = ['cuda'] 
     555    elif args.engine == "all": 
     556        loaders = ['dll'] 
     557        if use_opencl(): 
     558            loaders.append('opencl') 
     559        if use_cuda(): 
     560            loaders.append('cuda') 
     561    else: 
     562        print("unknown engine " + args.engine) 
     563        return 1 
     564 
     565    runner = TestRunner(verbosity=args.verbose, **test_args) 
     566    result = runner.run(make_suite(loaders, args.models)) 
     567    return 1 if result.failures or result.errors else 0 
     568 
     569 
    506570if __name__ == "__main__": 
    507     sys.exit(main(*sys.argv[1:])) 
     571    sys.exit(main()) 
  • sasmodels/modelinfo.py

    rbd547d0 rd8eaa3d  
    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 
     
    291290    example, might be used to set the value of a shape parameter. 
    292291 
    293     These values are set by :func:`make_parameter_table` and 
     292    Control parameters are used for variant models such as :ref:`rpa` which 
     293    have different cases with different parameters, as well as models 
     294    like :ref:`spherical_sld` with its user defined number of shells. 
     295    The control parameter should appear in the parameter table along with the 
     296    parameters it is is controlling.  For variant models, use *[CASES]* in 
     297    place of the parameter limits Within the parameter definition table, 
     298    with case names such as:: 
     299 
     300         CASES = ["diblock copolymer", "triblock copolymer", ...] 
     301 
     302    This should give *limits=[[case1, case2, ...]]*, but the model loader 
     303    translates it to *limits=[0, len(CASES)-1]*, and adds *choices=CASES* to 
     304    the :class:`Parameter` definition. Note that models can use a list of 
     305    cases as a parameter without it being a control parameter.  Either way, 
     306    the parameter is sent to the model evaluator as *float(choice_num)*, 
     307    where choices are numbered from 0. :meth:`ModelInfo.get_hidden_parameters` 
     308    will determine which parameers to display. 
     309 
     310    The class contructor should not be called directly, but instead the 
     311    parameter table is built using :func:`make_parameter_table` and 
    294312    :func:`parse_parameter` therein. 
    295313    """ 
     
    405423      parameters counted as n individual parameters p1, p2, ... 
    406424 
     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 
    407429    * *call_parameters* is the complete list of parameters to the kernel, 
    408430      including scale and background, with vector parameters recorded as 
     
    417439    parameters don't use vector notation, and instead use p1, p2, ... 
    418440    """ 
    419     # scale and background are implicit parameters 
    420     COMMON = [Parameter(*p) for p in COMMON_PARAMETERS] 
    421  
    422441    def __init__(self, parameters): 
    423442        # 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 
    424449        self.kernel_parameters = parameters 
    425450        self._set_vector_lengths() 
    426  
    427451        self.npars = sum(p.length for p in self.kernel_parameters) 
    428452        self.nmagnetic = sum(p.length for p in self.kernel_parameters 
     
    431455        if self.nmagnetic: 
    432456            self.nvalues += 3 + 3*self.nmagnetic 
    433  
    434457        self.call_parameters = self._get_call_parameters() 
    435458        self.defaults = self._get_defaults() 
     
    471494                         if p.polydisperse and p.type not in ('orientation', 'magnetic')) 
    472495        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() 
    473507 
    474508    def check_angles(self): 
     
    570604    def _get_call_parameters(self): 
    571605        # type: () -> List[Parameter] 
    572         full_list = self.COMMON[:] 
     606        full_list = self.common_parameters[:] 
    573607        for p in self.kernel_parameters: 
    574608            if p.length == 1: 
     
    673707 
    674708        # Gather the user parameters in order 
    675         result = control + self.COMMON 
     709        result = control + self.common_parameters 
    676710        for p in self.kernel_parameters: 
    677711            if not is2d and p.type in ('orientation', 'magnetic'): 
     
    722756 
    723757#: Set of variables defined in the model that might contain C code 
    724 C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'c_code'] 
     758C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'shell_volume', 'c_code'] 
    725759 
    726760def _find_source_lines(model_info, kernel_module): 
     
    771805        # Custom sum/multi models 
    772806        return kernel_module.model_info 
     807 
    773808    info = ModelInfo() 
     809 
     810    # Build the parameter table 
    774811    #print("make parameter table", kernel_module.parameters) 
    775812    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. 
    776823    demo = expand_pars(parameters, getattr(kernel_module, 'demo', None)) 
     824 
    777825    filename = abspath(kernel_module.__file__).replace('.pyc', '.py') 
    778826    kernel_id = splitext(basename(filename))[0] 
     
    792840    info.category = getattr(kernel_module, 'category', None) 
    793841    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) 
    794845    info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y']) 
    795846    # Note: custom.load_custom_kernel_module assumes the C sources are defined 
     
    797848    info.source = getattr(kernel_module, 'source', []) 
    798849    info.c_code = getattr(kernel_module, 'c_code', None) 
     850    info.effective_radius = getattr(kernel_module, 'effective_radius', None) 
    799851    # TODO: check the structure of the tests 
    800852    info.tests = getattr(kernel_module, 'tests', []) 
    801     info.ER = getattr(kernel_module, 'ER', None) # type: ignore 
    802     info.VR = getattr(kernel_module, 'VR', None) # type: ignore 
    803853    info.form_volume = getattr(kernel_module, 'form_volume', None) # type: ignore 
     854    info.shell_volume = getattr(kernel_module, 'shell_volume', None) # type: ignore 
    804855    info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore 
    805856    info.Iqxy = getattr(kernel_module, 'Iqxy', None) # type: ignore 
     
    813864    info.single = getattr(kernel_module, 'single', not callable(info.Iq)) 
    814865    info.random = getattr(kernel_module, 'random', None) 
    815  
    816     # multiplicity info 
    817     control_pars = [p.id for p in parameters.kernel_parameters if p.is_control] 
    818     default_control = control_pars[0] if control_pars else None 
    819     info.control = getattr(kernel_module, 'control', default_control) 
    820866    info.hidden = getattr(kernel_module, 'hidden', None) # type: ignore 
     867 
     868    # Set control flag for explicitly set parameters, e.g., in the RPA model. 
     869    control = getattr(kernel_module, 'control', None) 
     870    if control is not None: 
     871        parameters[control].is_control = True 
    821872 
    822873    if callable(info.Iq) and parameters.has_2d: 
     
    825876    info.lineno = {} 
    826877    _find_source_lines(info, kernel_module) 
    827  
    828878    return info 
    829879 
     
    872922    #: *sphere*hardsphere* or *cylinder+sphere*. 
    873923    composition = None      # type: Optional[Tuple[str, List[ModelInfo]]] 
    874     #: Name of the control parameter for a variant model such as :ref:`rpa`. 
    875     #: The *control* parameter should appear in the parameter table, with 
    876     #: limits defined as *[CASES]*, for case names such as 
    877     #: *CASES = ["diblock copolymer", "triblock copolymer", ...]*. 
    878     #: This should give *limits=[[case1, case2, ...]]*, but the 
    879     #: model loader translates this to *limits=[0, len(CASES)-1]*, and adds 
    880     #: *choices=CASES* to the :class:`Parameter` definition. Note that 
    881     #: models can use a list of cases as a parameter without it being a 
    882     #: control parameter.  Either way, the parameter is sent to the model 
    883     #: evaluator as *float(choice_num)*, where choices are numbered from 0. 
    884     #: See also :attr:`hidden`. 
    885     control = None          # type: str 
    886924    #: Different variants require different parameters.  In order to show 
    887925    #: just the parameters needed for the variant selected by :attr:`control`, 
     
    918956    #: provided in the file. 
    919957    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] 
    920964    #: List of C source files used to define the model.  The source files 
    921965    #: should define the *Iq* function, and possibly *Iqac* or *Iqabc* if the 
    922966    #: model defines orientation parameters. Files containing the most basic 
    923967    #: functions must appear first in the list, followed by the files that 
    924     #: use those functions.  Form factors are indicated by providing 
    925     #: an :attr:`ER` function. 
     968    #: use those functions. 
    926969    source = None           # type: List[str] 
    927     #: The set of tests that must pass.  The format of the tests is described 
    928     #: in :mod:`model_test`. 
    929     tests = None            # type: List[TestCondition] 
    930     #: Returns the effective radius of the model given its volume parameters. 
    931     #: The presence of *ER* indicates that the model is a form factor model 
    932     #: that may be used together with a structure factor to form an implicit 
    933     #: multiplication model. 
    934     #: 
    935     #: The parameters to the *ER* function must be marked with type *volume*. 
    936     #: in the parameter table.  They will appear in the same order as they 
    937     #: do in the table.  The values passed to *ER* will be vectors, with one 
    938     #: value for each polydispersity condition.  For example, if the model 
    939     #: is polydisperse over both length and radius, then both length and 
    940     #: radius will have the same number of values in the vector, with one 
    941     #: value for each *length X radius*.  If only *radius* is polydisperse, 
    942     #: then the value for *length* will be repeated once for each value of 
    943     #: *radius*.  The *ER* function should return one effective radius for 
    944     #: each parameter set.  Multiplicity parameters will be received as 
    945     #: arrays, with one row per polydispersity condition. 
    946     ER = None               # type: Optional[Callable[[np.ndarray], np.ndarray]] 
    947     #: Returns the occupied volume and the total volume for each parameter set. 
    948     #: See :attr:`ER` for details on the parameters. 
    949     VR = None               # type: Optional[Callable[[np.ndarray], Tuple[np.ndarray, np.ndarray]]] 
    950     #: Arbitrary C code containing supporting functions, etc., to be inserted 
    951     #: after everything in source.  This can include Iq and Iqxy functions with 
    952     #: the full function signature, including all parameters. 
    953     c_code = None 
     970    #: inline source code, added after all elements of source 
     971    c_code = None           # type: Optional[str] 
    954972    #: Returns the form volume for python-based models.  Form volume is needed 
    955973    #: for volume normalization in the polydispersity integral.  If no 
     
    959977    #: C code, either defined as a string, or in the sources. 
    960978    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]] 
     987    #: Computes the effective radius of the shape given the volume parameters. 
     988    #: Only needed for models defined in python that can be used for 
     989    #: monodisperse approximation for non-dilute solutions, P@S.  The first 
     990    #: argument is the integer effective radius mode, with default 0. 
     991    effective_radius = None  # type: Union[None, Callable[[int, np.ndarray], float]] 
    961992    #: Returns *I(q, a, b, ...)* for parameters *a*, *b*, etc. defined 
    962993    #: by the parameter table.  *Iq* can be defined as a python function, or 
     
    9701001    #: will return *I(q, a, b, ...)*.  Multiplicity parameters are sent as 
    9711002    #: pointers to doubles.  Constants in floating point expressions should 
    972     #: include the decimal point. See :mod:`generate` for more details. 
     1003    #: include the decimal point. See :mod:`generate` for more details. If 
     1004    #: *have_Fq* is True, then Iq should return an interleaved array of 
     1005    #: $[\sum F(q_1), \sum F^2(q_1), \ldots, \sum F(q_n), \sum F^2(q_n)]$. 
    9731006    Iq = None               # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
    9741007    #: Returns *I(qab, qc, a, b, ...)*.  The interface follows :attr:`Iq`. 
     
    9931026    #: Line numbers for symbols defining C code 
    9941027    lineno = None           # type: Dict[str, int] 
     1028    #: The set of tests that must pass.  The format of the tests is described 
     1029    #: in :mod:`model_test`. 
     1030    tests = None            # type: List[TestCondition] 
    9951031 
    9961032    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*