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 .modelinfo import Parameter, ParameterTable, ModelInfo |
---|
17 | |
---|
18 | def make_mixture_info(parts): |
---|
19 | """ |
---|
20 | Create info block for product model. |
---|
21 | """ |
---|
22 | flatten = [] |
---|
23 | for part in parts: |
---|
24 | if part['composition'] and part['composition'][0] == 'mixture': |
---|
25 | flatten.extend(part['compostion'][1]) |
---|
26 | else: |
---|
27 | flatten.append(part) |
---|
28 | parts = flatten |
---|
29 | |
---|
30 | # Build new parameter list |
---|
31 | pars = [] |
---|
32 | for k, part in enumerate(parts): |
---|
33 | # Parameter prefix per model, A_, B_, ... |
---|
34 | # Note that prefix must also be applied to id and length_control |
---|
35 | # to support vector parameters |
---|
36 | prefix = chr(ord('A')+k) + '_' |
---|
37 | pars.append(Parameter(prefix+'scale')) |
---|
38 | for p in part['parameters'].kernel_pars: |
---|
39 | p = copy(p) |
---|
40 | p.name = prefix+p.name |
---|
41 | p.id = prefix+p.id |
---|
42 | if p.length_control is not None: |
---|
43 | p.length_control = prefix+p.length_control |
---|
44 | pars.append(p) |
---|
45 | partable = ParameterTable(pars) |
---|
46 | |
---|
47 | model_info = ModelInfo() |
---|
48 | model_info.id = '+'.join(part['id']) |
---|
49 | model_info.name = ' + '.join(part['name']) |
---|
50 | model_info.filename = None |
---|
51 | model_info.title = 'Mixture model with ' + model_info.name |
---|
52 | model_info.description = model_info.title |
---|
53 | model_info.docs = model_info.title |
---|
54 | model_info.category = "custom" |
---|
55 | model_info.parameters = partable |
---|
56 | #model_info.single = any(part['single'] for part in parts) |
---|
57 | model_info.structure_factor = False |
---|
58 | model_info.variant_info = None |
---|
59 | #model_info.tests = [] |
---|
60 | #model_info.source = [] |
---|
61 | # Iq, Iqxy, form_volume, ER, VR and sesans |
---|
62 | # Remember the component info blocks so we can build the model |
---|
63 | model_info.composition = ('mixture', parts) |
---|
64 | |
---|
65 | |
---|
66 | class MixtureModel(object): |
---|
67 | def __init__(self, model_info, parts): |
---|
68 | self.info = model_info |
---|
69 | self.parts = parts |
---|
70 | |
---|
71 | def __call__(self, q_vectors): |
---|
72 | # Note: may be sending the q_vectors to the n times even though they |
---|
73 | # are only needed once. It would mess up modularity quite a bit to |
---|
74 | # handle this optimally, especially since there are many cases where |
---|
75 | # separate q vectors are needed (e.g., form in python and structure |
---|
76 | # in opencl; or both in opencl, but one in single precision and the |
---|
77 | # other in double precision). |
---|
78 | kernels = [part(q_vectors) for part in self.parts] |
---|
79 | return MixtureKernel(self.info, kernels) |
---|
80 | |
---|
81 | def release(self): |
---|
82 | """ |
---|
83 | Free resources associated with the model. |
---|
84 | """ |
---|
85 | for part in self.parts: |
---|
86 | part.release() |
---|
87 | |
---|
88 | |
---|
89 | class MixtureKernel(object): |
---|
90 | def __init__(self, model_info, kernels): |
---|
91 | dim = '2d' if kernels[0].q_input.is_2d else '1d' |
---|
92 | |
---|
93 | # fixed offsets starts at 2 for scale and background |
---|
94 | fixed_pars, pd_pars = [], [] |
---|
95 | offsets = [[2, 0]] |
---|
96 | #vol_index = [] |
---|
97 | def accumulate(fixed, pd, volume): |
---|
98 | # subtract 1 from fixed since we are removing background |
---|
99 | fixed_offset, pd_offset = offsets[-1] |
---|
100 | #vol_index.extend(k+pd_offset for k,v in pd if v in volume) |
---|
101 | offsets.append([fixed_offset + len(fixed) - 1, pd_offset + len(pd)]) |
---|
102 | pd_pars.append(pd) |
---|
103 | if dim == '2d': |
---|
104 | for p in kernels: |
---|
105 | partype = p.info.partype |
---|
106 | accumulate(partype['fixed-2d'], partype['pd-2d'], partype['volume']) |
---|
107 | else: |
---|
108 | for p in kernels: |
---|
109 | partype = p.info.partype |
---|
110 | accumulate(partype['fixed-1d'], partype['pd-1d'], partype['volume']) |
---|
111 | |
---|
112 | #self.vol_index = vol_index |
---|
113 | self.offsets = offsets |
---|
114 | self.fixed_pars = fixed_pars |
---|
115 | self.pd_pars = pd_pars |
---|
116 | self.info = model_info |
---|
117 | self.kernels = kernels |
---|
118 | self.results = None |
---|
119 | |
---|
120 | def __call__(self, fixed_pars, pd_pars, cutoff=1e-5): |
---|
121 | scale, background = fixed_pars[0:2] |
---|
122 | total = 0.0 |
---|
123 | self.results = [] # remember the parts for plotting later |
---|
124 | for k in range(len(self.offsets)-1): |
---|
125 | start_fixed, start_pd = self.offsets[k] |
---|
126 | end_fixed, end_pd = self.offsets[k+1] |
---|
127 | part_fixed = [fixed_pars[start_fixed], 0.0] + fixed_pars[start_fixed+1:end_fixed] |
---|
128 | part_pd = [pd_pars[start_pd], 0.0] + pd_pars[start_pd+1:end_pd] |
---|
129 | part_result = self.kernels[k](part_fixed, part_pd) |
---|
130 | total += part_result |
---|
131 | self.results.append(scale*sum+background) |
---|
132 | |
---|
133 | return scale*total + background |
---|
134 | |
---|
135 | def release(self): |
---|
136 | self.p_kernel.release() |
---|
137 | self.q_kernel.release() |
---|
138 | |
---|