[72a081d] | 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 | |
---|