source: sasmodels/sasmodels/convert.py @ aea515c

core_shell_microgelscostrafo411magnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since aea515c was aea515c, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

convert model parameters from SasView? 3.1 to sasmodels
  • Property mode set to 100644
File size: 19.8 KB
Line 
1"""
2Convert models to and from sasview.
3"""
4from __future__ import print_function, division
5
6import math
7import warnings
8
9from .conversion_table import CONVERSION_TABLE
10from .core import load_model_info
11
12# List of models which SasView versions don't contain the explicit 'scale' argument.
13# When converting such a model, please update this list.
14MODELS_WITHOUT_SCALE = [
15    'teubner_strey',
16    'broad_peak',
17    'two_lorentzian',
18    "two_power_law",
19    'gauss_lorentz_gel',
20    'be_polyelectrolyte',
21    'correlation_length',
22    'fractal_core_shell',
23    'binary_hard_sphere',
24    'raspberry'
25]
26
27# List of models which SasView versions don't contain the explicit 'background' argument.
28# When converting such a model, please update this list.
29MODELS_WITHOUT_BACKGROUND = [
30    'guinier',
31]
32
33MODELS_WITHOUT_VOLFRACTION = [
34    'fractal',
35    'vesicle',
36    'multilayer_vesicle',
37]
38
39MAGNETIC_SASVIEW_MODELS = [
40    'core_shell',
41    'core_multi_shell',
42    'cylinder',
43    'parallelepiped',
44    'sphere',
45]
46
47
48# Convert new style names for polydispersity info to old style names
49PD_DOT = [
50    ("", ""),
51    ("_pd", ".width"),
52    ("_pd_n", ".npts"),
53    ("_pd_nsigma", ".nsigmas"),
54    ("_pd_type", ".type"),
55    ]
56
57def _rescale(par, scale):
58    return [pk*scale for pk in par] if isinstance(par, list) else par*scale
59
60def _is_sld(model_info, id):
61    """
62    Return True if parameter is a magnetic magnitude or SLD parameter.
63    """
64    if id.startswith('M0:'):
65        return True
66    if '_pd' in id or '.' in id:
67        return False
68    for p in model_info.parameters.call_parameters:
69        if p.id == id:
70            return p.type == 'sld'
71    # check through kernel parameters in case it is a named as a vector
72    for p in model_info.parameters.kernel_parameters:
73        if p.id == id:
74            return p.type == 'sld'
75    raise ValueError("unknown parameter %r in conversion"%id)
76
77def _rescale_sld(model_info, pars, scale):
78    """
79    rescale all sld parameters in the new model definition by *scale* so the
80    numbers are nicer.  Relies on the fact that all sld parameters in the
81    new model definition end with sld.  For backward conversion use
82    *scale=1e-6*.  For forward conversion use *scale=1e6*.
83    """
84    return dict((id, (_rescale(v, scale) if _is_sld(model_info, id) else v))
85                for id, v in pars.items())
86
87
88def _get_translation_table(model_info):
89    _, translation = CONVERSION_TABLE.get(model_info.id, [None, {}])
90    translation = translation.copy()
91    for p in model_info.parameters.kernel_parameters:
92        if p.length > 1:
93            newid = p.id
94            oldid = translation.get(p.id, p.id)
95            translation.pop(newid, None)
96            for k in range(1, p.length+1):
97                if newid+str(k) not in translation:
98                    translation[newid+str(k)] = oldid+str(k)
99    # Remove control parameter from the result
100    if model_info.control:
101        translation[model_info.control] = "CONTROL"
102    return translation
103
104# ========= FORWARD CONVERSION sasview 3.x => sasmodels ===========
105def _dot_pd_to_underscore_pd(par):
106    if par.endswith(".width"):
107        return par[:-6]+"_pd"
108    elif par.endswith(".type"):
109        return par[:-5]+"_pd_type"
110    elif par.endswith(".nsigmas"):
111        return par[:-8]+"_pd_nsigma"
112    elif par.endswith(".npts"):
113        return par[:-5]+"_pd_n"
114    else:
115        return par
116
117def _pd_to_underscores(pars):
118    return dict((_dot_pd_to_underscore_pd(k), v) for k, v in pars.items())
119
120def _convert_pars(pars, mapping):
121    """
122    Rename the parameters and any associated polydispersity attributes.
123    """
124    newpars = pars.copy()
125    for new, old in mapping.items():
126        if old == new: continue
127        if old is None: continue
128        for underscore, dot in PD_DOT:
129            source = old+dot
130            if source in newpars:
131                if new is not None:
132                    target = new+dot
133                else:
134                    target = None
135                if source != target:
136                    if target:
137                        newpars[target] = pars[old+dot]
138                    del newpars[source]
139    return newpars
140
141
142def _conversion_target(model_name):
143    """
144    Find the sasmodel name which translates into the sasview name.
145
146    Note: *CoreShellEllipsoidModel* translates into *core_shell_ellipsoid:1*.
147    This is necessary since there is only one variant in sasmodels for the
148    two variants in sasview.
149    """
150    for sasmodels_name, [sasview_name, _] in CONVERSION_TABLE.items():
151        if sasview_name == model_name:
152            return sasmodels_name
153    return None
154
155
156def _hand_convert(name, oldpars):
157    if name == 'core_shell_parallelepiped':
158        # Make sure pd on rim parameters defaults to zero
159        # ... probably not necessary.
160        oldpars['rimA.width'] = 0.0
161        oldpars['rimB.width'] = 0.0
162        oldpars['rimC.width'] = 0.0
163    elif name == 'hollow_cylinder':
164        # now uses radius and thickness
165        thickness = oldpars['radius'] - oldpars['core_radius']
166        pd = oldpars['radius.width']*oldpars['radius']/thickness
167        oldpars['radius'] = thickness
168        oldpars['radius.width'] = pd
169    elif name == 'pearl_necklace':
170        pass
171        #_remove_pd(oldpars, 'num_pearls', name)
172        #_remove_pd(oldpars, 'thick_string', name)
173    elif name == 'polymer_micelle':
174        if 'ndensity' in oldpars:
175            oldpars['ndensity'] /= 1e15
176    elif name == 'rpa':
177        # convert scattering lengths from femtometers to centimeters
178        for p in "L1", "L2", "L3", "L4":
179            if p in oldpars: oldpars[p] /= 1e-13
180    elif name == 'spherical_sld':
181        oldpars["CONTROL"] += 1
182    elif name == 'teubner_strey':
183        # basically undoing the entire Teubner-Strey calculations here.
184        #    drho = (sld_a - sld_b)
185        #    k = 2.0*math.pi*xi/d
186        #    a2 = (1.0 + k**2)**2
187        #    c1 = 2.0 * xi**2 * (1.0 - k**2)
188        #    c2 = xi**4
189        #    prefactor = 8.0*math.pi*phi*(1.0-phi)*drho**2*c2/xi
190        #    scale = 1e-4*prefactor
191        #    oldpars['scale'] = a2/scale
192        #    oldpars['c1'] = c1/scale
193        #    oldpars['c2'] = c2/scale
194
195        # need xi, d, sld_a, sld_b, phi=volfraction_a
196        # assume contrast is 1.0e-6, scale=1, background=0
197        sld_a, sld_b = 1.0, 0.
198        drho = sld_a - sld_b
199
200        # find xi
201        p_scale = oldpars['scale']
202        p_c1 = oldpars['c1']
203        p_c2= oldpars['c2']
204        xi = math.sqrt(2/(math.sqrt(p_scale/p_c2) + 0.5*p_c1/p_c2))
205
206        # find d from xi
207        k = math.sqrt(1 - 0.5*p_c1/p_c2*xi**2)
208        d = 2*math.pi*xi/k
209
210        # solve quadratic phi (1-phi) = scale/(1e-4 8 pi drho^2 xi^3)
211        # favour volume fraction in [0, 0.5]
212        c = xi/(p_c2*1e-4*8.0*math.pi*drho**2)
213        phi = 0.5 - math.sqrt(0.25 - c)
214
215        # scale sld_a by 1e-6 because the translator will scale it back
216        oldpars.update(volfraction_a=phi, xi=xi, d=d, sld_a=sld_a*1e-6,
217                       sld_b=sld_b, scale=1.0)
218        oldpars.pop('c1')
219        oldpars.pop('c2')
220
221    return oldpars
222
223def convert_model(name, pars, use_underscore=False):
224    """
225    Convert model from old style parameter names to new style.
226    """
227    newname = _conversion_target(name)
228    if newname is None:
229        return name, pars
230    if newname.endswith(':1'):
231        model_info = load_model_info(newname[:-2])
232        # Know that the table exists and isn't multiplicity so grab it directly
233        translation = CONVERSION_TABLE[newname]
234    else:
235        model_info = load_model_info(newname)
236        translation = _get_translation_table(model_info)
237    newpars = _hand_convert(newname, pars.copy())
238    newpars = _convert_pars(newpars, translation)
239    newpars = _rescale_sld(model_info, newpars, 1e6)
240    newpars.setdefault('scale', 1.0)
241    newpars.setdefault('background', 0.0)
242    if use_underscore:
243        newpars = _pd_to_underscores(newpars)
244    return newname, newpars
245
246
247# ========= BACKWARD CONVERSION sasmodels => sasview 3.x ===========
248
249def _revert_pars(pars, mapping):
250    """
251    Rename the parameters and any associated polydispersity attributes.
252    """
253    newpars = pars.copy()
254
255    for new, old in mapping.items():
256        for underscore, dot in PD_DOT:
257            if old and old+underscore == new+dot:
258                continue
259            if new+underscore in newpars:
260                if old is not None:
261                    newpars[old+dot] = pars[new+underscore]
262                del newpars[new+underscore]
263    for k in list(newpars.keys()):
264        for underscore, dot in PD_DOT[1:]:  # skip "" => ""
265            if k.endswith(underscore):
266                newpars[k[:-len(underscore)]+dot] = newpars[k]
267                del newpars[k]
268    return newpars
269
270def revert_name(model_info):
271    oldname, _ = CONVERSION_TABLE.get(model_info.id, [None, {}])
272    return oldname
273
274def _remove_pd(pars, key, name):
275    """
276    Remove polydispersity from the parameter list.
277
278    Note: operates in place
279    """
280    # Bumps style parameter names
281    width = pars.pop(key+".width", 0.0)
282    n_points = pars.pop(key+".npts", 0)
283    if width != 0.0 and n_points != 0:
284        warnings.warn("parameter %s not polydisperse in sasview %s"%(key, name))
285    pars.pop(key+".nsigmas", None)
286    pars.pop(key+".type", None)
287    return pars
288
289def _trim_vectors(model_info, pars, oldpars):
290    _, translation = CONVERSION_TABLE.get(model_info.id, [None, {}])
291    for p in model_info.parameters.kernel_parameters:
292        if p.length_control is not None:
293            n = int(pars[p.length_control])
294            oldname = translation.get(p.id, p.id)
295            for k in range(n+1, p.length+1):
296                for _, old in PD_DOT:
297                    oldpars.pop(oldname+str(k)+old, None)
298    return oldpars
299
300def revert_pars(model_info, pars):
301    """
302    Convert model from new style parameter names to old style.
303    """
304    if model_info.composition is not None:
305        composition_type, parts = model_info.composition
306        if composition_type == 'product':
307            translation = {'scale':'scale_factor'}
308            translation.update(_get_translation_table(parts[0]))
309            translation.update(_get_translation_table(parts[1]))
310        else:
311            raise NotImplementedError("cannot convert to sasview sum")
312    else:
313        translation = _get_translation_table(model_info)
314    oldpars = _revert_pars(_rescale_sld(model_info, pars, 1e-6), translation)
315    oldpars = _trim_vectors(model_info, pars, oldpars)
316
317    # Make sure the control parameter is an integer
318    if "CONTROL" in oldpars:
319        oldpars["CONTROL"] = int(oldpars["CONTROL"])
320
321    # Note: update compare.constrain_pars to match
322    name = model_info.id
323    if name in MODELS_WITHOUT_SCALE or model_info.structure_factor:
324        if oldpars.pop('scale', 1.0) != 1.0:
325            warnings.warn("parameter scale not used in sasview %s"%name)
326    if name in MODELS_WITHOUT_BACKGROUND or model_info.structure_factor:
327        if oldpars.pop('background', 0.0) != 0.0:
328            warnings.warn("parameter background not used in sasview %s"%name)
329
330    # Remove magnetic parameters from non-magnetic sasview models
331    if name not in MAGNETIC_SASVIEW_MODELS:
332        oldpars = dict((k, v) for k, v in oldpars.items() if ':' not in k)
333
334    # If it is a product model P*S, then check the individual forms for special
335    # cases.  Note: despite the structure factor alone not having scale or
336    # background, the product model does, so this is below the test for
337    # models without scale or background.
338    namelist = name.split('*') if '*' in name else [name]
339    for name in namelist:
340        if name in MODELS_WITHOUT_VOLFRACTION:
341            del oldpars['volfraction']
342        elif name == 'core_multi_shell':
343            # kill extra shells
344            for k in range(5, 11):
345                oldpars.pop('sld_shell'+str(k), 0)
346                oldpars.pop('thick_shell'+str(k), 0)
347                oldpars.pop('mtheta:sld'+str(k), 0)
348                oldpars.pop('mphi:sld'+str(k), 0)
349                oldpars.pop('M0:sld'+str(k), 0)
350                _remove_pd(oldpars, 'sld_shell'+str(k), 'sld')
351                _remove_pd(oldpars, 'thick_shell'+str(k), 'thickness')
352        elif name == 'core_shell_parallelepiped':
353            _remove_pd(oldpars, 'rimA', name)
354            _remove_pd(oldpars, 'rimB', name)
355            _remove_pd(oldpars, 'rimC', name)
356        elif name == 'hollow_cylinder':
357            # now uses radius and thickness
358            thickness = oldpars['core_radius']
359            oldpars['radius'] += thickness
360            oldpars['radius.width'] *= thickness/oldpars['radius']
361        #elif name in ['mono_gauss_coil', 'poly_gauss_coil']:
362        #    del oldpars['i_zero']
363        elif name == 'onion':
364            oldpars.pop('n_shells', None)
365        elif name == 'pearl_necklace':
366            _remove_pd(oldpars, 'num_pearls', name)
367            _remove_pd(oldpars, 'thick_string', name)
368        elif name == 'polymer_micelle':
369            if 'ndensity' in oldpars:
370                oldpars['ndensity'] *= 1e15
371        elif name == 'rpa':
372            # convert scattering lengths from femtometers to centimeters
373            for p in "L1", "L2", "L3", "L4":
374                if p in oldpars: oldpars[p] *= 1e-13
375            if pars['case_num'] < 2:
376                for k in ("a", "b"):
377                    for p in ("L", "N", "Phi", "b", "v"):
378                        oldpars.pop(p+k, None)
379                for k in "Kab,Kac,Kad,Kbc,Kbd".split(','):
380                    oldpars.pop(k, None)
381            elif pars['case_num'] < 5:
382                for k in ("a",):
383                    for p in ("L", "N", "Phi", "b", "v"):
384                        oldpars.pop(p+k, None)
385                for k in "Kab,Kac,Kad".split(','):
386                    oldpars.pop(k, None)
387        elif name == 'spherical_sld':
388            oldpars["CONTROL"] -= 1
389            # remove polydispersity from shells
390            for k in range(1, 11):
391                _remove_pd(oldpars, 'thick_flat'+str(k), 'thickness')
392                _remove_pd(oldpars, 'thick_inter'+str(k), 'interface')
393            # remove extra shells
394            for k in range(int(pars['n_shells']), 11):
395                oldpars.pop('sld_flat'+str(k), 0)
396                oldpars.pop('thick_flat'+str(k), 0)
397                oldpars.pop('thick_inter'+str(k), 0)
398                oldpars.pop('func_inter'+str(k), 0)
399                oldpars.pop('nu_inter'+str(k), 0)
400        elif name == 'stacked_disks':
401            _remove_pd(oldpars, 'n_stacking', name)
402        elif name == 'teubner_strey':
403            # basically redoing the entire Teubner-Strey calculations here.
404            volfraction = oldpars.pop('volfraction_a')
405            xi = oldpars.pop('xi')
406            d = oldpars.pop('d')
407            sld_a = oldpars.pop('sld_a')
408            sld_b = oldpars.pop('sld_b')
409            drho = 1e6*(sld_a - sld_b)  # conversion autoscaled these
410            k = 2.0*math.pi*xi/d
411            a2 = (1.0 + k**2)**2
412            c1 = 2.0 * xi**2 * (1.0 - k**2)
413            c2 = xi**4
414            prefactor = 8.0*math.pi*volfraction*(1.0-volfraction)*drho**2*c2/xi
415            scale = 1e-4*prefactor
416            oldpars['scale'] = a2/scale
417            oldpars['c1'] = c1/scale
418            oldpars['c2'] = c2/scale
419
420    #print("convert from",list(sorted(pars)))
421    #print("convert to",list(sorted(oldpars.items())))
422    return oldpars
423
424def constrain_new_to_old(model_info, pars):
425    """
426    Restrict parameter values to those that will match sasview.
427    """
428    name = model_info.id
429    # Note: update convert.revert_model to match
430    if name in MODELS_WITHOUT_SCALE or model_info.structure_factor:
431        pars['scale'] = 1
432    if name in MODELS_WITHOUT_BACKGROUND or model_info.structure_factor:
433        pars['background'] = 0
434    # sasview multiplies background by structure factor
435    if '*' in name:
436        pars['background'] = 0
437
438    # Shut off magnetism when comparing non-magnetic sasview models
439    if name not in MAGNETIC_SASVIEW_MODELS:
440        suppress_magnetism = False
441        for key in pars.keys():
442            if key.startswith("M0:"):
443                suppress_magnetism = suppress_magnetism or (pars[key] != 0)
444                pars[key] = 0
445        if suppress_magnetism:
446            warnings.warn("suppressing magnetism for comparison with sasview")
447
448    # Shut off theta polydispersity since algorithm has changed
449    if 'theta_pd_n' in pars:
450        if pars['theta_pd_n'] != 0:
451            warnings.warn("suppressing theta polydispersity for comparison with sasview")
452        pars['theta_pd_n'] = 0
453
454    # If it is a product model P*S, then check the individual forms for special
455    # cases.  Note: despite the structure factor alone not having scale or
456    # background, the product model does, so this is below the test for
457    # models without scale or background.
458    namelist = name.split('*') if '*' in name else [name]
459    for name in namelist:
460        if name in MODELS_WITHOUT_VOLFRACTION:
461            pars['volfraction'] = 1
462        if name == 'core_multi_shell':
463            pars['n'] = min(math.ceil(pars['n']), 4)
464        elif name == 'gel_fit':
465            pars['scale'] = 1
466        elif name == 'line':
467            pars['scale'] = 1
468            pars['background'] = 0
469        elif name == 'mono_gauss_coil':
470            pars['scale'] = 1
471        elif name == 'onion':
472            pars['n_shells'] = math.ceil(pars['n_shells'])
473        elif name == 'pearl_necklace':
474            pars['string_thickness_pd_n'] = 0
475            pars['number_of_pearls_pd_n'] = 0
476        elif name == 'poly_gauss_coil':
477            pars['scale'] = 1
478        elif name == 'rpa':
479            pars['case_num'] = int(pars['case_num'])
480        elif name == 'spherical_sld':
481            pars['n_shells'] = math.ceil(pars['n_shells'])
482            pars['n_steps'] = math.ceil(pars['n_steps'])
483            for k in range(1, 11):
484                pars['shape%d'%k] = math.trunc(pars['shape%d'%k]+0.5)
485            for k in range(2, 11):
486                pars['thickness%d_pd_n'%k] = 0
487                pars['interface%d_pd_n'%k] = 0
488        elif name == 'teubner_strey':
489            pars['scale'] = 1
490            if pars['volfraction_a'] > 0.5:
491                pars['volfraction_a'] = 1.0 - pars['volfraction_a']
492        elif name == 'unified_power_Rg':
493            pars['level'] = int(pars['level'])
494
495def _check_one(name, seed=None):
496    """
497    Generate a random set of parameters for *name*, and check that they can
498    be converted back to SasView 3.x and forward again to sasmodels.  Raises
499    an error if the parameters are changed.
500    """
501    from . import compare
502
503    model_info = load_model_info(name)
504
505    old_name = revert_name(model_info)
506    if old_name is None:
507        return
508
509    pars = compare.get_pars(model_info, use_demo=False)
510    pars = compare.randomize_pars(model_info, pars, seed=seed)
511    if name == "teubner_strey":
512        # T-S model is underconstrained, so fix the assumptions.
513        pars['sld_a'], pars['sld_b'] = 1.0, 0.0
514    compare.constrain_pars(model_info, pars)
515    constrain_new_to_old(model_info, pars)
516    old_pars = revert_pars(model_info, pars)
517    new_name, new_pars = convert_model(old_name, old_pars, use_underscore=True)
518    if 1:
519        print("==== %s in ====="%name)
520        print(str(compare.parlist(model_info, pars, True)))
521        print("==== %s ====="%old_name)
522        for k, v in sorted(old_pars.items()):
523            print(k, v)
524        print("==== %s out ====="%new_name)
525        print(str(compare.parlist(model_info, new_pars, True)))
526    assert name==new_name, "%r != %r"%(name, new_name)
527    for k, v in new_pars.items():
528        assert k in pars, "%s: %r appeared from conversion"%(name, k)
529        if isinstance(v, float):
530            assert abs(v-pars[k])<=abs(1e-12*v), "%s: %r  %s != %s"%(name, k, v, pars[k])
531        else:
532            assert v == pars[k], "%s: %r  %s != %s"%(name, k, v, pars[k])
533    for k, v in pars.items():
534        assert k in pars, "%s: %r not converted"%(name, k)
535
536def test_backward_forward():
537    from .core import list_models
538    for name in list_models('all'):
539        L = lambda: _check_one(name, seed=1)
540        L.description = name
541        yield L
Note: See TracBrowser for help on using the repository browser.