Changeset 298d2d4 in sasmodels


Ignore:
Timestamp:
Mar 26, 2019 11:26:07 AM (7 months ago)
Author:
awashington
Branches:
master
Children:
d522352
Parents:
c9d052d (diff), 598a354 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' into log_sesans

Files:
52 added
1 deleted
167 edited

Legend:

Unmodified
Added
Removed
  • .travis.yml

    r335271e rf3767c2  
    66    env: 
    77    - PY=2.7 
    8     - DEPLOY=True 
     8    #- DEPLOY=True 
    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

    r59485a4 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 
     
    139180the compiler. 
    140181 
    141 On Windows, set *SASCOMPILER=tinycc* for the tinycc compiler, 
    142 *SASCOMPILER=msvc* for the Microsoft Visual C compiler, 
    143 or *SASCOMPILER=mingw* for the MinGW compiler. If TinyCC is available 
     182On Windows, set *SAS_COMPILER=tinycc* for the tinycc compiler, 
     183*SAS_COMPILER=msvc* for the Microsoft Visual C compiler, 
     184or *SAS_COMPILER=mingw* for the MinGW compiler. If TinyCC is available 
    144185on the python path (it is provided with SasView), that will be the 
    145186default. If you want one of the other compilers, be sure to have it 
     
    154195*Document History* 
    155196 
    156 | 2017-09-27 Paul Kienzle 
     197| 2018-10-15 Paul Kienzle 
  • doc/guide/magnetism/magnetism.rst

    rbefe905 rdf87acf  
    8989 
    9090===========   ================================================================ 
    91  M0:sld       $D_M M_0$ 
    92  mtheta:sld   $\theta_M$ 
    93  mphi:sld     $\phi_M$ 
    94  up:angle     $\theta_\mathrm{up}$ 
    95  up:frac_i    $u_i$ = (spin up)/(spin up + spin down) *before* the sample 
    96  up:frac_f    $u_f$ = (spin up)/(spin up + spin down) *after* the sample 
     91 sld_M0       $D_M M_0$ 
     92 sld_mtheta   $\theta_M$ 
     93 sld_mphi     $\phi_M$ 
     94 up_frac_i    $u_i$ = (spin up)/(spin up + spin down) *before* the sample 
     95 up_frac_f    $u_f$ = (spin up)/(spin up + spin down) *after* the sample 
     96 up_angle     $\theta_\mathrm{up}$ 
    9797===========   ================================================================ 
    9898 
    9999.. note:: 
    100     The values of the 'up:frac_i' and 'up:frac_f' must be in the range 0 to 1. 
     100    The values of the 'up_frac_i' and 'up_frac_f' must be in the range 0 to 1. 
    101101 
    102102*Document History* 
  • doc/guide/pd/polydispersity.rst

    r29afc50 rd089a00  
    88.. _polydispersityhelp: 
    99 
    10 Polydispersity Distributions 
    11 ---------------------------- 
    12  
    13 With some models in sasmodels we can calculate the average intensity for a 
    14 population of particles that exhibit size and/or orientational 
    15 polydispersity. The resultant intensity is normalized by the average 
    16 particle volume such that 
     10Polydispersity & Orientational Distributions 
     11-------------------------------------------- 
     12 
     13For some models we can calculate the average intensity for a population of  
     14particles that possess size and/or orientational (ie, angular) distributions.  
     15In SasView we call the former *polydispersity* but use the parameter *PD* to  
     16parameterise both. In other words, the meaning of *PD* in a model depends on  
     17the actual parameter it is being applied too. 
     18 
     19The resultant intensity is then normalized by the average particle volume such  
     20that 
    1721 
    1822.. math:: 
     
    2024  P(q) = \text{scale} \langle F^* F \rangle / V + \text{background} 
    2125 
    22 where $F$ is the scattering amplitude and $\langle\cdot\rangle$ denotes an 
    23 average over the size distribution. 
     26where $F$ is the scattering amplitude and $\langle\cdot\rangle$ denotes an  
     27average over the distribution $f(x; \bar x, \sigma)$, giving 
     28 
     29.. math:: 
     30 
     31  P(q) = \frac{\text{scale}}{V} \int_\mathbb{R}  
     32  f(x; \bar x, \sigma) F^2(q, x)\, dx + \text{background} 
    2433 
    2534Each distribution is characterized by a center value $\bar x$ or 
    2635$x_\text{med}$, a width parameter $\sigma$ (note this is *not necessarily* 
    27 the standard deviation, so read the description carefully), the number of 
    28 sigmas $N_\sigma$ to include from the tails of the distribution, and the 
    29 number of points used to compute the average. The center of the distribution 
    30 is set by the value of the model parameter. The meaning of a polydispersity  
    31 parameter *PD* (not to be confused with a molecular weight distributions  
    32 in polymer science) in a model depends on the type of parameter it is being  
    33 applied too. 
     36the standard deviation, so read the description of the distribution carefully),  
     37the number of sigmas $N_\sigma$ to include from the tails of the distribution,  
     38and the number of points used to compute the average. The center of the  
     39distribution is set by the value of the model parameter. 
    3440 
    3541The distribution width applied to *volume* (ie, shape-describing) parameters  
    3642is relative to the center value such that $\sigma = \mathrm{PD} \cdot \bar x$.  
    37 However, the distribution width applied to *orientation* (ie, angle-describing)  
    38 parameters is just $\sigma = \mathrm{PD}$. 
     43However, the distribution width applied to *orientation* parameters is just  
     44$\sigma = \mathrm{PD}$. 
    3945 
    4046$N_\sigma$ determines how far into the tails to evaluate the distribution, 
    4147with larger values of $N_\sigma$ required for heavier tailed distributions. 
    4248The scattering in general falls rapidly with $qr$ so the usual assumption 
    43 that $G(r - 3\sigma_r)$ is tiny and therefore $f(r - 3\sigma_r)G(r - 3\sigma_r)$ 
     49that $f(r - 3\sigma_r)$ is tiny and therefore $f(r - 3\sigma_r)f(r - 3\sigma_r)$ 
    4450will not contribute much to the average may not hold when particles are large. 
    4551This, too, will require increasing $N_\sigma$. 
    4652 
    4753Users should note that the averaging computation is very intensive. Applying 
    48 polydispersion to multiple parameters at the same time or increasing the 
    49 number of points in the distribution will require patience! However, the 
    50 calculations are generally more robust with more data points or more angles. 
     54polydispersion and/or orientational distributions to multiple parameters at  
     55the same time, or increasing the number of points in the distribution, will  
     56require patience! However, the calculations are generally more robust with  
     57more data points or more angles. 
    5158 
    5259The following distribution functions are provided: 
     
    6471Additional distributions are under consideration. 
    6572 
     73**Beware: when the Polydispersity & Orientational Distribution panel in SasView is** 
     74**first opened, the default distribution for all parameters is the Gaussian Distribution.** 
     75**This may not be suitable. See Suggested Applications below.** 
     76 
     77.. note:: In 2009 IUPAC decided to introduce the new term 'dispersity' to replace  
     78           the term 'polydispersity' (see `Pure Appl. Chem., (2009), 81(2),  
     79           351-353 <http://media.iupac.org/publications/pac/2009/pdf/8102x0351.pdf>`_  
     80           in order to make the terminology describing distributions of chemical  
     81           properties unambiguous. However, these terms are unrelated to the  
     82           proportional size distributions and orientational distributions used in  
     83           SasView models. 
     84 
    6685Suggested Applications 
    6786^^^^^^^^^^^^^^^^^^^^^^ 
    6887 
    69 If applying polydispersion to parameters describing particle sizes, use 
     88If applying polydispersion to parameters describing particle sizes, consider using 
    7089the Lognormal or Schulz distributions. 
    7190 
    7291If applying polydispersion to parameters describing interfacial thicknesses 
    73 or angular orientations, use the Gaussian or Boltzmann distributions. 
     92or angular orientations, consider using the Gaussian or Boltzmann distributions. 
    7493 
    7594If applying polydispersion to parameters describing angles, use the Uniform  
     
    318337^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 
    319338 
    320 Many commercial Dynamic Light Scattering (DLS) instruments produce a size 
    321 polydispersity parameter, sometimes even given the symbol $p$\ ! This 
    322 parameter is defined as the relative standard deviation coefficient of 
    323 variation of the size distribution and is NOT the same as the polydispersity 
    324 parameters in the Lognormal and Schulz distributions above (though they all 
    325 related) except when the DLS polydispersity parameter is <0.13. 
    326  
    327 .. math:: 
    328  
    329     p_{DLS} = \sqrt(\nu / \bar x^2) 
    330  
    331 where $\nu$ is the variance of the distribution and $\bar x$ is the mean 
    332 value of $x$. 
     339Several measures of polydispersity abound in Dynamic Light Scattering (DLS) and  
     340it should not be assumed that any of the following can be simply equated with  
     341the polydispersity *PD* parameter used in SasView. 
     342 
     343The dimensionless **Polydispersity Index (PI)** is a measure of the width of the  
     344distribution of autocorrelation function decay rates (*not* the distribution of  
     345particle sizes itself, though the two are inversely related) and is defined by  
     346ISO 22412:2017 as 
     347 
     348.. math:: 
     349 
     350    PI = \mu_{2} / \bar \Gamma^2 
     351 
     352where $\mu_\text{2}$ is the second cumulant, and $\bar \Gamma^2$ is the  
     353intensity-weighted average value, of the distribution of decay rates. 
     354 
     355*If the distribution of decay rates is Gaussian* then 
     356 
     357.. math:: 
     358 
     359    PI = \sigma^2 / 2\bar \Gamma^2 
     360 
     361where $\sigma$ is the standard deviation, allowing a **Relative Polydispersity (RP)**  
     362to be defined as 
     363 
     364.. math:: 
     365 
     366    RP = \sigma / \bar \Gamma = \sqrt{2 \cdot PI} 
     367 
     368PI values smaller than 0.05 indicate a highly monodisperse system. Values  
     369greater than 0.7 indicate significant polydispersity. 
     370 
     371The **size polydispersity P-parameter** is defined as the relative standard  
     372deviation coefficient of variation   
     373 
     374.. math:: 
     375 
     376    P = \sqrt\nu / \bar R 
     377 
     378where $\nu$ is the variance of the distribution and $\bar R$ is the mean 
     379value of $R$. Here, the product $P \bar R$ is *equal* to the standard  
     380deviation of the Lognormal distribution. 
     381 
     382P values smaller than 0.13 indicate a monodisperse system. 
    333383 
    334384For more information see: 
    335 S King, C Washington & R Heenan, *Phys Chem Chem Phys*, (2005), 7, 143 
     385 
     386`ISO 22412:2017, International Standards Organisation (2017) <https://www.iso.org/standard/65410.html>`_. 
     387 
     388`Polydispersity: What does it mean for DLS and Chromatography <http://www.materials-talks.com/blog/2014/10/23/polydispersity-what-does-it-mean-for-dls-and-chromatography/>`_. 
     389 
     390`Dynamic Light Scattering: Common Terms Defined, Whitepaper WP111214. Malvern Instruments (2011) <http://www.biophysics.bioc.cam.ac.uk/wp-content/uploads/2011/02/DLS_Terms_defined_Malvern.pdf>`_. 
     391 
     392S King, C Washington & R Heenan, *Phys Chem Chem Phys*, (2005), 7, 143. 
     393 
     394T Allen, in *Particle Size Measurement*, 4th Edition, Chapman & Hall, London (1990). 
    336395 
    337396.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
     
    343402| 2018-03-20 Steve King 
    344403| 2018-04-04 Steve King 
     404| 2018-08-09 Steve King 
  • doc/guide/plugin.rst

    r7e6bc45e rc94ab04  
    273273`Form_Factors`_ for more details. 
    274274 
     275**model_info = ...** lets you define a model directly, for example, by 
     276loading and modifying existing models.  This is done implicitly by 
     277:func:`sasmodels.core.load_model_info`, which can create a mixture model 
     278from a pair of existing models.  For example:: 
     279 
     280    from sasmodels.core import load_model_info 
     281    model_info = load_model_info('sphere+cylinder') 
     282 
     283See :class:`sasmodels.modelinfo.ModelInfo` for details about the model 
     284attributes that are defined. 
     285 
    275286Model Parameters 
    276287................ 
     
    291302 
    292303**Note: The order of the parameters in the definition will be the order of the 
    293 parameters in the user interface and the order of the parameters in Iq(), 
    294 Iqac(), Iqabc() and form_volume(). And** *scale* **and** *background* 
    295 **parameters are implicit to all models, so they do not need to be included 
    296 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.** 
    297308 
    298309- **"name"** is the name of the parameter shown on the FitPage. 
     
    363374    scattered intensity. 
    364375 
    365   - "volume" parameters are passed to Iq(), Iqac(), Iqabc() and form_volume(), 
    366     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. 
    367378 
    368379  - "orientation" parameters are not passed, but instead are combined with 
     
    424435appropriately smeared pattern. 
    425436 
     437Each .py file also contains a function:: 
     438 
     439        def random(): 
     440        ... 
     441 
     442This function provides a model-specific random parameter set which shows model 
     443features in the USANS to SANS range.  For example, core-shell sphere sets the 
     444outer radius of the sphere logarithmically in `[20, 20,000]`, which sets the Q 
     445value for the transition from flat to falling.  It then uses a beta distribution 
     446to set the percentage of the shape which is shell, giving a preference for very 
     447thin or very thick shells (but never 0% or 100%).  Using `-sets=10` in sascomp 
     448should show a reasonable variety of curves over the default sascomp q range. 
     449The parameter set is returned as a dictionary of `{parameter: value, ...}`. 
     450Any model parameters not included in the dictionary will default according to 
     451the code in the `_randomize_one()` function from sasmodels/compare.py. 
     452 
    426453Python Models 
    427454............. 
     
    476503used. 
    477504 
     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 
     514**Note: Pure python models do not yet support direct computation of the** 
     515**average of $F(q)$ and $F^2(q)$. Neither do they support orientational** 
     516**distributions or magnetism (use C models if these are required).** 
     517 
    478518Embedded C Models 
    479519................. 
     
    487527This expands into the equivalent C code:: 
    488528 
    489     #include <math.h> 
    490529    double Iq(double q, double par1, double par2, ...); 
    491530    double Iq(double q, double par1, double par2, ...) 
     
    496535*form_volume* defines the volume of the shape. As in python models, it 
    497536includes only the volume parameters. 
     537 
     538*form_volume* defines the volume of the shell for hollow shapes. As in 
     539python models, it includes only the volume parameters. 
    498540 
    499541**source=['fn.c', ...]** includes the listed C source files in the 
     
    533575listed in *source*. 
    534576 
     577Structure Factors 
     578................. 
     579 
     580Structure factor calculations may need the underlying $<F(q)>$ and $<F^2(q)>$ 
     581rather than $I(q)$.  This is used to compute $\beta = <F(q)>^2/<F^2(q)>$ in 
     582the decoupling approximation to the structure factor. 
     583 
     584Instead of defining the *Iq* function, models can define *Fq* as 
     585something like:: 
     586 
     587    double Fq(double q, double *F1, double *F2, double par1, double par2, ...); 
     588    double Fq(double q, double *F1, double *F2, double par1, double par2, ...) 
     589    { 
     590        // Polar integration loop over all orientations. 
     591        ... 
     592        *F1 = 1e-2 * total_F1 * contrast * volume; 
     593        *F2 = 1e-4 * total_F2 * square(contrast * volume); 
     594        return I(q, par1, par2, ...); 
     595    } 
     596 
     597If the volume fraction scale factor is built into the model (as occurs for 
     598the vesicle model, for example), then scale *F1* by $\surd V_f$ so that 
     599$\beta$ is computed correctly. 
     600 
     601Structure factor calculations are not yet supported for oriented shapes. 
     602 
     603Note: only available as a separate C file listed in *source*, or within 
     604a *c_code* block within the python model definition file. 
     605 
    535606Oriented Shapes 
    536607............... 
     
    544615laboratory frame and beam travelling along $-z$. 
    545616 
    546 The oriented C model is called using *Iqabc(qa, qb, qc, par1, par2, ...)* where 
     617The oriented C model (oriented pure Python models are not supported)  
     618is called using *Iqabc(qa, qb, qc, par1, par2, ...)* where 
    547619*par1*, etc. are the parameters to the model.  If the shape is rotationally 
    548620symmetric about *c* then *psi* is not needed, and the model is called 
     
    648720to compute the proper magnetism and orientation, which you can implement 
    649721using *Iqxy(qx, qy, par1, par2, ...)*. 
     722 
     723**Note: Magnetism is not supported in pure Python models.** 
    650724 
    651725Special Functions 
     
    685759    erf, erfc, tgamma, lgamma:  **do not use** 
    686760        Special functions that should be part of the standard, but are missing 
    687         or inaccurate on some platforms. Use sas_erf, sas_erfc and sas_gamma 
    688         instead (see below). Note: lgamma(x) has not yet been tested. 
     761        or inaccurate on some platforms. Use sas_erf, sas_erfc, sas_gamma 
     762        and sas_lgamma instead (see below). 
    689763 
    690764Some non-standard constants and functions are also provided: 
     
    753827        Gamma function sas_gamma\ $(x) = \Gamma(x)$. 
    754828 
    755         The standard math function, tgamma(x) is unstable for $x < 1$ 
     829        The standard math function, tgamma(x), is unstable for $x < 1$ 
    756830        on some platforms. 
    757831 
    758832        :code:`source = ["lib/sas_gamma.c", ...]` 
    759833        (`sas_gamma.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gamma.c>`_) 
     834 
     835    sas_gammaln(x): 
     836        log gamma function sas_gammaln\ $(x) = \log \Gamma(|x|)$. 
     837 
     838        The standard math function, lgamma(x), is incorrect for single 
     839        precision on some platforms. 
     840 
     841        :code:`source = ["lib/sas_gammainc.c", ...]` 
     842        (`sas_gammainc.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gammainc.c>`_) 
     843 
     844    sas_gammainc(a, x), sas_gammaincc(a, x): 
     845        Incomplete gamma function 
     846        sas_gammainc\ $(a, x) = \int_0^x t^{a-1}e^{-t}\,dt / \Gamma(a)$ 
     847        and complementary incomplete gamma function 
     848        sas_gammaincc\ $(a, x) = \int_x^\infty t^{a-1}e^{-t}\,dt / \Gamma(a)$ 
     849 
     850        :code:`source = ["lib/sas_gammainc.c", ...]` 
     851        (`sas_gammainc.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gammainc.c>`_) 
    760852 
    761853    sas_erf(x), sas_erfc(x): 
     
    795887        If $n$ = 0 or 1, it uses sas_J0($x$) or sas_J1($x$), respectively. 
    796888 
     889        Warning: JN(n,x) can be very inaccurate (0.1%) for x not in [0.1, 100]. 
     890 
    797891        The standard math function jn(n, x) is not available on all platforms. 
    798892 
     
    803897        Sine integral Si\ $(x) = \int_0^x \tfrac{\sin t}{t}\,dt$. 
    804898 
     899        Warning: Si(x) can be very inaccurate (0.1%) for x in [0.1, 100]. 
     900 
    805901        This function uses Taylor series for small and large arguments: 
    806902 
    807         For large arguments, 
     903        For large arguments use the following Taylor series, 
    808904 
    809905        .. math:: 
     
    822918 
    823919        :code:`source = ["lib/Si.c", ...]` 
    824         (`Si.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/Si.c>`_) 
     920        (`Si.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_Si.c>`_) 
    825921 
    826922    sas_3j1x_x(x): 
     
    9741070          "radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, 
    9751071         0.2, 0.228843], 
    976         [{"radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, "ER", 120.], 
    977         [{"radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, "VR", 1.], 
     1072        [{"radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, 
     1073         0.1, None, None, 120., None, 1.],  # q, F, F^2, R_eff, V, form:shell 
     1074        [{"@S": "hardsphere"}, 0.1, None], 
    9781075    ] 
    9791076 
    9801077 
    981 **tests=[[{parameters}, q, result], ...]** is a list of lists. 
     1078**tests=[[{parameters}, q, Iq], ...]** is a list of lists. 
    9821079Each list is one test and contains, in order: 
    9831080 
     
    9911088- input and output values can themselves be lists if you have several 
    9921089  $q$ values to test for the same model parameters. 
    993 - for testing *ER* and *VR*, give the inputs as "ER" and "VR" respectively; 
    994   the output for *VR* should be the sphere/shell ratio, not the individual 
    995   sphere and shell values. 
     1090- for testing effective radius, volume and form:shell volume ratio, use the 
     1091  extended form of the tests results, with *None, None, R_eff, V, V_r* 
     1092  instead of *Iq*.  This calls the kernel *Fq* function instead of *Iq*. 
     1093- for testing F and F^2 (used for beta approximation) do the same as the 
     1094  effective radius test, but include values for the first two elements, 
     1095  $<F(q)>$ and $<F^2(q)>$. 
     1096- for testing interaction between form factor and structure factor, specify 
     1097  the structure factor name in the parameters as *{"@S": "name", ...}* with 
     1098  the remaining list of parameters defined by the *P@S* product model. 
    9961099 
    9971100.. _Test_Your_New_Model: 
     
    10091112and a check that the model runs. 
    10101113 
    1011 If you are not using sasmodels from SasView, skip this step. 
    1012  
    10131114Recommended Testing 
    10141115................... 
     1116 
     1117**NB: For now, this more detailed testing is only possible if you have a  
     1118SasView build environment available!** 
    10151119 
    10161120If the model compiles and runs, you can next run the unit tests that 
  • doc/guide/scripting.rst

    r4aa5dce r23df833  
    1010The key functions are :func:`sasmodels.core.load_model` for loading the 
    1111model definition and compiling the kernel and 
    12 :func:`sasmodels.data.load_data` for calling sasview to load the data. Need 
    13 the data because that defines the resolution function and the q values to 
    14 evaluate. If there is no data, then use :func:`sasmodels.data.empty_data1D` 
    15 or :func:`sasmodels.data.empty_data2D` to create some data with a given $q$. 
    16  
    17 Using sasmodels through bumps 
    18 ============================= 
    19  
    20 With the data and the model, you can wrap it in a *bumps* model with 
     12:func:`sasmodels.data.load_data` for calling sasview to load the data. 
     13 
     14Preparing data 
     15============== 
     16 
     17Usually you will load data via the sasview loader, with the 
     18:func:`sasmodels.data.load_data` function.  For example:: 
     19 
     20    from sasmodels.data import load_data 
     21    data = load_data("sasmodels/example/093191_201.dat") 
     22 
     23You may want to apply a data mask, such a beam stop, and trim high $q$:: 
     24 
     25    from sasmodels.data import set_beam_stop 
     26    set_beam_stop(data, qmin, qmax) 
     27 
     28The :func:`sasmodels.data.set_beam_stop` method simply sets the *mask* 
     29attribute for the data. 
     30 
     31The data defines the resolution function and the q values to evaluate, so 
     32even if you simulating experiments prior to making measurements, you still 
     33need a data object for reference. Use :func:`sasmodels.data.empty_data1D` 
     34or :func:`sasmodels.data.empty_data2D` to create a container with a 
     35given $q$ and $\Delta q/q$.  For example:: 
     36 
     37    import numpy as np 
     38    from sasmodels.data import empty_data1D 
     39 
     40    # 120 points logarithmically spaced from 0.005 to 0.2, with dq/q = 5% 
     41    q = np.logspace(np.log10(5e-3), np.log10(2e-1), 120) 
     42    data = empty_data1D(q, resolution=0.05) 
     43 
     44To use a more realistic model of resolution, or to load data from a file 
     45format not understood by SasView, you can use :class:`sasmodels.data.Data1D` 
     46or :class:`sasmodels.data.Data2D` directly.  The 1D data uses 
     47*x*, *y*, *dx* and *dy* for $x = q$ and $y = I(q)$, and 2D data uses 
     48*x*, *y*, *z*, *dx*, *dy*, *dz* for $x, y = qx, qy$ and $z = I(qx, qy)$. 
     49[Note: internally, the Data2D object uses SasView conventions, 
     50*qx_data*, *qy_data*, *data*, *dqx_data*, *dqy_data*, and *err_data*.] 
     51 
     52For USANS data, use 1D data, but set *dxl* and *dxw* attributes to 
     53indicate slit resolution:: 
     54 
     55    data.dxl = 0.117 
     56 
     57See :func:`sasmodels.resolution.slit_resolution` for details. 
     58 
     59SESANS data is more complicated; if your SESANS format is not supported by 
     60SasView you need to define a number of attributes beyond *x*, *y*.  For 
     61example:: 
     62 
     63    SElength = np.linspace(0, 2400, 61) # [A] 
     64    data = np.ones_like(SElength) 
     65    err_data = np.ones_like(SElength)*0.03 
     66 
     67    class Source: 
     68        wavelength = 6 # [A] 
     69        wavelength_unit = "A" 
     70    class Sample: 
     71        zacceptance = 0.1 # [A^-1] 
     72        thickness = 0.2 # [cm] 
     73 
     74    class SESANSData1D: 
     75        #q_zmax = 0.23 # [A^-1] 
     76        lam = 0.2 # [nm] 
     77        x = SElength 
     78        y = data 
     79        dy = err_data 
     80        sample = Sample() 
     81    data = SESANSData1D() 
     82 
     83    x, y = ... # create or load sesans 
     84    data = smd.Data 
     85 
     86The *data* module defines various data plotters as well. 
     87 
     88Using sasmodels directly 
     89======================== 
     90 
     91Once you have a computational kernel and a data object, you can evaluate 
     92the model for various parameters using 
     93:class:`sasmodels.direct_model.DirectModel`.  The resulting object *f* 
     94will be callable as *f(par=value, ...)*, returning the $I(q)$ for the $q$ 
     95values in the data.  For example:: 
     96 
     97    import numpy as np 
     98    from sasmodels.data import empty_data1D 
     99    from sasmodels.core import load_model 
     100    from sasmodels.direct_model import DirectModel 
     101 
     102    # 120 points logarithmically spaced from 0.005 to 0.2, with dq/q = 5% 
     103    q = np.logspace(np.log10(5e-3), np.log10(2e-1), 120) 
     104    data = empty_data1D(q, resolution=0.05) 
     105    kernel = load_model("ellipsoid) 
     106    f = DirectModel(data, kernel) 
     107    Iq = f(radius_polar=100) 
     108 
     109Polydispersity information is set with special parameter names: 
     110 
     111    * *par_pd* for polydispersity width, $\Delta p/p$, 
     112    * *par_pd_n* for the number of points in the distribution, 
     113    * *par_pd_type* for the distribution type (as a string), and 
     114    * *par_pd_nsigmas* for the limits of the distribution. 
     115 
     116Using sasmodels through the bumps optimizer 
     117=========================================== 
     118 
     119Like DirectModel, you can wrap data and a kernel in a *bumps* model with 
    21120class:`sasmodels.bumps_model.Model` and create an 
    22 class:`sasmodels.bump_model.Experiment` that you can fit with the *bumps* 
     121class:`sasmodels.bumps_model.Experiment` that you can fit with the *bumps* 
    23122interface. Here is an example from the *example* directory such as 
    24123*example/model.py*:: 
     
    75174    SasViewCom bumps.cli example/model.py --preview 
    76175 
    77 Using sasmodels directly 
    78 ======================== 
    79  
    80 Bumps has a notion of parameter boxes in which you can set and retrieve 
    81 values.  Instead of using bumps, you can create a directly callable function 
    82 with :class:`sasmodels.direct_model.DirectModel`.  The resulting object *f* 
    83 will be callable as *f(par=value, ...)*, returning the $I(q)$ for the $q$ 
    84 values in the data.  Polydisperse parameters use the same naming conventions 
    85 as in the bumps model, with e.g., radius_pd being the polydispersity associated 
    86 with radius. 
     176Calling the computation kernel 
     177============================== 
    87178 
    88179Getting a simple function that you can call on a set of q values and return 
     
    97188python kernel.  Once the kernel is in hand, we can then marshal a set of 
    98189parameters into a :class:`sasmodels.details.CallDetails` object and ship it to 
    99 the kernel using the :func:`sansmodels.direct_model.call_kernel` function.  An 
    100 example should help, *example/cylinder_eval.py*:: 
    101  
    102     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 
    103198    from matplotlib import pyplot as plt 
    104199    from sasmodels.core import load_model 
    105     from sasmodels.direct_model import call_kernel 
     200    from sasmodels.direct_model import call_kernel, call_Fq 
    106201 
    107202    model = load_model('cylinder') 
    108203    q = logspace(-3, -1, 200) 
    109204    kernel = model.make_kernel([q]) 
    110     Iq = call_kernel(kernel, dict(radius=200.)) 
    111     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() 
    112216    plt.show() 
    113217 
    114 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:: 
    115233 
    116234    SasViewCom example/cylinder_eval.py 
  • doc/rst_prolog

    r30b60d2 r2c12061  
    99.. |Ang^-3| replace:: |Ang|\ :sup:`-3` 
    1010.. |Ang^-4| replace:: |Ang|\ :sup:`-4` 
     11.. |nm^-1| replace:: nm\ :sup:`-1` 
    1112.. |cm^-1| replace:: cm\ :sup:`-1` 
    1213.. |cm^2| replace:: cm\ :sup:`2` 
  • 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() 
  • example/model_ellipsoid_hayter_msa.py

    r8a5f021 r9a99993  
    1616 
    1717# DEFINE THE MODEL 
    18 kernel = load_model('ellipsoid*hayter_msa') 
     18kernel = load_model('ellipsoid@hayter_msa') 
    1919 
    2020pars = dict(scale=6.4, background=0.06, sld=0.33, sld_solvent=2.15, radius_polar=14.0, 
  • example/multiscatfit.py

    r49d1f8b8 r2c4a190  
    1515 
    1616    # Show the model without fitting 
    17     PYTHONPATH=..:../explore:../../bumps:../../sasview/src python multiscatfit.py 
     17    PYTHONPATH=..:../../bumps:../../sasview/src python multiscatfit.py 
    1818 
    1919    # Run the fit 
    20     PYTHONPATH=..:../explore:../../bumps:../../sasview/src ../../bumps/run.py \ 
     20    PYTHONPATH=..:../../bumps:../../sasview/src ../../bumps/run.py \ 
    2121    multiscatfit.py --store=/tmp/t1 
    2222 
     
    5555    ) 
    5656 
     57# Tie the model to the data 
     58M = Experiment(data=data, model=model) 
     59 
     60# Stack mulitple scattering on top of the existing resolution function. 
     61M.resolution = MultipleScattering(resolution=M.resolution, probability=0.) 
     62 
    5763# SET THE FITTING PARAMETERS 
    5864model.radius_polar.range(15, 3000) 
     
    6571model.scale.range(0, 0.1) 
    6672 
    67 # Mulitple scattering probability parameter 
    68 # HACK: the probability is stuffed in as an extra parameter to the experiment. 
    69 probability = Parameter(name="probability", value=0.0) 
    70 probability.range(0.0, 0.9) 
     73# The multiple scattering probability parameter is in the resolution function 
     74# instead of the scattering function, so access it through M.resolution 
     75M.scattering_probability.range(0.0, 0.9) 
    7176 
    72 M = Experiment(data=data, model=model, extra_pars={'probability': probability}) 
    73  
    74 # Stack mulitple scattering on top of the existing resolution function. 
    75 # Because resolution functions in sasview don't have fitting parameters, 
    76 # we instead allow the multiple scattering calculator to take a function 
    77 # instead of a probability.  This function returns the current value of 
    78 # the parameter. ** THIS IS TEMPORARY ** when multiple scattering is 
    79 # properly integrated into sasmodels and sasview, its fittable parameter 
    80 # will be treated like the model parameters. 
    81 M.resolution = MultipleScattering(resolution=M.resolution, 
    82                                   probability=lambda: probability.value, 
    83                                   ) 
    84 M._kernel_inputs = M.resolution.q_calc 
     77# Let bumps know that we are fitting this experiment 
    8578problem = FitProblem(M) 
    8679 
  • 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

    r2a7e20e rcd28947  
    9595            neg:    [-100,100] 
    9696 
     97        For arbitrary range use "start:stop:steps:scale" where scale is 
     98        one of log, lin, or linear. 
     99 
    97100        *diff* is "relative", "absolute" or "none" 
    98101 
     
    102105        linear = not xrange.startswith("log") 
    103106        if xrange == "zoom": 
    104             lin_min, lin_max, lin_steps = 1000, 1010, 2000 
     107            start, stop, steps = 1000, 1010, 2000 
    105108        elif xrange == "neg": 
    106             lin_min, lin_max, lin_steps = -100.1, 100.1, 2000 
     109            start, stop, steps = -100.1, 100.1, 2000 
    107110        elif xrange == "linear": 
    108             lin_min, lin_max, lin_steps = 1, 1000, 2000 
    109             lin_min, lin_max, lin_steps = 0.001, 2, 2000 
     111            start, stop, steps = 1, 1000, 2000 
     112            start, stop, steps = 0.001, 2, 2000 
    110113        elif xrange == "log": 
    111             log_min, log_max, log_steps = -3, 5, 400 
     114            start, stop, steps = -3, 5, 400 
    112115        elif xrange == "logq": 
    113             log_min, log_max, log_steps = -4, 1, 400 
     116            start, stop, steps = -4, 1, 400 
     117        elif ':' in xrange: 
     118            parts = xrange.split(':') 
     119            linear = parts[3] != "log" if len(parts) == 4 else True 
     120            steps = int(parts[2]) if len(parts) > 2 else 400 
     121            start = float(parts[0]) 
     122            stop = float(parts[1]) 
     123 
    114124        else: 
    115125            raise ValueError("unknown range "+xrange) 
     
    121131            # value to x in the given precision. 
    122132            if linear: 
    123                 lin_min = max(lin_min, self.limits[0]) 
    124                 lin_max = min(lin_max, self.limits[1]) 
    125                 qrf = np.linspace(lin_min, lin_max, lin_steps, dtype='single') 
    126                 #qrf = np.linspace(lin_min, lin_max, lin_steps, dtype='double') 
     133                start = max(start, self.limits[0]) 
     134                stop = min(stop, self.limits[1]) 
     135                qrf = np.linspace(start, stop, steps, dtype='single') 
     136                #qrf = np.linspace(start, stop, steps, dtype='double') 
    127137                qr = [mp.mpf(float(v)) for v in qrf] 
    128                 #qr = mp.linspace(lin_min, lin_max, lin_steps) 
     138                #qr = mp.linspace(start, stop, steps) 
    129139            else: 
    130                 log_min = np.log10(max(10**log_min, self.limits[0])) 
    131                 log_max = np.log10(min(10**log_max, self.limits[1])) 
    132                 qrf = np.logspace(log_min, log_max, log_steps, dtype='single') 
    133                 #qrf = np.logspace(log_min, log_max, log_steps, dtype='double') 
     140                start = np.log10(max(10**start, self.limits[0])) 
     141                stop = np.log10(min(10**stop, self.limits[1])) 
     142                qrf = np.logspace(start, stop, steps, dtype='single') 
     143                #qrf = np.logspace(start, stop, steps, dtype='double') 
    134144                qr = [mp.mpf(float(v)) for v in qrf] 
    135                 #qr = [10**v for v in mp.linspace(log_min, log_max, log_steps)] 
     145                #qr = [10**v for v in mp.linspace(start, stop, steps)] 
    136146 
    137147        target = self.call_mpmath(qr, bits=500) 
     
    176186    """ 
    177187    if diff == "relative": 
    178         err = np.array([abs((t-a)/t) for t, a in zip(target, actual)], 'd') 
     188        err = np.array([(abs((t-a)/t) if t != 0 else a) for t, a in zip(target, actual)], 'd') 
    179189        #err = np.clip(err, 0, 1) 
    180190        pylab.loglog(x, err, '-', label=label) 
     
    197207    return model_info 
    198208 
     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. 
     212A = 1 
     213def parse_extra_pars(): 
     214    """ 
     215    Parse the command line looking for the second parameter "A=..." for the 
     216    gammainc/gammaincc functions. 
     217    """ 
     218    global A 
     219 
     220    A_str = str(A) 
     221    pop = [] 
     222    for k, v in enumerate(sys.argv[1:]): 
     223        if v.startswith("A="): 
     224            A_str = v[2:] 
     225            pop.append(k+1) 
     226    if pop: 
     227        sys.argv = [v for k, v in enumerate(sys.argv) if k not in pop] 
     228        A = float(A_str) 
     229 
     230parse_extra_pars() 
     231 
    199232 
    200233# =============== FUNCTION DEFINITIONS ================ 
     
    299332) 
    300333add_function( 
     334    name="gammaln(x)", 
     335    mp_function=mp.loggamma, 
     336    np_function=scipy.special.gammaln, 
     337    ocl_function=make_ocl("return sas_gammaln(q);", "sas_gammaln", ["lib/sas_gammainc.c"]), 
     338    #ocl_function=make_ocl("return lgamma(q);", "sas_gammaln"), 
     339) 
     340add_function( 
     341    # Note: "a" is given as A=... on the command line via parse_extra_pars 
     342    name="gammainc(x)", 
     343    mp_function=lambda x, a=A: mp.gammainc(a, a=0, b=x)/mp.gamma(a), 
     344    np_function=lambda x, a=A: scipy.special.gammainc(a, x), 
     345    ocl_function=make_ocl("return sas_gammainc(%.15g,q);"%A, "sas_gammainc", ["lib/sas_gammainc.c"]), 
     346) 
     347add_function( 
     348    # Note: "a" is given as A=... on the command line via parse_extra_pars 
     349    name="gammaincc(x)", 
     350    mp_function=lambda x, a=A: mp.gammainc(a, a=x, b=mp.inf)/mp.gamma(a), 
     351    np_function=lambda x, a=A: scipy.special.gammaincc(a, x), 
     352    ocl_function=make_ocl("return sas_gammaincc(%.15g,q);"%A, "sas_gammaincc", ["lib/sas_gammainc.c"]), 
     353) 
     354add_function( 
    301355    name="erf(x)", 
    302356    mp_function=mp.erf, 
     
    357411) 
    358412add_function( 
    359     name="debye", 
     413    name="gauss_coil", 
    360414    mp_function=lambda x: 2*(mp.exp(-x**2) + x**2 - 1)/x**4, 
    361415    np_function=lambda x: 2*(np.expm1(-x**2) + x**2)/x**4, 
    362416    ocl_function=make_ocl(""" 
    363417    const double qsq = q*q; 
    364     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) { 
    365421        const double x = qsq; 
    366422        if (0) { // 0.36 single 
     
    372428            const double B1=3./8., B2=3./56., B3=1./336.; 
    373429            return (((A3*x + A2)*x + A1)*x + 1.)/(((B3*x + B2)*x + B1)*x + 1.); 
    374         } else if (1) { // 1.0 for single, 0.25 for double 
     430        } else if (0) { // 1.0 for single, 0.25 for double 
    375431            // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}] 
    376432            const double A1=1./15., A2=1./60, A3=0., A4=1./75600.; 
     
    385441                  /(((((B5*x + B4)*x + B3)*x + B2)*x + B1)*x + 1.); 
    386442        } 
    387     } else if (qsq < 1.) { // Taylor series; 0.9 for single, 0.25 for double 
     443    } else if (qsq < 0.8) { 
    388444        const double x = qsq; 
    389445        const double C0 = +1.; 
     
    463519lanczos_gamma = """\ 
    464520    const double coeff[] = { 
    465             76.18009172947146,     -86.50532032941677, 
    466             24.01409824083091,     -1.231739572450155, 
     521            76.18009172947146, -86.50532032941677, 
     522            24.01409824083091, -1.231739572450155, 
    467523            0.1208650973866179e-2,-0.5395239384953e-5 
    468524            }; 
     
    475531""" 
    476532add_function( 
    477     name="log gamma(x)", 
     533    name="loggamma(x)", 
    478534    mp_function=mp.loggamma, 
    479535    np_function=scipy.special.gammaln, 
     
    599655 
    600656ALL_FUNCTIONS = set(FUNCTIONS.keys()) 
    601 ALL_FUNCTIONS.discard("loggamma")  # OCL version not ready yet 
     657ALL_FUNCTIONS.discard("loggamma")  # use cephes-based gammaln instead 
    602658ALL_FUNCTIONS.discard("3j1/x:taylor") 
    603659ALL_FUNCTIONS.discard("3j1/x:trig") 
     
    615671    -r indicates that the relative error should be plotted (default), 
    616672    -x<range> indicates the steps in x, where <range> is one of the following 
    617       log indicates log stepping in [10^-3, 10^5] (default) 
    618       logq indicates log stepping in [10^-4, 10^1] 
    619       linear indicates linear stepping in [1, 1000] 
    620       zoom indicates linear stepping in [1000, 1010] 
    621       neg indicates linear stepping in [-100.1, 100.1] 
    622 and name is "all" or one of: 
     673        log indicates log stepping in [10^-3, 10^5] (default) 
     674        logq indicates log stepping in [10^-4, 10^1] 
     675        linear indicates linear stepping in [1, 1000] 
     676        zoom indicates linear stepping in [1000, 1010] 
     677        neg indicates linear stepping in [-100.1, 100.1] 
     678        start:stop:n[:stepping] indicates an n-step plot in [start, stop] 
     679            or [10^start, 10^stop] if stepping is "log" (default n=400) 
     680Some functions (notably gammainc/gammaincc) have an additional parameter A 
     681which can be set from the command line as A=value.  Default is A=1. 
     682 
     683Name is one of: 
    623684    """+names) 
    624685    sys.exit(1) 
  • extra/build_linux.sh

    rdb8756e r0d362b7  
    1414# Run tests 
    1515STATUS=0 
    16 python -m sasmodels.model_test opencl_and_dll all || STATUS=$? 
     16python -m sasmodels.model_test -v all || STATUS=$? 
    1717python -m sasmodels.resolution || STATUS=$? 
    1818 
  • extra/pylint.rc

    r6e68289 rdf20715  
    6767    locally-enabled, 
    6868    locally-disabled, 
     69    #no-else-return, 
     70    #invalid-name, 
     71    fixme, 
     72    too-many-lines, 
     73    too-many-locals, 
     74    too-many-branches, 
     75    too-many-statements, 
     76    #too-many-instance-attributes, 
     77    #too-many-return-statements, 
     78    #no-self-use, 
    6979 
    7080[REPORTS] 
     
    109119good-names=_, 
    110120    input, id, 
    111     h, i, j, k, n, p, v, w, x, y, z, 
    112     q, qx, qy, qz, Q, Qx, Qy, Qz, 
    113     dt, dx, dy, dz, dq, 
     121    h, i, j, k, n, p, u, v, w, x, y, z, r, dr, 
     122    q, qx, qy, qz, Q, Qx, Qy, Qz, nx, ny, nz, nq, nr, 
     123    dt, dx, dy, dz, dq, S, P, a, b, c, 
    114124    Iq, dIq, Iqxy, Iq_calc, 
    115125    ER, call_ER, VR, call_VR, 
     
    133143 
    134144# Regular expression matching correct variable names 
    135 variable-rgx=[a-z_][a-z0-9_{2D}{1D}]{2,30}$ 
     145variable-rgx=[a-z_][a-z0-9_{2D}{1D}]{1,30}$ 
    136146 
    137147# Naming hint for variable names 
     
    151161 
    152162# Regular expression matching correct argument names 
    153 argument-rgx=[a-z_][a-z0-9_]{2,30}$ 
     163argument-rgx=[a-z_][a-z0-9_]{1,30}$ 
    154164 
    155165# Naming hint for argument names 
     
    280290# (useful for modules/projects where namespaces are manipulated during runtime 
    281291# and thus existing member attributes cannot be deduced by static analysis 
    282 ignored-modules=bumps,sas,numpy,numpy.random,scipy,scipy.special 
     292ignored-modules=bumps,sas,numpy,numpy.random,scipy,scipy.special,pyopencl,pycuda,matplotlib.cm,ipyvolume,ipywidget 
    283293 
    284294# List of classes names for which member attributes should not be checked 
  • sasmodels/__init__.py

    re65c3ba r37f38ff  
    1414defining new models. 
    1515""" 
    16 __version__ = "0.97" 
     16__version__ = "0.99" 
    1717 
    1818def data_files(): 
  • sasmodels/bumps_model.py

    r49d1f8b8 rb297ba9  
    2626    from .kernel import KernelModel 
    2727    from .modelinfo import ModelInfo 
     28    from .resolution import Resolution 
    2829    Data = Union[Data1D, Data2D] 
    2930except ImportError: 
     
    3536    # when bumps is not on the path. 
    3637    from bumps.names import Parameter # type: ignore 
     38    from bumps.parameter import Reference # type: ignore 
    3739except ImportError: 
    3840    pass 
     
    139141    def __init__(self, data, model, cutoff=1e-5, name=None, extra_pars=None): 
    140142        # type: (Data, Model, float) -> None 
     143        # Allow resolution function to define fittable parameters.  We do this 
     144        # by creating reference parameters within the resolution object rather 
     145        # than modifying the object itself to use bumps parameters.  We need 
     146        # to reset the parameters each time the object has changed.  These 
     147        # additional parameters need to be returned from the fitting engine. 
     148        # To make them available to the user, they are added as top-level 
     149        # attributes to the experiment object.  The only change to the 
     150        # resolution function is that it needs an optional 'fittable' attribute 
     151        # which maps the internal name to the user visible name for the 
     152        # for the parameter. 
     153        self._resolution = None 
     154        self._resolution_pars = {} 
    141155        # remember inputs so we can inspect from outside 
    142156        self.name = data.filename if name is None else name 
     
    145159        self._interpret_data(data, model.sasmodel) 
    146160        self._cache = {} 
     161        # CRUFT: no longer need extra parameters 
     162        # Multiple scattering probability is now retrieved directly from the 
     163        # multiple scattering resolution function. 
    147164        self.extra_pars = extra_pars 
    148165 
     
    162179        return len(self.Iq) 
    163180 
     181    @property 
     182    def resolution(self): 
     183        # type: () -> Union[None, Resolution] 
     184        """ 
     185        :class:`sasmodels.Resolution` applied to the data, if any. 
     186        """ 
     187        return self._resolution 
     188 
     189    @resolution.setter 
     190    def resolution(self, value): 
     191        # type: (Resolution) -> None 
     192        """ 
     193        :class:`sasmodels.Resolution` applied to the data, if any. 
     194        """ 
     195        self._resolution = value 
     196 
     197        # Remove old resolution fitting parameters from experiment 
     198        for name in self._resolution_pars: 
     199            delattr(self, name) 
     200 
     201        # Create new resolution fitting parameters 
     202        res_pars = getattr(self._resolution, 'fittable', {}) 
     203        self._resolution_pars = { 
     204            name: Reference(self._resolution, refname, name=name) 
     205            for refname, name in res_pars.items() 
     206        } 
     207 
     208        # Add new resolution fitting parameters as experiment attributes 
     209        for name, ref in self._resolution_pars.items(): 
     210            setattr(self, name, ref) 
     211 
    164212    def parameters(self): 
    165213        # type: () -> Dict[str, Parameter] 
     
    168216        """ 
    169217        pars = self.model.parameters() 
    170         if self.extra_pars: 
     218        if self.extra_pars is not None: 
    171219            pars.update(self.extra_pars) 
     220        pars.update(self._resolution_pars) 
    172221        return pars 
    173222 
     
    230279        Save the model parameters and data into a file. 
    231280 
    232         Not Implemented. 
     281        Not Implemented except for sesans fits. 
    233282        """ 
    234283        if self.data_type == "sesans": 
  • sasmodels/compare.py

    r1fbadb2 rb297ba9  
    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 
     
    107108    -title="note" adds note to the plot title, after the model name 
    108109    -weights shows weights plots for the polydisperse parameters 
     110    -profile shows the sld profile if the model has a plottable sld profile 
    109111 
    110112    === output options === 
    111113    -edit starts the parameter explorer 
    112114    -help/-html shows the model docs instead of running the model 
     115 
     116    === environment variables === 
     117    -DSAS_MODELPATH=path sets directory containing custom models 
     118    -DSAS_OPENCL=vendor:device|cuda:device|none sets the target GPU device 
     119    -DXDG_CACHE_HOME=~/.cache sets the pyopencl cache root (linux only) 
     120    -DSAS_COMPILER=tinycc|msvc|mingw|unix sets the DLL compiler 
     121    -DSAS_OPENMP=1 turns on OpenMP for the DLLs 
     122    -DSAS_DLL_PATH=path sets the path to the compiled modules 
    113123 
    114124The interpretation of quad precision depends on architecture, and may 
     
    343353 
    344354    # If it is a list of choices, pick one at random with equal probability 
    345     # In practice, the model specific random generator will override. 
    346355    par = model_info.parameters[name] 
    347     if len(par.limits) > 2:  # choice list 
    348         return np.random.randint(len(par.limits)) 
     356    if par.choices:  # choice list 
     357        return np.random.randint(len(par.choices)) 
    349358 
    350359    # If it is a fixed range, pick from it with equal probability. 
     
    359368 
    360369    # Limit magnetic SLDs to a smaller range, from zero to iron=5/A^2 
    361     if par.name.startswith('M0:'): 
     370    if par.name.endswith('_M0'): 
    362371        return np.random.uniform(0, 5) 
    363372 
     
    529538    magnetic_pars = [] 
    530539    for p in parameters.user_parameters(pars, is2d): 
    531         if any(p.id.startswith(x) for x in ('M0:', 'mtheta:', 'mphi:')): 
     540        if any(p.id.endswith(x) for x in ('_M0', '_mtheta', '_mphi')): 
    532541            continue 
    533542        if p.id.startswith('up:'): 
     
    541550            pdtype=pars.get(p.id+"_pd_type", 'gaussian'), 
    542551            relative_pd=p.relative_pd, 
    543             M0=pars.get('M0:'+p.id, 0.), 
    544             mphi=pars.get('mphi:'+p.id, 0.), 
    545             mtheta=pars.get('mtheta:'+p.id, 0.), 
     552            M0=pars.get(p.id+'_M0', 0.), 
     553            mphi=pars.get(p.id+'_mphi', 0.), 
     554            mtheta=pars.get(p.id+'_mtheta', 0.), 
    546555        ) 
    547556        lines.append(_format_par(p.name, **fields)) 
     
    609618    if suppress: 
    610619        for p in pars: 
    611             if p.startswith("M0:"): 
     620            if p.endswith("_M0"): 
    612621                pars[p] = 0 
    613622    else: 
     
    615624        first_mag = None 
    616625        for p in pars: 
    617             if p.startswith("M0:"): 
     626            if p.endswith("_M0"): 
    618627                any_mag |= (pars[p] != 0) 
    619628                if first_mag is None: 
     
    645654 
    646655def make_data(opts): 
    647     # type: (Dict[str, Any]) -> Tuple[Data, np.ndarray] 
     656    # type: (Dict[str, Any], float) -> Tuple[Data, np.ndarray] 
    648657    """ 
    649658    Generate an empty dataset, used with the model to set Q points 
     
    667676        if opts['zero']: 
    668677            q = np.hstack((0, q)) 
    669         data = empty_data1D(q, resolution=res) 
     678        # TODO: provide command line control of lambda and Delta lambda/lambda 
     679        #L, dLoL = 5, 0.14/np.sqrt(6)  # wavelength and 14% triangular FWHM 
     680        L, dLoL = 0, 0 
     681        data = empty_data1D(q, resolution=res, L=L, dL=L*dLoL) 
    670682        index = slice(None, None) 
    671683    return data, index 
     
    713725        set_integration_size(model_info, ngauss) 
    714726 
    715     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())): 
    716729        raise RuntimeError("OpenCL not available " + kernelcl.OPENCL_ERROR) 
    717730 
     
    759772        opts['pars'] = parse_pars(opts, maxdim=maxdim) 
    760773        if opts['pars'] is None: 
    761             return 
     774            return limits 
    762775        result = run_models(opts, verbose=True) 
    763776        if opts['plot']: 
     
    772785            dim = base._kernel.dim 
    773786            plot_weights(model_info, get_mesh(model_info, base_pars, dim=dim)) 
     787        if opts['show_profile']: 
     788            import pylab 
     789            base, comp = opts['engines'] 
     790            base_pars, comp_pars = opts['pars'] 
     791            have_base = base._kernel.info.profile is not None 
     792            have_comp = ( 
     793                comp is not None 
     794                and comp._kernel.info.profile is not None 
     795                and base_pars != comp_pars 
     796            ) 
     797            if have_base or have_comp: 
     798                pylab.figure() 
     799            if have_base: 
     800                plot_profile(base._kernel.info, **base_pars) 
     801            if have_comp: 
     802                plot_profile(comp._kernel.info, label='comp', **comp_pars) 
     803                pylab.legend() 
    774804    if opts['plot']: 
    775805        import matplotlib.pyplot as plt 
     
    777807    return limits 
    778808 
     809def plot_profile(model_info, label='base', **args): 
     810    # type: (ModelInfo, List[Tuple[float, np.ndarray, np.ndarray]]) -> None 
     811    """ 
     812    Plot the profile returned by the model profile method. 
     813 
     814    *model_info* defines model parameters, etc. 
     815 
     816    *mesh* is a list of tuples containing (*value*, *dispersity*, *weights*) 
     817    for each parameter, where (*dispersity*, *weights*) pairs are the 
     818    distributions to be plotted. 
     819    """ 
     820    import pylab 
     821 
     822    args = dict((k, v) for k, v in args.items() 
     823                if "_pd" not in k 
     824                and ":" not in k 
     825                and k not in ("background", "scale", "theta", "phi", "psi")) 
     826    args = args.copy() 
     827 
     828    args.pop('scale', 1.) 
     829    args.pop('background', 0.) 
     830    z, rho = model_info.profile(**args) 
     831    #pylab.interactive(True) 
     832    pylab.plot(z, rho, '-', label=label) 
     833    pylab.grid(True) 
     834    #pylab.show() 
     835 
     836 
     837 
    779838def run_models(opts, verbose=False): 
    780839    # type: (Dict[str, Any]) -> Dict[str, Any] 
     
    786845    base_n, comp_n = opts['count'] 
    787846    base_pars, comp_pars = opts['pars'] 
    788     data = opts['data'] 
     847    base_data, comp_data = opts['data'] 
    789848 
    790849    comparison = comp is not None 
     
    800859            print("%s t=%.2f ms, intensity=%.0f" 
    801860                  % (base.engine, base_time, base_value.sum())) 
    802         _show_invalid(data, base_value) 
     861        _show_invalid(base_data, base_value) 
    803862    except ImportError: 
    804863        traceback.print_exc() 
     
    812871                print("%s t=%.2f ms, intensity=%.0f" 
    813872                      % (comp.engine, comp_time, comp_value.sum())) 
    814             _show_invalid(data, comp_value) 
     873            _show_invalid(base_data, comp_value) 
    815874        except ImportError: 
    816875            traceback.print_exc() 
     
    837896    # work with trimmed data, not the full set 
    838897    sorted_err = np.sort(abs(err.compressed())) 
    839     if len(sorted_err) == 0: 
     898    if sorted_err.size == 0: 
    840899        print(label + "  no valid values") 
    841900        return 
     
    866925    have_base, have_comp = (base_value is not None), (comp_value is not None) 
    867926    base, comp = opts['engines'] 
    868     data = opts['data'] 
     927    base_data, comp_data = opts['data'] 
    869928    use_data = (opts['datafile'] is not None) and (have_base ^ have_comp) 
    870929 
    871930    # Plot if requested 
    872931    view = opts['view'] 
     932    #view = 'log' 
    873933    if limits is None: 
    874934        vmin, vmax = np.inf, -np.inf 
     
    884944        if have_comp: 
    885945            plt.subplot(131) 
    886         plot_theory(data, base_value, view=view, use_data=use_data, limits=limits) 
     946        plot_theory(base_data, base_value, view=view, use_data=use_data, limits=limits) 
    887947        plt.title("%s t=%.2f ms"%(base.engine, base_time)) 
    888948        #cbar_title = "log I" 
     
    891951            plt.subplot(132) 
    892952        if not opts['is2d'] and have_base: 
    893             plot_theory(data, base_value, view=view, use_data=use_data, limits=limits) 
    894         plot_theory(data, comp_value, view=view, use_data=use_data, limits=limits) 
     953            plot_theory(comp_data, base_value, view=view, use_data=use_data, limits=limits) 
     954        plot_theory(comp_data, comp_value, view=view, use_data=use_data, limits=limits) 
    895955        plt.title("%s t=%.2f ms"%(comp.engine, comp_time)) 
    896956        #cbar_title = "log I" 
     
    908968            err[err > cutoff] = cutoff 
    909969        #err,errstr = base/comp,"ratio" 
    910         plot_theory(data, None, resid=err, view=errview, use_data=use_data) 
     970        # Note: base_data only since base and comp have same q values (though 
     971        # perhaps different resolution), and we are plotting the difference 
     972        # at each q 
     973        plot_theory(base_data, None, resid=err, view=errview, use_data=use_data) 
    911974        plt.xscale('log' if view == 'log' and not opts['is2d'] else 'linear') 
    912975        plt.legend(['P%d'%(k+1) for k in range(setnum+1)], loc='best') 
     
    9421005OPTIONS = [ 
    9431006    # Plotting 
    944     'plot', 'noplot', 'weights', 
     1007    'plot', 'noplot', 
     1008    'weights', 'profile', 
    9451009    'linear', 'log', 'q4', 
    9461010    'rel', 'abs', 
     
    10581122 
    10591123    invalid = [o[1:] for o in flags 
    1060                if o[1:] not in NAME_OPTIONS 
    1061                and not any(o.startswith('-%s='%t) for t in VALUE_OPTIONS)] 
     1124               if not (o[1:] in NAME_OPTIONS 
     1125                       or any(o.startswith('-%s='%t) for t in VALUE_OPTIONS) 
     1126                       or o.startswith('-D'))] 
    10621127    if invalid: 
    10631128        print("Invalid options: %s"%(", ".join(invalid))) 
     
    10751140        'qmax'      : 0.05, 
    10761141        'nq'        : 128, 
    1077         'res'       : 0.0, 
     1142        'res'       : '0.0', 
    10781143        'noise'     : 0.0, 
    10791144        'accuracy'  : 'Low', 
     
    10871152        'rel_err'   : True, 
    10881153        'explore'   : False, 
    1089         'use_demo'  : True, 
     1154        'use_demo'  : False, 
    10901155        'zero'      : False, 
    10911156        'html'      : False, 
     
    10961161        'count'     : '1', 
    10971162        'show_weights' : False, 
     1163        'show_profile' : False, 
    10981164        'sphere'    : 0, 
    10991165        'ngauss'    : '0', 
     
    11151181        elif arg.startswith('-q='): 
    11161182            opts['qmin'], opts['qmax'] = [float(v) for v in arg[3:].split(':')] 
    1117         elif arg.startswith('-res='):      opts['res'] = float(arg[5:]) 
     1183        elif arg.startswith('-res='):      opts['res'] = arg[5:] 
    11181184        elif arg.startswith('-noise='):    opts['noise'] = float(arg[7:]) 
    11191185        elif arg.startswith('-sets='):     opts['sets'] = int(arg[6:]) 
     
    11561222        elif arg == '-default': opts['use_demo'] = False 
    11571223        elif arg == '-weights': opts['show_weights'] = True 
     1224        elif arg == '-profile': opts['show_profile'] = True 
    11581225        elif arg == '-html':    opts['html'] = True 
    11591226        elif arg == '-help':    opts['html'] = True 
     1227        elif arg.startswith('-D'): 
     1228            var, val = arg[2:].split('=') 
     1229            os.environ[var] = val 
    11601230    # pylint: enable=bad-whitespace,C0321 
    11611231 
     
    11731243    if opts['qmin'] is None: 
    11741244        opts['qmin'] = 0.001*opts['qmax'] 
    1175     if opts['datafile'] is not None: 
    1176         data = load_data(os.path.expanduser(opts['datafile'])) 
    1177     else: 
    1178         data, _ = make_data(opts) 
    11791245 
    11801246    comparison = any(PAR_SPLIT in v for v in values) 
     
    12161282        opts['cutoff'] = [float(opts['cutoff'])]*2 
    12171283 
    1218     base = make_engine(model_info[0], data, opts['engine'][0], 
     1284    if PAR_SPLIT in opts['res']: 
     1285        opts['res'] = [float(k) for k in opts['res'].split(PAR_SPLIT, 2)] 
     1286        comparison = True 
     1287    else: 
     1288        opts['res'] = [float(opts['res'])]*2 
     1289 
     1290    if opts['datafile'] is not None: 
     1291        data0 = load_data(os.path.expanduser(opts['datafile'])) 
     1292        data = data0, data0 
     1293    else: 
     1294        # Hack around the fact that make_data doesn't take a pair of resolutions 
     1295        res = opts['res'] 
     1296        opts['res'] = res[0] 
     1297        data0, _ = make_data(opts) 
     1298        if res[0] != res[1]: 
     1299            opts['res'] = res[1] 
     1300            data1, _ = make_data(opts) 
     1301        else: 
     1302            data1 = data0 
     1303        opts['res'] = res 
     1304        data = data0, data1 
     1305 
     1306    base = make_engine(model_info[0], data[0], opts['engine'][0], 
    12191307                       opts['cutoff'][0], opts['ngauss'][0]) 
    12201308    if comparison: 
    1221         comp = make_engine(model_info[1], data, opts['engine'][1], 
     1309        comp = make_engine(model_info[1], data[1], opts['engine'][1], 
    12221310                           opts['cutoff'][1], opts['ngauss'][1]) 
    12231311    else: 
  • sasmodels/convert.py

    ra69d8cd rb297ba9  
    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 
     
    165168    if version == (3, 1, 2): 
    166169        oldpars = _hand_convert_3_1_2_to_4_1(name, oldpars) 
     170    if version < (4, 2, 0): 
     171        oldpars = _rename_magnetic_pars(oldpars) 
    167172    return oldpars 
     173 
     174def _rename_magnetic_pars(pars): 
     175    """ 
     176    Change from M0:par to par_M0, etc. 
     177    """ 
     178    keys = list(pars.items()) 
     179    for k in keys: 
     180        if k.startswith('M0:'): 
     181            pars[k[3:]+'_M0'] = pars.pop(k) 
     182        elif k.startswith('mtheta:'): 
     183            pars[k[7:]+'_mtheta'] = pars.pop(k) 
     184        elif k.startswith('mphi:'): 
     185            pars[k[5:]+'_mphi'] = pars.pop(k) 
     186        elif k.startswith('up:'): 
     187            pars['up_'+k[3:]] = pars.pop(k) 
     188    return pars 
    168189 
    169190def _hand_convert_3_1_2_to_4_1(name, oldpars): 
     
    363384 
    364385def revert_name(model_info): 
     386    """Translate model name back to the name used in SasView 3.x""" 
    365387    oldname, _ = CONVERSION_TABLE.get(model_info.id, [None, {}]) 
    366388    return oldname 
     
    632654 
    633655def test_backward_forward(): 
     656    """ 
     657    Test conversion of model parameters from 4.x to 3.x and back. 
     658    """ 
    634659    from .core import list_models 
    635660    L = lambda name: _check_one(name, seed=1) 
  • sasmodels/core.py

    r3221de0 rb297ba9  
    2121from . import mixture 
    2222from . import kernelpy 
     23from . import kernelcuda 
    2324from . import kernelcl 
    2425from . import kerneldll 
     
    3738if CUSTOM_MODEL_PATH == "": 
    3839    CUSTOM_MODEL_PATH = joinpath(os.path.expanduser("~"), ".sasmodels", "custom_models") 
    39     if not os.path.isdir(CUSTOM_MODEL_PATH): 
    40         os.makedirs(CUSTOM_MODEL_PATH) 
     40    #if not os.path.isdir(CUSTOM_MODEL_PATH): 
     41    #    os.makedirs(CUSTOM_MODEL_PATH) 
    4142 
    4243# TODO: refactor composite model support 
     
    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) 
     
    233240        if not callable(model_info.Iq): 
    234241            source = generate.make_source(model_info)['dll'] 
    235             old_path = kerneldll.DLL_PATH 
     242            old_path = kerneldll.SAS_DLL_PATH 
    236243            try: 
    237                 kerneldll.DLL_PATH = path 
     244                kerneldll.SAS_DLL_PATH = path 
    238245                dll = kerneldll.make_dll(source, model_info, dtype=numpy_dtype) 
    239246            finally: 
    240                 kerneldll.DLL_PATH = old_path 
     247                kerneldll.SAS_DLL_PATH = old_path 
    241248            compiled_dlls.append(dll) 
    242249    return compiled_dlls 
     
    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 
    303323 
    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))) 
    313  
    314324def test_composite_order(): 
     325    """ 
     326    Check that mixture models produce the same result independent of ordder. 
     327    """ 
    315328    def test_models(fst, snd): 
    316329        """Confirm that two models produce the same parameters""" 
     
    327340 
    328341    def build_test(first, second): 
     342        """Construct pair model test""" 
    329343        test = lambda description: test_models(first, second) 
    330344        description = first + " vs. " + second 
     
    365379    model = load_model("cylinder@hardsphere*sphere") 
    366380    actual = [p.name for p in model.info.parameters.kernel_parameters] 
    367     target = ("sld sld_solvent radius length theta phi volfraction" 
     381    target = ("sld sld_solvent radius length theta phi" 
     382              " radius_effective volfraction " 
     383              " structure_factor_mode radius_effective_mode" 
    368384              " A_sld A_sld_solvent A_radius").split() 
    369385    assert target == actual, "%s != %s"%(target, actual) 
    370386 
     387def list_models_main(): 
     388    # type: () -> int 
     389    """ 
     390    Run list_models as a main program.  See :func:`list_models` for the 
     391    kinds of models that can be requested on the command line. 
     392    """ 
     393    import sys 
     394    kind = sys.argv[1] if len(sys.argv) > 1 else "all" 
     395    try: 
     396        models = list_models(kind) 
     397        print("\n".join(models)) 
     398    except Exception: 
     399        print(list_models.__doc__) 
     400        return 1 
     401    return 0 
     402 
    371403if __name__ == "__main__": 
    372404    list_models_main() 
  • sasmodels/custom/__init__.py

    r0f48f1e rd321747  
    1212import sys 
    1313import os 
    14 from os.path import basename, splitext 
     14from os.path import basename, splitext, join as joinpath, exists, dirname 
    1515 
    1616try: 
     
    1818    from importlib.util import spec_from_file_location, module_from_spec  # type: ignore 
    1919    def load_module_from_path(fullname, path): 
     20        # type: (str, str) -> "module" 
    2021        """load module from *path* as *fullname*""" 
    2122        spec = spec_from_file_location(fullname, os.path.expanduser(path)) 
     
    2728    import imp 
    2829    def load_module_from_path(fullname, path): 
     30        # type: (str, str) -> "module" 
    2931        """load module from *path* as *fullname*""" 
    3032        # Clear out old definitions, if any 
     
    3739        return module 
    3840 
     41_MODULE_CACHE = {} # type: Dict[str, Tuple("module", int)] 
     42_MODULE_DEPENDS = {} # type: Dict[str, List[str]] 
     43_MODULE_DEPENDS_STACK = [] # type: List[str] 
    3944def load_custom_kernel_module(path): 
     45    # type: str -> "module" 
    4046    """load SAS kernel from *path* as *sasmodels.custom.modelname*""" 
    4147    # Pull off the last .ext if it exists; there may be others 
    4248    name = basename(splitext(path)[0]) 
    43     # Placing the model in the 'sasmodels.custom' name space. 
    44     kernel_module = load_module_from_path('sasmodels.custom.'+name, 
    45                                           os.path.expanduser(path)) 
    46     return kernel_module 
     49    path = os.path.expanduser(path) 
     50 
     51    # Reload module if necessary. 
     52    if need_reload(path): 
     53        # Assume the module file is the only dependency 
     54        _MODULE_DEPENDS[path] = set([path]) 
     55 
     56        # Load the module while pushing it onto the dependency stack.  If 
     57        # this triggers any submodules, then they will add their dependencies 
     58        # to this module as the "working_on" parent.  Pop the stack when the 
     59        # module is loaded. 
     60        _MODULE_DEPENDS_STACK.append(path) 
     61        module = load_module_from_path('sasmodels.custom.'+name, path) 
     62        _MODULE_DEPENDS_STACK.pop() 
     63 
     64        # Include external C code in the dependencies.  We are looking 
     65        # for module.source and assuming that it is a list of C source files 
     66        # relative to the module itself.  Any files that do not exist, 
     67        # such as those in the standard libraries, will be ignored. 
     68        # TODO: look in builtin module path for standard c sources 
     69        # TODO: share code with generate.model_sources 
     70        c_sources = getattr(module, 'source', None) 
     71        if isinstance(c_sources, (list, tuple)): 
     72            _MODULE_DEPENDS[path].update(_find_sources(path, c_sources)) 
     73 
     74        # Cache the module, and tag it with the newest timestamp 
     75        timestamp = max(os.path.getmtime(f) for f in _MODULE_DEPENDS[path]) 
     76        _MODULE_CACHE[path] = module, timestamp 
     77 
     78        #print("loading", os.path.basename(path), _MODULE_CACHE[path][1], 
     79        #    [os.path.basename(p) for p in _MODULE_DEPENDS[path]]) 
     80 
     81    # Add path and all its dependence to the parent module, if there is one. 
     82    if _MODULE_DEPENDS_STACK: 
     83        working_on = _MODULE_DEPENDS_STACK[-1] 
     84        _MODULE_DEPENDS[working_on].update(_MODULE_DEPENDS[path]) 
     85 
     86    return _MODULE_CACHE[path][0] 
     87 
     88def need_reload(path): 
     89    # type: str -> bool 
     90    """ 
     91    Return True if any path dependencies have a timestamp newer than the time 
     92    when the path was most recently loaded. 
     93    """ 
     94    # TODO: fails if a dependency has a modification time in the future 
     95    # If the newest dependency has a time stamp in the future, then this 
     96    # will be recorded as the cached time.  When a second dependency 
     97    # is updated to the current time stamp, it will still be considered 
     98    # older than the current build and the reload will not be triggered. 
     99    # Could instead treat all future times as 0 here and in the code above 
     100    # which records the newest timestamp.  This will force a reload when 
     101    # the future time is reached, but other than that should perform 
     102    # correctly.  Probably not worth the extra code... 
     103    _, cache_time = _MODULE_CACHE.get(path, (None, -1)) 
     104    depends = _MODULE_DEPENDS.get(path, [path]) 
     105    #print("reload", any(cache_time < os.path.getmtime(p) for p in depends)) 
     106    #for f in depends: print(">>>  ", f, os.path.getmtime(f)) 
     107    return any(cache_time < os.path.getmtime(p) for p in depends) 
     108 
     109def _find_sources(path, source_list): 
     110    # type: (str, List[str]) -> List[str] 
     111    """ 
     112    Return a list of the sources relative to base file; ignore any that 
     113    are not found. 
     114    """ 
     115    root = dirname(path) 
     116    found = [] 
     117    for source_name in source_list: 
     118        source_path = joinpath(root, source_name) 
     119        if exists(source_path): 
     120            found.append(source_path) 
     121    return found 
  • sasmodels/data.py

    rd86f0fc rb297ba9  
    3636 
    3737import numpy as np  # type: ignore 
     38from numpy import sqrt, sin, cos, pi 
    3839 
    3940# pylint: disable=unused-import 
     
    182183    *x* is spin echo length and *y* is polarization (P/P0). 
    183184    """ 
     185    isSesans = True 
    184186    def __init__(self, **kw): 
    185187        Data1D.__init__(self, **kw) 
     
    300302        self.wavelength_unit = "A" 
    301303 
    302  
    303 def empty_data1D(q, resolution=0.0): 
     304class Sample(object): 
     305    """ 
     306    Sample attributes. 
     307    """ 
     308    def __init__(self): 
     309        # type: () -> None 
     310        pass 
     311 
     312def empty_data1D(q, resolution=0.0, L=0., dL=0.): 
    304313    # type: (np.ndarray, float) -> Data1D 
    305     """ 
     314    r""" 
    306315    Create empty 1D data using the given *q* as the x value. 
    307316 
    308     *resolution* dq/q defaults to 5%. 
     317    rms *resolution* $\Delta q/q$ defaults to 0%.  If wavelength *L* and rms 
     318    wavelength divergence *dL* are defined, then *resolution* defines 
     319    rms $\Delta \theta/\theta$ for the lowest *q*, with $\theta$ derived from 
     320    $q = 4\pi/\lambda \sin(\theta)$. 
    309321    """ 
    310322 
     
    313325    Iq, dIq = None, None 
    314326    q = np.asarray(q) 
    315     data = Data1D(q, Iq, dx=resolution * q, dy=dIq) 
     327    if L != 0 and resolution != 0: 
     328        theta = np.arcsin(q*L/(4*pi)) 
     329        dtheta = theta[0]*resolution 
     330        ## Solving Gaussian error propagation from 
     331        ##   Dq^2 = (dq/dL)^2 DL^2 + (dq/dtheta)^2 Dtheta^2 
     332        ## gives 
     333        ##   (Dq/q)^2 = (DL/L)**2 + (Dtheta/tan(theta))**2 
     334        ## Take the square root and multiply by q, giving 
     335        ##   Dq = (4*pi/L) * sqrt((sin(theta)*dL/L)**2 + (cos(theta)*dtheta)**2) 
     336        dq = (4*pi/L) * sqrt((sin(theta)*dL/L)**2 + (cos(theta)*dtheta)**2) 
     337    else: 
     338        dq = resolution * q 
     339    data = Data1D(q, Iq, dx=dq, dy=dIq) 
    316340    data.filename = "fake data" 
    317341    return data 
     
    476500            mdata = masked_array(data.y, data.mask.copy()) 
    477501            mdata[~np.isfinite(mdata)] = masked 
    478             if view is 'log': 
     502            if view == 'log': 
    479503                mdata[mdata <= 0] = masked 
    480504            plt.errorbar(data.x, scale*mdata, yerr=data.dy, fmt='.') 
     
    486510            # Note: masks merge, so any masked theory points will stay masked, 
    487511            # and the data mask will be added to it. 
    488             mtheory = masked_array(theory, data.mask.copy()) 
     512            #mtheory = masked_array(theory, data.mask.copy()) 
     513            theory_x = data.x[data.mask == 0] 
     514            mtheory = masked_array(theory) 
    489515            mtheory[~np.isfinite(mtheory)] = masked 
    490             if view is 'log': 
     516            if view == 'log': 
    491517                mtheory[mtheory <= 0] = masked 
    492             plt.plot(data.x, scale*mtheory, '-') 
     518            plt.plot(theory_x, scale*mtheory, '-') 
    493519            all_positive = all_positive and (mtheory > 0).all() 
    494520            some_present = some_present or (mtheory.count() > 0) 
     
    526552 
    527553    if use_resid: 
    528         mresid = masked_array(resid, data.mask.copy()) 
     554        theory_x = data.x[data.mask == 0] 
     555        mresid = masked_array(resid) 
    529556        mresid[~np.isfinite(mresid)] = masked 
    530557        some_present = (mresid.count() > 0) 
     
    532559        if num_plots > 1: 
    533560            plt.subplot(1, num_plots, use_calc + 2) 
    534         plt.plot(data.x, mresid, '.') 
     561        plt.plot(theory_x, mresid, '.') 
    535562        plt.xlabel("$q$/A$^{-1}$") 
    536563        plt.ylabel('residuals') 
  • sasmodels/details.py

    r885753a rb297ba9  
    236236    npars = kernel.info.parameters.npars 
    237237    nvalues = kernel.info.parameters.nvalues 
    238     scalars = [value for value, _dispersity, _weight in mesh] 
     238    scalars = [value for value, dispersity, weight in mesh] 
    239239    # skipping scale and background when building values and weights 
    240     _values, dispersity, weights = zip(*mesh[2:npars+2]) if npars else ((), (), ()) 
    241     #weights = correct_theta_weights(kernel.info.parameters, dispersity, weights) 
    242     length = np.array([len(w) for w in weights]) 
     240    value, dispersity, weight = zip(*mesh[2:npars+2]) if npars else ((), (), ()) 
     241    #weight = correct_theta_weights(kernel.info.parameters, dispersity, weight) 
     242    length = np.array([len(w) for w in weight]) 
    243243    offset = np.cumsum(np.hstack((0, length))) 
    244244    call_details = make_details(kernel.info, length, offset[:-1], offset[-1]) 
     
    246246    data_len = nvalues + 2*sum(len(v) for v in dispersity) 
    247247    extra = (32 - data_len%32)%32 
    248     data = np.hstack((scalars,) + dispersity + weights + ZEROS[:extra]) 
     248    data = np.hstack((scalars,) + dispersity + weight + ZEROS[:extra]) 
    249249    data = data.astype(kernel.dtype) 
    250250    is_magnetic = convert_magnetism(kernel.info.parameters, data) 
  • sasmodels/direct_model.py

    rd18d6dd rb297ba9  
    3131from . import resolution2d 
    3232from .details import make_kernel_args, dispersion_mesh 
     33from .product import RADIUS_MODE_ID 
    3334 
    3435# pylint: disable=unused-import 
     
    6364    return calculator(call_details, values, cutoff, is_magnetic) 
    6465 
    65 def call_ER(model_info, pars): 
    66     # type: (ModelInfo, ParameterSet) -> float 
    67     """ 
    68     Call the model ER function using *values*. 
    69  
    70     *model_info* is either *model.info* if you have a loaded model, 
    71     or *kernel.info* if you have a model kernel prepared for evaluation. 
    72     """ 
    73     if model_info.ER is None: 
    74         return 1.0 
    75     elif not model_info.parameters.form_volume_parameters: 
    76         # handle the case where ER is provided but model is not polydisperse 
    77         return model_info.ER() 
    78     else: 
    79         value, weight = _vol_pars(model_info, pars) 
    80         individual_radii = model_info.ER(*value) 
    81         return np.sum(weight*individual_radii) / np.sum(weight) 
    82  
    83  
    84 def call_VR(model_info, pars): 
    85     # type: (ModelInfo, ParameterSet) -> float 
    86     """ 
    87     Call the model VR function using *pars*. 
    88  
    89     *model_info* is either *model.info* if you have a loaded model, 
    90     or *kernel.info* if you have a model kernel prepared for evaluation. 
    91     """ 
    92     if model_info.VR is None: 
    93         return 1.0 
    94     elif not model_info.parameters.form_volume_parameters: 
    95         # handle the case where ER is provided but model is not polydisperse 
    96         return model_info.VR() 
    97     else: 
    98         value, weight = _vol_pars(model_info, pars) 
    99         whole, part = model_info.VR(*value) 
    100         return np.sum(weight*part)/np.sum(weight*whole) 
    101  
    102  
    103 def call_profile(model_info, **pars): 
    104     # type: (ModelInfo, ...) -> Tuple[np.ndarray, np.ndarray, Tuple[str, str]] 
     66def call_Fq(calculator, pars, cutoff=0., mono=False): 
     67    # type: (Kernel, ParameterSet, float, bool) -> np.ndarray 
     68    """ 
     69    Like :func:`call_kernel`, but returning F, F^2, R_eff, V_shell, V_form/V_shell. 
     70 
     71    For solid objects V_shell is equal to V_form and the volume ratio is 1. 
     72 
     73    Use parameter *radius_effective_mode* to select the effective radius 
     74    calculation. 
     75    """ 
     76    R_eff_type = int(pars.pop(RADIUS_MODE_ID, 1.0)) 
     77    mesh = get_mesh(calculator.info, pars, dim=calculator.dim, mono=mono) 
     78    #print("pars", list(zip(*mesh))[0]) 
     79    call_details, values, is_magnetic = make_kernel_args(calculator, mesh) 
     80    #print("values:", values) 
     81    return calculator.Fq(call_details, values, cutoff, is_magnetic, R_eff_type) 
     82 
     83def call_profile(model_info, pars=None): 
     84    # type: (ModelInfo, ParameterSet) -> Tuple[np.ndarray, np.ndarray, Tuple[str, str]] 
    10585    """ 
    10686    Returns the profile *x, y, (xlabel, ylabel)* representing the model. 
    10787    """ 
     88    if pars is None: 
     89        pars = {} 
    10890    args = {} 
    10991    for p in model_info.parameters.kernel_parameters: 
     
    241223            else: 
    242224                Iq, dIq = None, None 
    243             #self._theory = np.zeros_like(q) 
    244             q_vectors = [res.q_calc] 
    245225        elif self.data_type == 'Iqxy': 
    246226            #if not model.info.parameters.has_2d: 
     
    250230            qmax = getattr(data, 'qmax', np.inf) 
    251231            accuracy = getattr(data, 'accuracy', 'Low') 
    252             index = ~data.mask & (q >= qmin) & (q <= qmax) 
     232            index = (data.mask == 0) & (q >= qmin) & (q <= qmax) 
    253233            if data.data is not None: 
    254234                index &= ~np.isnan(data.data) 
     
    259239            res = resolution2d.Pinhole2D(data=data, index=index, 
    260240                                         nsigma=3.0, accuracy=accuracy) 
    261             #self._theory = np.zeros_like(self.Iq) 
    262             q_vectors = res.q_calc 
    263241        elif self.data_type == 'Iq': 
    264242            index = (data.x >= data.qmin) & (data.x <= data.qmax) 
     243            mask = getattr(data, 'mask', None) 
     244            if mask is not None: 
     245                index &= (mask == 0) 
    265246            if data.y is not None: 
    266247                index &= ~np.isnan(data.y) 
     
    282263            else: 
    283264                res = resolution.Perfect1D(data.x[index]) 
    284  
    285             #self._theory = np.zeros_like(self.Iq) 
    286             q_vectors = [res.q_calc] 
    287265        elif self.data_type == 'Iq-oriented': 
    288266            index = (data.x >= data.qmin) & (data.x <= data.qmax) 
     
    300278                                      qx_width=data.dxw[index], 
    301279                                      qy_width=data.dxl[index]) 
    302             q_vectors = res.q_calc 
    303280        else: 
    304281            raise ValueError("Unknown data type") # never gets here 
     
    306283        # Remember function inputs so we can delay loading the function and 
    307284        # so we can save/restore state 
    308         self._kernel_inputs = q_vectors 
    309285        self._kernel = None 
    310286        self.Iq, self.dIq, self.index = Iq, dIq, index 
     
    343319        # type: (ParameterSet, float) -> np.ndarray 
    344320        if self._kernel is None: 
    345             self._kernel = self._model.make_kernel(self._kernel_inputs) 
     321            # TODO: change interfaces so that resolution returns kernel inputs 
     322            # Maybe have resolution always return a tuple, or maybe have 
     323            # make_kernel accept either an ndarray or a pair of ndarrays. 
     324            kernel_inputs = self.resolution.q_calc 
     325            if isinstance(kernel_inputs, np.ndarray): 
     326                kernel_inputs = (kernel_inputs,) 
     327            self._kernel = self._model.make_kernel(kernel_inputs) 
    346328 
    347329        # Need to pull background out of resolution for multiple scattering 
    348         background = pars.get('background', 0.) 
     330        default_background = self._model.info.parameters.common_parameters[1].default 
     331        background = pars.get('background', default_background) 
    349332        pars = pars.copy() 
    350333        pars['background'] = 0. 
     
    399382        Generate a plottable profile. 
    400383        """ 
    401         return call_profile(self.model.info, **pars) 
     384        return call_profile(self.model.info, pars) 
    402385 
    403386def main(): 
     
    415398    model_name = sys.argv[1] 
    416399    call = sys.argv[2].upper() 
    417     if call != "ER_VR": 
    418         try: 
    419             values = [float(v) for v in call.split(',')] 
    420         except ValueError: 
    421             values = [] 
    422         if len(values) == 1: 
    423             q, = values 
    424             data = empty_data1D([q]) 
    425         elif len(values) == 2: 
    426             qx, qy = values 
    427             data = empty_data2D([qx], [qy]) 
    428         else: 
    429             print("use q or qx,qy or ER or VR") 
    430             sys.exit(1) 
     400    try: 
     401        values = [float(v) for v in call.split(',')] 
     402    except ValueError: 
     403        values = [] 
     404    if len(values) == 1: 
     405        q, = values 
     406        data = empty_data1D([q]) 
     407    elif len(values) == 2: 
     408        qx, qy = values 
     409        data = empty_data2D([qx], [qy]) 
    431410    else: 
    432         data = empty_data1D([0.001])  # Data not used in ER/VR 
     411        print("use q or qx,qy") 
     412        sys.exit(1) 
    433413 
    434414    model_info = load_model_info(model_name) 
     
    438418                for pair in sys.argv[3:] 
    439419                for k, v in [pair.split('=')]) 
    440     if call == "ER_VR": 
    441         ER = call_ER(model_info, pars) 
    442         VR = call_VR(model_info, pars) 
    443         print(ER, VR) 
    444     else: 
    445         Iq = calculator(**pars) 
    446         print(Iq[0]) 
     420    Iq = calculator(**pars) 
     421    print(Iq[0]) 
    447422 
    448423if __name__ == "__main__": 
  • sasmodels/generate.py

    rd86f0fc rb297ba9  
    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 
     
    197190 
    198191def get_data_path(external_dir, target_file): 
     192    """ 
     193    Search for the target file relative in the installed application. 
     194 
     195    Search first in the location of the generate module in case we are 
     196    running directly from the distribution.  Search next to the python 
     197    executable for windows installs.  Search in the ../Resources directory 
     198    next to the executable for Mac OS/X installs. 
     199    """ 
    199200    path = abspath(dirname(__file__)) 
    200201    if exists(joinpath(path, target_file)): 
     
    232233# 
    233234# NOTE: there is an RST_PROLOG at the end of this file which is NOT 
    234 # used for the bundled documentation. Still as long as we are defining the macros 
    235 # in two places any new addition should define the macro in both places. 
     235# used for the bundled documentation. Still as long as we are defining the 
     236# macros in two places any new addition should define the macro in both places. 
    236237RST_UNITS = { 
    237238    "Ang": "|Ang|", 
     
    537538 
    538539def test_tag_float(): 
    539     """check that floating point constants are properly identified and tagged with 'f'""" 
     540    """Check that floating point constants are identified and tagged with 'f'""" 
    540541 
    541542    cases = """ 
     
    671672 
    672673 
    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*[(]", 
     674_IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ac|abc|xy))\s*[(]", 
    675675                           flags=re.MULTILINE) 
    676676def find_xy_mode(source): 
     
    701701    return 'qa' 
    702702 
     703# Note: not presently used.  Need to know whether Fq is available before 
     704# trying to compile the source.  Leave the code here in case we decide that 
     705# define have_Fq for each form factor is too tedious and error prone. 
     706_FQ_PATTERN = re.compile(r"(^|\s)void\s+Fq[(]", flags=re.MULTILINE) 
     707def contains_Fq(source): 
     708    # type: (List[str]) -> bool 
     709    """ 
     710    Return True if C source defines "void Fq(". 
     711    """ 
     712    for code in source: 
     713        if _FQ_PATTERN.search(code) is not None: 
     714            return True 
     715    return False 
     716 
     717_SHELL_VOLUME_PATTERN = re.compile(r"(^|\s)double\s+shell_volume[(]", 
     718                                   flags=re.MULTILINE) 
     719def contains_shell_volume(source): 
     720    # type: (List[str]) -> bool 
     721    """ 
     722    Return True if C source defines "double shell_volume(". 
     723    """ 
     724    for code in source: 
     725        if _SHELL_VOLUME_PATTERN.search(code) is not None: 
     726            return True 
     727    return False 
    703728 
    704729def _add_source(source, code, path, lineno=1): 
     
    730755    # dispersion.  Need to be careful that necessary parameters are available 
    731756    # for computing volume even if we allow non-disperse volume parameters. 
    732  
    733757    partable = model_info.parameters 
    734758 
     
    743767    for path, code in user_code: 
    744768        _add_source(source, code, path) 
    745  
    746769    if model_info.c_code: 
    747770        _add_source(source, model_info.c_code, model_info.filename, 
     
    751774    q, qx, qy, qab, qa, qb, qc \ 
    752775        = [Parameter(name=v) for v in 'q qx qy qab qa qb qc'.split()] 
     776 
    753777    # Generate form_volume function, etc. from body only 
    754778    if isinstance(model_info.form_volume, str): 
    755779        pars = partable.form_volume_parameters 
    756780        source.append(_gen_fn(model_info, 'form_volume', pars)) 
     781    if isinstance(model_info.shell_volume, str): 
     782        pars = partable.form_volume_parameters 
     783        source.append(_gen_fn(model_info, 'shell_volume', pars)) 
    757784    if isinstance(model_info.Iq, str): 
    758785        pars = [q] + partable.iq_parameters 
     
    768795        source.append(_gen_fn(model_info, 'Iqabc', pars)) 
    769796 
     797    # Check for shell_volume in source 
     798    is_hollow = contains_shell_volume(source) 
     799 
    770800    # What kind of 2D model do we need?  Is it consistent with the parameters? 
    771801    xy_mode = find_xy_mode(source) 
     
    780810            logger.warn("oriented shapes should define Iqac or Iqabc") 
    781811        else: 
    782             raise ValueError("Expected function Iqac or Iqabc for oriented shape") 
     812            raise ValueError("Expected Iqac or Iqabc for oriented shape") 
    783813 
    784814    # Define the parameter table 
     
    789819    source.append("\\\n".join(p.as_definition() 
    790820                              for p in partable.kernel_parameters)) 
    791  
    792821    # Define the function calls 
     822    call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) 0.0" 
    793823    if partable.form_volume_parameters: 
    794824        refs = _call_pars("_v.", partable.form_volume_parameters) 
    795         call_volume = "#define CALL_VOLUME(_v) form_volume(%s)"%(",".join(refs)) 
     825        if is_hollow: 
     826            call_volume = ( 
     827                "#define CALL_VOLUME(_form, _shell, _v) " 
     828                "do { _form = form_volume(%s); _shell = shell_volume(%s); } " 
     829                "while (0)") % ((",".join(refs),)*2) 
     830        else: 
     831            call_volume = ( 
     832                "#define CALL_VOLUME(_form, _shell, _v) " 
     833                "do { _form = _shell = form_volume(%s); } " 
     834                "while (0)") % (",".join(refs)) 
     835        if model_info.effective_radius_type: 
     836            call_effective_radius = ( 
     837                "#define CALL_EFFECTIVE_RADIUS(_mode, _v) " 
     838                "effective_radius(_mode, %s)") % (",".join(refs)) 
    796839    else: 
    797840        # Model doesn't have volume.  We could make the kernel run a little 
    798841        # faster by not using/transferring the volume normalizations, but 
    799842        # the ifdef's reduce readability more than is worthwhile. 
    800         call_volume = "#define CALL_VOLUME(v) 1.0" 
     843        call_volume = ( 
     844            "#define CALL_VOLUME(_form, _shell, _v) " 
     845            "do { _form = _shell = 1.0; } while (0)") 
    801846    source.append(call_volume) 
    802  
     847    source.append(call_effective_radius) 
    803848    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 
     849 
     850    if model_info.have_Fq: 
     851        pars = ",".join(["_q", "&_F1", "&_F2",] + model_refs) 
     852        call_iq = "#define CALL_FQ(_q, _F1, _F2, _v) Fq(%s)" % pars 
     853        clear_iq = "#undef CALL_FQ" 
     854    else: 
     855        pars = ",".join(["_q"] + model_refs) 
     856        call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars 
     857        clear_iq = "#undef CALL_IQ" 
    806858    if xy_mode == 'qabc': 
    807859        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs) 
     
    812864        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqac(%s)" % pars 
    813865        clear_iqxy = "#undef CALL_IQ_AC" 
    814     elif xy_mode == 'qa': 
     866    elif xy_mode == 'qa' and not model_info.have_Fq: 
    815867        pars = ",".join(["_qa"] + model_refs) 
    816868        call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars 
    817869        clear_iqxy = "#undef CALL_IQ_A" 
     870    elif xy_mode == 'qa' and model_info.have_Fq: 
     871        pars = ",".join(["_qa", "&_F1", "&_F2",] + model_refs) 
     872        # Note: uses rare C construction (expr1, expr2) which computes 
     873        # expr1 then expr2 and evaluates to expr2.  This allows us to 
     874        # leave it looking like a function even though it is returning 
     875        # its values by reference. 
     876        call_iqxy = "#define CALL_FQ_A(_qa,_F1,_F2,_v) (Fq(%s),_F2)" % pars 
     877        clear_iqxy = "#undef CALL_FQ_A" 
    818878    elif xy_mode == 'qxy': 
    819879        orientation_refs = _call_pars("_v.", partable.orientation_parameters) 
     
    831891    magpars = [k-2 for k, p in enumerate(partable.call_parameters) 
    832892               if p.type == 'sld'] 
    833  
    834893    # Fill in definitions for numbers of parameters 
    835894    source.append("#define MAX_PD %s"%partable.max_pd) 
     
    839898    source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars)) 
    840899    source.append("#define PROJECTION %d"%PROJECTION) 
    841  
    842     # 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) 
    846     result = { 
    847         'dll': '\n'.join(source+dll[0]+dll[1]+dll[2]), 
    848         'opencl': '\n'.join(source+ocl[0]+ocl[1]+ocl[2]), 
    849     } 
    850  
     900    wrappers = _kernels(kernel_code, call_iq, clear_iq, 
     901                        call_iqxy, clear_iqxy, model_info.name) 
     902    code = '\n'.join(source + wrappers[0] + wrappers[1] + wrappers[2]) 
     903 
     904    # Note: Identical code for dll and opencl.  This may not always be the case 
     905    # so leave them as separate entries in the returned value. 
     906    result = {'dll': code, 'opencl': code} 
    851907    return result 
    852908 
    853909 
    854 def _kernels(kernel, call_iq, call_iqxy, clear_iqxy, name): 
     910def _kernels(kernel, call_iq, clear_iq, call_iqxy, clear_iqxy, name): 
    855911    # type: ([str,str], str, str, str) -> List[str] 
    856912    code = kernel[0] 
     
    862918        '#line 1 "%s Iq"' % path, 
    863919        code, 
    864         "#undef CALL_IQ", 
     920        clear_iq, 
    865921        "#undef KERNEL_NAME", 
    866922        ] 
     
    9651021    docs = model_info.docs if model_info.docs is not None else "" 
    9661022    docs = convert_section_titles_to_boldface(docs) 
    967     pars = make_partable(model_info.parameters.COMMON 
    968                          + model_info.parameters.kernel_parameters) 
     1023    if model_info.structure_factor: 
     1024        pars = model_info.parameters.kernel_parameters 
     1025    else: 
     1026        pars = (model_info.parameters.common_parameters 
     1027                + model_info.parameters.kernel_parameters) 
     1028    partable = make_partable(pars) 
    9691029    subst = dict(id=model_info.id.replace('_', '-'), 
    9701030                 name=model_info.name, 
    9711031                 title=model_info.title, 
    972                  parameters=pars, 
     1032                 parameters=partable, 
    9731033                 returns=Sq_units if model_info.structure_factor else Iq_units, 
    9741034                 docs=docs) 
  • sasmodels/gengauss.py

    rff31782 rb297ba9  
    11#!/usr/bin/env python 
     2""" 
     3Generate the Gauss-Legendre integration points and save them as a C file. 
     4""" 
    25from __future__ import division, print_function 
    3  
    4 import sys 
    56 
    67import numpy as np 
     
    89 
    910def gengauss(n, path): 
     11    """ 
     12    Save the Gauss-Legendre integration points for length *n* into file *path*. 
     13    """ 
    1014    z, w = leggauss(n) 
    1115 
  • sasmodels/guyou.py

    rb3703f5 rb297ba9  
    6868 
    6969def ellipticJi(u, v, m): 
     70    """Returns [sn, cn, dn](u + iv|m).""" 
    7071    scalar = np.isscalar(u) and np.isscalar(v) and np.isscalar(m) 
    7172    u, v, m = np.broadcast_arrays(u, v, m) 
     
    100101    ] 
    101102 
    102 # calculate F(phi+ipsi|m). 
    103 # see Abramowitz and Stegun, 17.4.11. 
    104103def ellipticFi(phi, psi, m): 
     104    """Returns F(phi+ipsi|m). See Abramowitz and Stegun, 17.4.11.""" 
    105105    if np.any(phi == 0): 
    106106        scalar = np.isscalar(phi) and np.isscalar(psi) and np.isscalar(m) 
     
    207207        plt.plot(x, y, 'g') 
    208208 
    209     for long in range(-limit, limit+1, step): 
    210         x, y = guyou_invert(scale*long, scale*long_line) 
     209    for longitude in range(-limit, limit+1, step): 
     210        x, y = guyou_invert(scale*longitude, scale*long_line) 
    211211        plt.plot(x, y, 'b') 
    212212    plt.xlabel('longitude') 
  • sasmodels/jitter.py

    rb3703f5 r4e28511  
    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 
    20 from numpy import pi, cos, sin, sqrt, exp, degrees, radians 
    21  
    22 def draw_beam(axes, view=(0, 0)): 
     53from numpy import pi, cos, sin, sqrt, exp, log, degrees, radians, arccos, arctan2 
     54 
     55# Too many complaints about variable names from pylint: 
     56#    a, b, c, u, v, x, y, z, dx, dy, dz, px, py, pz, R, Rx, Ry, Rz, ... 
     57# pylint: disable=invalid-name 
     58 
     59def draw_beam(axes, view=(0, 0), alpha=0.5, steps=25): 
    2360    """ 
    2461    Draw the beam going from source at (0, 0, 1) to detector at (0, 0, -1) 
    2562    """ 
    2663    #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 
    30     u = np.linspace(0, 2 * np.pi, steps) 
    31     v = np.linspace(-1, 1, steps) 
     64    #axes.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=alpha) 
     65 
     66    u = np.linspace(0, 2 * pi, steps) 
     67    v = np.linspace(-1, 1, 2) 
    3268 
    3369    r = 0.02 
    34     x = r*np.outer(np.cos(u), np.ones_like(v)) 
    35     y = r*np.outer(np.sin(u), np.ones_like(v)) 
     70    x = r*np.outer(cos(u), np.ones_like(v)) 
     71    y = r*np.outer(sin(u), np.ones_like(v)) 
    3672    z = 1.3*np.outer(np.ones_like(u), v) 
    3773 
     
    4177    points = Rz(phi)*Ry(theta)*points 
    4278    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) 
     79    axes.plot_surface(x, y, z, color='yellow', alpha=alpha) 
     80 
     81    # TODO: draw endcaps on beam 
     82    ## Drawing tiny balls on the end will work 
     83    #draw_sphere(axes, radius=0.02, center=(0, 0, 1.3), color='yellow', alpha=alpha) 
     84    #draw_sphere(axes, radius=0.02, center=(0, 0, -1.3), color='yellow', alpha=alpha) 
     85    ## The following does not work 
     86    #triangles = [(0, i+1, i+2) for i in range(steps-2)] 
     87    #x_cap, y_cap = x[:, 0], y[:, 0] 
     88    #for z_cap in z[:, 0], z[:, -1]: 
     89    #    axes.plot_trisurf(x_cap, y_cap, z_cap, triangles, 
     90    #                      color='yellow', alpha=alpha) 
     91 
    4592 
    4693def draw_ellipsoid(axes, size, view, jitter, steps=25, alpha=1): 
    4794    """Draw an ellipsoid.""" 
    4895    a, b, c = size 
    49     u = np.linspace(0, 2 * np.pi, steps) 
    50     v = np.linspace(0, np.pi, steps) 
    51     x = a*np.outer(np.cos(u), np.sin(v)) 
    52     y = b*np.outer(np.sin(u), np.sin(v)) 
    53     z = c*np.outer(np.ones_like(u), np.cos(v)) 
     96    u = np.linspace(0, 2 * pi, steps) 
     97    v = np.linspace(0, pi, steps) 
     98    x = a*np.outer(cos(u), sin(v)) 
     99    y = b*np.outer(sin(u), sin(v)) 
     100    z = c*np.outer(np.ones_like(u), cos(v)) 
    54101    x, y, z = transform_xyz(view, jitter, x, y, z) 
    55102 
    56     axes.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha) 
     103    axes.plot_surface(x, y, z, color='w', alpha=alpha) 
    57104 
    58105    draw_labels(axes, view, jitter, [ 
     
    67114def draw_sc(axes, size, view, jitter, steps=None, alpha=1): 
    68115    """Draw points for simple cubic paracrystal""" 
     116    # pylint: disable=unused-argument 
    69117    atoms = _build_sc() 
    70118    _draw_crystal(axes, size, view, jitter, atoms=atoms) 
     
    72120def draw_fcc(axes, size, view, jitter, steps=None, alpha=1): 
    73121    """Draw points for face-centered cubic paracrystal""" 
     122    # pylint: disable=unused-argument 
    74123    # Build the simple cubic crystal 
    75124    atoms = _build_sc() 
     
    87136def draw_bcc(axes, size, view, jitter, steps=None, alpha=1): 
    88137    """Draw points for body-centered cubic paracrystal""" 
     138    # pylint: disable=unused-argument 
    89139    # Build the simple cubic crystal 
    90140    atoms = _build_sc() 
     
    123173    return atoms 
    124174 
    125 def draw_parallelepiped(axes, size, view, jitter, steps=None, alpha=1): 
    126     """Draw a parallelepiped.""" 
     175def draw_box(axes, size, view): 
     176    """Draw a wireframe box at a particular view.""" 
     177    a, b, c = size 
     178    x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1]) 
     179    y = b*np.array([+1, +1, -1, -1, +1, +1, -1, -1]) 
     180    z = c*np.array([+1, +1, +1, +1, -1, -1, -1, -1]) 
     181    x, y, z = transform_xyz(view, None, x, y, z) 
     182    def _draw(i, j): 
     183        axes.plot([x[i], x[j]], [y[i], y[j]], [z[i], z[j]], color='black') 
     184    _draw(0, 1) 
     185    _draw(0, 2) 
     186    _draw(0, 3) 
     187    _draw(7, 4) 
     188    _draw(7, 5) 
     189    _draw(7, 6) 
     190 
     191def draw_parallelepiped(axes, size, view, jitter, steps=None, 
     192                        color=(0.6, 1.0, 0.6), alpha=1): 
     193    """Draw a parallelepiped surface, with view and jitter.""" 
     194    # pylint: disable=unused-argument 
    127195    a, b, c = size 
    128196    x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1]) 
     
    141209 
    142210    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. 
     211    axes.plot_trisurf(x, y, triangles=tri, Z=z, color=color, alpha=alpha, 
     212                      linewidth=0) 
     213 
     214    # Colour the c+ face of the box. 
    146215    # Since I can't control face color, instead draw a thin box situated just 
    147216    # in front of the "c+" face.  Use the c face so that rotations about psi 
    148217    # rotate that face. 
    149     if 1: 
     218    if 0: # pylint: disable=using-constant-test 
     219        color = (1, 0.6, 0.6)  # pink 
    150220        x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1]) 
    151221        y = b*np.array([+1, +1, -1, -1, +1, +1, -1, -1]) 
    152222        z = c*np.array([+1, +1, +1, +1, -1, -1, -1, -1]) 
    153223        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) 
     224        axes.plot_trisurf(x, y, triangles=tri, Z=z, color=color, alpha=alpha) 
    155225 
    156226    draw_labels(axes, view, jitter, [ 
     
    163233    ]) 
    164234 
    165 def draw_sphere(axes, radius=10., steps=100): 
     235def draw_sphere(axes, radius=1.0, steps=25, 
     236                center=(0, 0, 0), color='w', alpha=1.): 
    166237    """Draw a sphere""" 
    167     u = np.linspace(0, 2 * np.pi, steps) 
    168     v = np.linspace(0, np.pi, steps) 
    169  
    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): 
     238    u = np.linspace(0, 2 * pi, steps) 
     239    v = np.linspace(0, pi, steps) 
     240 
     241    x = radius * np.outer(cos(u), sin(v)) + center[0] 
     242    y = radius * np.outer(sin(u), sin(v)) + center[1] 
     243    z = radius * np.outer(np.ones(np.size(u)), cos(v)) + center[2] 
     244    axes.plot_surface(x, y, z, color=color, alpha=alpha) 
     245    #axes.plot_wireframe(x, y, z) 
     246 
     247def draw_axes(axes, origin=(-1, -1, -1), length=(2, 2, 2)): 
     248    """Draw wireframe axes lines, with given origin and length""" 
     249    x, y, z = origin 
     250    dx, dy, dz = length 
     251    axes.plot([x, x+dx], [y, y], [z, z], color='black') 
     252    axes.plot([x, x], [y, y+dy], [z, z], color='black') 
     253    axes.plot([x, x], [y, y], [z, z+dz], color='black') 
     254 
     255def draw_person_on_sphere(axes, view, height=0.5, radius=1.0): 
     256    """ 
     257    Draw a person on the surface of a sphere. 
     258 
     259    *view* indicates (latitude, longitude, orientation) 
     260    """ 
     261    limb_offset = height * 0.05 
     262    head_radius = height * 0.10 
     263    head_height = height - head_radius 
     264    neck_length = head_radius * 0.50 
     265    shoulder_height = height - 2*head_radius - neck_length 
     266    torso_length = shoulder_height * 0.55 
     267    torso_radius = torso_length * 0.30 
     268    leg_length = shoulder_height - torso_length 
     269    arm_length = torso_length * 0.90 
     270 
     271    def _draw_part(y, z): 
     272        x = np.zeros_like(y) 
     273        xp, yp, zp = transform_xyz(view, None, x, y, z + radius) 
     274        axes.plot(xp, yp, zp, color='k') 
     275 
     276    # circle for head 
     277    u = np.linspace(0, 2 * pi, 40) 
     278    y = head_radius * cos(u) 
     279    z = head_radius * sin(u) + head_height 
     280    _draw_part(y, z) 
     281 
     282    # rectangle for body 
     283    y = np.array([-torso_radius, torso_radius, torso_radius, -torso_radius, -torso_radius]) 
     284    z = np.array([0., 0, torso_length, torso_length, 0]) + leg_length 
     285    _draw_part(y, z) 
     286 
     287    # arms 
     288    y = np.array([-torso_radius - limb_offset, -torso_radius - limb_offset, -torso_radius]) 
     289    z = np.array([shoulder_height - arm_length, shoulder_height, shoulder_height]) 
     290    _draw_part(y, z) 
     291    _draw_part(-y, z)  # pylint: disable=invalid-unary-operand-type 
     292 
     293    # legs 
     294    y = np.array([-torso_radius + limb_offset, -torso_radius + limb_offset]) 
     295    z = np.array([0, leg_length]) 
     296    _draw_part(y, z) 
     297    _draw_part(-y, z)  # pylint: disable=invalid-unary-operand-type 
     298 
     299    limits = [-radius-height, radius+height] 
     300    axes.set_xlim(limits) 
     301    axes.set_ylim(limits) 
     302    axes.set_zlim(limits) 
     303    axes.set_axis_off() 
     304 
     305def draw_jitter(axes, view, jitter, dist='gaussian', 
     306                size=(0.1, 0.4, 1.0), 
     307                draw_shape=draw_parallelepiped, 
     308                projection='equirectangular', 
     309                alpha=0.8, 
     310                views=None): 
    177311    """ 
    178312    Represent jitter as a set of shapes at different orientations. 
    179313    """ 
     314    project, project_weight = get_projection(projection) 
     315 
    180316    # set max diagonal to 0.95 
    181317    scale = 0.95/sqrt(sum(v**2 for v in size)) 
    182318    size = tuple(scale*v for v in size) 
    183319 
    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     ] 
    215320    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) 
     321    base = {'gaussian':3, 'rectangle':sqrt(3), 'uniform':1}[dist] 
     322    def _steps(delta): 
     323        if views is None: 
     324            n = max(3, min(25, 2*int(base*delta/5))) 
     325        else: 
     326            n = views 
     327        return base*delta*np.linspace(-1, 1, n) if delta > 0 else [0.] 
     328    for theta in _steps(dtheta): 
     329        for phi in _steps(dphi): 
     330            for psi in _steps(dpsi): 
     331                w = project_weight(theta, phi, 1.0, 1.0) 
     332                if w > 0: 
     333                    dview = project(theta, phi, psi) 
     334                    draw_shape(axes, size, view, dview, alpha=alpha) 
    227335    for v in 'xyz': 
    228336        a, b, c = size 
    229         lim = np.sqrt(a**2 + b**2 + c**2) 
     337        lim = sqrt(a**2 + b**2 + c**2) 
    230338        getattr(axes, 'set_'+v+'lim')([-lim, lim]) 
    231         getattr(axes, v+'axis').label.set_text(v) 
     339        #getattr(axes, v+'axis').label.set_text(v) 
    232340 
    233341PROJECTIONS = [ 
     
    237345    'azimuthal_equal_area', 
    238346] 
    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  
     347def get_projection(projection): 
     348 
     349    """ 
    245350    jitter projections 
    246351    <https://en.wikipedia.org/wiki/List_of_map_projections> 
     
    296401        Should allow free movement in theta, but phi is distorted. 
    297402    """ 
     403    # pylint: disable=unused-argument 
    298404    # TODO: try Kent distribution instead of a gaussian warped by projection 
    299405 
    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  
    313406    if projection == 'equirectangular':  #define PROJECTION 1 
    314         def _rotate(theta_i, phi_j): 
    315             return Rx(phi_j)*Ry(theta_i) 
     407        def _project(theta_i, phi_j, psi): 
     408            latitude, longitude = theta_i, phi_j 
     409            return latitude, longitude, psi, 'xyz' 
     410            #return Rx(phi_j)*Ry(theta_i) 
    316411        def _weight(theta_i, phi_j, w_i, w_j): 
    317412            return w_i*w_j*abs(cos(radians(theta_i))) 
    318413    elif projection == 'sinusoidal':  #define PROJECTION 2 
    319         def _rotate(theta_i, phi_j): 
     414        def _project(theta_i, phi_j, psi): 
    320415            latitude = theta_i 
    321416            scale = cos(radians(latitude)) 
    322417            longitude = phi_j/scale if abs(phi_j) < abs(scale)*180 else 0 
    323418            #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): 
     419            return latitude, longitude, psi, 'xyz' 
     420            #return Rx(longitude)*Ry(latitude) 
     421        def _project(theta_i, phi_j, w_i, w_j): 
    326422            latitude = theta_i 
    327423            scale = cos(radians(latitude)) 
     
    329425            return active*w_i*w_j 
    330426    elif projection == 'guyou':  #define PROJECTION 3  (eventually?) 
    331         def _rotate(theta_i, phi_j): 
     427        def _project(theta_i, phi_j, psi): 
    332428            from .guyou import guyou_invert 
    333429            #latitude, longitude = guyou_invert([theta_i], [phi_j]) 
    334430            longitude, latitude = guyou_invert([phi_j], [theta_i]) 
    335             return Rx(longitude[0])*Ry(latitude[0]) 
     431            return latitude, longitude, psi, 'xyz' 
     432            #return Rx(longitude[0])*Ry(latitude[0]) 
    336433        def _weight(theta_i, phi_j, w_i, w_j): 
    337434            return w_i*w_j 
    338     elif projection == 'azimuthal_equidistance':  # Note: Rz Ry, not Rx Ry 
    339         def _rotate(theta_i, phi_j): 
     435    elif projection == 'azimuthal_equidistance': 
     436        # Note that calculates angles for Rz Ry rather than Rx Ry 
     437        def _project(theta_i, phi_j, psi): 
    340438            latitude = sqrt(theta_i**2 + phi_j**2) 
    341             longitude = degrees(np.arctan2(phi_j, theta_i)) 
     439            longitude = degrees(arctan2(phi_j, theta_i)) 
    342440            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
    343             return Rz(longitude)*Ry(latitude) 
     441            return latitude, longitude, psi-longitude, 'zyz' 
     442            #R = Rz(longitude)*Ry(latitude)*Rz(psi) 
     443            #return R_to_xyz(R) 
     444            #return Rz(longitude)*Ry(latitude) 
    344445        def _weight(theta_i, phi_j, w_i, w_j): 
    345446            # Weighting for each point comes from the integral: 
     
    375476            return weight*w_i*w_j if latitude < 180 else 0 
    376477    elif projection == 'azimuthal_equal_area': 
    377         def _rotate(theta_i, phi_j): 
     478        # Note that calculates angles for Rz Ry rather than Rx Ry 
     479        def _project(theta_i, phi_j, psi): 
    378480            radius = min(1, sqrt(theta_i**2 + phi_j**2)/180) 
    379             latitude = 180-degrees(2*np.arccos(radius)) 
    380             longitude = degrees(np.arctan2(phi_j, theta_i)) 
     481            latitude = 180-degrees(2*arccos(radius)) 
     482            longitude = degrees(arctan2(phi_j, theta_i)) 
    381483            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
    382             return Rz(longitude)*Ry(latitude) 
     484            return latitude, longitude, psi, 'zyz' 
     485            #R = Rz(longitude)*Ry(latitude)*Rz(psi) 
     486            #return R_to_xyz(R) 
     487            #return Rz(longitude)*Ry(latitude) 
    383488        def _weight(theta_i, phi_j, w_i, w_j): 
    384489            latitude = sqrt(theta_i**2 + phi_j**2) 
     
    388493        raise ValueError("unknown projection %r"%projection) 
    389494 
     495    return _project, _weight 
     496 
     497def R_to_xyz(R): 
     498    """ 
     499    Return phi, theta, psi Tait-Bryan angles corresponding to the given rotation matrix. 
     500 
     501    Extracting Euler Angles from a Rotation Matrix 
     502    Mike Day, Insomniac Games 
     503    https://d3cw3dd2w32x2b.cloudfront.net/wp-content/uploads/2012/07/euler-angles1.pdf 
     504    Based on: Shoemake’s "Euler Angle Conversion", Graphics Gems IV, pp.  222-229 
     505    """ 
     506    phi = arctan2(R[1, 2], R[2, 2]) 
     507    theta = arctan2(-R[0, 2], sqrt(R[0, 0]**2 + R[0, 1]**2)) 
     508    psi = arctan2(R[0, 1], R[0, 0]) 
     509    return degrees(phi), degrees(theta), degrees(psi) 
     510 
     511def draw_mesh(axes, view, jitter, radius=1.2, n=11, dist='gaussian', 
     512              projection='equirectangular'): 
     513    """ 
     514    Draw the dispersion mesh showing the theta-phi orientations at which 
     515    the model will be evaluated. 
     516    """ 
     517 
     518    _project, _weight = get_projection(projection) 
     519    def _rotate(theta, phi, z): 
     520        dview = _project(theta, phi, 0.) 
     521        if dview[3] == 'zyz': 
     522            return Rz(dview[1])*Ry(dview[0])*z 
     523        else:  # dview[3] == 'xyz': 
     524            return Rx(dview[1])*Ry(dview[0])*z 
     525 
     526 
     527    dist_x = np.linspace(-1, 1, n) 
     528    weights = np.ones_like(dist_x) 
     529    if dist == 'gaussian': 
     530        dist_x *= 3 
     531        weights = exp(-0.5*dist_x**2) 
     532    elif dist == 'rectangle': 
     533        # Note: uses sasmodels ridiculous definition of rectangle width 
     534        dist_x *= sqrt(3) 
     535    elif dist == 'uniform': 
     536        pass 
     537    else: 
     538        raise ValueError("expected dist to be gaussian, rectangle or uniform") 
     539 
    390540    # mesh in theta, phi formed by rotating z 
    391     dtheta, dphi, dpsi = jitter 
     541    dtheta, dphi, dpsi = jitter  # pylint: disable=unused-variable 
    392542    z = np.matrix([[0], [0], [radius]]) 
    393     points = np.hstack([_rotate(theta_i, phi_j)*z 
     543    points = np.hstack([_rotate(theta_i, phi_j, z) 
    394544                        for theta_i in dtheta*dist_x 
    395545                        for phi_j in dphi*dist_x]) 
     
    469619    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple. 
    470620    """ 
    471     dtheta, dphi, dpsi = jitter 
    472     points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points 
     621    if jitter is None: 
     622        return points 
     623    # Hack to deal with the fact that azimuthal_equidistance uses euler angles 
     624    if len(jitter) == 4: 
     625        dtheta, dphi, dpsi, _ = jitter 
     626        points = Rz(dphi)*Ry(dtheta)*Rz(dpsi)*points 
     627    else: 
     628        dtheta, dphi, dpsi = jitter 
     629        points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points 
    473630    return points 
    474631 
     
    480637    """ 
    481638    theta, phi, psi = view 
    482     points = Rz(phi)*Ry(theta)*Rz(psi)*points 
     639    points = Rz(phi)*Ry(theta)*Rz(psi)*points # viewing angle 
     640    #points = Rz(psi)*Ry(pi/2-theta)*Rz(phi)*points # 1-D integration angles 
     641    #points = Rx(phi)*Ry(theta)*Rz(psi)*points  # angular dispersion angle 
    483642    return points 
     643 
     644def orient_relative_to_beam_quaternion(view, points): 
     645    """ 
     646    Apply the view transform to a set of points. 
     647 
     648    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple. 
     649 
     650    This variant uses quaternions rather than rotation matrices for the 
     651    computation.  It works but it is not used because it doesn't solve 
     652    any problems.  The challenge of mapping theta/phi/psi to SO(3) does 
     653    not disappear by calculating the transform differently. 
     654    """ 
     655    theta, phi, psi = view 
     656    x, y, z = [1, 0, 0], [0, 1, 0], [0, 0, 1] 
     657    q = Quaternion(1, [0, 0, 0]) 
     658    ## Compose a rotation about the three axes by rotating 
     659    ## the unit vectors before applying the rotation. 
     660    #q = Quaternion.from_angle_axis(theta, q.rot(x)) * q 
     661    #q = Quaternion.from_angle_axis(phi, q.rot(y)) * q 
     662    #q = Quaternion.from_angle_axis(psi, q.rot(z)) * q 
     663    ## The above turns out to be equivalent to reversing 
     664    ## the order of application, so ignore it and use below. 
     665    q = q * Quaternion.from_angle_axis(theta, x) 
     666    q = q * Quaternion.from_angle_axis(phi, y) 
     667    q = q * Quaternion.from_angle_axis(psi, z) 
     668    ## Reverse the order by post-multiply rather than pre-multiply 
     669    #q = Quaternion.from_angle_axis(theta, x) * q 
     670    #q = Quaternion.from_angle_axis(phi, y) * q 
     671    #q = Quaternion.from_angle_axis(psi, z) * q 
     672    #print("axes psi", q.rot(np.matrix([x, y, z]).T)) 
     673    return q.rot(points) 
     674#orient_relative_to_beam = orient_relative_to_beam_quaternion 
     675 
     676# === Quaterion class definition === BEGIN 
     677# Simple stand-alone quaternion class 
     678 
     679# Note: this code works but isn't unused since quaternions didn't solve the 
     680# representation problem.  Leave it here in case we want to revisit this later. 
     681 
     682#import numpy as np 
     683class Quaternion(object): 
     684    r""" 
     685    Quaternion(w, r) = w + ir[0] + jr[1] + kr[2] 
     686 
     687    Quaternion.from_angle_axis(theta, r) for a rotation of angle theta about 
     688    an axis oriented toward the direction r.  This defines a unit quaternion, 
     689    normalizing $r$ to the unit vector $\hat r$, and setting quaternion 
     690    $Q = \cos \theta + \sin \theta \hat r$ 
     691 
     692    Quaternion objects can be multiplied, which applies a rotation about the 
     693    given axis, allowing composition of rotations without risk of gimbal lock. 
     694    The resulting quaternion is applied to a set of points using *Q.rot(v)*. 
     695    """ 
     696    def __init__(self, w, r): 
     697        self.w = w 
     698        self.r = np.asarray(r, dtype='d') 
     699 
     700    @staticmethod 
     701    def from_angle_axis(theta, r): 
     702        """Build quaternion as rotation theta about axis r""" 
     703        theta = np.radians(theta)/2 
     704        r = np.asarray(r) 
     705        w = np.cos(theta) 
     706        r = np.sin(theta)*r/np.dot(r, r) 
     707        return Quaternion(w, r) 
     708 
     709    def __mul__(self, other): 
     710        """Multiply quaterions""" 
     711        if isinstance(other, Quaternion): 
     712            w = self.w*other.w - np.dot(self.r, other.r) 
     713            r = self.w*other.r + other.w*self.r + np.cross(self.r, other.r) 
     714            return Quaternion(w, r) 
     715        raise NotImplementedError("Quaternion * non-quaternion not implemented") 
     716 
     717    def rot(self, v): 
     718        """Transform point *v* by quaternion""" 
     719        v = np.asarray(v).T 
     720        use_transpose = (v.shape[-1] != 3) 
     721        if use_transpose: 
     722            v = v.T 
     723        v = v + np.cross(2*self.r, np.cross(self.r, v) + self.w*v) 
     724        #v = v + 2*self.w*np.cross(self.r, v) + np.cross(2*self.r, np.cross(self.r, v)) 
     725        if use_transpose: 
     726            v = v.T 
     727        return v.T 
     728 
     729    def conj(self): 
     730        """Conjugate quaternion""" 
     731        return Quaternion(self.w, -self.r) 
     732 
     733    def inv(self): 
     734        """Inverse quaternion""" 
     735        return self.conj()/self.norm()**2 
     736 
     737    def norm(self): 
     738        """Quaternion length""" 
     739        return np.sqrt(self.w**2 + np.sum(self.r**2)) 
     740 
     741    def __str__(self): 
     742        return "%g%+gi%+gj%+gk"%(self.w, self.r[0], self.r[1], self.r[2]) 
     743 
     744def test_qrot(): 
     745    """Quaternion checks""" 
     746    # Define rotation of 60 degrees around an axis in y-z that is 60 degrees 
     747    # from y.  The rotation axis is determined by rotating the point [0, 1, 0] 
     748    # about x. 
     749    ax = Quaternion.from_angle_axis(60, [1, 0, 0]).rot([0, 1, 0]) 
     750    q = Quaternion.from_angle_axis(60, ax) 
     751    # Set the point to be rotated, and its expected rotated position. 
     752    p = [1, -1, 2] 
     753    target = [(10+4*np.sqrt(3))/8, (1+2*np.sqrt(3))/8, (14-3*np.sqrt(3))/8] 
     754    #print(q, q.rot(p) - target) 
     755    assert max(abs(q.rot(p) - target)) < 1e-14 
     756#test_qrot() 
     757#import sys; sys.exit() 
     758# === Quaterion class definition === END 
    484759 
    485760# translate between number of dimension of dispersity and the number of 
     
    503778    or the top of the range, depending on whether *mode* is 'central' or 'top'. 
    504779    """ 
    505     if portion == 1.0: 
    506         return data.min(), data.max() 
    507     elif mode == 'central': 
    508         data = np.sort(data.flatten()) 
    509         offset = int(portion*len(data)/2 + 0.5) 
    510         return data[offset], data[-offset] 
    511     elif mode == 'top': 
    512         data = np.sort(data.flatten()) 
    513         offset = int(portion*len(data) + 0.5) 
    514         return data[offset], data[-1] 
     780    if portion < 1.0: 
     781        if mode == 'central': 
     782            data = np.sort(data.flatten()) 
     783            offset = int(portion*len(data)/2 + 0.5) 
     784            return data[offset], data[-offset] 
     785        if mode == 'top': 
     786            data = np.sort(data.flatten()) 
     787            offset = int(portion*len(data) + 0.5) 
     788            return data[offset], data[-1] 
     789    # Default: full range 
     790    return data.min(), data.max() 
    515791 
    516792def draw_scattering(calculator, axes, view, jitter, dist='gaussian'): 
     
    543819 
    544820    # compute the pattern 
    545     qx, qy = calculator._data.x_bins, calculator._data.y_bins 
     821    qx, qy = calculator.qxqy 
    546822    Iqxy = calculator(**pars).reshape(len(qx), len(qy)) 
    547823 
    548824    # scale it and draw it 
    549     Iqxy = np.log(Iqxy) 
     825    Iqxy = log(Iqxy) 
    550826    if calculator.limits: 
    551827        # use limits from orientation (0,0,0) 
     
    555831        vmin = vmax*10**-7 
    556832        #vmin, vmax = clipped_range(Iqxy, portion=portion, mode='top') 
     833    #vmin, vmax = Iqxy.min(), Iqxy.max() 
    557834    #print("range",(vmin,vmax)) 
    558835    #qx, qy = np.meshgrid(qx, qy) 
    559     if 0: 
     836    if 0:  # pylint: disable=using-constant-test 
    560837        level = np.asarray(255*(Iqxy - vmin)/(vmax - vmin), 'i') 
    561838        level[level < 0] = 0 
     839        from matplotlib import pylab as plt 
    562840        colors = plt.get_cmap()(level) 
    563         axes.plot_surface(qx, qy, -1.1, rstride=1, cstride=1, facecolors=colors) 
    564     elif 1: 
     841        #from matplotlib import cm 
     842        #colors = cm.coolwarm(level) 
     843        #colors = cm.gist_yarg(level) 
     844        #colors = cm.Wistia(level) 
     845        colors[level <= 0, 3] = 0.  # set floor to transparent 
     846        x, y = np.meshgrid(qx/qx.max(), qy/qy.max()) 
     847        axes.plot_surface(x, y, -1.1*np.ones_like(x), facecolors=colors) 
     848    elif 1:  # pylint: disable=using-constant-test 
    565849        axes.contourf(qx/qx.max(), qy/qy.max(), Iqxy, zdir='z', offset=-1.1, 
    566850                      levels=np.linspace(vmin, vmax, 24)) 
     
    593877    calculator = DirectModel(data, model) 
    594878 
     879    # Remember the data axes so we can plot the results 
     880    calculator.qxqy = (q, q) 
     881 
    595882    # stuff the values for non-orientation parameters into the calculator 
    596883    calculator.pars = pars.copy() 
     
    600887    # under rotation or angular dispersion 
    601888    Iqxy = calculator(theta=0, phi=0, psi=0, **calculator.pars) 
    602     Iqxy = np.log(Iqxy) 
     889    Iqxy = log(Iqxy) 
    603890    vmin, vmax = clipped_range(Iqxy, 0.95, mode='top') 
    604891    calculator.limits = vmin, vmax+1 
     
    691978} 
    692979 
     980 
    693981def run(model_name='parallelepiped', size=(10, 40, 100), 
     982        view=(0, 0, 0), jitter=(0, 0, 0), 
    694983        dist='gaussian', mesh=30, 
    695984        projection='equirectangular'): 
     
    701990 
    702991    *size* gives the dimensions (a, b, c) of the shape. 
     992 
     993    *view* gives the initial view (theta, phi, psi) of the shape. 
     994 
     995    *view* gives the initial jitter (dtheta, dphi, dpsi) of the shape. 
    703996 
    704997    *dist* is the type of dispersition: gaussian, rectangle, or uniform. 
     
    7201013    calculator, size = select_calculator(model_name, n=150, size=size) 
    7211014    draw_shape = DRAW_SHAPES.get(model_name, draw_parallelepiped) 
     1015    #draw_shape = draw_fcc 
    7221016 
    7231017    ## uncomment to set an independent the colour range for every view 
     
    7251019    calculator.limits = None 
    7261020 
    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 
     1021    PLOT_ENGINE(calculator, draw_shape, size, view, jitter, dist, mesh, projection) 
     1022 
     1023def _mpl_plot(calculator, draw_shape, size, view, jitter, dist, mesh, projection): 
     1024    # Note: travis-ci does not support mpl_toolkits.mplot3d, but this shouldn't be 
     1025    # an issue since we are lazy-loading the package on a path that isn't tested. 
     1026    # Importing mplot3d adds projection='3d' option to subplot 
     1027    import mpl_toolkits.mplot3d  # pylint: disable=unused-variable 
     1028    import matplotlib as mpl 
     1029    import matplotlib.pyplot as plt 
     1030    from matplotlib.widgets import Slider 
    7331031 
    7341032    ## create the plot window 
     
    7461044        pass 
    7471045 
    748     axcolor = 'lightgoldenrodyellow' 
     1046    # CRUFT: use axisbg instead of facecolor for matplotlib<2 
     1047    facecolor_prop = 'facecolor' if mpl.__version__ > '2' else 'axisbg' 
     1048    props = {facecolor_prop: 'lightgoldenrodyellow'} 
    7491049 
    7501050    ## 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) 
     1051    axes_theta = plt.axes([0.05, 0.15, 0.50, 0.04], **props) 
     1052    axes_phi = plt.axes([0.05, 0.10, 0.50, 0.04], **props) 
     1053    axes_psi = plt.axes([0.05, 0.05, 0.50, 0.04], **props) 
     1054    stheta = Slider(axes_theta, u'Ξ', -90, 90, valinit=0) 
     1055    sphi = Slider(axes_phi, u'φ', -180, 180, valinit=0) 
     1056    spsi = Slider(axes_psi, u'ψ', -180, 180, valinit=0) 
     1057 
     1058    axes_dtheta = plt.axes([0.70, 0.15, 0.20, 0.04], **props) 
     1059    axes_dphi = plt.axes([0.70, 0.1, 0.20, 0.04], **props) 
     1060    axes_dpsi = plt.axes([0.70, 0.05, 0.20, 0.04], **props) 
     1061 
    7611062    # Note: using ridiculous definition of rectangle distribution, whose width 
    7621063    # in sasmodels is sqrt(3) times the given width.  Divide by sqrt(3) to keep 
    7631064    # the maximum width to 90. 
    7641065    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  
     1066    sdtheta = Slider(axes_dtheta, u'Δξ', 0, 2*dlimit, valinit=0) 
     1067    sdphi = Slider(axes_dphi, u'Δφ', 0, 2*dlimit, valinit=0) 
     1068    sdpsi = Slider(axes_dpsi, u'Δψ', 0, 2*dlimit, valinit=0) 
     1069 
     1070    ## initial view and jitter 
     1071    theta, phi, psi = view 
     1072    stheta.set_val(theta) 
     1073    sphi.set_val(phi) 
     1074    spsi.set_val(psi) 
     1075    dtheta, dphi, dpsi = jitter 
     1076    sdtheta.set_val(dtheta) 
     1077    sdphi.set_val(dphi) 
     1078    sdpsi.set_val(dpsi) 
    7691079 
    7701080    ## callback to draw the new view 
    771     def update(val, axis=None): 
     1081    def _update(val, axis=None): 
     1082        # pylint: disable=unused-argument 
    7721083        view = stheta.val, sphi.val, spsi.val 
    7731084        jitter = sdtheta.val, sdphi.val, sdpsi.val 
    7741085        # set small jitter as 0 if multiple pd dims 
    7751086        dims = sum(v > 0 for v in jitter) 
    776         limit = [0, 0, 0.5, 5][dims] 
     1087        limit = [0, 0.5, 5, 5][dims] 
    7771088        jitter = [0 if v < limit else v for v in jitter] 
    7781089        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)) 
     1090 
     1091        ## Visualize as person on globe 
     1092        #draw_sphere(axes, radius=0.5) 
     1093        #draw_person_on_sphere(axes, view, radius=0.5) 
     1094 
     1095        ## Move beam instead of shape 
     1096        #draw_beam(axes, -view[:2]) 
     1097        #draw_jitter(axes, (0,0,0), (0,0,0), views=3) 
     1098 
     1099        ## Move shape and draw scattering 
     1100        draw_beam(axes, (0, 0), alpha=1.) 
     1101        #draw_person_on_sphere(axes, view, radius=1.2, height=0.5) 
     1102        draw_jitter(axes, view, jitter, dist=dist, size=size, alpha=1., 
     1103                    draw_shape=draw_shape, projection=projection, views=3) 
    7821104        draw_mesh(axes, view, jitter, dist=dist, n=mesh, projection=projection) 
    7831105        draw_scattering(calculator, axes, view, jitter, dist=dist) 
     1106 
    7841107        plt.gcf().canvas.draw() 
    7851108 
    7861109    ## bind control widgets to view updater 
    787     stheta.on_changed(lambda v: update(v, 'theta')) 
    788     sphi.on_changed(lambda v: update(v, 'phi')) 
    789     spsi.on_changed(lambda v: update(v, 'psi')) 
    790     sdtheta.on_changed(lambda v: update(v, 'dtheta')) 
    791     sdphi.on_changed(lambda v: update(v, 'dphi')) 
    792     sdpsi.on_changed(lambda v: update(v, 'dpsi')) 
     1110    stheta.on_changed(lambda v: _update(v, 'theta')) 
     1111    sphi.on_changed(lambda v: _update(v, 'phi')) 
     1112    spsi.on_changed(lambda v: _update(v, 'psi')) 
     1113    sdtheta.on_changed(lambda v: _update(v, 'dtheta')) 
     1114    sdphi.on_changed(lambda v: _update(v, 'dphi')) 
     1115    sdpsi.on_changed(lambda v: _update(v, 'dpsi')) 
    7931116 
    7941117    ## initialize view 
    795     update(None, 'phi') 
     1118    _update(None, 'phi') 
    7961119 
    7971120    ## go interactive 
    7981121    plt.show() 
    7991122 
     1123 
     1124def map_colors(z, kw): 
     1125    """ 
     1126    Process matplotlib-style colour arguments. 
     1127 
     1128    Pulls 'cmap', 'alpha', 'vmin', and 'vmax' from th *kw* dictionary, setting 
     1129    the *kw['color']* to an RGB array.  These are ignored if 'c' or 'color' are 
     1130    set inside *kw*. 
     1131    """ 
     1132    from matplotlib import cm 
     1133 
     1134    cmap = kw.pop('cmap', cm.coolwarm) 
     1135    alpha = kw.pop('alpha', None) 
     1136    vmin = kw.pop('vmin', z.min()) 
     1137    vmax = kw.pop('vmax', z.max()) 
     1138    c = kw.pop('c', None) 
     1139    color = kw.pop('color', c) 
     1140    if color is None: 
     1141        znorm = ((z - vmin) / (vmax - vmin)).clip(0, 1) 
     1142        color = cmap(znorm) 
     1143    elif isinstance(color, np.ndarray) and color.shape == z.shape: 
     1144        color = cmap(color) 
     1145    if alpha is None: 
     1146        if isinstance(color, np.ndarray): 
     1147            color = color[..., :3] 
     1148    else: 
     1149        color[..., 3] = alpha 
     1150    kw['color'] = color 
     1151 
     1152def make_vec(*args): 
     1153    """Turn all elements of *args* into numpy arrays""" 
     1154    #return [np.asarray(v, 'd').flatten() for v in args] 
     1155    return [np.asarray(v, 'd') for v in args] 
     1156 
     1157def make_image(z, kw): 
     1158    """Convert numpy array *z* into a *PIL* RGB image.""" 
     1159    import PIL.Image 
     1160    from matplotlib import cm 
     1161 
     1162    cmap = kw.pop('cmap', cm.coolwarm) 
     1163 
     1164    znorm = (z-z.min())/z.ptp() 
     1165    c = cmap(znorm) 
     1166    c = c[..., :3] 
     1167    rgb = np.asarray(c*255, 'u1') 
     1168    image = PIL.Image.fromarray(rgb, mode='RGB') 
     1169    return image 
     1170 
     1171 
     1172_IPV_MARKERS = { 
     1173    'o': 'sphere', 
     1174} 
     1175_IPV_COLORS = { 
     1176    'w': 'white', 
     1177    'k': 'black', 
     1178    'c': 'cyan', 
     1179    'm': 'magenta', 
     1180    'y': 'yellow', 
     1181    'r': 'red', 
     1182    'g': 'green', 
     1183    'b': 'blue', 
     1184} 
     1185def _ipv_fix_color(kw): 
     1186    alpha = kw.pop('alpha', None) 
     1187    color = kw.get('color', None) 
     1188    if isinstance(color, str): 
     1189        color = _IPV_COLORS.get(color, color) 
     1190        kw['color'] = color 
     1191    if alpha is not None: 
     1192        color = kw['color'] 
     1193        #TODO: convert color to [r, g, b, a] if not already 
     1194        if isinstance(color, (tuple, list)): 
     1195            if len(color) == 3: 
     1196                color = (color[0], color[1], color[2], alpha) 
     1197            else: 
     1198                color = (color[0], color[1], color[2], alpha*color[3]) 
     1199            color = np.array(color) 
     1200        if isinstance(color, np.ndarray) and color.shape[-1] == 4: 
     1201            color[..., 3] = alpha 
     1202            kw['color'] = color 
     1203 
     1204def _ipv_set_transparency(kw, obj): 
     1205    color = kw.get('color', None) 
     1206    if (isinstance(color, np.ndarray) 
     1207            and color.shape[-1] == 4 
     1208            and (color[..., 3] != 1.0).any()): 
     1209        obj.material.transparent = True 
     1210        obj.material.side = "FrontSide" 
     1211 
     1212def ipv_axes(): 
     1213    """ 
     1214    Build a matplotlib style Axes interface for ipyvolume 
     1215    """ 
     1216    import ipyvolume as ipv 
     1217 
     1218    class Axes(object): 
     1219        """ 
     1220        Matplotlib Axes3D style interface to ipyvolume renderer. 
     1221        """ 
     1222        # pylint: disable=no-self-use,no-init 
     1223        # transparency can be achieved by setting the following: 
     1224        #    mesh.color = [r, g, b, alpha] 
     1225        #    mesh.material.transparent = True 
     1226        #    mesh.material.side = "FrontSide" 
     1227        # smooth(ish) rotation can be achieved by setting: 
     1228        #    slide.continuous_update = True 
     1229        #    figure.animation = 0. 
     1230        #    mesh.material.x = x 
     1231        #    mesh.material.y = y 
     1232        #    mesh.material.z = z 
     1233        # maybe need to synchronize update of x/y/z to avoid shimmy when moving 
     1234        def plot(self, x, y, z, **kw): 
     1235            """mpl style plot interface for ipyvolume""" 
     1236            _ipv_fix_color(kw) 
     1237            x, y, z = make_vec(x, y, z) 
     1238            ipv.plot(x, y, z, **kw) 
     1239        def plot_surface(self, x, y, z, **kw): 
     1240            """mpl style plot_surface interface for ipyvolume""" 
     1241            facecolors = kw.pop('facecolors', None) 
     1242            if facecolors is not None: 
     1243                kw['color'] = facecolors 
     1244            _ipv_fix_color(kw) 
     1245            x, y, z = make_vec(x, y, z) 
     1246            h = ipv.plot_surface(x, y, z, **kw) 
     1247            _ipv_set_transparency(kw, h) 
     1248            #h.material.side = "DoubleSide" 
     1249            return h 
     1250        def plot_trisurf(self, x, y, triangles=None, Z=None, **kw): 
     1251            """mpl style plot_trisurf interface for ipyvolume""" 
     1252            kw.pop('linewidth', None) 
     1253            _ipv_fix_color(kw) 
     1254            x, y, z = make_vec(x, y, Z) 
     1255            if triangles is not None: 
     1256                triangles = np.asarray(triangles) 
     1257            h = ipv.plot_trisurf(x, y, z, triangles=triangles, **kw) 
     1258            _ipv_set_transparency(kw, h) 
     1259            return h 
     1260        def scatter(self, x, y, z, **kw): 
     1261            """mpl style scatter interface for ipyvolume""" 
     1262            x, y, z = make_vec(x, y, z) 
     1263            map_colors(z, kw) 
     1264            marker = kw.get('marker', None) 
     1265            kw['marker'] = _IPV_MARKERS.get(marker, marker) 
     1266            h = ipv.scatter(x, y, z, **kw) 
     1267            _ipv_set_transparency(kw, h) 
     1268            return h 
     1269        def contourf(self, x, y, v, zdir='z', offset=0, levels=None, **kw): 
     1270            """mpl style contourf interface for ipyvolume""" 
     1271            # pylint: disable=unused-argument 
     1272            # Don't use contour for now (although we might want to later) 
     1273            self.pcolor(x, y, v, zdir='z', offset=offset, **kw) 
     1274        def pcolor(self, x, y, v, zdir='z', offset=0, **kw): 
     1275            """mpl style pcolor interface for ipyvolume""" 
     1276            # pylint: disable=unused-argument 
     1277            x, y, v = make_vec(x, y, v) 
     1278            image = make_image(v, kw) 
     1279            xmin, xmax = x.min(), x.max() 
     1280            ymin, ymax = y.min(), y.max() 
     1281            x = np.array([[xmin, xmax], [xmin, xmax]]) 
     1282            y = np.array([[ymin, ymin], [ymax, ymax]]) 
     1283            z = x*0 + offset 
     1284            u = np.array([[0., 1], [0, 1]]) 
     1285            v = np.array([[0., 0], [1, 1]]) 
     1286            h = ipv.plot_mesh(x, y, z, u=u, v=v, texture=image, wireframe=False) 
     1287            _ipv_set_transparency(kw, h) 
     1288            h.material.side = "DoubleSide" 
     1289            return h 
     1290        def text(self, *args, **kw): 
     1291            """mpl style text interface for ipyvolume""" 
     1292            pass 
     1293        def set_xlim(self, limits): 
     1294            """mpl style set_xlim interface for ipyvolume""" 
     1295            ipv.xlim(*limits) 
     1296        def set_ylim(self, limits): 
     1297            """mpl style set_ylim interface for ipyvolume""" 
     1298            ipv.ylim(*limits) 
     1299        def set_zlim(self, limits): 
     1300            """mpl style set_zlim interface for ipyvolume""" 
     1301            ipv.zlim(*limits) 
     1302        def set_axes_on(self): 
     1303            """mpl style set_axes_on interface for ipyvolume""" 
     1304            ipv.style.axis_on() 
     1305        def set_axis_off(self): 
     1306            """mpl style set_axes_off interface for ipyvolume""" 
     1307            ipv.style.axes_off() 
     1308    return Axes() 
     1309 
     1310def _ipv_plot(calculator, draw_shape, size, view, jitter, dist, mesh, projection): 
     1311    from IPython.display import display 
     1312    import ipywidgets as widgets 
     1313    import ipyvolume as ipv 
     1314 
     1315    axes = ipv_axes() 
     1316 
     1317    def _draw(view, jitter): 
     1318        camera = ipv.gcf().camera 
     1319        #print(ipv.gcf().__dict__.keys()) 
     1320        #print(dir(ipv.gcf())) 
     1321        ipv.figure(animation=0.)  # no animation when updating object mesh 
     1322 
     1323        # set small jitter as 0 if multiple pd dims 
     1324        dims = sum(v > 0 for v in jitter) 
     1325        limit = [0, 0.5, 5, 5][dims] 
     1326        jitter = [0 if v < limit else v for v in jitter] 
     1327 
     1328        ## Visualize as person on globe 
     1329        #draw_beam(axes, (0, 0)) 
     1330        #draw_sphere(axes, radius=0.5) 
     1331        #draw_person_on_sphere(axes, view, radius=0.5) 
     1332 
     1333        ## Move beam instead of shape 
     1334        #draw_beam(axes, view=(-view[0], -view[1])) 
     1335        #draw_jitter(axes, view=(0,0,0), jitter=(0,0,0)) 
     1336 
     1337        ## Move shape and draw scattering 
     1338        draw_beam(axes, (0, 0), steps=25) 
     1339        draw_jitter(axes, view, jitter, dist=dist, size=size, alpha=1.0, 
     1340                    draw_shape=draw_shape, projection=projection) 
     1341        draw_mesh(axes, view, jitter, dist=dist, n=mesh, radius=0.95, 
     1342                  projection=projection) 
     1343        draw_scattering(calculator, axes, view, jitter, dist=dist) 
     1344 
     1345        draw_axes(axes, origin=(-1, -1, -1.1)) 
     1346        ipv.style.box_off() 
     1347        ipv.style.axes_off() 
     1348        ipv.xyzlabel(" ", " ", " ") 
     1349 
     1350        ipv.gcf().camera = camera 
     1351        ipv.show() 
     1352 
     1353 
     1354    trange, prange = (-180., 180., 1.), (-180., 180., 1.) 
     1355    dtrange, dprange = (0., 180., 1.), (0., 180., 1.) 
     1356 
     1357    ## Super simple interfaca, but uses non-ascii variable namese 
     1358    # Ξ φ ψ Δξ Δφ Δψ 
     1359    #def update(**kw): 
     1360    #    view = kw['Ξ'], kw['φ'], kw['ψ'] 
     1361    #    jitter = kw['Δξ'], kw['Δφ'], kw['Δψ'] 
     1362    #    draw(view, jitter) 
     1363    #widgets.interact(update, Ξ=trange, φ=prange, ψ=prange, Δξ=dtrange, Δφ=dprange, Δψ=dprange) 
     1364 
     1365    def _update(theta, phi, psi, dtheta, dphi, dpsi): 
     1366        _draw(view=(theta, phi, psi), jitter=(dtheta, dphi, dpsi)) 
     1367 
     1368    def _slider(name, slice, init=0.): 
     1369        return widgets.FloatSlider( 
     1370            value=init, 
     1371            min=slice[0], 
     1372            max=slice[1], 
     1373            step=slice[2], 
     1374            description=name, 
     1375            disabled=False, 
     1376            #continuous_update=True, 
     1377            continuous_update=False, 
     1378            orientation='horizontal', 
     1379            readout=True, 
     1380            readout_format='.1f', 
     1381            ) 
     1382    theta = _slider(u'Ξ', trange, view[0]) 
     1383    phi = _slider(u'φ', prange, view[1]) 
     1384    psi = _slider(u'ψ', prange, view[2]) 
     1385    dtheta = _slider(u'Δξ', dtrange, jitter[0]) 
     1386    dphi = _slider(u'Δφ', dprange, jitter[1]) 
     1387    dpsi = _slider(u'Δψ', dprange, jitter[2]) 
     1388    fields = { 
     1389        'theta': theta, 'phi': phi, 'psi': psi, 
     1390        'dtheta': dtheta, 'dphi': dphi, 'dpsi': dpsi, 
     1391    } 
     1392    ui = widgets.HBox([ 
     1393        widgets.VBox([theta, phi, psi]), 
     1394        widgets.VBox([dtheta, dphi, dpsi]) 
     1395    ]) 
     1396 
     1397    out = widgets.interactive_output(_update, fields) 
     1398    display(ui, out) 
     1399 
     1400 
     1401_ENGINES = { 
     1402    "matplotlib": _mpl_plot, 
     1403    "mpl": _mpl_plot, 
     1404    #"plotly": _plotly_plot, 
     1405    "ipvolume": _ipv_plot, 
     1406    "ipv": _ipv_plot, 
     1407} 
     1408PLOT_ENGINE = _ENGINES["matplotlib"] 
     1409def set_plotter(name): 
     1410    """ 
     1411    Setting the plotting engine to matplotlib/ipyvolume or equivalently mpl/ipv. 
     1412    """ 
     1413    global PLOT_ENGINE 
     1414    PLOT_ENGINE = _ENGINES[name] 
     1415 
    8001416def main(): 
     1417    """ 
     1418    Command line interface to the jitter viewer. 
     1419    """ 
    8011420    parser = argparse.ArgumentParser( 
    8021421        description="Display jitter", 
     
    8081427    parser.add_argument('-s', '--size', type=str, default='10,40,100', 
    8091428                        help='a,b,c lengths') 
     1429    parser.add_argument('-v', '--view', type=str, default='0,0,0', 
     1430                        help='initial view angles') 
     1431    parser.add_argument('-j', '--jitter', type=str, default='0,0,0', 
     1432                        help='initial angular dispersion') 
    8101433    parser.add_argument('-d', '--distribution', choices=DISTRIBUTIONS, 
    8111434                        default=DISTRIBUTIONS[0], 
     
    8161439                        help='oriented shape') 
    8171440    opts = parser.parse_args() 
    818     size = tuple(int(v) for v in opts.size.split(',')) 
    819     run(opts.shape, size=size, 
     1441    size = tuple(float(v) for v in opts.size.split(',')) 
     1442    view = tuple(float(v) for v in opts.view.split(',')) 
     1443    jitter = tuple(float(v) for v in opts.jitter.split(',')) 
     1444    run(opts.shape, size=size, view=view, jitter=jitter, 
    8201445        mesh=opts.mesh, dist=opts.distribution, 
    8211446        projection=opts.projection) 
  • sasmodels/kernel.py

    r2d81cfe rb297ba9  
    2323# pylint: enable=unused-import 
    2424 
     25 
    2526class KernelModel(object): 
     27    """ 
     28    Model definition for the compute engine. 
     29    """ 
    2630    info = None  # type: ModelInfo 
    2731    dtype = None # type: np.dtype 
    2832    def make_kernel(self, q_vectors): 
    2933        # type: (List[np.ndarray]) -> "Kernel" 
     34        """ 
     35        Instantiate a kernel for evaluating the model at *q_vectors*. 
     36        """ 
    3037        raise NotImplementedError("need to implement make_kernel") 
    3138 
    3239    def release(self): 
    3340        # type: () -> None 
     41        """ 
     42        Free resources associated with the kernel. 
     43        """ 
    3444        pass 
    3545 
     46 
    3647class Kernel(object): 
    37     #: kernel dimension, either "1d" or "2d" 
     48    """ 
     49    Instantiated model for the compute engine, applied to a particular *q*. 
     50 
     51    Subclasses should define :meth:`_call_kernel` to evaluate the kernel over 
     52    its inputs. 
     53    """ 
     54    #: Kernel dimension, either "1d" or "2d". 
    3855    dim = None  # type: str 
     56    #: Model info. 
    3957    info = None  # type: ModelInfo 
    40     results = None # type: List[np.ndarray] 
     58    #: Numerical precision for the computation. 
    4159    dtype = None  # type: np.dtype 
     60    #: q values at which the kernel is to be evaluated 
     61    q_input = None  # type: Any 
     62    #: Place to hold result of :meth:`_call_kernel` for subclass. 
     63    result = None # type: np.ndarray 
    4264 
    43     def __call__(self, call_details, values, cutoff, magnetic): 
     65    def Iq(self, call_details, values, cutoff, magnetic): 
    4466        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    45         raise NotImplementedError("need to implement __call__") 
     67        r""" 
     68        Returns I(q) from the polydisperse average scattering. 
     69 
     70        .. math:: 
     71 
     72            I(q) = \text{scale} \cdot P(q) + \text{background} 
     73 
     74        With the correct choice of model and contrast, setting *scale* to 
     75        the volume fraction $V_f$ of particles should match the measured 
     76        absolute scattering.  Some models (e.g., vesicle) have volume fraction 
     77        built into the model, and do not need an additional scale. 
     78        """ 
     79        _, F2, _, shell_volume, _ = self.Fq(call_details, values, cutoff, 
     80                                            magnetic, effective_radius_type=0) 
     81        combined_scale = values[0]/shell_volume 
     82        background = values[1] 
     83        return combined_scale*F2 + background 
     84    __call__ = Iq 
     85 
     86    def Fq(self, call_details, values, cutoff, magnetic, 
     87           effective_radius_type=0): 
     88        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
     89        r""" 
     90        Returns <F(q)>, <F(q)^2>, effective radius, shell volume and 
     91        form:shell volume ratio.  The <F(q)> term may be None if the 
     92        form factor does not support direct computation of $F(q)$ 
     93 
     94        $P(q) = <F^2(q)>/<V>$ is used for structure factor calculations, 
     95 
     96        .. math:: 
     97 
     98            I(q) = \text{scale} \cdot P(q) \cdot S(q) + \text{background} 
     99 
     100        For the beta approximation, this becomes 
     101 
     102        .. math:: 
     103 
     104            I(q) = \text{scale} P (1 + <F>^2/<F^2> (S - 1)) + \text{background} 
     105                 = \text{scale}/<V> (<F^2> + <F>^2 (S - 1)) + \text{background} 
     106 
     107        $<F(q)>$ and $<F^2(q)>$ are averaged by polydispersity in shape 
     108        and orientation, with each configuration $x_k$ having form factor 
     109        $F(q, x_k)$, weight $w_k$ and volume $V_k$.  The result is: 
     110 
     111        .. math:: 
     112 
     113            P(q)=\frac{\sum w_k F^2(q, x_k) / \sum w_k}{\sum w_k V_k / \sum w_k} 
     114 
     115        The form factor itself is scaled by volume and contrast to compute the 
     116        total scattering.  This is then squared, and the volume weighted 
     117        F^2 is then normalized by volume F.  For a given density, the number 
     118        of scattering centers is assumed to scale linearly with volume.  Later 
     119        scaling the resulting $P(q)$ by the volume fraction of particles 
     120        gives the total scattering on an absolute scale. Most models 
     121        incorporate the volume fraction into the overall scale parameter.  An 
     122        exception is vesicle, which includes the volume fraction parameter in 
     123        the model itself, scaling $F$ by $\surd V_f$ so that the math for 
     124        the beta approximation works out. 
     125 
     126        By scaling $P(q)$ by total weight $\sum w_k$, there is no need to make 
     127        sure that the polydisperisity distributions normalize to one.  In 
     128        particular, any distibution values $x_k$ outside the valid domain 
     129        of $F$ will not be included, and the distribution will be implicitly 
     130        truncated.  This is controlled by the parameter limits defined in the 
     131        model (which truncate the distribution before calling the kernel) as 
     132        well as any region excluded using the *INVALID* macro defined within 
     133        the model itself. 
     134 
     135        The volume used in the polydispersity calculation is the form volume 
     136        for solid objects or the shell volume for hollow objects.  Shell 
     137        volume should be used within $F$ so that the normalizing scale 
     138        represents the volume fraction of the shell rather than the entire 
     139        form.  This corresponds to the volume fraction of shell-forming 
     140        material added to the solvent. 
     141 
     142        The calculation of $S$ requires the effective radius and the 
     143        volume fraction of the particles.  The model can have several 
     144        different ways to compute effective radius, with the 
     145        *effective_radius_type* parameter used to select amongst them.  The 
     146        volume fraction of particles should be determined from the total 
     147        volume fraction of the form, not just the shell volume fraction. 
     148        This makes a difference for hollow shapes, which need to scale 
     149        the volume fraction by the returned volume ratio when computing $S$. 
     150        For solid objects, the shell volume is set to the form volume so 
     151        this scale factor evaluates to one and so can be used for both 
     152        hollow and solid shapes. 
     153        """ 
     154        self._call_kernel(call_details, values, cutoff, magnetic, 
     155                          effective_radius_type) 
     156        #print("returned",self.q_input.q, self.result) 
     157        nout = 2 if self.info.have_Fq and self.dim == '1d' else 1 
     158        total_weight = self.result[nout*self.q_input.nq + 0] 
     159        # Note: total_weight = sum(weight > cutoff), with cutoff >= 0, so it 
     160        # is okay to test directly against zero.  If weight is zero then I(q), 
     161        # etc. must also be zero. 
     162        if total_weight == 0.: 
     163            total_weight = 1. 
     164        # Note: shell_volume == form_volume for solid objects 
     165        form_volume = self.result[nout*self.q_input.nq + 1]/total_weight 
     166        shell_volume = self.result[nout*self.q_input.nq + 2]/total_weight 
     167        effective_radius = self.result[nout*self.q_input.nq + 3]/total_weight 
     168        if shell_volume == 0.: 
     169            shell_volume = 1. 
     170        F1 = (self.result[1:nout*self.q_input.nq:nout]/total_weight 
     171              if nout == 2 else None) 
     172        F2 = self.result[0:nout*self.q_input.nq:nout]/total_weight 
     173        return F1, F2, effective_radius, shell_volume, form_volume/shell_volume 
    46174 
    47175    def release(self): 
    48176        # type: () -> None 
     177        """ 
     178        Free resources associated with the kernel instance. 
     179        """ 
    49180        pass 
     181 
     182    def _call_kernel(self, call_details, values, cutoff, magnetic, 
     183                     effective_radius_type): 
     184        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
     185        """ 
     186        Call the kernel.  Subclasses defining kernels for particular execution 
     187        engines need to provide an implementation for this. 
     188        """ 
     189        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

    r7c35fda 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. 
     
    8488  out_spin = clip(out_spin, 0.0, 1.0); 
    8589  // Previous version of this function took the square root of the weights, 
    86   // under the assumption that  
     90  // under the assumption that 
    8791  // 
    8892  //     w*I(q, rho1, rho2, ...) = I(q, sqrt(w)*rho1, sqrt(w)*rho2, ...) 
     
    188192    QACRotation *rotation, 
    189193    double qx, double qy, 
    190     double *qa_out, double *qc_out) 
     194    double *qab_out, double *qc_out) 
    191195{ 
     196    // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2 
    192197    const double dqc = rotation->R31*qx + rotation->R32*qy; 
    193     // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2 
    194     const double dqa = sqrt(-dqc*dqc + qx*qx + qy*qy); 
    195  
    196     *qa_out = dqa; 
     198    const double dqab_sq = -dqc*dqc + qx*qx + qy*qy; 
     199    //*qab_out = sqrt(fabs(dqab_sq)); 
     200    *qab_out = dqab_sq > 0.0 ? sqrt(dqab_sq) : 0.0; 
    197201    *qc_out = dqc; 
    198202} 
     
    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 r069743a  
    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 
     
    5153from __future__ import print_function 
    5254 
     55import sys 
    5356import os 
    5457import warnings 
     
    5659import time 
    5760 
     61try: 
     62    from time import perf_counter as clock 
     63except ImportError: # CRUFT: python < 3.3 
     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): 
     101    # type: (str) -> str 
    89102    """ 
    90103    Quote the path if it is not already quoted. 
     
    96109    return '"'+v+'"' if v and ' ' in v and not v[0] in "\"'-" else v 
    97110 
     111 
    98112def fix_pyopencl_include(): 
     113    # type: (None) -> None 
    99114    """ 
    100115    Monkey patch pyopencl to allow spaces in include file path. 
    101116    """ 
    102     import pyopencl as cl 
    103     if hasattr(cl, '_DEFAULT_INCLUDE_OPTIONS'): 
    104         cl._DEFAULT_INCLUDE_OPTIONS = [quote_path(v) for v in cl._DEFAULT_INCLUDE_OPTIONS] 
     117    # pylint: disable=protected-access 
     118    import pyopencl 
     119    if hasattr(pyopencl, '_DEFAULT_INCLUDE_OPTIONS'): 
     120        pyopencl._DEFAULT_INCLUDE_OPTIONS = [ 
     121            quote_path(v) for v in pyopencl._DEFAULT_INCLUDE_OPTIONS 
     122            ] 
     123 
    105124 
    106125if HAVE_OPENCL: 
     
    115134MAX_LOOPS = 2048 
    116135 
    117  
    118136# Pragmas for enable OpenCL features.  Be sure to pro