source: sasmodels/sasmodels/convert.py @ e85aa2e

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

support scale ↔ scale_factor conversion in P*S models
  • Property mode set to 100644
File size: 20.0 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 == 'core_shell_ellipsoid:1':
164    elif name == 'hollow_cylinder':
165        # now uses radius and thickness
166        thickness = oldpars['radius'] - oldpars['core_radius']
167        pd = oldpars['radius.width']*oldpars['radius']/thickness
168        oldpars['radius'] = thickness
169        oldpars['radius.width'] = pd
170    elif name == 'pearl_necklace':
171        pass
172        #_remove_pd(oldpars, 'num_pearls', name)
173        #_remove_pd(oldpars, 'thick_string', name)
174    elif name == 'polymer_micelle':
175        if 'ndensity' in oldpars:
176            oldpars['ndensity'] /= 1e15
177    elif name == 'rpa':
178        # convert scattering lengths from femtometers to centimeters
179        for p in "L1", "L2", "L3", "L4":
180            if p in oldpars: oldpars[p] /= 1e-13
181    elif name == 'spherical_sld':
182        oldpars["CONTROL"] += 1
183    elif name == 'teubner_strey':
184        # basically undoing the entire Teubner-Strey calculations here.
185        #    drho = (sld_a - sld_b)
186        #    k = 2.0*math.pi*xi/d
187        #    a2 = (1.0 + k**2)**2
188        #    c1 = 2.0 * xi**2 * (1.0 - k**2)
189        #    c2 = xi**4
190        #    prefactor = 8.0*math.pi*phi*(1.0-phi)*drho**2*c2/xi
191        #    scale = 1e-4*prefactor
192        #    oldpars['scale'] = a2/scale
193        #    oldpars['c1'] = c1/scale
194        #    oldpars['c2'] = c2/scale
195
196        # need xi, d, sld_a, sld_b, phi=volfraction_a
197        # assume contrast is 1.0e-6, scale=1, background=0
198        sld_a, sld_b = 1.0, 0.
199        drho = sld_a - sld_b
200
201        # find xi
202        p_scale = oldpars['scale']
203        p_c1 = oldpars['c1']
204        p_c2= oldpars['c2']
205        xi = math.sqrt(2/(math.sqrt(p_scale/p_c2) + 0.5*p_c1/p_c2))
206
207        # find d from xi
208        k = math.sqrt(1 - 0.5*p_c1/p_c2*xi**2)
209        d = 2*math.pi*xi/k
210
211        # solve quadratic phi (1-phi) = xi/(1e-4 8 pi drho^2 c2)
212        # favour volume fraction in [0, 0.5]
213        c = xi / (1e-4 * 8.0 * math.pi * drho**2 * p_c2)
214        phi = 0.5 - math.sqrt(0.25 - c)
215
216        # scale sld_a by 1e-6 because the translator will scale it back
217        oldpars.update(volfraction_a=phi, xi=xi, d=d, sld_a=sld_a*1e-6,
218                       sld_b=sld_b, scale=1.0)
219        oldpars.pop('c1')
220        oldpars.pop('c2')
221
222    return oldpars
223
224def convert_model(name, pars, use_underscore=False):
225    """
226    Convert model from old style parameter names to new style.
227    """
228    newname = _conversion_target(name)
229    if newname is None:
230        return name, pars
231    if ':' in newname:   # core_shell_ellipsoid:1
232        model_info = load_model_info(newname[:-2])
233        # Know that the table exists and isn't multiplicity so grab it directly
234        # Can't use _get_translation_table since that will return the 'bare'
235        # version.
236        translation = CONVERSION_TABLE[newname]
237    else:
238        model_info = load_model_info(newname)
239        translation = _get_translation_table(model_info)
240    newpars = _hand_convert(newname, pars.copy())
241    newpars = _convert_pars(newpars, translation)
242    newpars = _rescale_sld(model_info, newpars, 1e6)
243    newpars.setdefault('scale', 1.0)
244    newpars.setdefault('background', 0.0)
245    if use_underscore:
246        newpars = _pd_to_underscores(newpars)
247    return newname, newpars
248
249
250# ========= BACKWARD CONVERSION sasmodels => sasview 3.x ===========
251
252def _revert_pars(pars, mapping):
253    """
254    Rename the parameters and any associated polydispersity attributes.
255    """
256    newpars = pars.copy()
257
258    for new, old in mapping.items():
259        for underscore, dot in PD_DOT:
260            if old and old+underscore == new+dot:
261                continue
262            if new+underscore in newpars:
263                if old is not None:
264                    newpars[old+dot] = pars[new+underscore]
265                del newpars[new+underscore]
266    for k in list(newpars.keys()):
267        for underscore, dot in PD_DOT[1:]:  # skip "" => ""
268            if k.endswith(underscore):
269                newpars[k[:-len(underscore)]+dot] = newpars[k]
270                del newpars[k]
271    return newpars
272
273def revert_name(model_info):
274    oldname, _ = CONVERSION_TABLE.get(model_info.id, [None, {}])
275    return oldname
276
277def _remove_pd(pars, key, name):
278    """
279    Remove polydispersity from the parameter list.
280
281    Note: operates in place
282    """
283    # Bumps style parameter names
284    width = pars.pop(key+".width", 0.0)
285    n_points = pars.pop(key+".npts", 0)
286    if width != 0.0 and n_points != 0:
287        warnings.warn("parameter %s not polydisperse in sasview %s"%(key, name))
288    pars.pop(key+".nsigmas", None)
289    pars.pop(key+".type", None)
290    return pars
291
292def _trim_vectors(model_info, pars, oldpars):
293    _, translation = CONVERSION_TABLE.get(model_info.id, [None, {}])
294    for p in model_info.parameters.kernel_parameters:
295        if p.length_control is not None:
296            n = int(pars[p.length_control])
297            oldname = translation.get(p.id, p.id)
298            for k in range(n+1, p.length+1):
299                for _, old in PD_DOT:
300                    oldpars.pop(oldname+str(k)+old, None)
301    return oldpars
302
303def revert_pars(model_info, pars):
304    """
305    Convert model from new style parameter names to old style.
306    """
307    if model_info.composition is not None:
308        composition_type, parts = model_info.composition
309        if composition_type == 'product':
310            translation = _get_translation_table(parts[0])
311            # structure factor models include scale:scale_factor mapping
312            translation.update(_get_translation_table(parts[1]))
313        else:
314            raise NotImplementedError("cannot convert to sasview sum")
315    else:
316        translation = _get_translation_table(model_info)
317    oldpars = _revert_pars(_rescale_sld(model_info, pars, 1e-6), translation)
318    oldpars = _trim_vectors(model_info, pars, oldpars)
319
320    # Make sure the control parameter is an integer
321    if "CONTROL" in oldpars:
322        oldpars["CONTROL"] = int(oldpars["CONTROL"])
323
324    # Note: update compare.constrain_pars to match
325    name = model_info.id
326    if name in MODELS_WITHOUT_SCALE or model_info.structure_factor:
327        if oldpars.pop('scale', 1.0) != 1.0:
328            warnings.warn("parameter scale not used in sasview %s"%name)
329    if name in MODELS_WITHOUT_BACKGROUND or model_info.structure_factor:
330        if oldpars.pop('background', 0.0) != 0.0:
331            warnings.warn("parameter background not used in sasview %s"%name)
332
333    # Remove magnetic parameters from non-magnetic sasview models
334    if name not in MAGNETIC_SASVIEW_MODELS:
335        oldpars = dict((k, v) for k, v in oldpars.items() if ':' not in k)
336
337    # If it is a product model P*S, then check the individual forms for special
338    # cases.  Note: despite the structure factor alone not having scale or
339    # background, the product model does, so this is below the test for
340    # models without scale or background.
341    namelist = name.split('*') if '*' in name else [name]
342    for name in namelist:
343        if name in MODELS_WITHOUT_VOLFRACTION:
344            del oldpars['volfraction']
345        elif name == 'core_multi_shell':
346            # kill extra shells
347            for k in range(5, 11):
348                oldpars.pop('sld_shell'+str(k), 0)
349                oldpars.pop('thick_shell'+str(k), 0)
350                oldpars.pop('mtheta:sld'+str(k), 0)
351                oldpars.pop('mphi:sld'+str(k), 0)
352                oldpars.pop('M0:sld'+str(k), 0)
353                _remove_pd(oldpars, 'sld_shell'+str(k), 'sld')
354                _remove_pd(oldpars, 'thick_shell'+str(k), 'thickness')
355        elif name == 'core_shell_parallelepiped':
356            _remove_pd(oldpars, 'rimA', name)
357            _remove_pd(oldpars, 'rimB', name)
358            _remove_pd(oldpars, 'rimC', name)
359        elif name == 'hollow_cylinder':
360            # now uses radius and thickness
361            thickness = oldpars['core_radius']
362            oldpars['radius'] += thickness
363            oldpars['radius.width'] *= thickness/oldpars['radius']
364        #elif name in ['mono_gauss_coil', 'poly_gauss_coil']:
365        #    del oldpars['i_zero']
366        elif name == 'onion':
367            oldpars.pop('n_shells', None)
368        elif name == 'pearl_necklace':
369            _remove_pd(oldpars, 'num_pearls', name)
370            _remove_pd(oldpars, 'thick_string', name)
371        elif name == 'polymer_micelle':
372            if 'ndensity' in oldpars:
373                oldpars['ndensity'] *= 1e15
374        elif name == 'rpa':
375            # convert scattering lengths from femtometers to centimeters
376            for p in "L1", "L2", "L3", "L4":
377                if p in oldpars: oldpars[p] *= 1e-13
378            if pars['case_num'] < 2:
379                for k in ("a", "b"):
380                    for p in ("L", "N", "Phi", "b", "v"):
381                        oldpars.pop(p+k, None)
382                for k in "Kab,Kac,Kad,Kbc,Kbd".split(','):
383                    oldpars.pop(k, None)
384            elif pars['case_num'] < 5:
385                for k in ("a",):
386                    for p in ("L", "N", "Phi", "b", "v"):
387                        oldpars.pop(p+k, None)
388                for k in "Kab,Kac,Kad".split(','):
389                    oldpars.pop(k, None)
390        elif name == 'spherical_sld':
391            oldpars["CONTROL"] -= 1
392            # remove polydispersity from shells
393            for k in range(1, 11):
394                _remove_pd(oldpars, 'thick_flat'+str(k), 'thickness')
395                _remove_pd(oldpars, 'thick_inter'+str(k), 'interface')
396            # remove extra shells
397            for k in range(int(pars['n_shells']), 11):
398                oldpars.pop('sld_flat'+str(k), 0)
399                oldpars.pop('thick_flat'+str(k), 0)
400                oldpars.pop('thick_inter'+str(k), 0)
401                oldpars.pop('func_inter'+str(k), 0)
402                oldpars.pop('nu_inter'+str(k), 0)
403        elif name == 'stacked_disks':
404            _remove_pd(oldpars, 'n_stacking', name)
405        elif name == 'teubner_strey':
406            # basically redoing the entire Teubner-Strey calculations here.
407            volfraction = oldpars.pop('volfraction_a')
408            xi = oldpars.pop('xi')
409            d = oldpars.pop('d')
410            sld_a = oldpars.pop('sld_a')
411            sld_b = oldpars.pop('sld_b')
412            drho = 1e6*(sld_a - sld_b)  # conversion autoscaled these
413            k = 2.0*math.pi*xi/d
414            a2 = (1.0 + k**2)**2
415            c1 = 2.0 * xi**2 * (1.0 - k**2)
416            c2 = xi**4
417            prefactor = 8.0*math.pi*volfraction*(1.0-volfraction)*drho**2*c2/xi
418            scale = 1e-4*prefactor
419            oldpars['scale'] = a2/scale
420            oldpars['c1'] = c1/scale
421            oldpars['c2'] = c2/scale
422
423    #print("convert from",list(sorted(pars)))
424    #print("convert to",list(sorted(oldpars.items())))
425    return oldpars
426
427def constrain_new_to_old(model_info, pars):
428    """
429    Restrict parameter values to those that will match sasview.
430    """
431    name = model_info.id
432    # Note: update convert.revert_model to match
433    if name in MODELS_WITHOUT_SCALE or model_info.structure_factor:
434        pars['scale'] = 1
435    if name in MODELS_WITHOUT_BACKGROUND or model_info.structure_factor:
436        pars['background'] = 0
437    # sasview multiplies background by structure factor
438    if '*' in name:
439        pars['background'] = 0
440
441    # Shut off magnetism when comparing non-magnetic sasview models
442    if name not in MAGNETIC_SASVIEW_MODELS:
443        suppress_magnetism = False
444        for key in pars.keys():
445            if key.startswith("M0:"):
446                suppress_magnetism = suppress_magnetism or (pars[key] != 0)
447                pars[key] = 0
448        if suppress_magnetism:
449            warnings.warn("suppressing magnetism for comparison with sasview")
450
451    # Shut off theta polydispersity since algorithm has changed
452    if 'theta_pd_n' in pars:
453        if pars['theta_pd_n'] != 0:
454            warnings.warn("suppressing theta polydispersity for comparison with sasview")
455        pars['theta_pd_n'] = 0
456
457    # If it is a product model P*S, then check the individual forms for special
458    # cases.  Note: despite the structure factor alone not having scale or
459    # background, the product model does, so this is below the test for
460    # models without scale or background.
461    namelist = name.split('*') if '*' in name else [name]
462    for name in namelist:
463        if name in MODELS_WITHOUT_VOLFRACTION:
464            pars['volfraction'] = 1
465        if name == 'core_multi_shell':
466            pars['n'] = min(math.ceil(pars['n']), 4)
467        elif name == 'gel_fit':
468            pars['scale'] = 1
469        elif name == 'line':
470            pars['scale'] = 1
471            pars['background'] = 0
472        elif name == 'mono_gauss_coil':
473            pars['scale'] = 1
474        elif name == 'onion':
475            pars['n_shells'] = math.ceil(pars['n_shells'])
476        elif name == 'pearl_necklace':
477            pars['string_thickness_pd_n'] = 0
478            pars['number_of_pearls_pd_n'] = 0
479        elif name == 'poly_gauss_coil':
480            pars['scale'] = 1
481        elif name == 'rpa':
482            pars['case_num'] = int(pars['case_num'])
483        elif name == 'spherical_sld':
484            pars['n_shells'] = math.ceil(pars['n_shells'])
485            pars['n_steps'] = math.ceil(pars['n_steps'])
486            for k in range(1, 11):
487                pars['shape%d'%k] = math.trunc(pars['shape%d'%k]+0.5)
488            for k in range(2, 11):
489                pars['thickness%d_pd_n'%k] = 0
490                pars['interface%d_pd_n'%k] = 0
491        elif name == 'teubner_strey':
492            pars['scale'] = 1
493            if pars['volfraction_a'] > 0.5:
494                pars['volfraction_a'] = 1.0 - pars['volfraction_a']
495        elif name == 'unified_power_Rg':
496            pars['level'] = int(pars['level'])
497
498def _check_one(name, seed=None):
499    """
500    Generate a random set of parameters for *name*, and check that they can
501    be converted back to SasView 3.x and forward again to sasmodels.  Raises
502    an error if the parameters are changed.
503    """
504    from . import compare
505
506    model_info = load_model_info(name)
507
508    old_name = revert_name(model_info)
509    if old_name is None:
510        return
511
512    pars = compare.get_pars(model_info, use_demo=False)
513    pars = compare.randomize_pars(model_info, pars, seed=seed)
514    if name == "teubner_strey":
515        # T-S model is underconstrained, so fix the assumptions.
516        pars['sld_a'], pars['sld_b'] = 1.0, 0.0
517    compare.constrain_pars(model_info, pars)
518    constrain_new_to_old(model_info, pars)
519    old_pars = revert_pars(model_info, pars)
520    new_name, new_pars = convert_model(old_name, old_pars, use_underscore=True)
521    if 1:
522        print("==== %s in ====="%name)
523        print(str(compare.parlist(model_info, pars, True)))
524        print("==== %s ====="%old_name)
525        for k, v in sorted(old_pars.items()):
526            print(k, v)
527        print("==== %s out ====="%new_name)
528        print(str(compare.parlist(model_info, new_pars, True)))
529    assert name==new_name, "%r != %r"%(name, new_name)
530    for k, v in new_pars.items():
531        assert k in pars, "%s: %r appeared from conversion"%(name, k)
532        if isinstance(v, float):
533            assert abs(v-pars[k])<=abs(1e-12*v), "%s: %r  %s != %s"%(name, k, v, pars[k])
534        else:
535            assert v == pars[k], "%s: %r  %s != %s"%(name, k, v, pars[k])
536    for k, v in pars.items():
537        assert k in pars, "%s: %r not converted"%(name, k)
538
539def test_backward_forward():
540    from .core import list_models
541    for name in list_models('all'):
542        L = lambda: _check_one(name, seed=1)
543        L.description = name
544        yield L
Note: See TracBrowser for help on using the repository browser.