Changeset cc8b183 in sasmodels


Ignore:
Timestamp:
Nov 9, 2018 2:12:05 PM (4 months ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
ticket_1156
Children:
f752b9b
Parents:
b3f4831 (diff), cf3d0ce (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'beta_approx' into ticket_1156

Files:
50 added
114 edited

Legend:

Unmodified
Added
Removed
  • doc/guide/gpu_setup.rst

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

    r57c609b 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 
     
    492492used. 
    493493 
     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 
    494506Embedded C Models 
    495507................. 
     
    503515This expands into the equivalent C code:: 
    504516 
    505     #include <math.h> 
    506517    double Iq(double q, double par1, double par2, ...); 
    507518    double Iq(double q, double par1, double par2, ...) 
     
    512523*form_volume* defines the volume of the shape. As in python models, it 
    513524includes 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. 
    514528 
    515529**source=['fn.c', ...]** includes the listed C source files in the 
     
    548562The INVALID define can go into *Iq*, or *c_code*, or an external C file 
    549563listed 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. 
    550593 
    551594Oriented Shapes 
     
    10121055          "radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, 
    10131056         0.2, 0.228843], 
    1014         [{"radius": 120., "radius_pd": 0.2, "radius_pd_n":45}, "ER", 120.], 
    1015         [{"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], 
    10161060    ] 
    10171061 
    10181062 
    1019 **tests=[[{parameters}, q, result], ...]** is a list of lists. 
     1063**tests=[[{parameters}, q, Iq], ...]** is a list of lists. 
    10201064Each list is one test and contains, in order: 
    10211065 
     
    10291073- input and output values can themselves be lists if you have several 
    10301074  $q$ values to test for the same model parameters. 
    1031 - for testing *ER* and *VR*, give the inputs as "ER" and "VR" respectively; 
    1032   the output for *VR* should be the sphere/shell ratio, not the individual 
    1033   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. 
    10341084 
    10351085.. _Test_Your_New_Model: 
  • doc/guide/scripting.rst

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

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

    rfba9ca0 raa8c6e0  
    403403) 
    404404add_function( 
    405     name="debye", 
     405    name="gauss_coil", 
    406406    mp_function=lambda x: 2*(mp.exp(-x**2) + x**2 - 1)/x**4, 
    407407    np_function=lambda x: 2*(np.expm1(-x**2) + x**2)/x**4, 
    408408    ocl_function=make_ocl(""" 
    409409    const double qsq = q*q; 
    410     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) { 
    411413        const double x = qsq; 
    412414        if (0) { // 0.36 single 
     
    418420            const double B1=3./8., B2=3./56., B3=1./336.; 
    419421            return (((A3*x + A2)*x + A1)*x + 1.)/(((B3*x + B2)*x + B1)*x + 1.); 
    420         } else if (1) { // 1.0 for single, 0.25 for double 
     422        } else if (0) { // 1.0 for single, 0.25 for double 
    421423            // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}] 
    422424            const double A1=1./15., A2=1./60, A3=0., A4=1./75600.; 
     
    431433                  /(((((B5*x + B4)*x + B3)*x + B2)*x + B1)*x + 1.); 
    432434        } 
    433     } else if (qsq < 1.) { // Taylor series; 0.9 for single, 0.25 for double 
     435    } else if (qsq < 0.8) { 
    434436        const double x = qsq; 
    435437        const double C0 = +1.; 
  • sasmodels/compare.py

    r610ef23 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 
     
    115116    === environment variables === 
    116117    -DSAS_MODELPATH=path sets directory containing custom models 
    117     -DSAS_OPENCL=vendor:device|none sets the target OpenCL device 
     118    -DSAS_OPENCL=vendor:device|cuda:device|none sets the target GPU device 
    118119    -DXDG_CACHE_HOME=~/.cache sets the pyopencl cache root (linux only) 
    119120    -DSAS_COMPILER=tinycc|msvc|mingw|unix sets the DLL compiler 
     
    352353 
    353354    # If it is a list of choices, pick one at random with equal probability 
    354     # In practice, the model specific random generator will override. 
    355355    par = model_info.parameters[name] 
    356     if len(par.limits) > 2:  # choice list 
    357         return np.random.randint(len(par.limits)) 
     356    if par.choices:  # choice list 
     357        return np.random.randint(len(par.choices)) 
    358358 
    359359    # If it is a fixed range, pick from it with equal probability. 
     
    725725        set_integration_size(model_info, ngauss) 
    726726 
    727     if dtype != "default" and not dtype.endswith('!') and not kernelcl.use_opencl(): 
     727    if (dtype != "default" and not dtype.endswith('!')  
     728            and not (kernelcl.use_opencl() or kernelcuda.use_cuda())): 
    728729        raise RuntimeError("OpenCL not available " + kernelcl.OPENCL_ERROR) 
    729730 
  • sasmodels/core.py

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

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

    r7b9e4dd r5024a56  
    3232from .details import make_kernel_args, dispersion_mesh 
    3333from .modelinfo import DEFAULT_BACKGROUND 
     34from .product import RADIUS_MODE_ID 
    3435 
    3536# pylint: disable=unused-import 
     
    6465    return calculator(call_details, values, cutoff, is_magnetic) 
    6566 
    66 def call_ER(model_info, pars): 
    67     # type: (ModelInfo, ParameterSet) -> float 
    68     """ 
    69     Call the model ER function using *values*. 
    70  
    71     *model_info* is either *model.info* if you have a loaded model, 
    72     or *kernel.info* if you have a model kernel prepared for evaluation. 
    73     """ 
    74     if model_info.ER is None: 
    75         return 1.0 
    76     elif not model_info.parameters.form_volume_parameters: 
    77         # handle the case where ER is provided but model is not polydisperse 
    78         return model_info.ER() 
    79     else: 
    80         value, weight = _vol_pars(model_info, pars) 
    81         individual_radii = model_info.ER(*value) 
    82         return np.sum(weight*individual_radii) / np.sum(weight) 
    83  
    84  
    85 def call_VR(model_info, pars): 
    86     # type: (ModelInfo, ParameterSet) -> float 
    87     """ 
    88     Call the model VR function using *pars*. 
    89  
    90     *model_info* is either *model.info* if you have a loaded model, 
    91     or *kernel.info* if you have a model kernel prepared for evaluation. 
    92     """ 
    93     if model_info.VR is None: 
    94         return 1.0 
    95     elif not model_info.parameters.form_volume_parameters: 
    96         # handle the case where ER is provided but model is not polydisperse 
    97         return model_info.VR() 
    98     else: 
    99         value, weight = _vol_pars(model_info, pars) 
    100         whole, part = model_info.VR(*value) 
    101         return np.sum(weight*part)/np.sum(weight*whole) 
    102  
    103  
    104 def call_profile(model_info, **pars): 
    105     # type: (ModelInfo, ...) -> Tuple[np.ndarray, np.ndarray, Tuple[str, str]] 
     67def call_Fq(calculator, pars, cutoff=0., mono=False): 
     68    # type: (Kernel, ParameterSet, float, bool) -> np.ndarray 
     69    """ 
     70    Like :func:`call_kernel`, but returning F, F^2, R_eff, V_shell, V_form/V_shell. 
     71 
     72    For solid objects V_shell is equal to V_form and the volume ratio is 1. 
     73 
     74    Use parameter *radius_effective_mode* to select the effective radius 
     75    calculation. 
     76    """ 
     77    R_eff_type = int(pars.pop(RADIUS_MODE_ID, 1.0)) 
     78    mesh = get_mesh(calculator.info, pars, dim=calculator.dim, mono=mono) 
     79    #print("pars", list(zip(*mesh))[0]) 
     80    call_details, values, is_magnetic = make_kernel_args(calculator, mesh) 
     81    #print("values:", values) 
     82    return calculator.Fq(call_details, values, cutoff, is_magnetic, R_eff_type) 
     83 
     84def call_profile(model_info, pars=None): 
     85    # type: (ModelInfo, ParameterSet) -> Tuple[np.ndarray, np.ndarray, Tuple[str, str]] 
    10686    """ 
    10787    Returns the profile *x, y, (xlabel, ylabel)* representing the model. 
    10888    """ 
     89    if pars is None: 
     90        pars = {} 
    10991    args = {} 
    11092    for p in model_info.parameters.kernel_parameters: 
     
    403385        Generate a plottable profile. 
    404386        """ 
    405         return call_profile(self.model.info, **pars) 
     387        return call_profile(self.model.info, pars) 
    406388 
    407389def main(): 
     
    419401    model_name = sys.argv[1] 
    420402    call = sys.argv[2].upper() 
    421     if call != "ER_VR": 
    422         try: 
    423             values = [float(v) for v in call.split(',')] 
    424         except ValueError: 
    425             values = [] 
    426         if len(values) == 1: 
    427             q, = values 
    428             data = empty_data1D([q]) 
    429         elif len(values) == 2: 
    430             qx, qy = values 
    431             data = empty_data2D([qx], [qy]) 
    432         else: 
    433             print("use q or qx,qy or ER or VR") 
    434             sys.exit(1) 
     403    try: 
     404        values = [float(v) for v in call.split(',')] 
     405    except ValueError: 
     406        values = [] 
     407    if len(values) == 1: 
     408        q, = values 
     409        data = empty_data1D([q]) 
     410    elif len(values) == 2: 
     411        qx, qy = values 
     412        data = empty_data2D([qx], [qy]) 
    435413    else: 
    436         data = empty_data1D([0.001])  # Data not used in ER/VR 
     414        print("use q or qx,qy") 
     415        sys.exit(1) 
    437416 
    438417    model_info = load_model_info(model_name) 
     
    442421                for pair in sys.argv[3:] 
    443422                for k, v in [pair.split('=')]) 
    444     if call == "ER_VR": 
    445         ER = call_ER(model_info, pars) 
    446         VR = call_VR(model_info, pars) 
    447         print(ER, VR) 
    448     else: 
    449         Iq = calculator(**pars) 
    450         print(Iq[0]) 
     423    Iq = calculator(**pars) 
     424    print(Iq[0]) 
    451425 
    452426if __name__ == "__main__": 
  • sasmodels/generate.py

    r6e45516 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        ] 
  • 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

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

    rd86f0fc 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

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

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

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

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

    rbd547d0 r39a06c9  
    164164    parameter.length = length 
    165165    parameter.length_control = control 
    166  
    167166    return parameter 
    168167 
     
    265264    will have a magnitude and a direction, which may be different from 
    266265    other sld parameters. The volume parameters are used for calls 
    267     to form_volume within the kernel (required for volume normalization) 
    268     and for calls to ER and VR for effective radius and volume ratio 
    269     respectively. 
     266    to form_volume within the kernel (required for volume normalization), 
     267    to shell_volume (for hollow shapes), and to effective_radius (for 
     268    structure factor interactions) respectively. 
    270269 
    271270    *description* is a short description of the parameter.  This will 
     
    424423        self.kernel_parameters = parameters 
    425424        self._set_vector_lengths() 
    426  
    427425        self.npars = sum(p.length for p in self.kernel_parameters) 
    428426        self.nmagnetic = sum(p.length for p in self.kernel_parameters 
     
    431429        if self.nmagnetic: 
    432430            self.nvalues += 3 + 3*self.nmagnetic 
    433  
    434431        self.call_parameters = self._get_call_parameters() 
    435432        self.defaults = self._get_defaults() 
     
    722719 
    723720#: Set of variables defined in the model that might contain C code 
    724 C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'c_code'] 
     721C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'shell_volume', 'c_code'] 
    725722 
    726723def _find_source_lines(model_info, kernel_module): 
     
    771768        # Custom sum/multi models 
    772769        return kernel_module.model_info 
     770 
    773771    info = ModelInfo() 
    774772    #print("make parameter table", kernel_module.parameters) 
     
    792790    info.category = getattr(kernel_module, 'category', None) 
    793791    info.structure_factor = getattr(kernel_module, 'structure_factor', False) 
     792    # TODO: find Fq by inspection 
     793    info.effective_radius_type = getattr(kernel_module, 'effective_radius_type', None) 
     794    info.have_Fq = getattr(kernel_module, 'have_Fq', False) 
    794795    info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y']) 
    795796    # Note: custom.load_custom_kernel_module assumes the C sources are defined 
     
    797798    info.source = getattr(kernel_module, 'source', []) 
    798799    info.c_code = getattr(kernel_module, 'c_code', None) 
     800    info.effective_radius = getattr(kernel_module, 'effective_radius', None) 
    799801    # TODO: check the structure of the tests 
    800802    info.tests = getattr(kernel_module, 'tests', []) 
    801     info.ER = getattr(kernel_module, 'ER', None) # type: ignore 
    802     info.VR = getattr(kernel_module, 'VR', None) # type: ignore 
    803803    info.form_volume = getattr(kernel_module, 'form_volume', None) # type: ignore 
     804    info.shell_volume = getattr(kernel_module, 'shell_volume', None) # type: ignore 
    804805    info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore 
    805806    info.Iqxy = getattr(kernel_module, 'Iqxy', None) # type: ignore 
     
    825826    info.lineno = {} 
    826827    _find_source_lines(info, kernel_module) 
    827  
    828828    return info 
    829829 
     
    918918    #: provided in the file. 
    919919    structure_factor = None # type: bool 
     920    #: True if the model defines an Fq function with signature 
     921    #: void Fq(double q, double *F1, double *F2, ...) 
     922    have_Fq = False 
     923    #: List of options for computing the effective radius of the shape, 
     924    #: or None if the model is not usable as a form factor model. 
     925    effective_radius_type = None   # type: List[str] 
    920926    #: List of C source files used to define the model.  The source files 
    921927    #: should define the *Iq* function, and possibly *Iqac* or *Iqabc* if the 
    922928    #: model defines orientation parameters. Files containing the most basic 
    923929    #: functions must appear first in the list, followed by the files that 
    924     #: use those functions.  Form factors are indicated by providing 
    925     #: an :attr:`ER` function. 
     930    #: use those functions. 
    926931    source = None           # type: List[str] 
    927     #: The set of tests that must pass.  The format of the tests is described 
    928     #: in :mod:`model_test`. 
    929     tests = None            # type: List[TestCondition] 
    930     #: Returns the effective radius of the model given its volume parameters. 
    931     #: The presence of *ER* indicates that the model is a form factor model 
    932     #: that may be used together with a structure factor to form an implicit 
    933     #: multiplication model. 
    934     #: 
    935     #: The parameters to the *ER* function must be marked with type *volume*. 
    936     #: in the parameter table.  They will appear in the same order as they 
    937     #: do in the table.  The values passed to *ER* will be vectors, with one 
    938     #: value for each polydispersity condition.  For example, if the model 
    939     #: is polydisperse over both length and radius, then both length and 
    940     #: radius will have the same number of values in the vector, with one 
    941     #: value for each *length X radius*.  If only *radius* is polydisperse, 
    942     #: then the value for *length* will be repeated once for each value of 
    943     #: *radius*.  The *ER* function should return one effective radius for 
    944     #: each parameter set.  Multiplicity parameters will be received as 
    945     #: arrays, with one row per polydispersity condition. 
    946     ER = None               # type: Optional[Callable[[np.ndarray], np.ndarray]] 
    947     #: Returns the occupied volume and the total volume for each parameter set. 
    948     #: See :attr:`ER` for details on the parameters. 
    949     VR = None               # type: Optional[Callable[[np.ndarray], Tuple[np.ndarray, np.ndarray]]] 
    950     #: Arbitrary C code containing supporting functions, etc., to be inserted 
    951     #: after everything in source.  This can include Iq and Iqxy functions with 
    952     #: the full function signature, including all parameters. 
    953     c_code = None 
     932    #: inline source code, added after all elements of source 
     933    c_code = None           # type: Optional[str] 
    954934    #: Returns the form volume for python-based models.  Form volume is needed 
    955935    #: for volume normalization in the polydispersity integral.  If no 
     
    959939    #: C code, either defined as a string, or in the sources. 
    960940    form_volume = None      # type: Union[None, str, Callable[[np.ndarray], float]] 
     941    #: Returns the shell volume for python-based models.  Form volume and 
     942    #: shell volume are needed for volume normalization in the polydispersity 
     943    #: integral and structure interactions for hollow shapes.  If no 
     944    #: parameters are *volume* parameters, then shell volume is not needed. 
     945    #: For C-based models, (with :attr:`sources` defined, or with :attr:`Iq` 
     946    #: defined using a string containing C code), shell_volume must also be 
     947    #: C code, either defined as a string, or in the sources. 
     948    shell_volume = None      # type: Union[None, str, Callable[[np.ndarray], float]] 
    961949    #: Returns *I(q, a, b, ...)* for parameters *a*, *b*, etc. defined 
    962950    #: by the parameter table.  *Iq* can be defined as a python function, or 
     
    993981    #: Line numbers for symbols defining C code 
    994982    lineno = None           # type: Dict[str, int] 
     983    #: The set of tests that must pass.  The format of the tests is described 
     984    #: in :mod:`model_test`. 
     985    tests = None            # type: List[TestCondition] 
    995986 
    996987    def __init__(self): 
  • sasmodels/models/_spherepy.py

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    r642046e rcc8b183  
    7979} 
    8080 
    81  
    8281static double Iqabc(double qa, double qb, double qc, 
    8382    double lattice_spacing, double lattice_distortion, double radius, 
  • sasmodels/models/flexible_cylinder_elliptical.c

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

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

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

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

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

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

    r2d81cfe r304c775  
    6262 
    6363#             ["name", "units", default, [lower, upper], "type","description"], 
    64 parameters = [["radius_effective", "Ang", 50.0, [0, inf], "volume", 
     64parameters = [["radius_effective", "Ang", 50.0, [0, inf], "", 
    6565               "effective radius of hard sphere"], 
    6666              ["volfraction", "", 0.2, [0, 0.74], "", 
    6767               "volume fraction of hard spheres"], 
    6868             ] 
    69  
    70 # No volume normalization despite having a volume parameter 
    71 # This should perhaps be volume normalized? 
    72 form_volume = """ 
    73     return 1.0; 
    74     """ 
    7569 
    7670Iq = r""" 
     
    168162    return pars 
    169163 
    170 # ER defaults to 0.0 
    171 # VR defaults to 1.0 
    172  
    173164demo = dict(radius_effective=200, volfraction=0.2, 
    174165            radius_effective_pd=0.1, radius_effective_pd_n=40) 
  • sasmodels/models/hayter_msa.py

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

    r108e70e r99658f6  
    11//#define INVALID(v) (v.radius_core >= v.radius) 
    22 
    3 // From Igor library 
    43static double 
    5 _hollow_cylinder_scaling(double integrand, double delrho, double volume) 
     4shell_volume(double radius, double thickness, double length) 
    65{ 
    7     return 1.0e-4 * square(volume * delrho) * integrand; 
     6    return M_PI*length*(square(radius+thickness) - radius*radius); 
     7} 
     8 
     9static double 
     10form_volume(double radius, double thickness, double length) 
     11{ 
     12    return M_PI*length*square(radius+thickness); 
     13} 
     14 
     15static double 
     16radius_from_excluded_volume(double radius, double thickness, double length) 
     17{ 
     18    const double radius_tot = radius + thickness; 
     19    return 0.5*cbrt(0.75*radius_tot*(2.0*radius_tot*length + (radius_tot + length)*(M_PI*radius_tot + length))); 
     20} 
     21 
     22static double 
     23radius_from_volume(double radius, double thickness, double length) 
     24{ 
     25    const double volume_outer_cyl = M_PI*square(radius + thickness)*length; 
     26    return cbrt(volume_outer_cyl/M_4PI_3); 
     27} 
     28 
     29static double 
     30radius_from_diagonal(double radius, double thickness, double length) 
     31{ 
     32    return sqrt(square(radius + thickness) + 0.25*square(length)); 
     33} 
     34 
     35static double 
     36effective_radius(int mode, double radius, double thickness, double length) 
     37{ 
     38    switch (mode) { 
     39    default: 
     40    case 1: // excluded volume 
     41        return radius_from_excluded_volume(radius, thickness, length); 
     42    case 2: // equivalent volume sphere 
     43        return radius_from_volume(radius, thickness, length); 
     44    case 3: // outer radius 
     45        return radius + thickness; 
     46    case 4: // half length 
     47        return 0.5*length; 
     48    case 5: // half outer min dimension 
     49        return (radius + thickness < 0.5*length ? radius + thickness : 0.5*length); 
     50    case 6: // half outer max dimension 
     51        return (radius + thickness > 0.5*length ? radius + thickness : 0.5*length); 
     52    case 7: // half outer diagonal 
     53        return radius_from_diagonal(radius,thickness,length); 
     54    } 
    855} 
    956 
     
    2269} 
    2370 
    24 static double 
    25 form_volume(double radius, double thickness, double length) 
    26 { 
    27     double v_shell = M_PI*length*(square(radius+thickness) - radius*radius); 
    28     return v_shell; 
    29 } 
    30  
    31  
    32 static double 
    33 Iq(double q, double radius, double thickness, double length, 
     71static void 
     72Fq(double q, double *F1, double *F2, double radius, double thickness, double length, 
    3473    double sld, double solvent_sld) 
    3574{ 
     
    3776    const double upper = 1.0;        //limits of numerical integral 
    3877 
    39     double summ = 0.0;            //initialize intergral 
     78    double total_F1 = 0.0;            //initialize intergral 
     79    double total_F2 = 0.0; 
    4080    for (int i=0;i<GAUSS_N;i++) { 
    4181        const double cos_theta = 0.5*( GAUSS_Z[i] * (upper-lower) + lower + upper ); 
     
    4383        const double form = _fq(q*sin_theta, q*cos_theta, 
    4484                                radius, thickness, length); 
    45         summ += GAUSS_W[i] * form * form; 
     85        total_F1 += GAUSS_W[i] * form; 
     86        total_F2 += GAUSS_W[i] * form * form; 
    4687    } 
     88    total_F1 *= 0.5*(upper-lower); 
     89    total_F2 *= 0.5*(upper-lower); 
     90    const double s = (sld - solvent_sld) * shell_volume(radius, thickness, length); 
     91    *F1 = 1e-2 * s * total_F1; 
     92    *F2 = 1e-4 * s*s * total_F2; 
     93} 
    4794 
    48     const double Aq = 0.5*summ*(upper-lower); 
    49     const double volume = form_volume(radius, thickness, length); 
    50     return _hollow_cylinder_scaling(Aq, solvent_sld - sld, volume); 
    51 } 
    5295 
    5396static double 
     
    57100{ 
    58101    const double form = _fq(qab, qc, radius, thickness, length); 
    59  
    60     const double vol = form_volume(radius, thickness, length); 
    61     return _hollow_cylinder_scaling(form*form, solvent_sld-sld, vol); 
     102    const double s = (sld - solvent_sld) * shell_volume(radius, thickness, length); 
     103    return 1.0e-4*square(s * form); 
    62104} 
  • sasmodels/models/hollow_cylinder.py

    r455aaa1 r99658f6  
    6060.. [#] L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and 
    6161   Neutron Scattering*, Plenum Press, New York, (1987) 
     62L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
    6263 
    6364Authorship and Verification 
     
    6970* **Last Reviewed by:** Paul Butler **Date:** September 06, 2018 
    7071""" 
     72from __future__ import division 
    7173 
    7274import numpy as np 
     
    99101 
    100102source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "hollow_cylinder.c"] 
    101  
    102 # pylint: disable=W0613 
    103 def ER(radius, thickness, length): 
    104     """ 
    105     :param radius:      Cylinder core radius 
    106     :param thickness:   Cylinder wall thickness 
    107     :param length:      Cylinder length 
    108     :return:            Effective radius 
    109     """ 
    110     router = radius + thickness 
    111     if router == 0 or length == 0: 
    112         return 0.0 
    113     len1 = router 
    114     len2 = length/2.0 
    115     term1 = len1*len1*2.0*len2/2.0 
    116     term2 = 1.0 + (len2/len1)*(1.0 + 1/len2/2.0)*(1.0 + pi*len1/len2/2.0) 
    117     ddd = 3.0*term1*term2 
    118     diam = pow(ddd, (1.0/3.0)) 
    119     return diam 
    120  
    121 def VR(radius, thickness, length): 
    122     """ 
    123     :param radius:      Cylinder radius 
    124     :param thickness:   Cylinder wall thickness 
    125     :param length:      Cylinder length 
    126     :return:            Volf ratio for P(q)*S(q) 
    127     """ 
    128     router = radius + thickness 
    129     vol_core = pi*radius*radius*length 
    130     vol_total = pi*router*router*length 
    131     vol_shell = vol_total - vol_core 
    132     return vol_total, vol_shell 
     103have_Fq = True 
     104effective_radius_type = [ 
     105    "excluded volume", "equivalent outer volume sphere", "outer radius", "half length", 
     106    "half outer min dimension", "half outer max dimension", 
     107    "half outer diagonal", 
     108    ] 
    133109 
    134110def random(): 
     
    158134qx = q*cos(pi/6.0) 
    159135qy = q*sin(pi/6.0) 
     136radius = parameters[0][2] 
     137thickness = parameters[1][2] 
     138length = parameters[2][2] 
    160139# Parameters for unit tests 
    161140tests = [ 
    162141    [{}, 0.00005, 1764.926], 
    163     [{}, 'VR', 0.55555556], 
     142    [{}, 0.1, None, None, 
     143     0.5*(0.75*(radius+thickness)*(2.0*(radius+thickness)*length + ((radius+thickness) + length)*(pi*(radius+thickness) + length)))**(1./3.),  # R_eff from excluded volume 
     144     pi*((radius+thickness)**2-radius**2)*length,  # shell volume 
     145     (radius+thickness)**2/((radius+thickness)**2 - radius**2), # form:shell ratio 
     146    ], 
    164147    [{}, 0.001, 1756.76], 
    165148    [{}, (qx, qy), 2.36885476192], 
  • sasmodels/models/hollow_rectangular_prism.c

    rd86f0fc r99658f6  
     1static double 
     2shell_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
     3{ 
     4    const double length_b = length_a * b2a_ratio; 
     5    const double length_c = length_a * c2a_ratio; 
     6    const double form_volume = length_a * length_b * length_c; 
     7    const double a_core = length_a - 2.0*thickness; 
     8    const double b_core = length_b - 2.0*thickness; 
     9    const double c_core = length_c - 2.0*thickness; 
     10    const double core_volume = a_core * b_core * c_core; 
     11    return form_volume - core_volume; 
     12} 
     13 
    114static double 
    215form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
    316{ 
    4     double length_b = length_a * b2a_ratio; 
    5     double length_c = length_a * c2a_ratio; 
    6     double a_core = length_a - 2.0*thickness; 
    7     double b_core = length_b - 2.0*thickness; 
    8     double c_core = length_c - 2.0*thickness; 
    9     double vol_core = a_core * b_core * c_core; 
    10     double vol_total = length_a * length_b * length_c; 
    11     double vol_shell = vol_total - vol_core; 
    12     return vol_shell; 
     17    const double length_b = length_a * b2a_ratio; 
     18    const double length_c = length_a * c2a_ratio; 
     19    const double form_volume = length_a * length_b * length_c; 
     20    return form_volume; 
    1321} 
    1422 
    1523static double 
    16 Iq(double q, 
     24radius_from_excluded_volume(double length_a, double b2a_ratio, double c2a_ratio) 
     25{ 
     26    const double r_equiv = sqrt(length_a*length_a*b2a_ratio/M_PI); 
     27    const double length_c = length_a*c2a_ratio; 
     28    return 0.5*cbrt(0.75*r_equiv*(2.0*r_equiv*length_c + (r_equiv + length_c)*(M_PI*r_equiv + length_c))); 
     29} 
     30 
     31static double 
     32effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
     33// NOTE length_a is external dimension! 
     34{ 
     35    switch (mode) { 
     36    default: 
     37    case 1: // equivalent cylinder excluded volume 
     38        return radius_from_excluded_volume(length_a, b2a_ratio, c2a_ratio); 
     39    case 2: // equivalent outer volume sphere 
     40        return cbrt(cube(length_a)*b2a_ratio*c2a_ratio/M_4PI_3); 
     41    case 3: // half length_a 
     42        return 0.5 * length_a; 
     43    case 4: // half length_b 
     44        return 0.5 * length_a*b2a_ratio; 
     45    case 5: // half length_c 
     46        return 0.5 * length_a*c2a_ratio; 
     47    case 6: // equivalent outer circular cross-section 
     48        return length_a*sqrt(b2a_ratio/M_PI); 
     49    case 7: // half ab diagonal 
     50        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio))); 
     51    case 8: // half diagonal 
     52        return 0.5*sqrt(square(length_a) * (1.0 + square(b2a_ratio) + square(c2a_ratio))); 
     53    } 
     54} 
     55 
     56static void 
     57Fq(double q, 
     58    double *F1, 
     59    double *F2, 
    1760    double sld, 
    1861    double solvent_sld, 
     
    3679    const double v2b = M_PI_2;  //phi integration limits 
    3780 
    38     double outer_sum = 0.0; 
     81    double outer_sum_F1 = 0.0; 
     82    double outer_sum_F2 = 0.0; 
    3983    for(int i=0; i<GAUSS_N; i++) { 
    4084 
     
    4690        const double termC2 = sas_sinx_x(q * (c_half-thickness)*cos(theta)); 
    4791 
    48         double inner_sum = 0.0; 
     92        double inner_sum_F1 = 0.0; 
     93        double inner_sum_F2 = 0.0; 
    4994        for(int j=0; j<GAUSS_N; j++) { 
    5095 
     
    64109            const double AP2 = vol_core * termA2 * termB2 * termC2; 
    65110 
    66             inner_sum += GAUSS_W[j] * square(AP1-AP2); 
     111            inner_sum_F1 += GAUSS_W[j] * (AP1-AP2); 
     112            inner_sum_F2 += GAUSS_W[j] * square(AP1-AP2); 
    67113        } 
    68         inner_sum *= 0.5 * (v2b-v2a); 
     114        inner_sum_F1 *= 0.5 * (v2b-v2a); 
     115        inner_sum_F2 *= 0.5 * (v2b-v2a); 
    69116 
    70         outer_sum += GAUSS_W[i] * inner_sum * sin(theta); 
     117        outer_sum_F1 += GAUSS_W[i] * inner_sum_F1 * sin(theta); 
     118        outer_sum_F2 += GAUSS_W[i] * inner_sum_F2 * sin(theta); 
    71119    } 
    72     outer_sum *= 0.5*(v1b-v1a); 
     120    outer_sum_F1 *= 0.5*(v1b-v1a); 
     121    outer_sum_F2 *= 0.5*(v1b-v1a); 
    73122 
    74123    // Normalize as in Eqn. (15) without the volume factor (as cancels with (V*DelRho)^2 normalization) 
    75124    // The factor 2 is due to the different theta integration limit (pi/2 instead of pi) 
    76     const double form = outer_sum/M_PI_2; 
     125    const double form_avg = outer_sum_F1/M_PI_2; 
     126    const double form_squared_avg = outer_sum_F2/M_PI_2; 
    77127 
    78128    // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization. 
    79     const double delrho = sld - solvent_sld; 
     129    const double contrast = sld - solvent_sld; 
    80130 
    81131    // Convert from [1e-12 A-1] to [cm-1] 
    82     return 1.0e-4 * delrho * delrho * form; 
     132    *F1 = 1.0e-2 * contrast * form_avg; 
     133    *F2 = 1.0e-4 * contrast * contrast * form_squared_avg; 
    83134} 
    84135 
  • sasmodels/models/hollow_rectangular_prism.py

    r455aaa1 r99658f6  
    9898 
    9999.. [#Nayuk2012] R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 
     100L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949).  
    100101 
    101102 
     
    132133               "Solvent scattering length density"], 
    133134              ["length_a", "Ang", 35, [0, inf], "volume", 
    134                "Shorter side of the parallelepiped"], 
     135               "Shortest, external, size of the parallelepiped"], 
    135136              ["b2a_ratio", "Ang", 1, [0, inf], "volume", 
    136137               "Ratio sides b/a"], 
     
    148149 
    149150source = ["lib/gauss76.c", "hollow_rectangular_prism.c"] 
    150  
    151 def ER(length_a, b2a_ratio, c2a_ratio, thickness): 
    152     """ 
    153     Return equivalent radius (ER) 
    154     thickness parameter not used 
    155     """ 
    156     b_side = length_a * b2a_ratio 
    157     c_side = length_a * c2a_ratio 
    158  
    159     # surface average radius (rough approximation) 
    160     surf_rad = sqrt(length_a * b_side / pi) 
    161  
    162     ddd = 0.75 * surf_rad * (2 * surf_rad * c_side + (c_side + surf_rad) * (c_side + pi * surf_rad)) 
    163     return 0.5 * (ddd) ** (1. / 3.) 
    164  
    165 def VR(length_a, b2a_ratio, c2a_ratio, thickness): 
    166     """ 
    167     Return shell volume and total volume 
    168     """ 
    169     b_side = length_a * b2a_ratio 
    170     c_side = length_a * c2a_ratio 
    171     a_core = length_a - 2.0*thickness 
    172     b_core = b_side - 2.0*thickness 
    173     c_core = c_side - 2.0*thickness 
    174     vol_core = a_core * b_core * c_core 
    175     vol_total = length_a * b_side * c_side 
    176     vol_shell = vol_total - vol_core 
    177     return vol_total, vol_shell 
    178  
     151have_Fq = True 
     152effective_radius_type = [ 
     153    "equivalent cylinder excluded volume", "equivalent outer volume sphere",  
     154    "half length_a", "half length_b", "half length_c", 
     155    "equivalent outer circular cross-section", 
     156    "half ab diagonal", "half diagonal", 
     157    ] 
    179158 
    180159def random(): 
  • sasmodels/models/hollow_rectangular_prism_thin_walls.c

    rd86f0fc r99658f6  
     1static double 
     2shell_volume(double length_a, double b2a_ratio, double c2a_ratio) 
     3{ 
     4    const double length_b = length_a * b2a_ratio; 
     5    const double length_c = length_a * c2a_ratio; 
     6    const double shell_volume = 2.0 * (length_a*length_b + length_a*length_c + length_b*length_c); 
     7    return shell_volume; 
     8} 
     9 
    110static double 
    211form_volume(double length_a, double b2a_ratio, double c2a_ratio) 
    312{ 
    4     double length_b = length_a * b2a_ratio; 
    5     double length_c = length_a * c2a_ratio; 
    6     double vol_shell = 2.0 * (length_a*length_b + length_a*length_c + length_b*length_c); 
    7     return vol_shell; 
     13    const double length_b = length_a * b2a_ratio; 
     14    const double length_c = length_a * c2a_ratio; 
     15    const double form_volume = length_a * length_b * length_c;