source: sasmodels/sasmodels/product.py @ 17bbadd

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

refactor so all model defintion queries use model_info; better documentation of model_info structure; initial implementation of product model (broken)

  • Property mode set to 100644
File size: 6.4 KB
Line 
1"""
2Product 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"""
13import numpy as np
14
15from .core import call_ER_VR
16from .generate import process_parameters
17
18SCALE=0
19BACKGROUND=1
20EFFECT_RADIUS=2
21VOLFRACTION=3
22
23def make_product_info(p_info, s_info):
24    """
25    Create info block for product model.
26    """
27    p_id, p_name, p_pars = p_info['id'], p_info['name'], p_info['parameters']
28    s_id, s_name, s_pars = s_info['id'], s_info['name'], s_info['parameters']
29    # We require models to start with scale and background
30    assert s_pars[SCALE][0] == 'scale'
31    assert s_pars[BACKGROUND][0] == 'background'
32    # We require structure factors to start with effect radius and volfraction
33    assert s_pars[EFFECT_RADIUS][0] == 'effect_radius'
34    assert s_pars[VOLFRACTION][0] == 'volfraction'
35    # Combine the parameter sets.  We are skipping the first three
36    # parameters of S since scale, background are defined in P and
37    # effect_radius is set from P.ER().
38    pars = p_pars + s_pars[3:]
39    # check for duplicates; can't use assertion since they may never be checked
40    if len(set(p[0] for p in pars)) != len(pars):
41        raise ValueError("Duplicate parameters in %s and %s"%(p_id))
42    # For comparison with sasview, determine the old parameters.
43    oldname = [p_info['oldname'], s_info['oldname']]
44    oldpars = {'scale':'scale_factor'}
45    oldpars.update(p_info['oldpars'])
46    oldpars.update(s_info['oldpars'])
47
48    model_info = {}
49    model_info['id'] = '*'.join((p_id, s_id))
50    model_info['name'] = ' X '.join((p_name, s_name))
51    model_info['filename'] = None
52    model_info['title'] = 'Product of %s and structure factor %s'%(p_name, s_name)
53    model_info['description'] = model_info['title']
54    model_info['category'] = "custom"
55    model_info['parameters'] = pars
56    # Remember the component info blocks so we can build the product model
57    model_info['composition'] = ('product', [p_info, s_info])
58    model_info['oldname'] = oldname
59    model_info['oldpars'] = oldpars
60    process_parameters(model_info)
61    return model_info
62
63class ProductModel(object):
64    def __init__(self, P, S):
65        self.P = P
66        self.S = S
67        self.info = make_product_info(P.info, S.info)
68
69    def __call__(self, q_vectors):
70        # Note: may be sending the q_vectors to the GPU twice even though they
71        # are only needed once.  It would mess up modularity quite a bit to
72        # handle this optimally, especially since there are many cases where
73        # separate q vectors are needed (e.g., form in python and structure
74        # in opencl; or both in opencl, but one in single precision and the
75        # other in double precision).
76        p_kernel = self.P(q_vectors)
77        s_kernel = self.S(q_vectors)
78        return ProductKernel(self.info, p_kernel, s_kernel)
79
80    def release(self):
81        """
82        Free resources associated with the model.
83        """
84        self.P.release()
85        self.S.release()
86
87
88class ProductKernel(object):
89    def __init__(self, model_info, p_kernel, s_kernel):
90        dim = '2d' if p_kernel.q_input.is_2d else '1d'
91
92        # Need to know if we want 2D and magnetic parameters when constructing
93        # a parameter map.
94        par_map = {}
95        p_info = p_kernel.info['partype']
96        s_info = s_kernel.info['partype']
97        vol_pars = set(p_info['volume'])
98        if dim == '2d':
99            num_p_fixed = len(p_info['fixed-2d'])
100            num_p_pd = len(p_info['pd-2d'])
101            num_s_fixed = len(s_info['fixed-2d'])
102            num_s_pd = len(s_info['pd-2d'])
103            # volume parameters are amongst the pd pars for P, not S
104            vol_par_idx = [k for k,v in enumerate(p_info['pd-2d'])
105                           if v in vol_pars]
106        else:
107            num_p_fixed = len(p_info['fixed-1d'])
108            num_p_pd = len(p_info['pd-1d'])
109            num_s_fixed = len(s_info['fixed-1d'])
110            num_s_pd = len(s_info['pd-1d'])
111            # volume parameters are amongst the pd pars for P, not S
112            vol_par_idx = [k for k,v in enumerate(p_info['pd-1d'])
113                           if v in vol_pars]
114
115        start = 0
116        par_map['p_fixed'] = np.arange(start, start+num_p_fixed)
117        # User doesn't set scale, background or effect_radius for S in P*S,
118        # so borrow values from end of p_fixed.  This makes volfraction the
119        # first S parameter.
120        start += num_p_fixed - 3
121        par_map['s_fixed'] = np.arange(start, start+num_s_fixed)
122        par_map['volfraction'] = num_p_fixed
123        start += num_s_fixed
124        # vol pars offset from the start of pd pars
125        par_map['vol_pars'] = [start+k for k in vol_par_idx]
126        par_map['p_pd'] = np.arange(start, start+num_p_pd)
127        start += num_p_pd
128        par_map['s_pd'] = np.arange(start, start+num_s_pd)  # should be empty...
129
130
131        self.fixed_pars = model_info['partype']['fixed-' + dim]
132        self.pd_pars = model_info['partype']['pd-' + dim]
133        self.info = model_info
134        self.p_kernel = p_kernel
135        self.s_kernel = s_kernel
136        self.par_map = par_map
137
138    def __call__(self, fixed_pars, pd_pars, cutoff=1e-5):
139        pars = fixed_pars + pd_pars
140        scale = pars[SCALE]
141        background = pars[BACKGROUND]
142        s_volfraction = pars[self.par_map['volfraction']]
143        p_fixed = [pars[k] for k in self.par_map['p_fixed']]
144        s_fixed = [pars[k] for k in self.par_map['s_fixed']]
145        p_pd = [pars[k] for k in self.par_map['p_pd']]
146        s_pd = [pars[k] for k in self.par_map['s_pd']]
147        vol_pars = [pars[k] for k in self.par_map['vol_pars']]
148        p_pars = dict(list(zip(self.fixed_pars, p_fixed))
149                      +list(zip(self.pd_pars, p_pd)))
150
151        effect_radius, vol_ratio = call_ER_VR(self.p_kernel.info, vol_pars)
152
153        p_fixed[SCALE] = s_volfraction
154        p_fixed[BACKGROUND] = 0.0
155        s_fixed[SCALE] = scale
156        s_fixed[BACKGROUND] = 0.0
157        s_fixed[EFFECT_RADIUS] = effect_radius
158        s_fixed[VOLFRACTION] = s_volfraction/vol_ratio
159
160        p_res = self.p_kernel(p_fixed, p_pd)
161        s_res = self.s_kernel(s_fixed, s_pd)
162
163        return p_res*s_res + background
164
165    def release(self):
166        self.p_kernel.release()
167        self.q_kernel.release()
168
Note: See TracBrowser for help on using the repository browser.