source: sasmodels/sasmodels/mixture.py @ a557a99

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since a557a99 was 6dc78e4, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

re-enable structure factor product

  • Property mode set to 100644
File size: 8.2 KB
Line 
1"""
2Mixture model
3-------------
4
5The product model multiplies the structure factor by the form factor,
6modulated by the effective radius of the form.  The resulting model
7has a attributes of both the model description (with parameters, etc.)
8and the module evaluator (with call, release, etc.).
9
10To use it, first load form factor P and structure factor S, then create
11*ProductModel(P, S)*.
12"""
13from __future__ import print_function
14
15from copy import copy
16import numpy as np  # type: ignore
17
18from .modelinfo import Parameter, ParameterTable, ModelInfo
19from .kernel import KernelModel, Kernel
20from .details import make_details
21
22try:
23    from typing import List
24except ImportError:
25    pass
26
27def make_mixture_info(parts):
28    # type: (List[ModelInfo]) -> ModelInfo
29    """
30    Create info block for mixture model.
31    """
32    flatten = []
33    for part in parts:
34        if part.composition and part.composition[0] == 'mixture':
35            flatten.extend(part.composition[1])
36        else:
37            flatten.append(part)
38    parts = flatten
39
40    # Build new parameter list
41    combined_pars = []
42    demo = {}
43    for k, part in enumerate(parts):
44        # Parameter prefix per model, A_, B_, ...
45        # Note that prefix must also be applied to id and length_control
46        # to support vector parameters
47        prefix = chr(ord('A')+k) + '_'
48        scale =  Parameter(prefix+'scale', default=1.0,
49                           description="model intensity for " + part.name)
50        combined_pars.append(scale)
51        for p in part.parameters.kernel_parameters:
52            p = copy(p)
53            p.name = prefix + p.name
54            p.id = prefix + p.id
55            if p.length_control is not None:
56                p.length_control = prefix + p.length_control
57            combined_pars.append(p)
58        demo.update((prefix+k, v) for k, v in part.demo.items()
59                    if k != "background")
60    #print("pars",combined_pars)
61    parameters = ParameterTable(combined_pars)
62    parameters.max_pd = sum(part.parameters.max_pd for part in parts)
63
64    model_info = ModelInfo()
65    model_info.id = '+'.join(part.id for part in parts)
66    model_info.name = ' + '.join(part.name for part in parts)
67    model_info.filename = None
68    model_info.title = 'Mixture model with ' + model_info.name
69    model_info.description = model_info.title
70    model_info.docs = model_info.title
71    model_info.category = "custom"
72    model_info.parameters = parameters
73    #model_info.single = any(part['single'] for part in parts)
74    model_info.structure_factor = False
75    model_info.variant_info = None
76    #model_info.tests = []
77    #model_info.source = []
78    # Iq, Iqxy, form_volume, ER, VR and sesans
79    # Remember the component info blocks so we can build the model
80    model_info.composition = ('mixture', parts)
81    model_info.demo = demo
82    return model_info
83
84
85class MixtureModel(KernelModel):
86    def __init__(self, model_info, parts):
87        # type: (ModelInfo, List[KernelModel]) -> None
88        self.info = model_info
89        self.parts = parts
90
91    def make_kernel(self, q_vectors):
92        # type: (List[np.ndarray]) -> MixtureKernel
93        # Note: may be sending the q_vectors to the n times even though they
94        # are only needed once.  It would mess up modularity quite a bit to
95        # handle this optimally, especially since there are many cases where
96        # separate q vectors are needed (e.g., form in python and structure
97        # in opencl; or both in opencl, but one in single precision and the
98        # other in double precision).
99        kernels = [part.make_kernel(q_vectors) for part in self.parts]
100        return MixtureKernel(self.info, kernels)
101
102    def release(self):
103        # type: () -> None
104        """
105        Free resources associated with the model.
106        """
107        for part in self.parts:
108            part.release()
109
110
111class MixtureKernel(Kernel):
112    def __init__(self, model_info, kernels):
113        # type: (ModelInfo, List[Kernel]) -> None
114        self.dim = kernels[0].dim
115        self.info =  model_info
116        self.kernels = kernels
117        self.dtype = self.kernels[0].dtype
118        self.results = []  # type: List[np.ndarray]
119
120    def __call__(self, call_details, values, cutoff, magnetic):
121        # type: (CallDetails, np.ndarray, np.ndarry, float, bool) -> np.ndarray
122        scale, background = values[0:2]
123        total = 0.0
124        # remember the parts for plotting later
125        self.results = []  # type: List[np.ndarray]
126        offset = 2 # skip scale & background
127        parts = MixtureParts(self.info, self.kernels, call_details, values)
128        for kernel, kernel_details, kernel_values in parts:
129            #print("calling kernel", kernel.info.name)
130            result = kernel(kernel_details, kernel_values, cutoff, magnetic)
131            #print(kernel.info.name, result)
132            total += result
133            self.results.append(result)
134
135        return scale*total + background
136
137    def release(self):
138        # type: () -> None
139        for k in self.kernels:
140            k.release()
141
142
143class MixtureParts(object):
144    def __init__(self, model_info, kernels, call_details, values):
145        # type: (ModelInfo, List[Kernel], CallDetails, np.ndarray) -> None
146        self.model_info = model_info
147        self.parts = model_info.composition[1]
148        self.kernels = kernels
149        self.call_details = call_details
150        self.values = values
151        self.spin_index = model_info.parameters.npars + 2
152        #call_details.show(values)
153
154    def __iter__(self):
155        # type: () -> PartIterable
156        self.part_num = 0
157        self.par_index = 2
158        self.mag_index = self.spin_index + 3
159        return self
160
161    def next(self):
162        # type: () -> Tuple[List[Callable], CallDetails, np.ndarray]
163        if self.part_num >= len(self.parts):
164            raise StopIteration()
165        info = self.parts[self.part_num]
166        kernel = self.kernels[self.part_num]
167        call_details = self._part_details(info, self.par_index)
168        values = self._part_values(info, self.par_index, self.mag_index)
169        values = values.astype(kernel.dtype)
170        #call_details.show(values)
171
172        self.part_num += 1
173        self.par_index += info.parameters.npars + 1
174        self.mag_index += 3 * len(info.parameters.magnetism_index)
175
176        return kernel, call_details, values
177
178    def _part_details(self, info, par_index):
179        # type: (ModelInfo, int) -> CallDetails
180        full = self.call_details
181        # par_index is index into values array of the current parameter,
182        # which includes the initial scale and background parameters.
183        # We want the index into the weight length/offset for each parameter.
184        # Exclude the initial scale and background, so subtract two, but each
185        # component has its own scale factor which we need to skip when
186        # constructing the details for the kernel, so add one, giving a
187        # net subtract one.
188        index = slice(par_index - 1, par_index - 1 + info.parameters.npars)
189        length = full.length[index]
190        offset = full.offset[index]
191        # The complete weight vector is being sent to each part so that
192        # offsets don't need to be adjusted.
193        part = make_details(info, length, offset, full.num_weights)
194        return part
195
196    def _part_values(self, info, par_index, mag_index):
197        # type: (ModelInfo, int, int) -> np.ndarray
198        #print(info.name, par_index, self.values[par_index:par_index + info.parameters.npars + 1])
199        scale = self.values[par_index]
200        pars = self.values[par_index + 1:par_index + info.parameters.npars + 1]
201        nmagnetic = len(info.parameters.magnetism_index)
202        if nmagnetic:
203            spin_state = self.values[self.spin_index:self.spin_index + 3]
204            mag_index = self.values[mag_index:mag_index + 3 * nmagnetic]
205        else:
206            spin_state = []
207            mag_index = []
208        nvalues = self.model_info.parameters.nvalues
209        nweights = self.call_details.num_weights
210        weights = self.values[nvalues:nvalues+2*nweights]
211        zero = self.values.dtype.type(0.)
212        values = [[scale, zero], pars, spin_state, mag_index, weights]
213        # Pad value array to a 32 value boundary
214        spacer = (32 - sum(len(v) for v in values)%32)%32
215        values.append([zero]*spacer)
216        values = np.hstack(values).astype(self.kernels[0].dtype)
217        return values
Note: See TracBrowser for help on using the repository browser.