Changeset 9150036 in sasmodels


Ignore:
Timestamp:
Mar 6, 2019 6:05:28 PM (2 months ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
f64b154, bd91e8f, c11d09f, 02226a2
Parents:
e589e9a (diff), 37f38ff (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' into beta_approx

Files:
50 added
116 edited

Legend:

Unmodified
Added
Removed
  • doc/guide/plugin.rst

    raa8c6e0 r9150036  
    272272structure factor to account for interactions between particles.  See 
    273273`Form_Factors`_ for more details. 
     274 
     275**model_info = ...** lets you define a model directly, for example, by 
     276loading and modifying existing models.  This is done implicitly by 
     277:func:`sasmodels.core.load_model_info`, which can create a mixture model 
     278from a pair of existing models.  For example:: 
     279 
     280    from sasmodels.core import load_model_info 
     281    model_info = load_model_info('sphere+cylinder') 
     282 
     283See :class:`sasmodels.modelinfo.ModelInfo` for details about the model 
     284attributes that are defined. 
    274285 
    275286Model Parameters 
     
    894905             - \frac{\sin(x)}{x}\left(\frac{1}{x} - \frac{3!}{x^3} + \frac{5!}{x^5} - \frac{7!}{x^7}\right) 
    895906 
    896         For small arguments , 
     907        For small arguments, 
    897908 
    898909        .. math:: 
  • example/multiscatfit.py

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

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

    r49d1f8b8 r2c4a190  
    3535    # when bumps is not on the path. 
    3636    from bumps.names import Parameter # type: ignore 
     37    from bumps.parameter import Reference # type: ignore 
    3738except ImportError: 
    3839    pass 
     
    139140    def __init__(self, data, model, cutoff=1e-5, name=None, extra_pars=None): 
    140141        # type: (Data, Model, float) -> None 
     142        # Allow resolution function to define fittable parameters.  We do this 
     143        # by creating reference parameters within the resolution object rather 
     144        # than modifying the object itself to use bumps parameters.  We need 
     145        # to reset the parameters each time the object has changed.  These 
     146        # additional parameters need to be returned from the fitting engine. 
     147        # To make them available to the user, they are added as top-level 
     148        # attributes to the experiment object.  The only change to the 
     149        # resolution function is that it needs an optional 'fittable' attribute 
     150        # which maps the internal name to the user visible name for the 
     151        # for the parameter. 
     152        self._resolution = None 
     153        self._resolution_pars = {} 
    141154        # remember inputs so we can inspect from outside 
    142155        self.name = data.filename if name is None else name 
     
    145158        self._interpret_data(data, model.sasmodel) 
    146159        self._cache = {} 
     160        # CRUFT: no longer need extra parameters 
     161        # Multiple scattering probability is now retrieved directly from the 
     162        # multiple scattering resolution function. 
    147163        self.extra_pars = extra_pars 
    148164 
     
    162178        return len(self.Iq) 
    163179 
     180    @property 
     181    def resolution(self): 
     182        return self._resolution 
     183 
     184    @resolution.setter 
     185    def resolution(self, value): 
     186        self._resolution = value 
     187 
     188        # Remove old resolution fitting parameters from experiment 
     189        for name in self._resolution_pars: 
     190            delattr(self, name) 
     191 
     192        # Create new resolution fitting parameters 
     193        res_pars = getattr(self._resolution, 'fittable', {}) 
     194        self._resolution_pars = { 
     195            name: Reference(self._resolution, refname, name=name) 
     196            for refname, name in res_pars.items() 
     197        } 
     198 
     199        # Add new resolution fitting parameters as experiment attributes 
     200        for name, ref in self._resolution_pars.items(): 
     201            setattr(self, name, ref) 
     202 
    164203    def parameters(self): 
    165204        # type: () -> Dict[str, Parameter] 
     
    168207        """ 
    169208        pars = self.model.parameters() 
    170         if self.extra_pars: 
     209        if self.extra_pars is not None: 
    171210            pars.update(self.extra_pars) 
     211        pars.update(self._resolution_pars) 
    172212        return pars 
    173213 
  • sasmodels/direct_model.py

    rc1799d3 r9150036  
    224224            else: 
    225225                Iq, dIq = None, None 
    226             #self._theory = np.zeros_like(q) 
    227             q_vectors = [res.q_calc] 
    228226        elif self.data_type == 'Iqxy': 
    229227            #if not model.info.parameters.has_2d: 
     
    242240            res = resolution2d.Pinhole2D(data=data, index=index, 
    243241                                         nsigma=3.0, accuracy=accuracy) 
    244             #self._theory = np.zeros_like(self.Iq) 
    245             q_vectors = res.q_calc 
    246242        elif self.data_type == 'Iq': 
    247243            index = (data.x >= data.qmin) & (data.x <= data.qmax) 
     
    268264            else: 
    269265                res = resolution.Perfect1D(data.x[index]) 
    270  
    271             #self._theory = np.zeros_like(self.Iq) 
    272             q_vectors = [res.q_calc] 
    273266        elif self.data_type == 'Iq-oriented': 
    274267            index = (data.x >= data.qmin) & (data.x <= data.qmax) 
     
    286279                                      qx_width=data.dxw[index], 
    287280                                      qy_width=data.dxl[index]) 
    288             q_vectors = res.q_calc 
    289281        else: 
    290282            raise ValueError("Unknown data type") # never gets here 
     
    292284        # Remember function inputs so we can delay loading the function and 
    293285        # so we can save/restore state 
    294         self._kernel_inputs = q_vectors 
    295286        self._kernel = None 
    296287        self.Iq, self.dIq, self.index = Iq, dIq, index 
     
    329320        # type: (ParameterSet, float) -> np.ndarray 
    330321        if self._kernel is None: 
    331             self._kernel = self._model.make_kernel(self._kernel_inputs) 
     322            # TODO: change interfaces so that resolution returns kernel inputs 
     323            # Maybe have resolution always return a tuple, or maybe have 
     324            # make_kernel accept either an ndarray or a pair of ndarrays. 
     325            kernel_inputs = self.resolution.q_calc 
     326            if isinstance(kernel_inputs, np.ndarray): 
     327                kernel_inputs = (kernel_inputs,) 
     328            self._kernel = self._model.make_kernel(kernel_inputs) 
    332329 
    333330        # Need to pull background out of resolution for multiple scattering 
  • sasmodels/multiscat.py

    rb3703f5 r2c4a190  
    342342 
    343343    *probability* is related to the expected number of scattering 
    344     events in the sample $\lambda$ as $p = 1 = e^{-\lambda}$.  As a 
    345     hack to allow probability to be a fitted parameter, the "value" 
    346     can be a function that takes no parameters and returns the current 
    347     value of the probability.  *coverage* determines how many scattering 
    348     steps to consider.  The default is 0.99, which sets $n$ such that 
    349     $1 \ldots n$ covers 99% of the Poisson probability mass function. 
     344    events in the sample $\lambda$ as $p = 1 - e^{-\lambda}$. 
     345    *coverage* determines how many scattering steps to consider.  The 
     346    default is 0.99, which sets $n$ such that $1 \ldots n$ covers 99% 
     347    of the Poisson probability mass function. 
    350348 
    351349    *is2d* is True then 2D scattering is used, otherwise it accepts 
     
    399397        self.qmin = qmin 
    400398        self.nq = nq 
    401         self.probability = probability 
     399        self.probability = 0. if probability is None else probability 
    402400        self.coverage = coverage 
    403401        self.is2d = is2d 
     
    456454        self.Iqxy = None # type: np.ndarray 
    457455 
     456        # Label probability as a fittable parameter, and give its external name 
     457        # Note that the external name must be a valid python identifier, since 
     458        # is will be set as an experiment attribute. 
     459        self.fittable = {'probability': 'scattering_probability'} 
     460 
    458461    def apply(self, theory): 
    459462        if self.is2d: 
     
    463466        Iq_calc = Iq_calc.reshape(self.nq, self.nq) 
    464467 
     468        # CRUFT: don't need probability as a function anymore 
    465469        probability = self.probability() if callable(self.probability) else self.probability 
    466470        coverage = self.coverage 
  • sasmodels/sasview_model.py

    ra8a1d48 r9150036  
    2525from . import core 
    2626from . import custom 
     27from . import kernelcl 
    2728from . import product 
    2829from . import generate 
     
    3031from . import modelinfo 
    3132from .details import make_kernel_args, dispersion_mesh 
     33from .kernelcl import reset_environment 
    3234 
    3335# pylint: disable=unused-import 
     
    6870#: has changed since we last reloaded. 
    6971_CACHED_MODULE = {}  # type: Dict[str, "module"] 
     72 
     73def reset_environment(): 
     74    # type: () -> None 
     75    """ 
     76    Clear the compute engine context so that the GUI can change devices. 
     77 
     78    This removes all compiled kernels, even those that are active on fit 
     79    pages, but they will be restored the next time they are needed. 
     80    """ 
     81    kernelcl.reset_environment() 
     82    for model in MODELS.values(): 
     83        model._model = None 
    7084 
    7185def find_model(modelname): 
     
    696710    def _calculate_Iq(self, qx, qy=None): 
    697711        if self._model is None: 
    698             self._model = core.build_model(self._model_info) 
     712            # Only need one copy of the compiled kernel regardless of how many 
     713            # times it is used, so store it in the class.  Also, to reset the 
     714            # compute engine, need to clear out all existing compiled kernels, 
     715            # which is much easier to do if we store them in the class. 
     716            self.__class__._model = core.build_model(self._model_info) 
    699717        if qy is not None: 
    700718            q_vectors = [np.asarray(qx), np.asarray(qy)] 
  • README.rst

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    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 rc1799d3  
    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 
     
    405404      parameters counted as n individual parameters p1, p2, ... 
    406405 
     406    * *common_parameters* is the list of common parameters, with a unique 
     407      copy for each model so that structure factors can have a default 
     408      background of 0.0. 
     409 
    407410    * *call_parameters* is the complete list of parameters to the kernel, 
    408411      including scale and background, with vector parameters recorded as 
     
    417420    parameters don't use vector notation, and instead use p1, p2, ... 
    418421    """ 
    419     # scale and background are implicit parameters 
    420     COMMON = [Parameter(*p) for p in COMMON_PARAMETERS] 
    421  
    422422    def __init__(self, parameters): 
    423423        # type: (List[Parameter]) -> None 
     424 
     425        # scale and background are implicit parameters 
     426        # Need them to be unique to each model in case they have different 
     427        # properties, such as default=0.0 for structure factor backgrounds. 
     428        self.common_parameters = [Parameter(*p) for p in COMMON_PARAMETERS] 
     429 
    424430        self.kernel_parameters = parameters 
    425431        self._set_vector_lengths() 
    426  
    427432        self.npars = sum(p.length for p in self.kernel_parameters) 
    428433        self.nmagnetic = sum(p.length for p in self.kernel_parameters 
     
    431436        if self.nmagnetic: 
    432437            self.nvalues += 3 + 3*self.nmagnetic 
    433  
    434438        self.call_parameters = self._get_call_parameters() 
    435439        self.defaults = self._get_defaults() 
     
    471475                         if p.polydisperse and p.type not in ('orientation', 'magnetic')) 
    472476        self.pd_2d = set(p.name for p in self.call_parameters if p.polydisperse) 
     477 
     478    def set_zero_background(self): 
     479        """ 
     480        Set the default background to zero for this model.  This is done for 
     481        structure factor models. 
     482        """ 
     483        # type: () -> None 
     484        # Make sure background is the second common parameter. 
     485        assert self.common_parameters[1].id == "background" 
     486        self.common_parameters[1].default = 0.0 
     487        self.defaults = self._get_defaults() 
    473488 
    474489    def check_angles(self): 
     
    570585    def _get_call_parameters(self): 
    571586        # type: () -> List[Parameter] 
    572         full_list = self.COMMON[:] 
     587        full_list = self.common_parameters[:] 
    573588        for p in self.kernel_parameters: 
    574589            if p.length == 1: 
     
    673688 
    674689        # Gather the user parameters in order 
    675         result = control + self.COMMON 
     690        result = control + self.common_parameters 
    676691        for p in self.kernel_parameters: 
    677692            if not is2d and p.type in ('orientation', 'magnetic'): 
     
    722737 
    723738#: Set of variables defined in the model that might contain C code 
    724 C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'c_code'] 
     739C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'shell_volume', 'c_code'] 
    725740 
    726741def _find_source_lines(model_info, kernel_module): 
     
    771786        # Custom sum/multi models 
    772787        return kernel_module.model_info 
     788 
    773789    info = ModelInfo() 
     790 
     791    # Build the parameter table 
    774792    #print("make parameter table", kernel_module.parameters) 
    775793    parameters = make_parameter_table(getattr(kernel_module, 'parameters', [])) 
     794 
     795    # background defaults to zero for structure factor models 
     796    structure_factor = getattr(kernel_module, 'structure_factor', False) 
     797    if structure_factor: 
     798        parameters.set_zero_background() 
     799 
     800    # TODO: remove demo parameters 
     801    # The plots in the docs are generated from the model default values. 
     802    # Sascomp set parameters from the command line, and so doesn't need 
     803    # demo values for testing. 
    776804    demo = expand_pars(parameters, getattr(kernel_module, 'demo', None)) 
     805 
    777806    filename = abspath(kernel_module.__file__).replace('.pyc', '.py') 
    778807    kernel_id = splitext(basename(filename))[0] 
     
    792821    info.category = getattr(kernel_module, 'category', None) 
    793822    info.structure_factor = getattr(kernel_module, 'structure_factor', False) 
     823    # TODO: find Fq by inspection 
     824    info.effective_radius_type = getattr(kernel_module, 'effective_radius_type', None) 
     825    info.have_Fq = getattr(kernel_module, 'have_Fq', False) 
    794826    info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y']) 
    795827    # Note: custom.load_custom_kernel_module assumes the C sources are defined 
     
    797829    info.source = getattr(kernel_module, 'source', []) 
    798830    info.c_code = getattr(kernel_module, 'c_code', None) 
     831    info.effective_radius = getattr(kernel_module, 'effective_radius', None) 
    799832    # TODO: check the structure of the tests 
    800833    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 
    803834    info.form_volume = getattr(kernel_module, 'form_volume', None) # type: ignore 
     835    info.shell_volume = getattr(kernel_module, 'shell_volume', None) # type: ignore 
    804836    info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore 
    805837    info.Iqxy = getattr(kernel_module, 'Iqxy', None) # type: ignore 
     
    825857    info.lineno = {} 
    826858    _find_source_lines(info, kernel_module) 
    827  
    828859    return info 
    829860 
     
    918949    #: provided in the file. 
    919950    structure_factor = None # type: bool 
     951    #: True if the model defines an Fq function with signature 
     952    #: void Fq(double q, double *F1, double *F2, ...) 
     953    have_Fq = False 
     954    #: List of options for computing the effective radius of the shape, 
     955    #: or None if the model is not usable as a form factor model. 
     956    effective_radius_type = None   # type: List[str] 
    920957    #: List of C source files used to define the model.  The source files 
    921958    #: should define the *Iq* function, and possibly *Iqac* or *Iqabc* if the 
    922959    #: model defines orientation parameters. Files containing the most basic 
    923960    #: 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. 
     961    #: use those functions. 
    926962    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 
     963    #: inline source code, added after all elements of source 
     964    c_code = None           # type: Optional[str] 
    954965    #: Returns the form volume for python-based models.  Form volume is needed 
    955966    #: for volume normalization in the polydispersity integral.  If no 
     
    959970    #: C code, either defined as a string, or in the sources. 
    960971    form_volume = None      # type: Union[None, str, Callable[[np.ndarray], float]] 
     972    #: Returns the shell volume for python-based models.  Form volume and 
     973    #: shell volume are needed for volume normalization in the polydispersity 
     974    #: integral and structure interactions for hollow shapes.  If no 
     975    #: parameters are *volume* parameters, then shell volume is not needed. 
     976    #: For C-based models, (with :attr:`sources` defined, or with :attr:`Iq` 
     977    #: defined using a string containing C code), shell_volume must also be 
     978    #: C code, either defined as a string, or in the sources. 
     979    shell_volume = None      # type: Union[None, str, Callable[[np.ndarray], float]] 
    961980    #: Returns *I(q, a, b, ...)* for parameters *a*, *b*, etc. defined 
    962981    #: by the parameter table.  *Iq* can be defined as a python function, or 
     
    9931012    #: Line numbers for symbols defining C code 
    9941013    lineno = None           # type: Dict[str, int] 
     1014    #: The set of tests that must pass.  The format of the tests is described 
     1015    #: in :mod:`model_test`. 
     1016    tests = None            # type: List[TestCondition] 
    9951017 
    9961018    def __init__(self): 
  • sasmodels/models/_spherepy.py

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    <
    r108e70e r99658f6  
    11//#define INVALID(v) (v.radius_core >= v.radius) 
    22 
    3 // From Igor library 
    43static double 
    5 _hollow_cylinder_scaling(double integrand, double delrho, double volume) 
     4shell_volume(double radius, double thickness, double length) 
    65{ 
    7     return 1.0e-4 * square(volume * delrho) * integrand; 
     6    return M_PI*length*(square(radius+thickness) - radius*radius); 
     7} 
     8 
     9static double 
     10form_volume(double radius, double thickness, double length) 
     11{ 
     12    return M_PI*length*square(radius+thickness); 
     13} 
     14 
     15static double 
     16radius_from_excluded_volume(double radius, double thickness, double length) 
     17{