source: sasmodels/sasmodels/kernel.py @ 36a2418

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 36a2418 was 36a2418, checked in by Paul Kienzle <pkienzle@…>, 5 years ago

Merge branch 'beta_approx' into ticket-1015-gpu-mem-error

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