Changeset 765d025 in sasmodels


Ignore:
Timestamp:
Oct 30, 2018 12:05:27 PM (21 months ago)
Author:
Paul Kienzle <pkienzle@…>
Children:
646eeaa
Parents:
1662ebe (diff), aa8c6e0 (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 remote-tracking branch 'upstream/beta_approx'

Files:
51 added
1 deleted
127 edited

Legend:

Unmodified
Added
Removed
  • .travis.yml

    r335271e r5c36bf1  
    66    env: 
    77    - PY=2.7 
    8     - DEPLOY=True 
     8    #- DEPLOY=True 
    99  - os: linux 
    1010    env: 
  • 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

    r4f5afc9 rdf87acf  
    3939 
    4040.. math:: 
    41     -- &= ((1-u_i)(1-u_f))^{1/4} \\ 
    42     -+ &= ((1-u_i)(u_f))^{1/4} \\ 
    43     +- &= ((u_i)(1-u_f))^{1/4} \\ 
    44     ++ &= ((u_i)(u_f))^{1/4} 
     41    -- &= (1-u_i)(1-u_f) \\ 
     42    -+ &= (1-u_i)(u_f) \\ 
     43    +- &= (u_i)(1-u_f) \\ 
     44    ++ &= (u_i)(u_f) 
    4545 
    4646Ideally the experiment would measure the pure spin states independently and 
     
    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* 
     
    104104| 2015-05-02 Steve King 
    105105| 2017-11-15 Paul Kienzle 
     106| 2018-06-02 Adam Washington 
  • 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 raa8c6e0  
    291291 
    292292**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.** 
     293parameters in the user interface and the order of the parameters in Fq(), Iq(), 
     294Iqac(), Iqabc(), form_volume() and shell_volume(). 
     295And** *scale* **and** *background* **parameters are implicit to all models, 
     296so they do not need to be included in the parameter table.** 
    297297 
    298298- **"name"** is the name of the parameter shown on the FitPage. 
     
    363363    scattered intensity. 
    364364 
    365   - "volume" parameters are passed to Iq(), Iqac(), Iqabc() and form_volume(), 
    366     and have polydispersity loops generated automatically. 
     365  - "volume" parameters are passed to Fq(), Iq(), Iqac(), Iqabc(), form_volume() 
     366    and shell_volume(), and have polydispersity loops generated automatically. 
    367367 
    368368  - "orientation" parameters are not passed, but instead are combined with 
     
    424424appropriately smeared pattern. 
    425425 
     426Each .py file also contains a function:: 
     427 
     428        def random(): 
     429        ... 
     430 
     431This function provides a model-specific random parameter set which shows model 
     432features in the USANS to SANS range.  For example, core-shell sphere sets the 
     433outer radius of the sphere logarithmically in `[20, 20,000]`, which sets the Q 
     434value for the transition from flat to falling.  It then uses a beta distribution 
     435to set the percentage of the shape which is shell, giving a preference for very 
     436thin or very thick shells (but never 0% or 100%).  Using `-sets=10` in sascomp 
     437should show a reasonable variety of curves over the default sascomp q range. 
     438The parameter set is returned as a dictionary of `{parameter: value, ...}`. 
     439Any model parameters not included in the dictionary will default according to 
     440the code in the `_randomize_one()` function from sasmodels/compare.py. 
     441 
    426442Python Models 
    427443............. 
     
    476492used. 
    477493 
     494Hollow shapes, where the volume fraction of particle corresponds to the 
     495material in the shell rather than the volume enclosed by the shape, must 
     496also define a *shell_volume(par1, par2, ...)* function.  The parameters 
     497are the same as for *form_volume*.  The *I(q)* calculation should use 
     498*shell_volume* squared as its scale factor for the volume normalization. 
     499The structure factor calculation needs *form_volume* in order to properly 
     500scale the volume fraction parameter, so both functions are required for 
     501hollow shapes. 
     502 
     503Note: Pure python models do not yet support direct computation of the 
     504average of $F(q)$ and $F^2(q)$. 
     505 
    478506Embedded C Models 
    479507................. 
     
    487515This expands into the equivalent C code:: 
    488516 
    489     #include <math.h> 
    490517    double Iq(double q, double par1, double par2, ...); 
    491518    double Iq(double q, double par1, double par2, ...) 
     
    496523*form_volume* defines the volume of the shape. As in python models, it 
    497524includes only the volume parameters. 
     525 
     526*form_volume* defines the volume of the shell for hollow shapes. As in 
     527python models, it includes only the volume parameters. 
    498528 
    499529**source=['fn.c', ...]** includes the listed C source files in the 
     
    532562The INVALID define can go into *Iq*, or *c_code*, or an external C file 
    533563listed in *source*. 
     564 
     565Structure Factors 
     566................. 
     567 
     568Structure factor calculations may need the underlying $<F(q)>$ and $<F^2(q)>$ 
     569rather than $I(q)$.  This is used to compute $\beta = <F(q)>^2/<F^2(q)>$ in 
     570the decoupling approximation to the structure factor. 
     571 
     572Instead of defining the *Iq* function, models can define *Fq* as 
     573something like:: 
     574 
     575    double Fq(double q, double *F1, double *F2, double par1, double par2, ...); 
     576    double Fq(double q, double *F1, double *F2, double par1, double par2, ...) 
     577    { 
     578        // Polar integration loop over all orientations. 
     579        ... 
     580        *F1 = 1e-2 * total_F1 * contrast * volume; 
     581        *F2 = 1e-4 * total_F2 * square(contrast * volume); 
     582        return I(q, par1, par2, ...); 
     583    } 
     584 
     585If the volume fraction scale factor is built into the model (as occurs for 
     586the vesicle model, for example), then scale *F1* by $\surd V_f$ so that 
     587$\beta$ is computed correctly. 
     588 
     589Structure factor calculations are not yet supported for oriented shapes. 
     590 
     591Note: only available as a separate C file listed in *source*, or within 
     592a *c_code* block within the python model definition file. 
    534593 
    535594Oriented Shapes 
     
    685744    erf, erfc, tgamma, lgamma:  **do not use** 
    686745        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. 
     746        or inaccurate on some platforms. Use sas_erf, sas_erfc, sas_gamma 
     747        and sas_lgamma instead (see below). 
    689748 
    690749Some non-standard constants and functions are also provided: 
     
    753812        Gamma function sas_gamma\ $(x) = \Gamma(x)$. 
    754813 
    755         The standard math function, tgamma(x) is unstable for $x < 1$ 
     814        The standard math function, tgamma(x), is unstable for $x < 1$ 
    756815        on some platforms. 
    757816 
    758817        :code:`source = ["lib/sas_gamma.c", ...]` 
    759818        (`sas_gamma.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gamma.c>`_) 
     819 
     820    sas_gammaln(x): 
     821        log gamma function sas_gammaln\ $(x) = \log \Gamma(|x|)$. 
     822 
     823        The standard math function, lgamma(x), is incorrect for single 
     824        precision on some platforms. 
     825 
     826        :code:`source = ["lib/sas_gammainc.c", ...]` 
     827        (`sas_gammainc.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gammainc.c>`_) 
     828 
     829    sas_gammainc(a, x), sas_gammaincc(a, x): 
     830        Incomplete gamma function 
     831        sas_gammainc\ $(a, x) = \int_0^x t^{a-1}e^{-t}\,dt / \Gamma(a)$ 
     832        and complementary incomplete gamma function 
     833        sas_gammaincc\ $(a, x) = \int_x^\infty t^{a-1}e^{-t}\,dt / \Gamma(a)$ 
     834 
     835        :code:`source = ["lib/sas_gammainc.c", ...]` 
     836        (`sas_gammainc.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_gammainc.c>`_) 
    760837 
    761838    sas_erf(x), sas_erfc(x): 
     
    795872        If $n$ = 0 or 1, it uses sas_J0($x$) or sas_J1($x$), respectively. 
    796873 
     874        Warning: JN(n,x) can be very inaccurate (0.1%) for x not in [0.1, 100]. 
     875 
    797876        The standard math function jn(n, x) is not available on all platforms. 
    798877 
     
    803882        Sine integral Si\ $(x) = \int_0^x \tfrac{\sin t}{t}\,dt$. 
    804883 
     884        Warning: Si(x) can be very inaccurate (0.1%) for x in [0.1, 100]. 
     885 
    805886        This function uses Taylor series for small and large arguments: 
    806887 
    807         For large arguments, 
     888        For large arguments use the following Taylor series, 
    808889 
    809890        .. math:: 
     
    813894             - \frac{\sin(x)}{x}\left(\frac{1}{x} - \frac{3!}{x^3} + \frac{5!}{x^5} - \frac{7!}{x^7}\right) 
    814895 
    815         For small arguments, 
     896        For small arguments , 
    816897 
    817898        .. math:: 
     
    822903 
    823904        :code:`source = ["lib/Si.c", ...]` 
    824         (`Si.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/Si.c>`_) 
     905        (`Si.c <https://github.com/SasView/sasmodels/tree/master/sasmodels/models/lib/sas_Si.c>`_) 
    825906 
    826907    sas_3j1x_x(x): 
     
    9741055          "radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, 
    9751056         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.], 
     1057        [{"radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, 
     1058         0.1, None, None, 120., None, 1.],  # q, F, F^2, R_eff, V, form:shell 
     1059        [{"@S": "hardsphere"}, 0.1, None], 
    9781060    ] 
    9791061 
    9801062 
    981 **tests=[[{parameters}, q, result], ...]** is a list of lists. 
     1063**tests=[[{parameters}, q, Iq], ...]** is a list of lists. 
    9821064Each list is one test and contains, in order: 
    9831065 
     
    9911073- input and output values can themselves be lists if you have several 
    9921074  $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. 
     1075- for testing effective radius, volume and form:shell volume ratio, use the 
     1076  extended form of the tests results, with *None, None, R_eff, V, V_r* 
     1077  instead of *Iq*.  This calls the kernel *Fq* function instead of *Iq*. 
     1078- for testing F and F^2 (used for beta approximation) do the same as the 
     1079  effective radius test, but include values for the first two elements, 
     1080  $<F(q)>$ and $<F^2(q)>$. 
     1081- for testing interaction between form factor and structure factor, specify 
     1082  the structure factor name in the parameters as *{"@S": "name", ...}* with 
     1083  the remaining list of parameters defined by the *P@S* product model. 
    9961084 
    9971085.. _Test_Your_New_Model: 
  • doc/guide/scripting.rst

    r4aa5dce rbd7630d  
    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 
  • 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/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, 
  • explore/precision.py

    r2a7e20e raa8c6e0  
    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 two parameter functions 
     210A = 1 
     211def parse_extra_pars(): 
     212    global A 
     213 
     214    A_str = str(A) 
     215    pop = [] 
     216    for k, v in enumerate(sys.argv[1:]): 
     217        if v.startswith("A="): 
     218            A_str = v[2:] 
     219            pop.append(k+1) 
     220    if pop: 
     221        sys.argv = [v for k, v in enumerate(sys.argv) if k not in pop] 
     222        A = float(A_str) 
     223 
     224parse_extra_pars() 
     225 
    199226 
    200227# =============== FUNCTION DEFINITIONS ================ 
     
    299326) 
    300327add_function( 
     328    name="gammaln(x)", 
     329    mp_function=mp.loggamma, 
     330    np_function=scipy.special.gammaln, 
     331    ocl_function=make_ocl("return sas_gammaln(q);", "sas_gammaln", ["lib/sas_gammainc.c"]), 
     332    #ocl_function=make_ocl("return lgamma(q);", "sas_gammaln"), 
     333) 
     334add_function( 
     335    name="gammainc(x)", 
     336    mp_function=lambda x, a=A: mp.gammainc(a, a=0, b=x)/mp.gamma(a), 
     337    np_function=lambda x, a=A: scipy.special.gammainc(a, x), 
     338    ocl_function=make_ocl("return sas_gammainc(%.15g,q);"%A, "sas_gammainc", ["lib/sas_gammainc.c"]), 
     339) 
     340add_function( 
     341    name="gammaincc(x)", 
     342    mp_function=lambda x, a=A: mp.gammainc(a, a=x, b=mp.inf)/mp.gamma(a), 
     343    np_function=lambda x, a=A: scipy.special.gammaincc(a, x), 
     344    ocl_function=make_ocl("return sas_gammaincc(%.15g,q);"%A, "sas_gammaincc", ["lib/sas_gammainc.c"]), 
     345) 
     346add_function( 
    301347    name="erf(x)", 
    302348    mp_function=mp.erf, 
     
    357403) 
    358404add_function( 
    359     name="debye", 
     405    name="gauss_coil", 
    360406    mp_function=lambda x: 2*(mp.exp(-x**2) + x**2 - 1)/x**4, 
    361407    np_function=lambda x: 2*(np.expm1(-x**2) + x**2)/x**4, 
    362408    ocl_function=make_ocl(""" 
    363409    const double qsq = q*q; 
    364     if (qsq < 1.0) { // Pade approximation 
     410    // For double: use O(5) Pade with 0.5 cutoff (10 mad + 1 divide) 
     411    // For single: use O(7) Taylor with 0.8 cutoff (7 mad) 
     412    if (qsq < 0.0) { 
    365413        const double x = qsq; 
    366414        if (0) { // 0.36 single 
     
    372420            const double B1=3./8., B2=3./56., B3=1./336.; 
    373421            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 
     422        } else if (0) { // 1.0 for single, 0.25 for double 
    375423            // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}] 
    376424            const double A1=1./15., A2=1./60, A3=0., A4=1./75600.; 
     
    385433                  /(((((B5*x + B4)*x + B3)*x + B2)*x + B1)*x + 1.); 
    386434        } 
    387     } else if (qsq < 1.) { // Taylor series; 0.9 for single, 0.25 for double 
     435    } else if (qsq < 0.8) { 
    388436        const double x = qsq; 
    389437        const double C0 = +1.; 
     
    463511lanczos_gamma = """\ 
    464512    const double coeff[] = { 
    465             76.18009172947146,     -86.50532032941677, 
    466             24.01409824083091,     -1.231739572450155, 
     513            76.18009172947146, -86.50532032941677, 
     514            24.01409824083091, -1.231739572450155, 
    467515            0.1208650973866179e-2,-0.5395239384953e-5 
    468516            }; 
     
    475523""" 
    476524add_function( 
    477     name="log gamma(x)", 
     525    name="loggamma(x)", 
    478526    mp_function=mp.loggamma, 
    479527    np_function=scipy.special.gammaln, 
     
    599647 
    600648ALL_FUNCTIONS = set(FUNCTIONS.keys()) 
    601 ALL_FUNCTIONS.discard("loggamma")  # OCL version not ready yet 
     649ALL_FUNCTIONS.discard("loggamma")  # use cephes-based gammaln instead 
    602650ALL_FUNCTIONS.discard("3j1/x:taylor") 
    603651ALL_FUNCTIONS.discard("3j1/x:trig") 
     
    615663    -r indicates that the relative error should be plotted (default), 
    616664    -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: 
     665        log indicates log stepping in [10^-3, 10^5] (default) 
     666        logq indicates log stepping in [10^-4, 10^1] 
     667        linear indicates linear stepping in [1, 1000] 
     668        zoom indicates linear stepping in [1000, 1010] 
     669        neg indicates linear stepping in [-100.1, 100.1] 
     670        start:stop:n[:stepping] indicates an n-step plot in [start, stop] 
     671            or [10^start, 10^stop] if stepping is "log" (default n=400) 
     672Some functions (notably gammainc/gammaincc) have an additional parameter A 
     673which can be set from the command line as A=value.  Default is A=1. 
     674 
     675Name is one of: 
    623676    """+names) 
    624677    sys.exit(1) 
  • sasmodels/__init__.py

    re65c3ba ra1ec908  
    1414defining new models. 
    1515""" 
    16 __version__ = "0.97" 
     16__version__ = "0.98" 
    1717 
    1818def data_files(): 
  • sasmodels/autoc.py

    r15be191 r765d025  
    3030    'sas_erfc': ['lib/polevl.c', 'lib/sas_erf.c'], 
    3131    'sas_gamma': ['lib/sas_gamma.c'], 
     32    'sas_gammaln': ['lib/sas_gammainc.c'], 
     33    'sas_gammainc': ['lib/sas_gammainc.c'], 
     34    'sas_gammaincc': ['lib/sas_gammainc.c'], 
    3235    'sas_J0': ['lib/polevl.c', 'lib/sas_J0.c'], 
    3336    'sas_J1': ['lib/polevl.c', 'lib/sas_J1.c'], 
     
    4144    # type: ("ModelInfo", ModuleType) -> bool 
    4245    """ 
    43     convert Iq, Iqxy and form_volume to c 
     46    Convert Iq, Iqxy, form_volume, etc. to c 
     47 
     48    Returns list of warnings 
    4449    """ 
    4550    # Check if there is already C code 
     
    4752        return 
    4853 
    49     public_methods = "Iq", "Iqac", "Iqabc", "Iqxy", "form_volume" 
     54    public_methods = ("Iq", "Iqac", "Iqabc", "Iqxy",  
     55            "form_volume", "shell_volume", "effective_radius") 
    5056 
    5157    tagged = [] # type: List[str] 
     
    120126    # translate source 
    121127    ordered_code = [code[name] for name in py2c.ordered_dag(depends) if name in code] 
    122     functions = py2c.translate(ordered_code, constants) 
     128    functions, warnings = py2c.translate(ordered_code, constants) 
    123129    snippets.extend(functions) 
    124130 
     
    127133    info.c_code = "".join(snippets) 
    128134    info.Iq = info.Iqac = info.Iqabc = info.Iqxy = info.form_volume = None 
     135 
     136    return warnings 
     137 
  • sasmodels/compare.py

    r1fbadb2 r07646b6  
    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 
     
    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() 
     
    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', 
     
    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 r610ef23  
    165165    if version == (3, 1, 2): 
    166166        oldpars = _hand_convert_3_1_2_to_4_1(name, oldpars) 
     167    if version < (4, 2, 0): 
     168        oldpars = _rename_magnetic_pars(oldpars) 
    167169    return oldpars 
     170 
     171def _rename_magnetic_pars(pars): 
     172    """ 
     173    Change from M0:par to par_M0, etc. 
     174    """ 
     175    keys = list(pars.items()) 
     176    for k in keys: 
     177        if k.startswith('M0:'): 
     178            pars[k[3:]+'_M0'] = pars.pop(k) 
     179        elif k.startswith('mtheta:'): 
     180            pars[k[7:]+'_mtheta'] = pars.pop(k) 
     181        elif k.startswith('mphi:'): 
     182            pars[k[5:]+'_mphi'] = pars.pop(k) 
     183        elif k.startswith('up:'): 
     184            pars['up_'+k[3:]] = pars.pop(k) 
     185    return pars 
    168186 
    169187def _hand_convert_3_1_2_to_4_1(name, oldpars): 
  • sasmodels/core.py

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

    rd18d6dd r304c775  
    3131from . import resolution2d 
    3232from .details import make_kernel_args, dispersion_mesh 
     33from .modelinfo import DEFAULT_BACKGROUND 
    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  
     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, V_form/V_shell. 
     70    """ 
     71    R_eff_type = int(pars.pop('radius_effective_type', 1.0)) 
     72    mesh = get_mesh(calculator.info, pars, dim=calculator.dim, mono=mono) 
     73    #print("pars", list(zip(*mesh))[0]) 
     74    call_details, values, is_magnetic = make_kernel_args(calculator, mesh) 
     75    #print("values:", values) 
     76    return calculator.Fq(call_details, values, cutoff, is_magnetic, R_eff_type) 
    10277 
    10378def call_profile(model_info, **pars): 
     
    250225            qmax = getattr(data, 'qmax', np.inf) 
    251226            accuracy = getattr(data, 'accuracy', 'Low') 
    252             index = ~data.mask & (q >= qmin) & (q <= qmax) 
     227            index = (data.mask == 0) & (q >= qmin) & (q <= qmax) 
    253228            if data.data is not None: 
    254229                index &= ~np.isnan(data.data) 
     
    263238        elif self.data_type == 'Iq': 
    264239            index = (data.x >= data.qmin) & (data.x <= data.qmax) 
     240            mask = getattr(data, 'mask', None) 
     241            if mask is not None: 
     242                index &= (mask == 0) 
    265243            if data.y is not None: 
    266244                index &= ~np.isnan(data.y) 
     
    346324 
    347325        # Need to pull background out of resolution for multiple scattering 
    348         background = pars.get('background', 0.) 
     326        background = pars.get('background', DEFAULT_BACKGROUND) 
    349327        pars = pars.copy() 
    350328        pars['background'] = 0. 
     
    415393    model_name = sys.argv[1] 
    416394    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) 
     395    try: 
     396        values = [float(v) for v in call.split(',')] 
     397    except ValueError: 
     398        values = [] 
     399    if len(values) == 1: 
     400        q, = values 
     401        data = empty_data1D([q]) 
     402    elif len(values) == 2: 
     403        qx, qy = values 
     404        data = empty_data2D([qx], [qy]) 
    431405    else: 
    432         data = empty_data1D([0.001])  # Data not used in ER/VR 
     406        print("use q or qx,qy") 
     407        sys.exit(1) 
    433408 
    434409    model_info = load_model_info(model_name) 
     
    438413                for pair in sys.argv[3:] 
    439414                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]) 
     415    Iq = calculator(**pars) 
     416    print(Iq[0]) 
    447417 
    448418if __name__ == "__main__": 
  • sasmodels/generate.py

    rd86f0fc r39a06c9  
    2424    dimension, or 1.0 if no volume normalization is required. 
    2525 
    26     *ER(p1, p2, ...)* returns the effective radius of the form with 
    27     particular dimensions. 
    28  
    29     *VR(p1, p2, ...)* returns the volume ratio for core-shell style forms. 
     26    *shell_volume(p1, p2, ...)* returns the volume of the shell for forms 
     27    which are hollow. 
     28 
     29    *effective_radius(mode, p1, p2, ...)* returns the effective radius of 
     30    the form with particular dimensions.  Mode determines the type of 
     31    effective radius returned, with mode=1 for equivalent volume. 
    3032 
    3133    #define INVALID(v) (expr)  returns False if v.parameter is invalid 
     
    7274background.  This will work correctly even when polydispersity is off. 
    7375 
    74 *ER* and *VR* are python functions which operate on parameter vectors. 
    75 The constructor code will generate the necessary vectors for computing 
    76 them with the desired polydispersity. 
    7776The kernel module must set variables defining the kernel meta data: 
    7877 
     
    106105    create the OpenCL kernel functions.  The files defining the functions 
    107106    need to be listed before the files which use the functions. 
    108  
    109     *ER* is a python function defining the effective radius.  If it is 
    110     not present, the effective radius is 0. 
    111  
    112     *VR* is a python function defining the volume ratio.  If it is not 
    113     present, the volume ratio is 1. 
    114107 
    115108    *form_volume*, *Iq*, *Iqac*, *Iqabc* are strings containing 
     
    671664 
    672665 
    673 # type in IQXY pattern could be single, float, double, long double, ... 
    674 _IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ab?c|xy))\s*[(]", 
     666_IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ac|abc|xy))\s*[(]", 
    675667                           flags=re.MULTILINE) 
    676668def find_xy_mode(source): 
     
    701693    return 'qa' 
    702694 
     695# Note: not presently used.  Need to know whether Fq is available before 
     696# trying to compile the source.  Leave the code here in case we decide that 
     697# define have_Fq for each form factor is too tedious and error prone. 
     698_FQ_PATTERN = re.compile(r"(^|\s)void\s+Fq[(]", flags=re.MULTILINE) 
     699def contains_Fq(source): 
     700    # type: (List[str]) -> bool 
     701    """ 
     702    Return True if C source defines "void Fq(". 
     703    """ 
     704    for code in source: 
     705        m = _FQ_PATTERN.search(code) 
     706        if m is not None: 
     707            return True 
     708    return False 
     709 
     710_SHELL_VOLUME_PATTERN = re.compile(r"(^|\s)double\s+shell_volume[(]", flags=re.MULTILINE) 
     711def contains_shell_volume(source): 
     712    # type: (List[str]) -> bool 
     713    """ 
     714    Return True if C source defines "void Fq(". 
     715    """ 
     716    for code in source: 
     717        m = _SHELL_VOLUME_PATTERN.search(code) 
     718        if m is not None: 
     719            return True 
     720    return False 
    703721 
    704722def _add_source(source, code, path, lineno=1): 
     
    730748    # dispersion.  Need to be careful that necessary parameters are available 
    731749    # for computing volume even if we allow non-disperse volume parameters. 
    732  
    733750    partable = model_info.parameters 
    734751 
     
    743760    for path, code in user_code: 
    744761        _add_source(source, code, path) 
    745  
    746762    if model_info.c_code: 
    747763        _add_source(source, model_info.c_code, model_info.filename, 
     
    751767    q, qx, qy, qab, qa, qb, qc \ 
    752768        = [Parameter(name=v) for v in 'q qx qy qab qa qb qc'.split()] 
     769 
    753770    # Generate form_volume function, etc. from body only 
    754771    if isinstance(model_info.form_volume, str): 
    755772        pars = partable.form_volume_parameters 
    756773        source.append(_gen_fn(model_info, 'form_volume', pars)) 
     774    if isinstance(model_info.shell_volume, str): 
     775        pars = partable.form_volume_parameters 
     776        source.append(_gen_fn(model_info, 'shell_volume', pars)) 
    757777    if isinstance(model_info.Iq, str): 
    758778        pars = [q] + partable.iq_parameters 
     
    767787        pars = [qa, qb, qc] + partable.iq_parameters 
    768788        source.append(_gen_fn(model_info, 'Iqabc', pars)) 
     789 
     790    # Check for shell_volume in source 
     791    is_hollow = contains_shell_volume(source) 
    769792 
    770793    # What kind of 2D model do we need?  Is it consistent with the parameters? 
     
    789812    source.append("\\\n".join(p.as_definition() 
    790813                              for p in partable.kernel_parameters)) 
    791  
    792814    # Define the function calls 
     815    call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) 0.0" 
    793816    if partable.form_volume_parameters: 
    794817        refs = _call_pars("_v.", partable.form_volume_parameters) 
    795         call_volume = "#define CALL_VOLUME(_v) form_volume(%s)"%(",".join(refs)) 
     818        if is_hollow: 
     819            call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = form_volume(%s); _shell = shell_volume(%s); } while (0)"%((",".join(refs),)*2) 
     820        else: 
     821            call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = _shell = form_volume(%s); } while (0)"%(",".join(refs)) 
     822        if model_info.effective_radius_type: 
     823            call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) effective_radius(_mode, %s)"%(",".join(refs)) 
    796824    else: 
    797825        # Model doesn't have volume.  We could make the kernel run a little 
    798826        # faster by not using/transferring the volume normalizations, but 
    799827        # the ifdef's reduce readability more than is worthwhile. 
    800         call_volume = "#define CALL_VOLUME(v) 1.0" 
     828        call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = _shell = 1.0; } while (0)" 
    801829    source.append(call_volume) 
    802  
     830    source.append(call_effective_radius) 
    803831    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 
     832 
     833    if model_info.have_Fq: 
     834        pars = ",".join(["_q", "&_F1", "&_F2",] + model_refs) 
     835        call_iq = "#define CALL_FQ(_q, _F1, _F2, _v) Fq(%s)" % pars 
     836        clear_iq = "#undef CALL_FQ" 
     837    else: 
     838        pars = ",".join(["_q"] + model_refs) 
     839        call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars 
     840        clear_iq = "#undef CALL_IQ" 
    806841    if xy_mode == 'qabc': 
    807842        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs) 
     
    812847        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqac(%s)" % pars 
    813848        clear_iqxy = "#undef CALL_IQ_AC" 
    814     elif xy_mode == 'qa': 
     849    elif xy_mode == 'qa' and not model_info.have_Fq: 
    815850        pars = ",".join(["_qa"] + model_refs) 
    816851        call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars 
    817852        clear_iqxy = "#undef CALL_IQ_A" 
     853    elif xy_mode == 'qa' and model_info.have_Fq: 
     854        pars = ",".join(["_qa", "&_F1", "&_F2",] + model_refs) 
     855        # Note: uses rare C construction (expr1, expr2) which computes 
     856        # expr1 then expr2 and evaluates to expr2.  This allows us to 
     857        # leave it looking like a function even though it is returning 
     858        # its values by reference. 
     859        call_iqxy = "#define CALL_FQ_A(_qa,_F1,_F2,_v) (Fq(%s),_F2)" % pars 
     860        clear_iqxy = "#undef CALL_FQ_A" 
    818861    elif xy_mode == 'qxy': 
    819862        orientation_refs = _call_pars("_v.", partable.orientation_parameters) 
     
    831874    magpars = [k-2 for k, p in enumerate(partable.call_parameters) 
    832875               if p.type == 'sld'] 
    833  
    834876    # Fill in definitions for numbers of parameters 
    835877    source.append("#define MAX_PD %s"%partable.max_pd) 
     
    839881    source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars)) 
    840882    source.append("#define PROJECTION %d"%PROJECTION) 
    841  
    842883    # 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) 
     884    ocl = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name) 
     885    dll = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name) 
     886 
    846887    result = { 
    847888        'dll': '\n'.join(source+dll[0]+dll[1]+dll[2]), 
    848889        'opencl': '\n'.join(source+ocl[0]+ocl[1]+ocl[2]), 
    849890    } 
    850  
    851891    return result 
    852892 
    853893 
    854 def _kernels(kernel, call_iq, call_iqxy, clear_iqxy, name): 
     894def _kernels(kernel, call_iq, clear_iq, call_iqxy, clear_iqxy, name): 
    855895    # type: ([str,str], str, str, str) -> List[str] 
    856896    code = kernel[0] 
     
    862902        '#line 1 "%s Iq"' % path, 
    863903        code, 
    864         "#undef CALL_IQ", 
     904        clear_iq, 
    865905        "#undef KERNEL_NAME", 
    866906        ] 
     
    9651005    docs = model_info.docs if model_info.docs is not None else "" 
    9661006    docs = convert_section_titles_to_boldface(docs) 
    967     pars = make_partable(model_info.parameters.COMMON 
    968                          + model_info.parameters.kernel_parameters) 
     1007    if model_info.structure_factor: 
     1008        pars = model_info.parameters.kernel_parameters 
     1009    else: 
     1010        pars = model_info.parameters.COMMON + model_info.parameters.kernel_parameters 
     1011    partable = make_partable(pars) 
    9691012    subst = dict(id=model_info.id.replace('_', '-'), 
    9701013                 name=model_info.name, 
    9711014                 title=model_info.title, 
    972                  parameters=pars, 
     1015                 parameters=partable, 
    9731016                 returns=Sq_units if model_info.structure_factor else Iq_units, 
    9741017                 docs=docs) 
  • sasmodels/jitter.py

    rb3703f5 r1198f90  
    774774        # set small jitter as 0 if multiple pd dims 
    775775        dims = sum(v > 0 for v in jitter) 
    776         limit = [0, 0, 0.5, 5][dims] 
     776        limit = [0, 0.5, 5][dims] 
    777777        jitter = [0 if v < limit else v for v in jitter] 
    778778        axes.cla() 
  • sasmodels/kernel.py

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

    rdc6f601 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. 
     
    8387  in_spin = clip(in_spin, 0.0, 1.0); 
    8488  out_spin = clip(out_spin, 0.0, 1.0); 
    85   // Note: sasview 3.1 scaled all slds by sqrt(weight) and assumed that 
     89  // Previous version of this function took the square root of the weights, 
     90  // under the assumption that 
     91  // 
    8692  //     w*I(q, rho1, rho2, ...) = I(q, sqrt(w)*rho1, sqrt(w)*rho2, ...) 
    87   // which is likely to be the case for simple models. 
    88   weight[0] = sqrt((1.0-in_spin) * (1.0-out_spin)); // dd 
    89   weight[1] = sqrt((1.0-in_spin) * out_spin);       // du.real 
    90   weight[2] = sqrt(in_spin * (1.0-out_spin));       // ud.real 
    91   weight[3] = sqrt(in_spin * out_spin);             // uu 
     93  // 
     94  // However, since the weights are applied to the final intensity and 
     95  // are not interned inside the I(q) function, we want the full 
     96  // weight and not the square root.  Any function using 
     97  // set_spin_weights as part of calculating an amplitude will need to 
     98  // manually take that square root, but there is currently no such 
     99  // function. 
     100  weight[0] = (1.0-in_spin) * (1.0-out_spin); // dd 
     101  weight[1] = (1.0-in_spin) * out_spin;       // du 
     102  weight[2] = in_spin * (1.0-out_spin);       // ud 
     103  weight[3] = in_spin * out_spin;             // uu 
    92104  weight[4] = weight[1]; // du.imag 
    93105  weight[5] = weight[2]; // ud.imag 
     
    180192    QACRotation *rotation, 
    181193    double qx, double qy, 
    182     double *qa_out, double *qc_out) 
     194    double *qab_out, double *qc_out) 
    183195{ 
     196    // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2 
    184197    const double dqc = rotation->R31*qx + rotation->R32*qy; 
    185     // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2 
    186     const double dqa = sqrt(-dqc*dqc + qx*qx + qy*qy); 
    187  
    188     *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; 
    189201    *qc_out = dqc; 
    190202} 
     
    262274#endif // _QABC_SECTION 
    263275 
    264  
    265276// ==================== KERNEL CODE ======================== 
    266  
    267277kernel 
    268278void KERNEL_NAME( 
    269     int32_t nq,                 // number of q values 
    270     const int32_t pd_start,     // where we are in the dispersity loop 
    271     const int32_t pd_stop,      // where we are stopping in the dispersity loop 
    272     global const ProblemDetails *details, 
    273     global const double *values, 
    274     global const double *q, // nq q values, with padding to boundary 
    275     global double *result,  // nq+1 return values, again with padding 
    276     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 
    277288    ) 
    278289{ 
    279 #ifdef USE_OPENCL 
     290#if defined(USE_GPU) 
    280291  // who we are and what element we are working with 
     292  #if defined(USE_OPENCL) 
    281293  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 
    282297  if (q_index >= nq) return; 
    283298#else 
     
    330345  // 
    331346  // The code differs slightly between opencl and dll since opencl is only 
    332   // 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 
    333348  // version must loop over all q. 
    334   #ifdef USE_OPENCL 
    335     double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
    336     double this_result = (pd_start == 0 ? 0.0 : result[q_index]); 
    337   #else // !USE_OPENCL 
    338     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 
    339368    if (pd_start == 0) { 
    340369      #ifdef USE_OPENMP 
    341370      #pragma omp parallel for 
    342371      #endif 
    343       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 
    344378    } 
    345379    //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]); 
    346 #endif // !USE_OPENCL 
     380#endif // !USE_GPU 
    347381 
    348382 
     
    367401  const int n4 = pd_length[4]; 
    368402  const int p4 = pd_par[4]; 
    369   global const double *v4 = pd_value + pd_offset[4]; 
    370   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]; 
    371405  int i4 = (pd_start/pd_stride[4])%n4;  // position in level 4 at pd_start 
    372406 
     
    403437          FETCH_Q         // set qx,qy from the q input vector 
    404438          APPLY_ROTATION  // convert qx,qy to qa,qb,qc 
    405           CALL_KERNEL     // scattering = Iqxy(qa, qb, qc, p1, p2, ...) 
     439          CALL_KERNEL     // F2 = Iqxy(qa, qb, qc, p1, p2, ...) 
    406440 
    407441      ++step;  // increment counter representing position in dispersity mesh 
     
    434468// inner loop and defines the macros that use them. 
    435469 
    436 #if defined(CALL_IQ) 
    437   // 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> 
    438497  double qk; 
    439498  #define FETCH_Q() do { qk = q[q_index]; } while (0) 
    440499  #define BUILD_ROTATION() do {} while(0) 
    441500  #define APPLY_ROTATION() do {} while(0) 
    442   #define CALL_KERNEL() CALL_IQ(qk, local_values.table) 
     501  #define CALL_KERNEL() CALL_IQ(qk,local_values.table) 
    443502 
    444503#elif defined(CALL_IQ_A) 
     
    474533  #define APPLY_ROTATION() qabc_apply(&rotation, qx, qy, &qa, &qb, &qc) 
    475534  #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table) 
     535 
    476536#elif defined(CALL_IQ_XY) 
    477537  // direct call to qx,qy calculator 
     
    487547// logic in the IQ_AC and IQ_ABC branches.  This will become more important 
    488548// if we implement more projections, or more complicated projections. 
    489 #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 
    490551  #define APPLY_PROJECTION() const double weight=weight0 
    491552#elif defined(CALL_IQ_XY) // pass orientation to the model 
     
    554615  const int n##_LOOP = details->pd_length[_LOOP]; \ 
    555616  const int p##_LOOP = details->pd_par[_LOOP]; \ 
    556   global const double *v##_LOOP = pd_value + details->pd_offset[_LOOP]; \ 
    557   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]; \ 
    558619  int i##_LOOP = (pd_start/details->pd_stride[_LOOP])%n##_LOOP; 
    559620 
     
    579640// Pointers to the start of the dispersity and weight vectors, if needed. 
    580641#if MAX_PD>0 
    581   global const double *pd_value = values + NUM_VALUES; 
    582   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; 
    583644#endif 
    584645 
     
    637698    // Note: weight==0 must always be excluded 
    638699    if (weight > cutoff) { 
    639       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      } 
    640708      BUILD_ROTATION(); 
    641709 
    642 #ifndef USE_OPENCL 
     710#if !defined(USE_GPU) 
    643711      // DLL needs to explicitly loop over the q values. 
    644712      #ifdef USE_OPENMP 
     
    646714      #endif 
    647715      for (q_index=0; q_index<nq; q_index++) 
    648 #endif // !USE_OPENCL 
     716#endif // !USE_GPU 
    649717      { 
    650718 
     
    655723        #if defined(MAGNETIC) && NUM_MAGNETIC > 0 
    656724          // Compute the scattering from the magnetic cross sections. 
    657           double scattering = 0.0; 
     725          double F2 = 0.0; 
    658726          const double qsq = qx*qx + qy*qy; 
    659727          if (qsq > 1.e-16) { 
     
    680748//  q_index, qx, qy, xs, sk, local_values.vector[sld_index], px, py, mx, my, mz); 
    681749                } 
    682                 scattering += xs_weight * CALL_KERNEL(); 
     750                F2 += xs_weight * CALL_KERNEL(); 
    683751              } 
    684752            } 
    685753          } 
    686754        #else  // !MAGNETIC 
    687           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 
    688760        #endif // !MAGNETIC 
    689 //printf("q_index:%d %g %g %g %g\n", q_index, scattering, weight0); 
    690  
    691         #ifdef USE_OPENCL 
    692           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 
    693770        #else // !USE_OPENCL 
    694           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 
    695777        #endif // !USE_OPENCL 
    696778      } 
    697779    } 
    698780  } 
    699  
    700781// close nested loops 
    701782++step; 
     
    716797#endif 
    717798 
    718 // Remember the current result and the updated norm. 
    719 #ifdef USE_OPENCL 
    720   result[q_index] = this_result; 
    721   if (q_index == 0) result[nq] = pd_norm; 
    722 //if (q_index == 0) printf("res: %g/%g\n", result[0], pd_norm); 
    723 #else // !USE_OPENCL 
    724   result[nq] = pd_norm; 
    725 //printf("res: %g/%g\n", result[0], pd_norm); 
    726 #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  } 
    727822 
    728823// ** clear the macros in preparation for the next kernel ** 
  • sasmodels/kernelcl.py

    rd86f0fc rf872fd1  
    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 
     
    5961 
    6062 
    61 # Attempt to setup opencl. This may fail if the opencl package is not 
     63# Attempt to setup opencl. This may fail if the pyopencl package is not 
    6264# installed or if it is installed but there are no devices available. 
    6365try: 
     
    7476 
    7577from . import generate 
     78from .generate import F32, F64 
    7679from .kernel import KernelModel, Kernel 
    7780 
     
    131134 
    132135def use_opencl(): 
    133     return HAVE_OPENCL and os.environ.get("SAS_OPENCL", "").lower() != "none" 
     136    sas_opencl = os.environ.get("SAS_OPENCL", "OpenCL").lower() 
     137    return HAVE_OPENCL and sas_opencl != "none" and not sas_opencl.startswith("cuda") 
    134138 
    135139ENV = None 
     
    162166    Return true if device supports the requested precision. 
    163167    """ 
    164     if dtype == generate.F32: 
     168    if dtype == F32: 
    165169        return True 
    166     elif dtype == generate.F64: 
     170    elif dtype == F64: 
    167171        return "cl_khr_fp64" in device.extensions 
    168     elif dtype == generate.F16: 
    169         return "cl_khr_fp16" in device.extensions 
    170172    else: 
     173        # Not supporting F16 type since it isn't accurate enough 
    171174        return False 
    172175 
     
    179182        cl.kernel_work_group_info.PREFERRED_WORK_GROUP_SIZE_MULTIPLE, 
    180183        queue.device) 
    181  
    182 def _stretch_input(vector, dtype, extra=1e-3, boundary=32): 
    183     # type: (np.ndarray, np.dtype, float, int) -> np.ndarray 
    184     """ 
    185     Stretch an input vector to the correct boundary. 
    186  
    187     Performance on the kernels can drop by a factor of two or more if the 
    188     number of values to compute does not fall on a nice power of two 
    189     boundary.   The trailing additional vector elements are given a 
    190     value of *extra*, and so f(*extra*) will be computed for each of 
    191     them.  The returned array will thus be a subset of the computed array. 
    192  
    193     *boundary* should be a power of 2 which is at least 32 for good 
    194     performance on current platforms (as of Jan 2015).  It should 
    195     probably be the max of get_warp(kernel,queue) and 
    196     device.min_data_type_align_size//4. 
    197     """ 
    198     remainder = vector.size % boundary 
    199     if remainder != 0: 
    200         size = vector.size + (boundary - remainder) 
    201         vector = np.hstack((vector, [extra] * (size - vector.size))) 
    202     return np.ascontiguousarray(vector, dtype=dtype) 
    203  
    204184 
    205185def compile_model(context, source, dtype, fast=False): 
     
    239219    """ 
    240220    GPU context, with possibly many devices, and one queue per device. 
     221 
     222    Because the environment can be reset during a live program (e.g., if the 
     223    user changes the active GPU device in the GUI), everything associated 
     224    with the device context must be cached in the environment and recreated 
     225    if the environment changes.  The *cache* attribute is a simple dictionary 
     226    which holds keys and references to objects, such as compiled kernels and 
     227    allocated buffers.  The running program should check in the cache for 
     228    long lived objects and create them if they are not there.  The program 
     229    should not hold onto cached objects, but instead only keep them active 
     230    for the duration of a function call.  When the environment is destroyed 
     231    then the *release* method for each active cache item is called before 
     232    the environment is freed.  This means that each cl buffer should be 
     233    in its own cache entry. 
    241234    """ 
    242235    def __init__(self): 
    243236        # type: () -> None 
    244237        # find gpu context 
    245         #self.context = cl.create_some_context() 
    246  
    247         self.context = None 
    248         if 'SAS_OPENCL' in os.environ: 
    249             #Setting PYOPENCL_CTX as a SAS_OPENCL to create cl context 
    250             os.environ["PYOPENCL_CTX"] = os.environ["SAS_OPENCL"] 
    251         if 'PYOPENCL_CTX' in os.environ: 
    252             self._create_some_context() 
    253  
    254         if not self.context: 
    255             self.context = _get_default_context() 
     238        context_list = _create_some_context() 
     239 
     240        # Find a context for F32 and for F64 (maybe the same one). 
     241        # F16 isn't good enough. 
     242        self.context = {} 
     243        for dtype in (F32, F64): 
     244            for context in context_list: 
     245                if has_type(context.devices[0], dtype): 
     246                    self.context[dtype] = context 
     247                    break 
     248            else: 
     249                self.context[dtype] = None 
     250 
     251        # Build a queue for each context 
     252        self.queue = {} 
     253        context = self.context[F32] 
     254        self.queue[F32] = cl.CommandQueue(context, context.devices[0]) 
     255        if self.context[F64] == self.context[F32]: 
     256            self.queue[F64] = self.queue[F32] 
     257        else: 
     258            context = self.context[F64] 
     259            self.queue[F64] = cl.CommandQueue(context, context.devices[0]) 
    256260 
    257261        # Byte boundary for data alignment 
    258         #self.data_boundary = max(d.min_data_type_align_size 
    259         #                         for d in self.context.devices) 
    260         self.queues = [cl.CommandQueue(context, context.devices[0]) 
    261                        for context in self.context] 
     262        #self.data_boundary = max(context.devices[0].min_data_type_align_size 
     263        #                         for context in self.context.values()) 
     264 
     265        # Cache for compiled programs, and for items in context 
    262266        self.compiled = {} 
     267        self.cache = {} 
    263268 
    264269    def has_type(self, dtype): 
     
    267272        Return True if all devices support a given type. 
    268273        """ 
    269         return any(has_type(d, dtype) 
    270                    for context in self.context 
    271                    for d in context.devices) 
    272  
    273     def get_queue(self, dtype): 
    274         # type: (np.dtype) -> cl.CommandQueue 
    275         """ 
    276         Return a command queue for the kernels of type dtype. 
    277         """ 
    278         for context, queue in zip(self.context, self.queues): 
    279             if all(has_type(d, dtype) for d in context.devices): 
    280                 return queue 
    281  
    282     def get_context(self, dtype): 
    283         # type: (np.dtype) -> cl.Context 
    284         """ 
    285         Return a OpenCL context for the kernels of type dtype. 
    286         """ 
    287         for context in self.context: 
    288             if all(has_type(d, dtype) for d in context.devices): 
    289                 return context 
    290  
    291     def _create_some_context(self): 
    292         # type: () -> cl.Context 
    293         """ 
    294         Protected call to cl.create_some_context without interactivity.  Use 
    295         this if SAS_OPENCL is set in the environment.  Sets the *context* 
    296         attribute. 
    297         """ 
    298         try: 
    299             self.context = [cl.create_some_context(interactive=False)] 
    300         except Exception as exc: 
    301             warnings.warn(str(exc)) 
    302             warnings.warn("pyopencl.create_some_context() failed") 
    303             warnings.warn("the environment variable 'SAS_OPENCL' might not be set correctly") 
     274        return self.context.get(dtype, None) is not None 
    304275 
    305276    def compile_program(self, name, source, dtype, fast, timestamp): 
     
    318289            del self.compiled[key] 
    319290        if key not in self.compiled: 
    320             context = self.get_context(dtype) 
     291            context = self.context[dtype] 
    321292            logging.info("building %s for OpenCL %s", key, 
    322293                         context.devices[0].name.strip()) 
    323             program = compile_model(self.get_context(dtype), 
     294            program = compile_model(self.context[dtype], 
    324295                                    str(source), dtype, fast) 
    325296            self.compiled[key] = (program, timestamp) 
    326297        return program 
     298 
     299    def free_buffer(self, key): 
     300        if key in self.cache: 
     301            self.cache[key].release() 
     302            del self.cache[key] 
     303 
     304    def __del__(self): 
     305        for v in self.cache.values(): 
     306            release = getattr(v, 'release', lambda: None) 
     307            release() 
     308        self.cache = {} 
     309 
     310_CURRENT_ID = 0 
     311def unique_id(): 
     312    global _CURRENT_ID 
     313    _CURRENT_ID += 1 
     314    return _CURRENT_ID 
     315 
     316def _create_some_context(): 
     317    # type: () -> cl.Context 
     318    """ 
     319    Protected call to cl.create_some_context without interactivity. 
     320 
     321    Uses SAS_OPENCL or PYOPENCL_CTX if they are set in the environment, 
     322    otherwise scans for the most appropriate device using 
     323    :func:`_get_default_context`.  Ignore *SAS_OPENCL=OpenCL*, which 
     324    indicates that an OpenCL device should be used without specifying 
     325    which one (and not a CUDA device, or no GPU). 
     326    """ 
     327    # Assume we do not get here if SAS_OPENCL is None or CUDA 
     328    sas_opencl = os.environ.get('SAS_OPENCL', 'opencl') 
     329    if sas_opencl.lower() != 'opencl': 
     330        # Setting PYOPENCL_CTX as a SAS_OPENCL to create cl context 
     331        os.environ["PYOPENCL_CTX"] = sas_opencl 
     332 
     333    if 'PYOPENCL_CTX' in os.environ: 
     334        try: 
     335            return [cl.create_some_context(interactive=False)] 
     336        except Exception as exc: 
     337            warnings.warn(str(exc)) 
     338            warnings.warn("pyopencl.create_some_context() failed") 
     339            warnings.warn("the environment variable 'SAS_OPENCL' or 'PYOPENCL_CTX' might not be set correctly") 
     340 
     341    return _get_default_context() 
    327342 
    328343def _get_default_context(): 
     
    404419        self.dtype = dtype 
    405420        self.fast = fast 
    406         self.program = None # delay program creation 
    407         self._kernels = None 
     421        self.timestamp = generate.ocl_timestamp(self.info) 
     422        self._cache_key = unique_id() 
    408423 
    409424    def __getstate__(self): 
     
    414429        # type: (Tuple[ModelInfo, str, np.dtype, bool]) -> None 
    415430        self.info, self.source, self.dtype, self.fast = state 
    416         self.program = None 
    417431 
    418432    def make_kernel(self, q_vectors): 
    419433        # type: (List[np.ndarray]) -> "GpuKernel" 
    420         if self.program is None: 
    421             compile_program = environment().compile_program 
    422             timestamp = generate.ocl_timestamp(self.info) 
    423             self.program = compile_program( 
     434        return GpuKernel(self, q_vectors) 
     435 
     436    @property 
     437    def Iq(self): 
     438        return self._fetch_kernel('Iq') 
     439 
     440    def fetch_kernel(self, name): 
     441        # type: (str) -> cl.Kernel 
     442        """ 
     443        Fetch the kernel from the environment by name, compiling it if it 
     444        does not already exist. 
     445        """ 
     446        gpu = environment() 
     447        key = self._cache_key 
     448        if key not in gpu.cache: 
     449            program = gpu.compile_program( 
    424450                self.info.name, 
    425451                self.source['opencl'], 
    426452                self.dtype, 
    427453                self.fast, 
    428                 timestamp) 
     454                self.timestamp) 
    429455            variants = ['Iq', 'Iqxy', 'Imagnetic'] 
    430456            names = [generate.kernel_name(self.info, k) for k in variants] 
    431             kernels = [getattr(self.program, k) for k in names] 
    432             self._kernels = dict((k, v) for k, v in zip(variants, kernels)) 
    433         is_2d = len(q_vectors) == 2 
    434         if is_2d: 
    435             kernel = [self._kernels['Iqxy'], self._kernels['Imagnetic']] 
     457            kernels = [getattr(program, k) for k in names] 
     458            data = dict((k, v) for k, v in zip(variants, kernels)) 
     459            # keep a handle to program so GC doesn't collect 
     460            data['program'] = program 
     461            gpu.cache[key] = data 
    436462        else: 
    437             kernel = [self._kernels['Iq']]*2 
    438         return GpuKernel(kernel, self.dtype, self.info, q_vectors) 
    439  
    440     def release(self): 
    441         # type: () -> None 
    442         """ 
    443         Free the resources associated with the model. 
    444         """ 
    445         if self.program is not None: 
    446             self.program = None 
    447  
    448     def __del__(self): 
    449         # type: () -> None 
    450         self.release() 
     463            data = gpu.cache[key] 
     464        return data[name] 
    451465 
    452466# TODO: check that we don't need a destructor for buffers which go out of scope 
     
    473487        # type: (List[np.ndarray], np.dtype) -> None 
    474488        # TODO: do we ever need double precision q? 
    475         env = environment() 
    476489        self.nq = q_vectors[0].size 
    477490        self.dtype = np.dtype(dtype) 
     
    482495        # architectures tested so far. 
    483496        if self.is_2d: 
    484             # Note: 16 rather than 15 because result is 1 longer than input. 
    485             width = ((self.nq+16)//16)*16 
     497            width = ((self.nq+15)//16)*16 
    486498            self.q = np.empty((width, 2), dtype=dtype) 
    487499            self.q[:self.nq, 0] = q_vectors[0] 
    488500            self.q[:self.nq, 1] = q_vectors[1] 
    489501        else: 
    490             # Note: 32 rather than 31 because result is 1 longer than input. 
    491             width = ((self.nq+32)//32)*32 
     502            width = ((self.nq+31)//32)*32 
    492503            self.q = np.empty(width, dtype=dtype) 
    493504            self.q[:self.nq] = q_vectors[0] 
    494505        self.global_size = [self.q.shape[0]] 
    495         context = env.get_context(self.dtype) 
    496         #print("creating inputs of size", self.global_size) 
    497         self.q_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
    498                              hostbuf=self.q) 
     506        self._cache_key = unique_id() 
     507 
     508    @property 
     509    def q_b(self): 
     510        """Lazy creation of q buffer so it can survive context reset""" 
     511        env = environment() 
     512        key = self._cache_key 
     513        if key not in env.cache: 
     514            context = env.context[self.dtype] 
     515            #print("creating inputs of size", self.global_size) 
     516            buffer = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
     517                               hostbuf=self.q) 
     518            env.cache[key] = buffer 
     519        return env.cache[key] 
    499520 
    500521    def release(self): 
    501522        # type: () -> None 
    502523        """ 
    503         Free the memory. 
    504         """ 
    505         if self.q_b is not None: 
    506             self.q_b.release() 
    507             self.q_b = None 
     524        Free the buffer associated with the q value 
     525        """ 
     526        environment().free_buffer(id(self)) 
    508527 
    509528    def __del__(self): 
     
    515534    Callable SAS kernel. 
    516535 
    517     *kernel* is the GpuKernel object to call 
    518  
    519     *model_info* is the module information 
    520  
    521     *q_vectors* is the q vectors at which the kernel should be evaluated 
     536    *model* is the GpuModel object to call 
     537 
     538    The following attributes are defined: 
     539 
     540    *info* is the module information 
    522541 
    523542    *dtype* is the kernel precision 
     543 
     544    *dim* is '1d' or '2d' 
     545 
     546    *result* is a vector to contain the results of the call 
    524547 
    525548    The resulting call method takes the *pars*, a list of values for 
     
    531554    Call :meth:`release` when done with the kernel instance. 
    532555    """ 
    533     def __init__(self, kernel, dtype, model_info, q_vectors): 
     556    def __init__(self, model, q_vectors): 
    534557        # type: (cl.Kernel, np.dtype, ModelInfo, List[np.ndarray]) -> None 
    535         q_input = GpuInput(q_vectors, dtype) 
    536         self.kernel = kernel 
    537         self.info = model_info 
    538         self.dtype = dtype 
    539         self.dim = '2d' if q_input.is_2d else '1d' 
    540         # plus three for the normalization values 
    541         self.result = np.empty(q_input.nq+1, dtype) 
    542  
    543         # Inputs and outputs for each kernel call 
    544         # Note: res may be shorter than res_b if global_size != nq 
     558        dtype = model.dtype 
     559        self.q_input = GpuInput(q_vectors, dtype) 
     560        self._model = model 
     561        # F16 isn't sufficient, so don't support it 
     562        self._as_dtype = np.float64 if dtype == generate.F64 else np.float32 
     563        self._cache_key = unique_id() 
     564 
     565        # attributes accessed from the outside 
     566        self.dim = '2d' if self.q_input.is_2d else '1d' 
     567        self.info = model.info 
     568        self.dtype = model.dtype 
     569 
     570        # holding place for the returned value 
     571        nout = 2 if self.info.have_Fq and self.dim == '1d' else 1 
     572        extra_q = 4  # total weight, form volume, shell volume and R_eff 
     573        self.result = np.empty(self.q_input.nq*nout+extra_q, dtype) 
     574 
     575    @property 
     576    def _result_b(self): 
     577        """Lazy creation of result buffer so it can survive context reset""" 
    545578        env = environment() 
    546         self.queue = env.get_queue(dtype) 
    547  
    548         self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 
    549                                   q_input.global_size[0] * dtype.itemsize) 
    550         self.q_input = q_input # allocated by GpuInput above 
    551  
    552         self._need_release = [self.result_b, self.q_input] 
    553         self.real = (np.float32 if dtype == generate.F32 
    554                      else np.float64 if dtype == generate.F64 
    555                      else np.float16 if dtype == generate.F16 
    556                      else np.float32)  # will never get here, so use np.float32 
    557  
    558     def __call__(self, call_details, values, cutoff, magnetic): 
     579        key = self._cache_key 
     580        if key not in env.cache: 
     581            context = env.context[self.dtype] 
     582            width = ((self.result.size+31)//32)*32 * self.dtype.itemsize 
     583            buffer = cl.Buffer(context, mf.READ_WRITE, width) 
     584            env.cache[key] = buffer 
     585        return env.cache[key] 
     586 
     587    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
    559588        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    560         context = self.queue.context 
    561         # Arrange data transfer to card 
     589        env = environment() 
     590        queue = env.queue[self._model.dtype] 
     591        context = queue.context 
     592 
     593        # Arrange data transfer to/from card 
     594        q_b = self.q_input.q_b 
     595        result_b = self._result_b 
    562596        details_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
    563597                              hostbuf=call_details.buffer) 
     
    565599                             hostbuf=values) 
    566600 
    567         kernel = self.kernel[1 if magnetic else 0] 
    568         args = [ 
     601        name = 'Iq' if self.dim == '1d' else 'Imagnetic' if magnetic else 'Iqxy' 
     602        kernel = self._model.fetch_kernel(name) 
     603        kernel_args = [ 
    569604            np.uint32(self.q_input.nq), None, None, 
    570             details_b, values_b, self.q_input.q_b, self.result_b, 
    571             self.real(cutoff), 
     605            details_b, values_b, q_b, result_b, 
     606            self._as_dtype(cutoff), 
     607            np.uint32(effective_radius_type), 
    572608        ] 
    573609        #print("Calling OpenCL") 
    574610        #call_details.show(values) 
    575         # Call kernel and retrieve results 
     611        #Call kernel and retrieve results 
    576612        wait_for = None 
    577613        last_nap = time.clock() 
     
    580616            stop = min(start + step, call_details.num_eval) 
    581617            #print("queuing",start,stop) 
    582             args[1:3] = [np.int32(start), np.int32(stop)] 
    583             wait_for = [kernel(self.queue, self.q_input.global_size, None, 
    584                                *args, wait_for=wait_for)] 
     618            kernel_args[1:3] = [np.int32(start), np.int32(stop)] 
     619            wait_for = [kernel(queue, self.q_input.global_size, None, 
     620                               *kernel_args, wait_for=wait_for)] 
    585621            if stop < call_details.num_eval: 
    586622                # Allow other processes to run 
     
    588624                current_time = time.clock() 
    589625                if current_time - last_nap > 0.5: 
    590                     time.sleep(0.05) 
     626                    time.sleep(0.001) 
    591627                    last_nap = current_time 
    592         cl.enqueue_copy(self.queue, self.result, self.result_b) 
     628        cl.enqueue_copy(queue, self.result, result_b, wait_for=wait_for) 
    593629        #print("result", self.result) 
    594630 
     
    598634                v.release() 
    599635 
    600         pd_norm = self.result[self.q_input.nq] 
    601         scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
    602         background = values[1] 
    603         #print("scale",scale,values[0],self.result[self.q_input.nq],background) 
    604         return scale*self.result[:self.q_input.nq] + background 
    605         # return self.result[:self.q_input.nq] 
    606  
    607636    def release(self): 
    608637        # type: () -> None 
     
    610639        Release resources associated with the kernel. 
    611640        """ 
    612         for v in self._need_release: 
    613             v.release() 
    614         self._need_release = [] 
     641        environment().free_buffer(id(self)) 
     642        self.q_input.release() 
    615643 
    616644    def __del__(self): 
  • sasmodels/kerneldll.py

    r1662ebe re44432d  
    9999    pass 
    100100# pylint: enable=unused-import 
     101 
     102# Compiler output is a byte stream that needs to be decode in python 3 
     103decode = (lambda s: s) if sys.version_info[0] < 3 else (lambda s: s.decode('utf8')) 
     104 
     105if "SAS_DLL_PATH" in os.environ: 
     106    SAS_DLL_PATH = os.environ["SAS_DLL_PATH"] 
     107else: 
     108    # Assume the default location of module DLLs is in .sasmodels/compiled_models. 
     109    SAS_DLL_PATH = os.path.join(os.path.expanduser("~"), ".sasmodels", "compiled_models") 
    101110 
    102111if "SAS_COMPILER" in os.environ: 
     
    161170        return CC + [source, "-o", output, "-lm"] 
    162171 
    163 # Assume the default location of module DLLs is in .sasmodels/compiled_models. 
    164 DLL_PATH = os.path.join(os.path.expanduser("~"), ".sasmodels", "compiled_models") 
    165  
    166172ALLOW_SINGLE_PRECISION_DLLS = True 
    167173 
     
    181187        subprocess.check_output(command, shell=shell, stderr=subprocess.STDOUT) 
    182188    except subprocess.CalledProcessError as exc: 
    183         raise RuntimeError("compile failed.\n%s\n%s" 
    184                            % (command_str, exc.output.decode())) 
     189        output = decode(exc.output) 
     190        raise RuntimeError("compile failed.\n%s\n%s"%(command_str, output)) 
    185191    if not os.path.exists(output): 
    186192        raise RuntimeError("compile failed.  File is in %r"%source) 
     
    201207        return path 
    202208 
    203     return joinpath(DLL_PATH, basename) 
     209    return joinpath(SAS_DLL_PATH, basename) 
    204210 
    205211 
     
    210216    exist yet if it hasn't been compiled. 
    211217    """ 
    212     return os.path.join(DLL_PATH, dll_name(model_info, dtype)) 
     218    return os.path.join(SAS_DLL_PATH, dll_name(model_info, dtype)) 
    213219 
    214220 
     
    229235    models are not allowed as DLLs. 
    230236 
    231     Set *sasmodels.kerneldll.DLL_PATH* to the compiled dll output path. 
     237    Set *sasmodels.kerneldll.SAS_DLL_PATH* to the compiled dll output path. 
     238    Alternatively, set the environment variable *SAS_DLL_PATH*. 
    232239    The default is in ~/.sasmodels/compiled_models. 
    233240    """ 
     
    248255    if need_recompile: 
    249256        # Make sure the DLL path exists 
    250         if not os.path.exists(DLL_PATH): 
    251             os.makedirs(DLL_PATH) 
     257        if not os.path.exists(SAS_DLL_PATH): 
     258            os.makedirs(SAS_DLL_PATH) 
    252259        basename = splitext(os.path.basename(dll))[0] + "_" 
    253260        system_fd, filename = tempfile.mkstemp(suffix=".c", prefix=basename) 
     
    312319 
    313320        # int, int, int, int*, double*, double*, double*, double*, double 
    314         argtypes = [ct.c_int32]*3 + [ct.c_void_p]*4 + [float_type] 
     321        argtypes = [ct.c_int32]*3 + [ct.c_void_p]*4 + [float_type, ct.c_int32] 
    315322        names = [generate.kernel_name(self.info, variant) 
    316323                 for variant in ("Iq", "Iqxy", "Imagnetic")] 
     
    372379    def __init__(self, kernel, model_info, q_input): 
    373380        # type: (Callable[[], np.ndarray], ModelInfo, PyInput) -> None 
     381        #,model_info,q_input) 
    374382        self.kernel = kernel 
    375383        self.info = model_info 
     
    377385        self.dtype = q_input.dtype 
    378386        self.dim = '2d' if q_input.is_2d else '1d' 
    379         self.result = np.empty(q_input.nq+1, q_input.dtype) 
     387        # leave room for f1/f2 results in case we need to compute beta for 1d models 
     388        nout = 2 if self.info.have_Fq else 1 
     389        # +4 for total weight, shell volume, effective radius, form volume 
     390        self.result = np.empty(q_input.nq*nout + 4, self.dtype) 
    380391        self.real = (np.float32 if self.q_input.dtype == generate.F32 
    381392                     else np.float64 if self.q_input.dtype == generate.F64 
    382393                     else np.float128) 
    383394 
    384     def __call__(self, call_details, values, cutoff, magnetic): 
    385         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    386  
     395    def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
     396        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
    387397        kernel = self.kernel[1 if magnetic else 0] 
    388398        args = [ 
     
    391401            None, # pd_stop pd_stride[MAX_PD] 
    392402            call_details.buffer.ctypes.data, # problem 
    393             values.ctypes.data,  #pars 
    394             self.q_input.q.ctypes.data, #q 
     403            values.ctypes.data,  # pars 
     404            self.q_input.q.ctypes.data, # q 
    395405            self.result.ctypes.data,   # results 
    396406            self.real(cutoff), # cutoff 
     407            effective_radius_type, # cutoff 
    397408        ] 
    398409        #print("Calling DLL") 
     
    404415            kernel(*args) # type: ignore 
    405416 
    406         #print("returned",self.q_input.q, self.result) 
    407         pd_norm = self.result[self.q_input.nq] 
    408         scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
    409         background = values[1] 
    410         #print("scale",scale,background) 
    411         return scale*self.result[:self.q_input.nq] + background 
    412  
    413417    def release(self): 
    414418        # type: () -> None 
  • sasmodels/kernelpy.py

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

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

    • Property mode changed from 100644 to 100755
    r3221de0 raa8c6e0  
    55Usage:: 
    66 
    7     python -m sasmodels.model_test [opencl|dll|opencl_and_dll] model1 model2 ... 
     7    python -m sasmodels.model_test [opencl|cuda|dll] model1 model2 ... 
    88 
    99    if model1 is 'all', then all except the remaining models will be tested 
    1010 
    1111Each model is tested using the default parameters at q=0.1, (qx, qy)=(0.1, 0.1), 
    12 and the ER and VR are computed.  The return values at these points are not 
    13 considered.  The test is only to verify that the models run to completion, 
    14 and do not produce inf or NaN. 
     12and Fq is called to make sure R_eff, volume and volume ratio are computed. 
     13The return values at these points are not considered.  The test is only to 
     14verify that the models run to completion, and do not produce inf or NaN. 
    1515 
    1616Tests are defined with the *tests* attribute in the model.py file.  *tests* 
    1717is a list of individual tests to run, where each test consists of the 
    1818parameter values for the test, the q-values and the expected results.  For 
    19 the effective radius test, the q-value should be 'ER'.  For the VR test, 
    20 the q-value should be 'VR'.  For 1-D tests, either specify the q value or 
    21 a list of q-values, and the corresponding I(q) value, or list of I(q) values. 
     19the effective radius test and volume ratio tests, use the extended output 
     20form, which checks each output of kernel.Fq. For 1-D tests, either specify 
     21the q value or a list of q-values, and the corresponding I(q) value, or 
     22list of I(q) values. 
    2223 
    2324That is:: 
     
    3233                        [I(qx1, qy1), I(qx2, qy2), ...]], 
    3334 
    34         [ {parameters}, 'ER', ER(pars) ], 
    35         [ {parameters}, 'VR', VR(pars) ], 
     35        [ {parameters}, q, F(q), F^2(q), R_eff, V, V_r ], 
    3636        ... 
    3737    ] 
     
    4747import sys 
    4848import unittest 
     49import traceback 
    4950 
    5051try: 
     
    5960from . import core 
    6061from .core import list_models, load_model_info, build_model 
    61 from .direct_model import call_kernel, call_ER, call_VR 
     62from .direct_model import call_kernel, call_Fq 
    6263from .exception import annotate_exception 
    6364from .modelinfo import expand_pars 
    6465from .kernelcl import use_opencl 
     66from .kernelcuda import use_cuda 
     67from . import product 
    6568 
    6669# pylint: disable=unused-import 
     
    7477# pylint: enable=unused-import 
    7578 
    76  
    7779def make_suite(loaders, models): 
    7880    # type: (List[str], List[str]) -> unittest.TestSuite 
     
    8082    Construct the pyunit test suite. 
    8183 
    82     *loaders* is the list of kernel drivers to use, which is one of 
    83     *["dll", "opencl"]*, *["dll"]* or *["opencl"]*.  For python models, 
    84     the python driver is always used. 
     84    *loaders* is the list of kernel drivers to use (dll, opencl or cuda). 
     85    For python model the python driver is always used. 
    8586 
    8687    *models* is the list of models to test, or *["all"]* to test all models. 
    8788    """ 
    88     ModelTestCase = _hide_model_case_from_nose() 
    8989    suite = unittest.TestSuite() 
    9090 
     
    9595        skip = [] 
    9696    for model_name in models: 
    97         if model_name in skip: 
    98             continue 
    99         model_info = load_model_info(model_name) 
    100  
    101         #print('------') 
    102         #print('found tests in', model_name) 
    103         #print('------') 
    104  
    105         # if ispy then use the dll loader to call pykernel 
    106         # don't try to call cl kernel since it will not be 
    107         # available in some environmentes. 
    108         is_py = callable(model_info.Iq) 
    109  
    110         # Some OpenCL drivers seem to be flaky, and are not producing the 
    111         # expected result.  Since we don't have known test values yet for 
    112         # all of our models, we are instead going to compare the results 
    113         # for the 'smoke test' (that is, evaluation at q=0.1 for the default 
    114         # parameters just to see that the model runs to completion) between 
    115         # the OpenCL and the DLL.  To do this, we define a 'stash' which is 
    116         # shared between OpenCL and DLL tests.  This is just a list.  If the 
    117         # list is empty (which it will be when DLL runs, if the DLL runs 
    118         # first), then the results are appended to the list.  If the list 
    119         # is not empty (which it will be when OpenCL runs second), the results 
    120         # are compared to the results stored in the first element of the list. 
    121         # This is a horrible stateful hack which only makes sense because the 
    122         # test suite is thrown away after being run once. 
    123         stash = [] 
    124  
    125         if is_py:  # kernel implemented in python 
    126             test_name = "%s-python"%model_name 
    127             test_method_name = "test_%s_python" % model_info.id 
     97        if model_name not in skip: 
     98            model_info = load_model_info(model_name) 
     99            _add_model_to_suite(loaders, suite, model_info) 
     100 
     101    return suite 
     102 
     103def _add_model_to_suite(loaders, suite, model_info): 
     104    ModelTestCase = _hide_model_case_from_nose() 
     105 
     106    #print('------') 
     107    #print('found tests in', model_name) 
     108    #print('------') 
     109 
     110    # if ispy then use the dll loader to call pykernel 
     111    # don't try to call cl kernel since it will not be 
     112    # available in some environmentes. 
     113    is_py = callable(model_info.Iq) 
     114 
     115    # Some OpenCL drivers seem to be flaky, and are not producing the 
     116    # expected result.  Since we don't have known test values yet for 
     117    # all of our models, we are instead going to compare the results 
     118    # for the 'smoke test' (that is, evaluation at q=0.1 for the default 
     119    # parameters just to see that the model runs to completion) between 
     120    # the OpenCL and the DLL.  To do this, we define a 'stash' which is 
     121    # shared between OpenCL and DLL tests.  This is just a list.  If the 
     122    # list is empty (which it will be when DLL runs, if the DLL runs 
     123    # first), then the results are appended to the list.  If the list 
     124    # is not empty (which it will be when OpenCL runs second), the results 
     125    # are compared to the results stored in the first element of the list. 
     126    # This is a horrible stateful hack which only makes sense because the 
     127    # test suite is thrown away after being run once. 
     128    stash = [] 
     129 
     130    if is_py:  # kernel implemented in python 
     131        test_name = "%s-python"%model_info.name 
     132        test_method_name = "test_%s_python" % model_info.id 
     133        test = ModelTestCase(test_name, model_info, 
     134                                test_method_name, 
     135                                platform="dll",  # so that 
     136                                dtype="double", 
     137                                stash=stash) 
     138        suite.addTest(test) 
     139    else:   # kernel implemented in C 
     140 
     141        # test using dll if desired 
     142        if 'dll' in loaders or not use_opencl(): 
     143            test_name = "%s-dll"%model_info.name 
     144            test_method_name = "test_%s_dll" % model_info.id 
    128145            test = ModelTestCase(test_name, model_info, 
    129                                  test_method_name, 
    130                                  platform="dll",  # so that 
    131                                  dtype="double", 
    132                                  stash=stash) 
     146                                    test_method_name, 
     147                                    platform="dll", 
     148                                    dtype="double", 
     149                                    stash=stash) 
    133150            suite.addTest(test) 
    134         else:   # kernel implemented in C 
    135  
    136             # test using dll if desired 
    137             if 'dll' in loaders or not use_opencl(): 
    138                 test_name = "%s-dll"%model_name 
    139                 test_method_name = "test_%s_dll" % model_info.id 
    140                 test = ModelTestCase(test_name, model_info, 
    141                                      test_method_name, 
    142                                      platform="dll", 
    143                                      dtype="double", 
    144                                      stash=stash) 
    145                 suite.addTest(test) 
    146  
    147             # test using opencl if desired and available 
    148             if 'opencl' in loaders and use_opencl(): 
    149                 test_name = "%s-opencl"%model_name 
    150                 test_method_name = "test_%s_opencl" % model_info.id 
    151                 # Using dtype=None so that the models that are only 
    152                 # correct for double precision are not tested using 
    153                 # single precision.  The choice is determined by the 
    154                 # presence of *single=False* in the model file. 
    155                 test = ModelTestCase(test_name, model_info, 
    156                                      test_method_name, 
    157                                      platform="ocl", dtype=None, 
    158                                      stash=stash) 
    159                 #print("defining", test_name) 
    160                 suite.addTest(test) 
    161  
    162     return suite 
     151 
     152        # test using opencl if desired and available 
     153        if 'opencl' in loaders and use_opencl(): 
     154            test_name = "%s-opencl"%model_info.name 
     155            test_method_name = "test_%s_opencl" % model_info.id 
     156            # Using dtype=None so that the models that are only 
     157            # correct for double precision are not tested using 
     158            # single precision.  The choice is determined by the 
     159            # presence of *single=False* in the model file. 
     160            test = ModelTestCase(test_name, model_info, 
     161                                    test_method_name, 
     162                                    platform="ocl", dtype=None, 
     163                                    stash=stash) 
     164            #print("defining", test_name) 
     165            suite.addTest(test) 
     166 
     167        # test using cuda if desired and available 
     168        if 'cuda' in loaders and use_cuda(): 
     169            test_name = "%s-cuda"%model_name 
     170            test_method_name = "test_%s_cuda" % model_info.id 
     171            # Using dtype=None so that the models that are only 
     172            # correct for double precision are not tested using 
     173            # single precision.  The choice is determined by the 
     174            # presence of *single=False* in the model file. 
     175            test = ModelTestCase(test_name, model_info, 
     176                                    test_method_name, 
     177                                    platform="cuda", dtype=None, 
     178                                    stash=stash) 
     179            #print("defining", test_name) 
     180            suite.addTest(test) 
     181 
    163182 
    164183def _hide_model_case_from_nose(): 
     
    199218                ({}, [0.001, 0.01, 0.1], [None]*3), 
    200219                ({}, [(0.1, 0.1)]*2, [None]*2), 
    201                 # test that ER/VR will run if they exist 
    202                 ({}, 'ER', None), 
    203                 ({}, 'VR', None), 
     220                # test that Fq will run, and return R_eff, V, V_r 
     221                ({}, 0.1, None, None, None, None, None), 
    204222                ] 
    205223            tests = smoke_tests 
     
    207225            if self.info.tests is not None: 
    208226                tests += self.info.tests 
     227            S_tests = [test for test in tests if '@S' in test[0]] 
     228            P_tests = [test for test in tests if '@S' not in test[0]] 
    209229            try: 
    210230                model = build_model(self.info, dtype=self.dtype, 
    211231                                    platform=self.platform) 
    212                 results = [self.run_one(model, test) for test in tests] 
     232                results = [self.run_one(model, test) for test in P_tests] 
     233                for test in S_tests: 
     234                    # pull the S model name out of the test defn 
     235                    pars = test[0].copy() 
     236                    s_name = pars.pop('@S') 
     237                    ps_test = [pars] + list(test[1:]) 
     238                    # build the P@S model 
     239                    s_info = load_model_info(s_name) 
     240                    ps_info = product.make_product_info(self.info, s_info) 
     241                    ps_model = build_model(ps_info, dtype=self.dtype, 
     242                                           platform=self.platform) 
     243                    # run the tests 
     244                    results.append(self.run_one(ps_model, ps_test)) 
     245 
    213246                if self.stash: 
    214247                    for test, target, actual in zip(tests, self.stash[0], results): 
     
    220253 
    221254                # Check for missing tests.  Only do so for the "dll" tests 
    222                 # to reduce noise from both opencl and dll, and because 
     255                # to reduce noise from both opencl and cuda, and because 
    223256                # python kernels use platform="dll". 
    224257                if self.platform == "dll": 
     
    235268        def _find_missing_tests(self): 
    236269            # type: () -> None 
    237             """make sure there are 1D, 2D, ER and VR tests as appropriate""" 
    238             model_has_VR = callable(self.info.VR) 
    239             model_has_ER = callable(self.info.ER) 
     270            """make sure there are 1D and 2D tests as appropriate""" 
    240271            model_has_1D = True 
    241272            model_has_2D = any(p.type == 'orientation' 
     
    245276            single = [test for test in self.info.tests 
    246277                      if not isinstance(test[2], list) and test[2] is not None] 
    247             tests_has_VR = any(test[1] == 'VR' for test in single) 
    248             tests_has_ER = any(test[1] == 'ER' for test in single) 
    249278            tests_has_1D_single = any(isinstance(test[1], float) for test in single) 
    250279            tests_has_2D_single = any(isinstance(test[1], tuple) for test in single) 
     
    259288 
    260289            missing = [] 
    261             if model_has_VR and not tests_has_VR: 
    262                 missing.append("VR") 
    263             if model_has_ER and not tests_has_ER: 
    264                 missing.append("ER") 
    265290            if model_has_1D and not (tests_has_1D_single or tests_has_1D_multiple): 
    266291                missing.append("1D") 
     
    273298            # type: (KernelModel, TestCondition) -> None 
    274299            """Run a single test case.""" 
    275             user_pars, x, y = test 
     300            user_pars, x, y = test[:3] 
    276301            pars = expand_pars(self.info.parameters, user_pars) 
    277302            invalid = invalid_pars(self.info.parameters, pars) 
     
    286311            self.assertEqual(len(y), len(x)) 
    287312 
    288             if x[0] == 'ER': 
    289                 actual = np.array([call_ER(model.info, pars)]) 
    290             elif x[0] == 'VR': 
    291                 actual = np.array([call_VR(model.info, pars)]) 
    292             elif isinstance(x[0], tuple): 
     313            if isinstance(x[0], tuple): 
    293314                qx, qy = zip(*x) 
    294315                q_vectors = [np.array(qx), np.array(qy)] 
    295                 kernel = model.make_kernel(q_vectors) 
    296                 actual = call_kernel(kernel, pars) 
    297316            else: 
    298317                q_vectors = [np.array(x)] 
    299                 kernel = model.make_kernel(q_vectors) 
     318 
     319            kernel = model.make_kernel(q_vectors) 
     320            if len(test) == 3: 
    300321                actual = call_kernel(kernel, pars) 
    301  
    302             self.assertTrue(len(actual) > 0) 
    303             self.assertEqual(len(y), len(actual)) 
    304  
    305             for xi, yi, actual_yi in zip(x, y, actual): 
     322                self._check_vectors(x, y, actual, 'I') 
     323                return actual 
     324            else: 
     325                y1 = y 
     326                y2 = test[3] if not isinstance(test[3], list) else [test[3]] 
     327                F1, F2, R_eff, volume, volume_ratio = call_Fq(kernel, pars) 
     328                if F1 is not None:  # F1 is none for models with Iq instead of Fq 
     329                    self._check_vectors(x, y1, F1, 'F') 
     330                self._check_vectors(x, y2, F2, 'F^2') 
     331                self._check_scalar(test[4], R_eff, 'R_eff') 
     332                self._check_scalar(test[5], volume, 'volume') 
     333                self._check_scalar(test[6], volume_ratio, 'form:shell ratio') 
     334                return F2 
     335 
     336        def _check_scalar(self, target, actual, name): 
     337            if target is None: 
     338                # smoke test --- make sure it runs and produces a value 
     339                self.assertTrue(not np.isnan(actual), 
     340                                'invalid %s: %s' % (name, actual)) 
     341            elif np.isnan(target): 
     342                # make sure nans match 
     343                self.assertTrue(np.isnan(actual), 
     344                                '%s: expected:%s; actual:%s' 
     345                                % (name, target, actual)) 
     346            else: 
     347                # is_near does not work for infinite values, so also test 
     348                # for exact values. 
     349                self.assertTrue(target == actual or is_near(target, actual, 5), 
     350                                '%s: expected:%s; actual:%s' 
     351                                % (name, target, actual)) 
     352 
     353        def _check_vectors(self, x, target, actual, name='I'): 
     354            self.assertTrue(len(actual) > 0, 
     355                            '%s(...) expected return'%name) 
     356            if target is None: 
     357                return 
     358            self.assertEqual(len(target), len(actual), 
     359                             '%s(...) returned wrong length'%name) 
     360            for xi, yi, actual_yi in zip(x, target, actual): 
    306361                if yi is None: 
    307362                    # smoke test --- make sure it runs and produces a value 
    308363                    self.assertTrue(not np.isnan(actual_yi), 
    309                                     'invalid f(%s): %s' % (xi, actual_yi)) 
     364                                    'invalid %s(%s): %s' % (name, xi, actual_yi)) 
    310365                elif np.isnan(yi): 
     366                    # make sure nans match 
    311367                    self.assertTrue(np.isnan(actual_yi), 
    312                                     'f(%s): expected:%s; actual:%s' 
    313                                     % (xi, yi, actual_yi)) 
     368                                    '%s(%s): expected:%s; actual:%s' 
     369                                    % (name, xi, yi, actual_yi)) 
    314370                else: 
    315371                    # is_near does not work for infinite values, so also test 
    316                     # for exact values.  Note that this will not 
     372                    # for exact values. 
    317373                    self.assertTrue(yi == actual_yi or is_near(yi, actual_yi, 5), 
    318                                     'f(%s); expected:%s; actual:%s' 
    319                                     % (xi, yi, actual_yi)) 
    320             return actual 
     374                                    '%s(%s); expected:%s; actual:%s' 
     375                                    % (name, xi, yi, actual_yi)) 
    321376 
    322377    return ModelTestCase 
     
    330385    invalid = [] 
    331386    for par in sorted(pars.keys()): 
     387        # special handling of R_eff mode, which is not a usual parameter 
     388        if par == 'radius_effective_type': 
     389            continue 
    332390        parts = par.split('_pd') 
    333391        if len(parts) > 1 and parts[1] not in ("", "_n", "nsigma", "type"): 
     
    348406    return abs(target-actual)/shift < 1.5*10**-digits 
    349407 
    350 def run_one(model): 
    351     # type: (str) -> str 
    352     """ 
    353     Run the tests for a single model, printing the results to stdout. 
    354  
    355     *model* can by a python file, which is handy for checking user defined 
    356     plugin models. 
     408# CRUFT: old interface; should be deprecated and removed 
     409def run_one(model_name): 
     410    # msg = "use check_model(model_info) rather than run_one(model_name)" 
     411    # warnings.warn(msg, category=DeprecationWarning, stacklevel=2) 
     412    try: 
     413        model_info = load_model_info(model_name) 
     414    except Exception: 
     415        output = traceback.format_exc() 
     416        return output 
     417 
     418    success, output = check_model(model_info) 
     419    return output 
     420 
     421def check_model(model_info): 
     422    # type: (ModelInfo) -> str 
     423    """ 
     424    Run the tests for a single model, capturing the output. 
     425 
     426    Returns success status and the output string. 
    357427    """ 
    358428    # Note that running main() directly did not work from within the 
     
    368438 
    369439    # Build a test suite containing just the model 
    370     loaders = ['opencl'] if use_opencl() else ['dll'] 
    371     models = [model] 
    372     try: 
    373         suite = make_suite(loaders, models) 
    374     except Exception: 
    375         import traceback 
    376         stream.writeln(traceback.format_exc()) 
    377         return 
    378     # Run the test suite 
    379     suite.run(result) 
    380  
    381     # Print the failures and errors 
    382     for _, tb in result.errors: 
    383         stream.writeln(tb) 
    384     for _, tb in result.failures: 
    385         stream.writeln(tb) 
     440    loaders = ['opencl' if use_opencl() else 'cuda' if use_cuda() else 'dll'] 
     441    suite = unittest.TestSuite() 
     442    _add_model_to_suite(loaders, suite, model_info) 
    386443 
    387444    # Warn if there are no user defined tests. 
     
    393450    # iterator since we don't have direct access to the list of tests in the 
    394451    # test suite. 
     452    # In Qt5 suite.run() will clear all tests in the suite after running 
     453    # with no way of retaining them for the test below, so let's check 
     454    # for user tests before running the suite. 
    395455    for test in suite: 
    396456        if not test.info.tests: 
    397             stream.writeln("Note: %s has no user defined tests."%model) 
     457            stream.writeln("Note: %s has no user defined tests."%model_info.name) 
    398458        break 
    399459    else: 
    400460        stream.writeln("Note: no test suite created --- this should never happen") 
    401461 
     462    # Run the test suite 
     463    suite.run(result) 
     464 
     465    # Print the failures and errors 
     466    for _, tb in result.errors: 
     467        stream.writeln(tb) 
     468    for _, tb in result.failures: 
     469        stream.writeln(tb) 
     470 
    402471    output = stream.getvalue() 
    403472    stream.close() 
    404     return output 
     473    return result.wasSuccessful(), output 
    405474 
    406475 
     
    430499        loaders = ['opencl'] 
    431500        models = models[1:] 
     501    elif models and models[0] == 'cuda': 
     502        if not use_cuda(): 
     503            print("cuda is not available") 
     504            return 1 
     505        loaders = ['cuda'] 
     506        models = models[1:] 
    432507    elif models and models[0] == 'dll': 
    433508        # TODO: test if compiler is available? 
    434509        loaders = ['dll'] 
    435510        models = models[1:] 
    436     elif models and models[0] == 'opencl_and_dll': 
    437         loaders = ['opencl', 'dll'] if use_opencl() else ['dll'] 
    438         models = models[1:] 
    439511    else: 
    440         loaders = ['opencl', 'dll'] if use_opencl() else ['dll'] 
     512        loaders = ['dll'] 
     513        if use_opencl(): 
     514            loaders.append('opencl') 
     515        if use_cuda(): 
     516            loaders.append('cuda') 
    441517    if not models: 
    442518        print("""\ 
    443519usage: 
    444   python -m sasmodels.model_test [-v] [opencl|dll] model1 model2 ... 
     520  python -m sasmodels.model_test [-v] [opencl|cuda|dll] model1 model2 ... 
    445521 
    446522If -v is included on the command line, then use verbose output. 
    447523 
    448 If neither opencl nor dll is specified, then models will be tested with 
    449 both OpenCL and dll; the compute target is ignored for pure python models. 
     524If no platform is specified, then models will be tested with dll, and 
     525if available, OpenCL and CUDA; the compute target is ignored for pure python models. 
    450526 
    451527If model1 is 'all', then all except the remaining models will be tested. 
     
    467543    Run "nosetests sasmodels" on the command line to invoke it. 
    468544    """ 
    469     loaders = ['opencl', 'dll'] if use_opencl() else ['dll'] 
     545    loaders = ['dll'] 
     546    if use_opencl(): 
     547        loaders.append('opencl') 
     548    if use_cuda(): 
     549        loaders.append('cuda') 
    470550    tests = make_suite(loaders, ['all']) 
    471551    def build_test(test): 
  • sasmodels/modelinfo.py

    r1662ebe r765d025  
    5050# Note that scale and background cannot be coordinated parameters whose value 
    5151# depends on the some polydisperse parameter with the current implementation 
     52DEFAULT_BACKGROUND = 1e-3 
    5253COMMON_PARAMETERS = [ 
    5354    ("scale", "", 1, (0.0, np.inf), "", "Source intensity"), 
    54     ("background", "1/cm", 1e-3, (-np.inf, np.inf), "", "Source background"), 
     55    ("background", "1/cm", DEFAULT_BACKGROUND, (-np.inf, np.inf), "", "Source background"), 
    5556] 
    5657assert (len(COMMON_PARAMETERS) == 2 
     
    168169    parameter.length = length 
    169170    parameter.length_control = control 
    170  
    171171    return parameter 
    172172 
     
    269269    will have a magnitude and a direction, which may be different from 
    270270    other sld parameters. The volume parameters are used for calls 
    271     to form_volume within the kernel (required for volume normalization) 
    272     and for calls to ER and VR for effective radius and volume ratio 
    273     respectively. 
     271    to form_volume within the kernel (required for volume normalization), 
     272    to shell_volume (for hollow shapes), and to effective_radius (for 
     273    structure factor interactions) respectively. 
    274274 
    275275    *description* is a short description of the parameter.  This will 
     
    428428        self.kernel_parameters = parameters 
    429429        self._set_vector_lengths() 
    430  
    431430        self.npars = sum(p.length for p in self.kernel_parameters) 
    432431        self.nmagnetic = sum(p.length for p in self.kernel_parameters 
     
    435434        if self.nmagnetic: 
    436435            self.nvalues += 3 + 3*self.nmagnetic 
    437  
    438436        self.call_parameters = self._get_call_parameters() 
    439437        self.defaults = self._get_defaults() 
     
    470468        self.is_asymmetric = any(p.name == 'psi' for p in self.kernel_parameters) 
    471469        self.magnetism_index = [k for k, p in enumerate(self.call_parameters) 
    472                                 if p.id.startswith('M0:')] 
     470                                if p.id.endswith('_M0')] 
    473471 
    474472        self.pd_1d = set(p.name for p in self.call_parameters 
     
    590588        if self.nmagnetic > 0: 
    591589            full_list.extend([ 
    592                 Parameter('up:frac_i', '', 0., [0., 1.], 
     590                Parameter('up_frac_i', '', 0., [0., 1.], 
    593591                          'magnetic', 'fraction of spin up incident'), 
    594                 Parameter('up:frac_f', '', 0., [0., 1.], 
     592                Parameter('up_frac_f', '', 0., [0., 1.], 
    595593                          'magnetic', 'fraction of spin up final'), 
    596                 Parameter('up:angle', 'degress', 0., [0., 360.], 
     594                Parameter('up_angle', 'degrees', 0., [0., 360.], 
    597595                          'magnetic', 'spin up angle'), 
    598596            ]) 
     
    600598            for p in slds: 
    601599                full_list.extend([ 
    602                     Parameter('M0:'+p.id, '1e-6/Ang^2', 0., [-np.inf, np.inf], 
     600                    Parameter(p.id+'_M0', '1e-6/Ang^2', 0., [-np.inf, np.inf], 
    603601                              'magnetic', 'magnetic amplitude for '+p.description), 
    604                     Parameter('mtheta:'+p.id, 'degrees', 0., [-90., 90.], 
     602                    Parameter(p.id+'_mtheta', 'degrees', 0., [-90., 90.], 
    605603                              'magnetic', 'magnetic latitude for '+p.description), 
    606                     Parameter('mphi:'+p.id, 'degrees', 0., [-180., 180.], 
     604                    Parameter(p.id+'_mphi', 'degrees', 0., [-180., 180.], 
    607605                              'magnetic', 'magnetic longitude for '+p.description), 
    608606                ]) 
     
    644642 
    645643        Parameters marked as sld will automatically have a set of associated 
    646         magnetic parameters (m0:p, mtheta:p, mphi:p), as well as polarization 
    647         information (up:theta, up:frac_i, up:frac_f). 
     644        magnetic parameters (p_M0, p_mtheta, p_mphi), as well as polarization 
     645        information (up_theta, up_frac_i, up_frac_f). 
    648646        """ 
    649647        # control parameters go first 
     
    672670            result.append(expanded_pars[name]) 
    673671            if is2d: 
    674                 for tag in 'M0:', 'mtheta:', 'mphi:': 
    675                     if tag+name in expanded_pars: 
    676                         result.append(expanded_pars[tag+name]) 
     672                for tag in '_M0', '_mtheta', '_mphi': 
     673                    if name+tag in expanded_pars: 
     674                        result.append(expanded_pars[name+tag]) 
    677675 
    678676        # Gather the user parameters in order 
     
    707705                append_group(p.id) 
    708706 
    709         if is2d and 'up:angle' in expanded_pars: 
     707        if is2d and 'up_angle' in expanded_pars: 
    710708            result.extend([ 
    711                 expanded_pars['up:frac_i'], 
    712                 expanded_pars['up:frac_f'], 
    713                 expanded_pars['up:angle'], 
     709                expanded_pars['up_frac_i'], 
     710                expanded_pars['up_frac_f'], 
     711                expanded_pars['up_angle'], 
    714712            ]) 
    715713 
     
    726724 
    727725#: Set of variables defined in the model that might contain C code 
    728 C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'c_code'] 
     726C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'shell_volume', 'c_code'] 
    729727 
    730728def _find_source_lines(model_info, kernel_module): 
     
    775773        # Custom sum/multi models 
    776774        return kernel_module.model_info 
     775 
    777776    info = ModelInfo() 
    778777    #print("make parameter table", kernel_module.parameters) 
     
    796795    info.category = getattr(kernel_module, 'category', None) 
    797796    info.structure_factor = getattr(kernel_module, 'structure_factor', False) 
     797    # TODO: find Fq by inspection 
     798    info.effective_radius_type = getattr(kernel_module, 'effective_radius_type', None) 
     799    info.have_Fq = getattr(kernel_module, 'have_Fq', False) 
    798800    info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y']) 
     801    # Note: custom.load_custom_kernel_module assumes the C sources are defined 
     802    # by this attribute. 
    799803    info.source = getattr(kernel_module, 'source', []) 
    800804    info.c_code = getattr(kernel_module, 'c_code', None) 
     805    info.effective_radius = getattr(kernel_module, 'effective_radius', None) 
    801806    # TODO: check the structure of the tests 
    802807    info.tests = getattr(kernel_module, 'tests', []) 
    803     info.ER = getattr(kernel_module, 'ER', None) # type: ignore 
    804     info.VR = getattr(kernel_module, 'VR', None) # type: ignore 
    805808    info.form_volume = getattr(kernel_module, 'form_volume', None) # type: ignore 
     809    info.shell_volume = getattr(kernel_module, 'shell_volume', None) # type: ignore 
    806810    info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore 
    807811    info.Iqxy = getattr(kernel_module, 'Iqxy', None) # type: ignore 
     
    811815    info.profile = getattr(kernel_module, 'profile', None) # type: ignore 
    812816    info.sesans = getattr(kernel_module, 'sesans', None) # type: ignore 
    813     # Default single and opencl to True for C models.  Python models have callable Iq. 
    814817    info.random = getattr(kernel_module, 'random', None) 
    815818 
     
    824827    if getattr(kernel_module, 'py2c', False): 
    825828        try: 
    826             autoc.convert(info, kernel_module) 
     829            warnings = autoc.convert(info, kernel_module) 
    827830        except Exception as exc: 
    828             logger.warn(str(exc) + " while converting %s from C to python"%name) 
    829  
     831            warnings = [str(exc)] 
     832        if warnings: 
     833            warnings.append("while converting %s from C to python"%name) 
     834            if len(warnings) > 2: 
     835                warnings = "\n".join(warnings) 
     836            else: 
     837                warnings = " ".join(warnings) 
     838            logger.warn(warnings) 
     839 
     840    # Default single and opencl to True for C models.  Python models have callable Iq. 
    830841    # Needs to come after autoc.convert since the Iq symbol may have been 
    831842    # converted from python to C 
     
    928939    #: provided in the file. 
    929940    structure_factor = None # type: bool 
     941    #: True if the model defines an Fq function with signature 
     942    #: void Fq(double q, double *F1, double *F2, ...) 
     943    have_Fq = False 
     944    #: List of options for computing the effective radius of the shape, 
     945    #: or None if the model is not usable as a form factor model. 
     946    effective_radius_type = None   # type: List[str] 
    930947    #: List of C source files used to define the model.  The source files 
    931948    #: should define the *Iq* function, and possibly *Iqac* or *Iqabc* if the 
    932949    #: model defines orientation parameters. Files containing the most basic 
    933950    #: functions must appear first in the list, followed by the files that 
    934     #: use those functions.  Form factors are indicated by providing 
    935     #: an :attr:`ER` function. 
     951    #: use those functions. 
    936952    source = None           # type: List[str] 
    937     #: The set of tests that must pass.  The format of the tests is described 
    938     #: in :mod:`model_test`. 
    939     tests = None            # type: List[TestCondition] 
    940     #: Returns the effective radius of the model given its volume parameters. 
    941     #: The presence of *ER* indicates that the model is a form factor model 
    942     #: that may be used together with a structure factor to form an implicit 
    943     #: multiplication model. 
    944     #: 
    945     #: The parameters to the *ER* function must be marked with type *volume*. 
    946     #: in the parameter table.  They will appear in the same order as they 
    947     #: do in the table.  The values passed to *ER* will be vectors, with one 
    948     #: value for each polydispersity condition.  For example, if the model 
    949     #: is polydisperse over both length and radius, then both length and 
    950     #: radius will have the same number of values in the vector, with one 
    951     #: value for each *length X radius*.  If only *radius* is polydisperse, 
    952     #: then the value for *length* will be repeated once for each value of 
    953     #: *radius*.  The *ER* function should return one effective radius for 
    954     #: each parameter set.  Multiplicity parameters will be received as 
    955     #: arrays, with one row per polydispersity condition. 
    956     ER = None               # type: Optional[Callable[[np.ndarray], np.ndarray]] 
    957     #: Returns the occupied volume and the total volume for each parameter set. 
    958     #: See :attr:`ER` for details on the parameters. 
    959     VR = None               # type: Optional[Callable[[np.ndarray], Tuple[np.ndarray, np.ndarray]]] 
    960     #: Arbitrary C code containing supporting functions, etc., to be inserted 
    961     #: after everything in source.  This can include Iq and Iqxy functions with 
    962     #: the full function signature, including all parameters. 
    963     c_code = None 
     953    #: inline source code, added after all elements of source 
     954    c_code = None           # type: Optional[str] 
    964955    #: Returns the form volume for python-based models.  Form volume is needed 
    965956    #: for volume normalization in the polydispersity integral.  If no 
     
    969960    #: C code, either defined as a string, or in the sources. 
    970961    form_volume = None      # type: Union[None, str, Callable[[np.ndarray], float]] 
     962    #: Returns the shell volume for python-based models.  Form volume and 
     963    #: shell volume are needed for volume normalization in the polydispersity 
     964    #: integral and structure interactions for hollow shapes.  If no 
     965    #: parameters are *volume* parameters, then shell volume is not needed. 
     966    #: For C-based models, (with :attr:`sources` defined, or with :attr:`Iq` 
     967    #: defined using a string containing C code), shell_volume must also be 
     968    #: C code, either defined as a string, or in the sources. 
     969    shell_volume = None      # type: Union[None, str, Callable[[np.ndarray], float]] 
    971970    #: Returns *I(q, a, b, ...)* for parameters *a*, *b*, etc. defined 
    972971    #: by the parameter table.  *Iq* can be defined as a python function, or 
     
    10031002    #: Line numbers for symbols defining C code 
    10041003    lineno = None           # type: Dict[str, int] 
     1004    #: The set of tests that must pass.  The format of the tests is described 
     1005    #: in :mod:`model_test`. 
     1006    tests = None            # type: List[TestCondition] 
    10051007 
    10061008    def __init__(self): 
     
    10261028                         for k in range(control+1, p.length+1) 
    10271029                         if p.length > 1) 
     1030            for p in self.parameters.kernel_parameters: 
     1031                if p.length > 1 and p.type == "sld": 
     1032                    for k in range(control+1, p.length+1): 
     1033                        base = p.id+str(k) 
     1034                        hidden.update((base+"_M0", base+"_mtheta", base+"_mphi")) 
    10281035        return hidden 
  • sasmodels/models/_cylpy.py

    rc01ed3e r765d025  
    177177    return 1.0e-4 * square(s * form) 
    178178 
    179 def ER(radius, length): 
    180     """ 
    181     Return equivalent radius (ER) 
    182     """ 
    183     ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius)) 
    184     return 0.5 * (ddd) ** (1. / 3.) 
     179def effective_radius(mode, radius, length): 
     180    mode = int(mode) 
     181    if mode == 1: 
     182        ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius)) 
     183        return 0.5 * (ddd) ** (1. / 3.) 
     184    else: 
     185        return 0. 
    185186 
    186187def random(): 
  • sasmodels/models/_spherepy.py

    r108e70e r304c775  
    11r""" 
    22For information about polarised and magnetic scattering, see 
    3 the :doc:`magnetic help <../sasgui/perspectives/fitting/mag_help>` documentation. 
     3the :ref:`magnetism` documentation. 
    44 
    55Definition 
     
    7171    return 1.333333333333333 * pi * radius ** 3 
    7272 
     73def effective_radius(mode, radius): 
     74    return radius 
     75 
    7376def Iq(q, sld, sld_solvent, radius): 
    7477    #print "q",q 
     
    105108sesans.vectorized = True  # sesans accepts an array of z values 
    106109 
    107 def ER(radius): 
    108     return radius 
    109  
    110 # VR defaults to 1.0 
    111  
    112110demo = dict(scale=1, background=0, 
    113111            sld=6, sld_solvent=1, 
  • sasmodels/models/barbell.c

    r108e70e rd42dd4a  
    6363 
    6464static double 
    65 Iq(double q, double sld, double solvent_sld, 
     65radius_from_volume(double radius_bell, double radius, double length) 
     66{ 
     67    const double vol_barbell = form_volume(radius_bell,radius,length); 
     68    return cbrt(vol_barbell/M_4PI_3); 
     69} 
     70 
     71static double 
     72radius_from_totallength(double radius_bell, double radius, double length) 
     73{ 
     74    const double hdist = sqrt(square(radius_bell) - square(radius)); 
     75    return 0.5*length + hdist + radius_bell; 
     76} 
     77 
     78static double 
     79effective_radius(int mode, double radius_bell, double radius, double length) 
     80{ 
     81    switch (mode) { 
     82    default: 
     83    case 1: // equivalent sphere 
     84        return radius_from_volume(radius_bell, radius , length); 
     85    case 2: // radius 
     86        return radius; 
     87    case 3: // half length 
     88        return 0.5*length; 
     89    case 4: // half total length 
     90        return radius_from_totallength(radius_bell,radius,length); 
     91    } 
     92} 
     93 
     94static void 
     95Fq(double q,double *F1, double *F2, double sld, double solvent_sld, 
    6696    double radius_bell, double radius, double length) 
    6797{ 
     
    72102    const double zm = M_PI_4; 
    73103    const double zb = M_PI_4; 
    74     double total = 0.0; 
     104    double total_F1 = 0.0; 
     105    double total_F2 = 0.0; 
    75106    for (int i = 0; i < GAUSS_N; i++){ 
    76107        const double alpha= GAUSS_Z[i]*zm + zb; 
     
    78109        SINCOS(alpha, sin_alpha, cos_alpha); 
    79110        const double Aq = _fq(q*sin_alpha, q*cos_alpha, h, radius_bell, radius, half_length); 
    80         total += GAUSS_W[i] * Aq * Aq * sin_alpha; 
     111        total_F1 += GAUSS_W[i] * Aq * sin_alpha; 
     112        total_F2 += GAUSS_W[i] * Aq * Aq * sin_alpha; 
    81113    } 
    82114    // translate dx in [-1,1] to dx in [lower,upper] 
    83     const double form = total*zm; 
     115    const double form_avg = total_F1*zm; 
     116    const double form_squared_avg = total_F2*zm; 
    84117 
    85118    //Contrast 
    86119    const double s = (sld - solvent_sld); 
    87     return 1.0e-4 * s * s * form; 
     120    *F1 = 1.0e-2 * s * form_avg; 
     121    *F2 = 1.0e-4 * s * s * form_squared_avg; 
    88122} 
    89  
    90123 
    91124static double 
  • sasmodels/models/barbell.py

    r2d81cfe ree60aa7  
    116116 
    117117source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "barbell.c"] 
     118have_Fq = True 
     119effective_radius_type = [ 
     120    "equivalent sphere", "radius", "half length", "half total length", 
     121    ] 
    118122 
    119123def random(): 
  • sasmodels/models/bcc_paracrystal.py

    r2d81cfe rda7b26b  
    11r""" 
     2.. warning:: This model and this model description are under review following  
     3             concerns raised by SasView users. If you need to use this model,  
     4             please email help@sasview.org for the latest situation. *The  
     5             SasView Developers. September 2018.* 
     6 
    27Definition 
    38---------- 
     
    1318 
    1419    I(q) = \frac{\text{scale}}{V_p} V_\text{lattice} P(q) Z(q) 
    15  
    1620 
    1721where *scale* is the volume fraction of spheres, $V_p$ is the volume of the 
     
    97101 
    98102Authorship and Verification 
    99 ---------------------------- 
     103--------------------------- 
    100104 
    101105* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
  • sasmodels/models/be_polyelectrolyte.py

    ref07e95 rca77fc1  
    11r""" 
     2.. note:: Please read the Validation section below. 
     3 
    24Definition 
    35---------- 
     
    1113 
    1214    I(q) = K\frac{q^2+k^2}{4\pi L_b\alpha ^2} 
    13     \frac{1}{1+r_{0}^2(q^2+k^2)(q^2-12hC_a/b^2)} + background 
     15    \frac{1}{1+r_{0}^4(q^2+k^2)(q^2-12hC_a/b^2)} + background 
    1416 
    1517    k^2 = 4\pi L_b(2C_s + \alpha C_a) 
    1618 
    17     r_{0}^2 = \frac{1}{\alpha \sqrt{C_a} \left( b/\sqrt{48\pi L_b}\right)} 
     19    r_{0}^2 = \frac{b}{\alpha \sqrt{C_a 48\pi L_b}} 
    1820 
    1921where 
    2022 
    2123$K$ is the contrast factor for the polymer which is defined differently than in 
    22 other models and is given in barns where $1 barn = 10^{-24} cm^2$.  $K$ is 
     24other models and is given in barns where 1 $barn = 10^{-24}$ $cm^2$.  $K$ is 
    2325defined as: 
    2426 
     
    2931    a = b_p - (v_p/v_s) b_s 
    3032 
    31 where $b_p$ and $b_s$ are sum of the scattering lengths of the atoms 
    32 constituting the monomer of the polymer and the sum of the scattering lengths 
    33 of the atoms constituting the solvent molecules respectively, and $v_p$ and 
    34 $v_s$ are the partial molar volume of the polymer and the solvent respectively 
    35  
    36 $L_b$ is the Bjerrum length(|Ang|) - **Note:** This parameter needs to be 
    37 kept constant for a given solvent and temperature! 
    38  
    39 $h$ is the virial parameter (|Ang^3|/mol) - **Note:** See [#Borue]_ for the 
    40 correct interpretation of this parameter.  It incorporates second and third 
    41 virial coefficients and can be Negative. 
    42  
    43 $b$ is the monomer length(|Ang|), $C_s$ is the concentration of monovalent 
    44 salt(mol/L), $\alpha$ is the ionization degree (ionization degree : ratio of 
    45 charged monomers  to total number of monomers), $C_a$ is the polymer molar 
    46 concentration(mol/L), and $background$ is the incoherent background. 
     33where: 
     34 
     35- $b_p$ and $b_s$ are **sum of the scattering lengths of the atoms** 
     36  constituting the polymer monomer and the solvent molecules, respectively. 
     37 
     38- $v_p$ and $v_s$ are the partial molar volume of the polymer and the  
     39  solvent, respectively. 
     40 
     41- $L_b$ is the Bjerrum length (|Ang|) - **Note:** This parameter needs to be 
     42  kept constant for a given solvent and temperature! 
     43 
     44- $h$ is the virial parameter (|Ang^3|) - **Note:** See [#Borue]_ for the 
     45  correct interpretation of this parameter.  It incorporates second and third 
     46  virial coefficients and can be *negative*. 
     47 
     48- $b$ is the monomer length (|Ang|). 
     49 
     50- $C_s$ is the concentration of monovalent salt(1/|Ang^3| - internally converted from mol/L). 
     51 
     52- $\alpha$ is the degree of ionization (the ratio of charged monomers to the total  
     53  number of monomers) 
     54 
     55- $C_a$ is the polymer molar concentration (1/|Ang^3| - internally converted from mol/L) 
     56 
     57- $background$ is the incoherent background. 
    4758 
    4859For 2D data the scattering intensity is calculated in the same way as 1D, 
     
    5263 
    5364    q = \sqrt{q_x^2 + q_y^2} 
     65 
     66Validation 
     67---------- 
     68 
     69As of the last revision, this code is believed to be correct.  However it 
     70needs further validation and should be used with caution at this time.  The 
     71history of this code goes back to a 1998 implementation. It was recently noted 
     72that in that implementation, while both the polymer concentration and salt 
     73concentration were converted from experimental units of mol/L to more 
     74dimensionally useful units of 1/|Ang^3|, only the converted version of the 
     75polymer concentration was actually being used in the calculation while the 
     76unconverted salt concentration (still in apparent units of mol/L) was being  
     77used.  This was carried through to Sasmodels as used for SasView 4.1 (though  
     78the line of code converting the salt concentration to the new units was removed  
     79somewhere along the line). Simple dimensional analysis of the calculation shows  
     80that the converted salt concentration should be used, which the original code  
     81suggests was the intention, so this has now been corrected (for SasView 4.2).  
     82Once better validation has been performed this note will be removed. 
    5483 
    5584References 
     
    6695 
    6796* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    68 * **Last Modified by:** Paul Kienzle **Date:** July 24, 2016 
    69 * **Last Reviewed by:** Paul Butler and Richard Heenan **Date:** October 07, 2016 
     97* **Last Modified by:** Paul Butler **Date:** September 25, 2018 
     98* **Last Reviewed by:** Paul Butler **Date:** September 25, 2018 
    7099""" 
    71100 
     
    92121    ["contrast_factor",       "barns",   10.0,  [-inf, inf], "", "Contrast factor of the polymer"], 
    93122    ["bjerrum_length",        "Ang",      7.1,  [0, inf],    "", "Bjerrum length"], 
    94     ["virial_param",          "Ang^3/mol", 12.0,  [-inf, inf], "", "Virial parameter"], 
     123    ["virial_param",          "Ang^3", 12.0,  [-inf, inf], "", "Virial parameter"], 
    95124    ["monomer_length",        "Ang",     10.0,  [0, inf],    "", "Monomer length"], 
    96125    ["salt_concentration",    "mol/L",    0.0,  [-inf, inf], "", "Concentration of monovalent salt"], 
     
    102131 
    103132def Iq(q, 
    104        contrast_factor=10.0, 
    105        bjerrum_length=7.1, 
    106        virial_param=12.0, 
    107        monomer_length=10.0, 
    108        salt_concentration=0.0, 
    109        ionization_degree=0.05, 
    110        polymer_concentration=0.7): 
     133       contrast_factor, 
     134       bjerrum_length, 
     135       virial_param, 
     136       monomer_length, 
     137       salt_concentration, 
     138       ionization_degree, 
     139       polymer_concentration): 
    111140    """ 
    112     :param q:                     Input q-value 
    113     :param contrast_factor:       Contrast factor of the polymer 
    114     :param bjerrum_length:        Bjerrum length 
    115     :param virial_param:          Virial parameter 
    116     :param monomer_length:        Monomer length 
    117     :param salt_concentration:    Concentration of monovalent salt 
    118     :param ionization_degree:     Degree of ionization 
    119     :param polymer_concentration: Polymer molar concentration 
    120     :return:                      1-D intensity 
     141    :params: see parameter table 
     142    :return: 1-D form factor for polyelectrolytes in low salt 
     143     
     144    parameter names, units, default values, and behavior (volume, sld etc) are 
     145    defined in the parameter table.  The concentrations are converted from 
     146    experimental mol/L to dimensionaly useful 1/A3 in first two lines 
    121147    """ 
    122148 
    123     concentration = polymer_concentration * 6.022136e-4 
    124  
    125     k_square = 4.0 * pi * bjerrum_length * (2*salt_concentration + 
    126                                             ionization_degree * concentration) 
    127  
    128     r0_square = 1.0/ionization_degree/sqrt(concentration) * \ 
     149    concentration_pol = polymer_concentration * 6.022136e-4 
     150    concentration_salt = salt_concentration * 6.022136e-4 
     151 
     152    k_square = 4.0 * pi * bjerrum_length * (2*concentration_salt + 
     153                                            ionization_degree * concentration_pol) 
     154 
     155    r0_square = 1.0/ionization_degree/sqrt(concentration_pol) * \ 
    129156                (monomer_length/sqrt((48.0*pi*bjerrum_length))) 
    130157 
     
    133160 
    134161    term2 = 1.0 + r0_square**2 * (q**2 + k_square) * \ 
    135         (q**2 - (12.0 * virial_param * concentration/(monomer_length**2))) 
     162        (q**2 - (12.0 * virial_param * concentration_pol/(monomer_length**2))) 
    136163 
    137164    return term1/term2 
     
    174201 
    175202    # Accuracy tests based on content in test/utest_other_models.py 
     203    # Note that these should some day be validated beyond this self validation 
     204    # (circular reasoning). -- i.e. the "good value," at least for those with 
     205    # non zero salt concentrations, were obtained by running the current 
     206    # model in SasView and copying the appropriate result here. 
     207    #    PDB -- Sep 26, 2018 
    176208    [{'contrast_factor':       10.0, 
    177209      'bjerrum_length':         7.1, 
     
    184216     }, 0.001, 0.0948379], 
    185217 
    186     # Additional tests with larger range of parameters 
    187218    [{'contrast_factor':       10.0, 
    188219      'bjerrum_length':       100.0, 
    189220      'virial_param':           3.0, 
    190       'monomer_length':         1.0, 
    191       'salt_concentration':    10.0, 
    192       'ionization_degree':      2.0, 
    193       'polymer_concentration': 10.0, 
     221      'monomer_length':         5.0, 
     222      'salt_concentration':     1.0, 
     223      'ionization_degree':      0.1, 
     224      'polymer_concentration':  1.0, 
    194225      'background':             0.0, 
    195      }, 0.1, -3.75693800588], 
     226     }, 0.1, 0.253469484], 
    196227 
    197228    [{'contrast_factor':       10.0, 
    198229      'bjerrum_length':       100.0, 
    199230      'virial_param':           3.0, 
    200       'monomer_length':         1.0, 
    201       'salt_concentration':    10.0, 
    202       'ionization_degree':      2.0, 
    203       'polymer_concentration': 10.0, 
    204       'background':           100.0 
    205      }, 5.0, 100.029142149], 
     231      'monomer_length':         5.0, 
     232      'salt_concentration':     1.0, 
     233      'ionization_degree':      0.1, 
     234      'polymer_concentration':  1.0, 
     235      'background':             1.0, 
     236     }, 0.05, 1.738358122], 
    206237 
    207238    [{'contrast_factor':     100.0, 
    208239      'bjerrum_length':       10.0, 
    209       'virial_param':        180.0, 
    210       'monomer_length':        1.0, 
     240      'virial_param':         12.0, 
     241      'monomer_length':       10.0, 
    211242      'salt_concentration':    0.1, 
    212243      'ionization_degree':     0.5, 
    213244      'polymer_concentration': 0.1, 
    214       'background':             0.0, 
    215      }, 200., 1.80664667511e-06], 
     245      'background':           0.01, 
     246     }, 0.5, 0.012881893], 
    216247    ] 
  • sasmodels/models/binary_hard_sphere.c

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

    r108e70e rd42dd4a  
    8585 
    8686static double 
    87 Iq(double q, double sld, double solvent_sld, 
     87radius_from_volume(double radius, double radius_cap, double length) 
     88{ 
     89    const double vol_cappedcyl = form_volume(radius,radius_cap,length); 
     90    return cbrt(vol_cappedcyl/M_4PI_3); 
     91} 
     92 
     93static double 
     94radius_from_totallength(double radius, double radius_cap, double length) 
     95{ 
     96    const double hc = radius_cap - sqrt(radius_cap*radius_cap - radius*radius); 
     97    return 0.5*length + hc; 
     98} 
     99 
     100static double 
     101effective_radius(int mode, double radius, double radius_cap, double length) 
     102{ 
     103    switch (mode) { 
     104    default: 
     105    case 1: // equivalent sphere 
     106        return radius_from_volume(radius, radius_cap, length); 
     107    case 2: // radius 
     108        return radius; 
     109    case 3: // half length 
     110        return 0.5*length; 
     111    case 4: // half total length 
     112        return radius_from_totallength(radius, radius_cap,length); 
     113    } 
     114} 
     115 
     116static void 
     117Fq(double q,double *F1, double *F2, double sld, double solvent_sld, 
    88118    double radius, double radius_cap, double length) 
    89119{ 
     
    94124    const double zm = M_PI_4; 
    95125    const double zb = M_PI_4; 
    96     double total = 0.0; 
     126    double total_F1 = 0.0; 
     127    double total_F2 = 0.0; 
    97128    for (int i=0; i<GAUSS_N ;i++) { 
    98129        const double theta = GAUSS_Z[i]*zm + zb; 
     
    103134        const double Aq = _fq(qab, qc, h, radius_cap, radius, half_length); 
    104135        // scale by sin_theta for spherical coord integration 
    105         total += GAUSS_W[i] * Aq * Aq * sin_theta; 
     136        total_F1 += GAUSS_W[i] * Aq * sin_theta; 
     137        total_F2 += GAUSS_W[i] * Aq * Aq * sin_theta; 
    106138    } 
    107139    // translate dx in [-1,1] to dx in [lower,upper] 
    108     const double form = total * zm; 
     140    const double form_avg = total_F1 * zm; 
     141    const double form_squared_avg = total_F2 * zm; 
    109142 
    110143    // Contrast 
    111144    const double s = (sld - solvent_sld); 
    112     return 1.0e-4 * s * s * form; 
     145    *F1 = 1.0e-2 * s * form_avg; 
     146    *F2 = 1.0e-4 * s * s * form_squared_avg; 
    113147} 
    114148 
  • sasmodels/models/capped_cylinder.py

    r2d81cfe ree60aa7  
    136136 
    137137source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "capped_cylinder.c"] 
     138have_Fq = True 
     139effective_radius_type = [ 
     140    "equivalent sphere", "radius", "half length", "half total length", 
     141    ] 
    138142 
    139143def random(): 
  • sasmodels/models/core_multi_shell.c

    rc3ccaec rd42dd4a  
    99 
    1010static double 
    11 form_volume(double core_radius, double fp_n, double thickness[]) 
     11outer_radius(double core_radius, double fp_n, double thickness[]) 
    1212{ 
    1313  double r = core_radius;