source: sasmodels/sasmodels/kernel.py @ a34b811

ticket-1257-vesicle-productticket_1156ticket_822_more_unit_tests
Last change on this file since a34b811 was a34b811, checked in by Paul Kienzle <pkienzle@…>, 3 months ago

use radius_effective/radius_effective_mode/radius_effective_modes consistently throughout the code

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