1 | """ |
---|
2 | Product 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 | import numpy as np |
---|
14 | |
---|
15 | from .details import dispersion_mesh |
---|
16 | from .modelinfo import suffix_parameter, ParameterTable, Parameter, ModelInfo |
---|
17 | |
---|
18 | # TODO: make estimates available to constraints |
---|
19 | #ESTIMATED_PARAMETERS = [ |
---|
20 | # ["est_effect_radius", "A", 0.0, [0, np.inf], "", "Estimated effective radius"], |
---|
21 | # ["est_volume_ratio", "", 1.0, [0, np.inf], "", "Estimated volume ratio"], |
---|
22 | #] |
---|
23 | |
---|
24 | # TODO: core_shell_sphere model has suppressed the volume ratio calculation |
---|
25 | # revert it after making VR and ER available at run time as constraints. |
---|
26 | def make_product_info(p_info, s_info): |
---|
27 | """ |
---|
28 | Create info block for product model. |
---|
29 | """ |
---|
30 | p_id, p_name, p_partable = p_info.id, p_info.name, p_info.parameters |
---|
31 | s_id, s_name, s_partable = s_info.id, s_info.name, s_info.parameters |
---|
32 | p_set = set(p.id for p in p_partable) |
---|
33 | s_set = set(p.id for p in s_partable) |
---|
34 | |
---|
35 | if p_set & s_set: |
---|
36 | # there is some overlap between the parameter names; tag the |
---|
37 | # overlapping S parameters with name_S |
---|
38 | s_pars = [(suffix_parameter(par, "_S") if par.id in p_set else par) |
---|
39 | for par in s_partable.kernel_parameters] |
---|
40 | pars = p_partable.kernel_parameters + s_pars |
---|
41 | else: |
---|
42 | pars= p_partable.kernel_parameters + s_partable.kernel_parameters |
---|
43 | |
---|
44 | model_info = ModelInfo() |
---|
45 | model_info.id = '*'.join((p_id, s_id)) |
---|
46 | model_info.name = ' X '.join((p_name, s_name)) |
---|
47 | model_info.filename = None |
---|
48 | model_info.title = 'Product of %s and %s'%(p_name, s_name) |
---|
49 | model_info.description = model_info.title |
---|
50 | model_info.docs = model_info.title |
---|
51 | model_info.category = "custom" |
---|
52 | model_info.parameters = ParameterTable(pars) |
---|
53 | #model_info.single = p_info.single and s_info.single |
---|
54 | model_info.structure_factor = False |
---|
55 | model_info.variant_info = None |
---|
56 | #model_info.tests = [] |
---|
57 | #model_info.source = [] |
---|
58 | # Iq, Iqxy, form_volume, ER, VR and sesans |
---|
59 | model_info.composition = ('product', [p_info, s_info]) |
---|
60 | return model_info |
---|
61 | |
---|
62 | class ProductModel(object): |
---|
63 | def __init__(self, model_info, P, S): |
---|
64 | self.info = model_info |
---|
65 | self.P = P |
---|
66 | self.S = S |
---|
67 | |
---|
68 | def __call__(self, q_vectors): |
---|
69 | # Note: may be sending the q_vectors to the GPU twice even though they |
---|
70 | # are only needed once. It would mess up modularity quite a bit to |
---|
71 | # handle this optimally, especially since there are many cases where |
---|
72 | # separate q vectors are needed (e.g., form in python and structure |
---|
73 | # in opencl; or both in opencl, but one in single precision and the |
---|
74 | # other in double precision). |
---|
75 | p_kernel = self.P(q_vectors) |
---|
76 | s_kernel = self.S(q_vectors) |
---|
77 | return ProductKernel(self.info, p_kernel, s_kernel) |
---|
78 | |
---|
79 | def release(self): |
---|
80 | """ |
---|
81 | Free resources associated with the model. |
---|
82 | """ |
---|
83 | self.P.release() |
---|
84 | self.S.release() |
---|
85 | |
---|
86 | |
---|
87 | class ProductKernel(object): |
---|
88 | def __init__(self, model_info, p_kernel, s_kernel): |
---|
89 | self.info = model_info |
---|
90 | self.p_kernel = p_kernel |
---|
91 | self.s_kernel = s_kernel |
---|
92 | |
---|
93 | def __call__(self, details, weights, values, cutoff): |
---|
94 | effect_radius, vol_ratio = call_ER_VR(self.p_kernel.info, vol_pars) |
---|
95 | |
---|
96 | p_fixed[SCALE] = s_volfraction |
---|
97 | p_fixed[BACKGROUND] = 0.0 |
---|
98 | s_fixed[SCALE] = scale |
---|
99 | s_fixed[BACKGROUND] = 0.0 |
---|
100 | s_fixed[2] = s_volfraction/vol_ratio |
---|
101 | s_pd[0] = [effect_radius], [1.0] |
---|
102 | |
---|
103 | p_res = self.p_kernel(p_details, p_weights, p_values, cutoff) |
---|
104 | s_res = self.s_kernel(s_details, s_weights, s_values, cutoff) |
---|
105 | #print s_fixed, s_pd, p_fixed, p_pd |
---|
106 | |
---|
107 | return p_res*s_res + background |
---|
108 | |
---|
109 | def release(self): |
---|
110 | self.p_kernel.release() |
---|
111 | self.q_kernel.release() |
---|
112 | |
---|
113 | def call_ER_VR(model_info, vol_pars): |
---|
114 | """ |
---|
115 | Return effect radius and volume ratio for the model. |
---|
116 | """ |
---|
117 | value, weight = dispersion_mesh(vol_pars) |
---|
118 | |
---|
119 | individual_radii = model_info.ER(*value) if model_info.ER else 1.0 |
---|
120 | whole, part = model_info.VR(*value) if model_info.VR else (1.0, 1.0) |
---|
121 | |
---|
122 | effect_radius = np.sum(weight*individual_radii) / np.sum(weight) |
---|
123 | volume_ratio = np.sum(weight*part)/np.sum(weight*whole) |
---|
124 | return effect_radius, volume_ratio |
---|