1 | """ |
---|
2 | Mixture model |
---|
3 | ------------- |
---|
4 | |
---|
5 | The product model multiplies the structure factor by the form factor, |
---|
6 | modulated by the effective radius of the form. The resulting model |
---|
7 | has a attributes of both the model description (with parameters, etc.) |
---|
8 | and the module evaluator (with call, release, etc.). |
---|
9 | |
---|
10 | To use it, first load form factor P and structure factor S, then create |
---|
11 | *ProductModel(P, S)*. |
---|
12 | """ |
---|
13 | from copy import copy |
---|
14 | import numpy as np |
---|
15 | |
---|
16 | from .generate import process_parameters, COMMON_PARAMETERS, Parameter |
---|
17 | |
---|
18 | SCALE=0 |
---|
19 | BACKGROUND=1 |
---|
20 | EFFECT_RADIUS=2 |
---|
21 | VOLFRACTION=3 |
---|
22 | |
---|
23 | def make_mixture_info(parts): |
---|
24 | """ |
---|
25 | Create info block for product model. |
---|
26 | """ |
---|
27 | flatten = [] |
---|
28 | for part in parts: |
---|
29 | if part['composition'] and part['composition'][0] == 'mixture': |
---|
30 | flatten.extend(part['compostion'][1]) |
---|
31 | else: |
---|
32 | flatten.append(part) |
---|
33 | parts = flatten |
---|
34 | |
---|
35 | # Build new parameter list |
---|
36 | pars = COMMON_PARAMETERS + [] |
---|
37 | for k, part in enumerate(parts): |
---|
38 | # Parameter prefix per model, A_, B_, ... |
---|
39 | prefix = chr(ord('A')+k) + '_' |
---|
40 | for p in part['parameters']: |
---|
41 | # No background on the individual mixture elements |
---|
42 | if p.name == 'background': |
---|
43 | continue |
---|
44 | # TODO: promote Parameter to a full class |
---|
45 | # this code knows too much about the implementation! |
---|
46 | p = list(p) |
---|
47 | if p[0] == 'scale': # allow model subtraction |
---|
48 | p[3] = [-np.inf, np.inf] # limits |
---|
49 | p[0] = prefix+p[0] # relabel parameters with A_, ... |
---|
50 | pars.append(p) |
---|
51 | |
---|
52 | model_info = {} |
---|
53 | model_info['id'] = '+'.join(part['id']) |
---|
54 | model_info['name'] = ' + '.join(part['name']) |
---|
55 | model_info['filename'] = None |
---|
56 | model_info['title'] = 'Mixture model with ' + model_info['name'] |
---|
57 | model_info['description'] = model_info['title'] |
---|
58 | model_info['docs'] = model_info['title'] |
---|
59 | model_info['category'] = "custom" |
---|
60 | model_info['parameters'] = pars |
---|
61 | #model_info['single'] = any(part['single'] for part in parts) |
---|
62 | model_info['structure_factor'] = False |
---|
63 | model_info['variant_info'] = None |
---|
64 | #model_info['tests'] = [] |
---|
65 | #model_info['source'] = [] |
---|
66 | # Iq, Iqxy, form_volume, ER, VR and sesans |
---|
67 | # Remember the component info blocks so we can build the model |
---|
68 | model_info['composition'] = ('mixture', parts) |
---|
69 | process_parameters(model_info) |
---|
70 | return model_info |
---|
71 | |
---|
72 | |
---|
73 | class MixtureModel(object): |
---|
74 | def __init__(self, model_info, parts): |
---|
75 | self.info = model_info |
---|
76 | self.parts = parts |
---|
77 | |
---|
78 | def __call__(self, q_vectors): |
---|
79 | # Note: may be sending the q_vectors to the n times even though they |
---|
80 | # are only needed once. It would mess up modularity quite a bit to |
---|
81 | # handle this optimally, especially since there are many cases where |
---|
82 | # separate q vectors are needed (e.g., form in python and structure |
---|
83 | # in opencl; or both in opencl, but one in single precision and the |
---|
84 | # other in double precision). |
---|
85 | kernels = [part(q_vectors) for part in self.parts] |
---|
86 | return MixtureKernel(self.info, kernels) |
---|
87 | |
---|
88 | def release(self): |
---|
89 | """ |
---|
90 | Free resources associated with the model. |
---|
91 | """ |
---|
92 | for part in self.parts: |
---|
93 | part.release() |
---|
94 | |
---|
95 | |
---|
96 | class MixtureKernel(object): |
---|
97 | def __init__(self, model_info, kernels): |
---|
98 | dim = '2d' if kernels[0].q_input.is_2d else '1d' |
---|
99 | |
---|
100 | # fixed offsets starts at 2 for scale and background |
---|
101 | fixed_pars, pd_pars = [], [] |
---|
102 | offsets = [[2, 0]] |
---|
103 | #vol_index = [] |
---|
104 | def accumulate(fixed, pd, volume): |
---|
105 | # subtract 1 from fixed since we are removing background |
---|
106 | fixed_offset, pd_offset = offsets[-1] |
---|
107 | #vol_index.extend(k+pd_offset for k,v in pd if v in volume) |
---|
108 | offsets.append([fixed_offset + len(fixed) - 1, pd_offset + len(pd)]) |
---|
109 | pd_pars.append(pd) |
---|
110 | if dim == '2d': |
---|
111 | for p in kernels: |
---|
112 | partype = p.info['partype'] |
---|
113 | accumulate(partype['fixed-2d'], partype['pd-2d'], partype['volume']) |
---|
114 | else: |
---|
115 | for p in kernels: |
---|
116 | partype = p.info['partype'] |
---|
117 | accumulate(partype['fixed-1d'], partype['pd-1d'], partype['volume']) |
---|
118 | |
---|
119 | #self.vol_index = vol_index |
---|
120 | self.offsets = offsets |
---|
121 | self.fixed_pars = fixed_pars |
---|
122 | self.pd_pars = pd_pars |
---|
123 | self.info = model_info |
---|
124 | self.kernels = kernels |
---|
125 | self.results = None |
---|
126 | |
---|
127 | def __call__(self, fixed_pars, pd_pars, cutoff=1e-5): |
---|
128 | scale, background = fixed_pars[0:2] |
---|
129 | total = 0.0 |
---|
130 | self.results = [] # remember the parts for plotting later |
---|
131 | for k in range(len(self.offsets)-1): |
---|
132 | start_fixed, start_pd = self.offsets[k] |
---|
133 | end_fixed, end_pd = self.offsets[k+1] |
---|
134 | part_fixed = [fixed_pars[start_fixed], 0.0] + fixed_pars[start_fixed+1:end_fixed] |
---|
135 | part_pd = [pd_pars[start_pd], 0.0] + pd_pars[start_pd+1:end_pd] |
---|
136 | part_result = self.kernels[k](part_fixed, part_pd) |
---|
137 | total += part_result |
---|
138 | self.results.append(scale*sum+background) |
---|
139 | |
---|
140 | return scale*total + background |
---|
141 | |
---|
142 | def release(self): |
---|
143 | self.p_kernel.release() |
---|
144 | self.q_kernel.release() |
---|
145 | |
---|