source: sasmodels/sasmodels/convert.py @ ef476e6

ticket-1257-vesicle-productticket_1156ticket_822_more_unit_tests
Last change on this file since ef476e6 was ef476e6, checked in by smk78, 5 years ago

Revert some of the changes done for #1059 committed yesterday

  • Property mode set to 100644
File size: 24.9 KB
RevLine 
[3964f92]1"""
2Convert models to and from sasview.
3"""
[1b49bf8]4from __future__ import print_function, division
[256dfe1]5
[c5ac2b2]6import math
[3964f92]7import warnings
[32e3c9b]8
[32398dc]9import numpy as np
10
[32e3c9b]11from .conversion_table import CONVERSION_TABLE
[aea515c]12from .core import load_model_info
[3964f92]13
14# List of models which SasView versions don't contain the explicit 'scale' argument.
15# When converting such a model, please update this list.
[310ddcb]16MODELS_WITHOUT_SCALE = [
[3964f92]17    'teubner_strey',
18    'broad_peak',
19    'two_lorentzian',
[34d6cab]20    "two_power_law",
[3964f92]21    'gauss_lorentz_gel',
22    'be_polyelectrolyte',
23    'correlation_length',
[74f7238]24    'fractal_core_shell',
[17bbadd]25    'binary_hard_sphere',
[bad8b12]26    'raspberry'
[3964f92]27]
28
29# List of models which SasView versions don't contain the explicit 'background' argument.
30# When converting such a model, please update this list.
[310ddcb]31MODELS_WITHOUT_BACKGROUND = [
[3964f92]32    'guinier',
33]
34
[256dfe1]35MODELS_WITHOUT_VOLFRACTION = [
36    'fractal',
37    'vesicle',
38    'multilayer_vesicle',
39]
40
[f1765a2]41MAGNETIC_SASVIEW_MODELS = [
42    'core_shell',
43    'core_multi_shell',
44    'cylinder',
45    'parallelepiped',
46    'sphere',
47]
48
[256dfe1]49
[f247314]50# Convert new style names for polydispersity info to old style names
[3964f92]51PD_DOT = [
52    ("_pd", ".width"),
53    ("_pd_n", ".npts"),
54    ("_pd_nsigma", ".nsigmas"),
55    ("_pd_type", ".type"),
[e1ea6b5]56    (".lower", ".lower"),
57    (".upper", ".upper"),
58    (".fittable", ".fittable"),
59    (".std", ".std"),
60    (".units", ".units"),
61    ("", "")
[3964f92]62    ]
[f247314]63
[aea515c]64def _rescale(par, scale):
65    return [pk*scale for pk in par] if isinstance(par, list) else par*scale
[3964f92]66
[e65c3ba]67def _is_sld(model_info, par):
[3964f92]68    """
[aea515c]69    Return True if parameter is a magnetic magnitude or SLD parameter.
[3964f92]70    """
[e65c3ba]71    if par.startswith('M0:'):
[32e3c9b]72        return True
[e65c3ba]73    if '_pd' in par or '.' in par:
[32e3c9b]74        return False
[aea515c]75    for p in model_info.parameters.call_parameters:
[e65c3ba]76        if p.id == par:
[32e3c9b]77            return p.type == 'sld'
78    # check through kernel parameters in case it is a named as a vector
[aea515c]79    for p in model_info.parameters.kernel_parameters:
[e65c3ba]80        if p.id == par:
[32e3c9b]81            return p.type == 'sld'
[d2fb4af]82    return False
[32e3c9b]83
[aea515c]84def _rescale_sld(model_info, pars, scale):
[3964f92]85    """
[aea515c]86    rescale all sld parameters in the new model definition by *scale* so the
[3964f92]87    numbers are nicer.  Relies on the fact that all sld parameters in the
[aea515c]88    new model definition end with sld.  For backward conversion use
89    *scale=1e-6*.  For forward conversion use *scale=1e6*.
[3964f92]90    """
[e65c3ba]91    return dict((par, (_rescale(v, scale) if _is_sld(model_info, par) else v))
92                for par, v in pars.items())
[3964f92]93
[aea515c]94
[e65c3ba]95def _get_translation_table(model_info, version=(3, 1, 2)):
[07c8d46]96    conv_param = CONVERSION_TABLE.get(version, {}).get(model_info.id, [None, {}])
97    translation = conv_param[1].copy()
[aea515c]98    for p in model_info.parameters.kernel_parameters:
99        if p.length > 1:
100            newid = p.id
101            oldid = translation.get(p.id, p.id)
102            translation.pop(newid, None)
103            for k in range(1, p.length+1):
104                if newid+str(k) not in translation:
105                    translation[newid+str(k)] = oldid+str(k)
106    # Remove control parameter from the result
[21c93c3]107    control_pars = [p.id for p in model_info.parameters.kernel_parameters
108                    if p.is_control]
109    if control_pars:
110        control_id = control_pars[0]
111        translation[control_id] = "CONTROL"
[aea515c]112    return translation
113
114# ========= FORWARD CONVERSION sasview 3.x => sasmodels ===========
115def _dot_pd_to_underscore_pd(par):
116    if par.endswith(".width"):
117        return par[:-6]+"_pd"
118    elif par.endswith(".type"):
119        return par[:-5]+"_pd_type"
120    elif par.endswith(".nsigmas"):
121        return par[:-8]+"_pd_nsigma"
122    elif par.endswith(".npts"):
123        return par[:-5]+"_pd_n"
124    else:
125        return par
126
127def _pd_to_underscores(pars):
128    return dict((_dot_pd_to_underscore_pd(k), v) for k, v in pars.items())
129
130def _convert_pars(pars, mapping):
[3964f92]131    """
[aea515c]132    Rename the parameters and any associated polydispersity attributes.
133    """
134    newpars = pars.copy()
135    for new, old in mapping.items():
[e65c3ba]136        if old == new:
137            continue
138        if old is None:
139            continue
140        for _, dot in PD_DOT:
[aea515c]141            source = old+dot
142            if source in newpars:
143                if new is not None:
144                    target = new+dot
145                else:
146                    target = None
147                if source != target:
148                    if target:
149                        newpars[target] = pars[old+dot]
150                    del newpars[source]
151    return newpars
[3964f92]152
[32398dc]153def _conversion_target(model_name, version=(3, 1, 2)):
[3964f92]154    """
[aea515c]155    Find the sasmodel name which translates into the sasview name.
156
157    Note: *CoreShellEllipsoidModel* translates into *core_shell_ellipsoid:1*.
158    This is necessary since there is only one variant in sasmodels for the
159    two variants in sasview.
160    """
[119fa13]161    for sasmodels_name, sasview_dict in \
[a297255]162            CONVERSION_TABLE.get(version, {}).items():
[119fa13]163        if sasview_dict[0] == model_name:
[aea515c]164            return sasmodels_name
165    return None
166
[e65c3ba]167def _hand_convert(name, oldpars, version=(3, 1, 2)):
168    if version == (3, 1, 2):
[54fb5d8]169        oldpars = _hand_convert_3_1_2_to_4_1(name, oldpars)
[610ef23]170    if version < (4, 2, 0):
171        oldpars = _rename_magnetic_pars(oldpars)
[54fb5d8]172    return oldpars
173
[610ef23]174def _rename_magnetic_pars(pars):
175    """
176    Change from M0:par to par_M0, etc.
177    """
178    keys = list(pars.items())
179    for k in keys:
180        if k.startswith('M0:'):
181            pars[k[3:]+'_M0'] = pars.pop(k)
182        elif k.startswith('mtheta:'):
183            pars[k[7:]+'_mtheta'] = pars.pop(k)
184        elif k.startswith('mphi:'):
185            pars[k[5:]+'_mphi'] = pars.pop(k)
186        elif k.startswith('up:'):
187            pars['up_'+k[3:]] = pars.pop(k)
188    return pars
189
[54fb5d8]190def _hand_convert_3_1_2_to_4_1(name, oldpars):
191    if name == 'core_shell_parallelepiped':
[aea515c]192        # Make sure pd on rim parameters defaults to zero
193        # ... probably not necessary.
194        oldpars['rimA.width'] = 0.0
195        oldpars['rimB.width'] = 0.0
196        oldpars['rimC.width'] = 0.0
[54fb5d8]197    elif name == 'core_shell_ellipsoid:1':
[037aa3d]198        # Reverse translation (from new to old), from core_shell_ellipsoid.c
199        #    equat_shell = equat_core + thick_shell
200        #    polar_core = equat_core * x_core
201        #    polar_shell = equat_core * x_core + thick_shell*x_polar_shell
202        # Forward translation (from old to new), inverting reverse translation:
203        #    thick_shell = equat_shell - equat_core
204        #    x_core = polar_core / equat_core
205        #    x_polar_shell = (polar_shell - polar_core)/(equat_shell - equat_core)
206        # Auto translation (old <=> new) happens after hand_convert
207        #    equat_shell <=> thick_shell
208        #    polar_core <=> x_core
209        #    polar_shell <=> x_polar_shell
210        # So...
211        equat_core, equat_shell = oldpars['equat_core'], oldpars['equat_shell']
212        polar_core, polar_shell = oldpars['polar_core'], oldpars['polar_shell']
213        oldpars['equat_shell'] = equat_shell - equat_core
214        oldpars['polar_core'] = polar_core / equat_core
215        oldpars['polar_shell'] = (polar_shell-polar_core)/(equat_shell-equat_core)
[54fb5d8]216    elif name == 'hollow_cylinder':
[aea515c]217        # now uses radius and thickness
218        thickness = oldpars['radius'] - oldpars['core_radius']
219        oldpars['radius'] = thickness
[f8f5a37]220        if 'radius.width' in oldpars:
221            pd = oldpars['radius.width']*oldpars['radius']/thickness
222            oldpars['radius.width'] = pd
[54fb5d8]223    elif name == 'multilayer_vesicle':
[077666e]224        if 'scale' in oldpars:
225            oldpars['volfraction'] = oldpars['scale']
226            oldpars['scale'] = 1.0
227        if 'scale.lower' in oldpars:
228            oldpars['volfraction.lower'] = oldpars['scale.lower']
229        if 'scale.upper' in oldpars:
230            oldpars['volfraction.upper'] = oldpars['scale.upper']
231        if 'scale.fittable' in oldpars:
232            oldpars['volfraction.fittable'] = oldpars['scale.fittable']
233        if 'scale.std' in oldpars:
234            oldpars['volfraction.std'] = oldpars['scale.std']
235        if 'scale.units' in oldpars:
236            oldpars['volfraction.units'] = oldpars['scale.units']
[54fb5d8]237    elif name == 'pearl_necklace':
[aea515c]238        pass
239        #_remove_pd(oldpars, 'num_pearls', name)
240        #_remove_pd(oldpars, 'thick_string', name)
[54fb5d8]241    elif name == 'polymer_micelle':
[aea515c]242        if 'ndensity' in oldpars:
243            oldpars['ndensity'] /= 1e15
[8d06779]244        if 'ndensity.lower' in oldpars:
[88fafac]245            oldpars['ndensity.lower'] /= 1e15
[8d06779]246        if 'ndensity.upper' in oldpars:
[88fafac]247            oldpars['ndensity.upper'] /= 1e15
[54fb5d8]248    elif name == 'rpa':
[aea515c]249        # convert scattering lengths from femtometers to centimeters
250        for p in "L1", "L2", "L3", "L4":
[88fafac]251            if p in oldpars:
252                oldpars[p] /= 1e-13
[8d06779]253            if p + ".lower" in oldpars:
[88fafac]254                oldpars[p + ".lower"] /= 1e-13
[8d06779]255            if p + ".upper" in oldpars:
[88fafac]256                oldpars[p + ".upper"] /= 1e-13
[54fb5d8]257    elif name == 'spherical_sld':
[16fd9e5]258        j = 0
259        while "func_inter" + str(j) in oldpars:
260            name = "func_inter" + str(j)
261            new_name = "shape" + str(j + 1)
262            if oldpars[name] == 'Erf(|nu|*z)':
263                oldpars[new_name] = int(0)
264            elif oldpars[name] == 'RPower(z^|nu|)':
265                oldpars[new_name] = int(1)
266            elif oldpars[name] == 'LPower(z^|nu|)':
267                oldpars[new_name] = int(2)
268            elif oldpars[name] == 'RExp(-|nu|*z)':
269                oldpars[new_name] = int(3)
270            elif oldpars[name] == 'LExp(-|nu|*z)':
271                oldpars[new_name] = int(4)
272            else:
273                oldpars[new_name] = int(0)
274            oldpars.pop(name)
275            oldpars['n_shells'] = str(j + 1)
276            j += 1
[54fb5d8]277    elif name == 'teubner_strey':
[aea515c]278        # basically undoing the entire Teubner-Strey calculations here.
279        #    drho = (sld_a - sld_b)
280        #    k = 2.0*math.pi*xi/d
281        #    a2 = (1.0 + k**2)**2
282        #    c1 = 2.0 * xi**2 * (1.0 - k**2)
283        #    c2 = xi**4
284        #    prefactor = 8.0*math.pi*phi*(1.0-phi)*drho**2*c2/xi
285        #    scale = 1e-4*prefactor
286        #    oldpars['scale'] = a2/scale
287        #    oldpars['c1'] = c1/scale
288        #    oldpars['c2'] = c2/scale
289
290        # need xi, d, sld_a, sld_b, phi=volfraction_a
291        # assume contrast is 1.0e-6, scale=1, background=0
292        sld_a, sld_b = 1.0, 0.
293        drho = sld_a - sld_b
294
295        # find xi
296        p_scale = oldpars['scale']
297        p_c1 = oldpars['c1']
[2d81cfe]298        p_c2 = oldpars['c2']
[f8f5a37]299        i_1 = 0.5*p_c1/p_c2
300        i_2 = math.sqrt(math.fabs(p_scale/p_c2))
301        i_3 = 2/(i_1 + i_2)
302        xi = math.sqrt(math.fabs(i_3))
[aea515c]303
304        # find d from xi
[f8f5a37]305        k = math.sqrt(math.fabs(1 - 0.5*p_c1/p_c2*xi**2))
[aea515c]306        d = 2*math.pi*xi/k
307
[999c87f]308        # solve quadratic phi (1-phi) = xi/(1e-4 8 pi drho^2 c2)
[aea515c]309        # favour volume fraction in [0, 0.5]
[999c87f]310        c = xi / (1e-4 * 8.0 * math.pi * drho**2 * p_c2)
[aea515c]311        phi = 0.5 - math.sqrt(0.25 - c)
312
313        # scale sld_a by 1e-6 because the translator will scale it back
314        oldpars.update(volfraction_a=phi, xi=xi, d=d, sld_a=sld_a*1e-6,
315                       sld_b=sld_b, scale=1.0)
316        oldpars.pop('c1')
317        oldpars.pop('c2')
318
319    return oldpars
320
[e65c3ba]321def convert_model(name, pars, use_underscore=False, model_version=(3, 1, 2)):
[aea515c]322    """
323    Convert model from old style parameter names to new style.
324    """
[4b4eee1]325    newpars = pars
[0795293]326    keys = sorted(CONVERSION_TABLE.keys())
327    for i, version in enumerate(keys):
[a297255]328        # Don't allow indices outside list
329        next_i = i + 1
330        if next_i == len(keys):
331            next_i = i
[0795293]332        # If the save state is from a later version, skip the check
[a297255]333        if model_version <= keys[next_i]:
334            newname = _conversion_target(name, version)
335        else:
336            newname = None
[0795293]337        # If no conversion is found, move on
[4b4eee1]338        if newname is None:
339            newname = name
340            continue
341        if ':' in newname:   # core_shell_ellipsoid:1
342            model_info = load_model_info(newname[:-2])
[298621b]343            # Know the table exists and isn't multiplicity so grab it directly
[4b4eee1]344            # Can't use _get_translation_table since that will return the 'bare'
345            # version.
[a297255]346            translation = CONVERSION_TABLE.get(version, {})[newname][1]
[4b4eee1]347        else:
348            model_info = load_model_info(newname)
349            translation = _get_translation_table(model_info, version)
350        newpars = _hand_convert(newname, newpars, version)
351        newpars = _convert_pars(newpars, translation)
[86e8e66]352        # TODO: Still not convinced this is the best check
[e65c3ba]353        if not model_info.structure_factor and version == (3, 1, 2):
[4b4eee1]354            newpars = _rescale_sld(model_info, newpars, 1e6)
355        newpars.setdefault('scale', 1.0)
356        newpars.setdefault('background', 0.0)
357        if use_underscore:
358            newpars = _pd_to_underscores(newpars)
359        name = newname
[aea515c]360    return newname, newpars
361
362# ========= BACKWARD CONVERSION sasmodels => sasview 3.x ===========
[3964f92]363
364def _revert_pars(pars, mapping):
365    """
366    Rename the parameters and any associated polydispersity attributes.
367    """
368    newpars = pars.copy()
369
370    for new, old in mapping.items():
[40a87fa]371        for underscore, dot in PD_DOT:
372            if old and old+underscore == new+dot:
[3964f92]373                continue
[40a87fa]374            if new+underscore in newpars:
[3964f92]375                if old is not None:
[40a87fa]376                    newpars[old+dot] = pars[new+underscore]
377                del newpars[new+underscore]
[3964f92]378    for k in list(newpars.keys()):
[40a87fa]379        for underscore, dot in PD_DOT[1:]:  # skip "" => ""
380            if k.endswith(underscore):
381                newpars[k[:-len(underscore)]+dot] = newpars[k]
[3964f92]382                del newpars[k]
383    return newpars
384
[f247314]385def revert_name(model_info):
[b297ba9]386    """Translate model name back to the name used in SasView 3.x"""
[40a87fa]387    oldname, _ = CONVERSION_TABLE.get(model_info.id, [None, {}])
[f247314]388    return oldname
389
[aea515c]390def _remove_pd(pars, key, name):
391    """
392    Remove polydispersity from the parameter list.
393
394    Note: operates in place
395    """
396    # Bumps style parameter names
397    width = pars.pop(key+".width", 0.0)
398    n_points = pars.pop(key+".npts", 0)
399    if width != 0.0 and n_points != 0:
400        warnings.warn("parameter %s not polydisperse in sasview %s"%(key, name))
401    pars.pop(key+".nsigmas", None)
402    pars.pop(key+".type", None)
403    return pars
[256dfe1]404
405def _trim_vectors(model_info, pars, oldpars):
406    _, translation = CONVERSION_TABLE.get(model_info.id, [None, {}])
407    for p in model_info.parameters.kernel_parameters:
408        if p.length_control is not None:
409            n = int(pars[p.length_control])
410            oldname = translation.get(p.id, p.id)
411            for k in range(n+1, p.length+1):
412                for _, old in PD_DOT:
413                    oldpars.pop(oldname+str(k)+old, None)
[f247314]414    return oldpars
415
[17bbadd]416def revert_pars(model_info, pars):
[3964f92]417    """
418    Convert model from new style parameter names to old style.
419    """
[6d6508e]420    if model_info.composition is not None:
421        composition_type, parts = model_info.composition
[f247314]422        if composition_type == 'product':
[e85aa2e]423            translation = _get_translation_table(parts[0])
424            # structure factor models include scale:scale_factor mapping
[6dc78e4]425            translation.update(_get_translation_table(parts[1]))
[f247314]426        else:
427            raise NotImplementedError("cannot convert to sasview sum")
428    else:
[256dfe1]429        translation = _get_translation_table(model_info)
[aea515c]430    oldpars = _revert_pars(_rescale_sld(model_info, pars, 1e-6), translation)
[6dc78e4]431    oldpars = _trim_vectors(model_info, pars, oldpars)
[f247314]432
[bd49c79]433    # Make sure the control parameter is an integer
434    if "CONTROL" in oldpars:
435        oldpars["CONTROL"] = int(oldpars["CONTROL"])
[3964f92]436
437    # Note: update compare.constrain_pars to match
[6d6508e]438    name = model_info.id
439    if name in MODELS_WITHOUT_SCALE or model_info.structure_factor:
[3964f92]440        if oldpars.pop('scale', 1.0) != 1.0:
441            warnings.warn("parameter scale not used in sasview %s"%name)
[6d6508e]442    if name in MODELS_WITHOUT_BACKGROUND or model_info.structure_factor:
[3964f92]443        if oldpars.pop('background', 0.0) != 0.0:
444            warnings.warn("parameter background not used in sasview %s"%name)
[17bbadd]445
[f1765a2]446    # Remove magnetic parameters from non-magnetic sasview models
447    if name not in MAGNETIC_SASVIEW_MODELS:
[40a87fa]448        oldpars = dict((k, v) for k, v in oldpars.items() if ':' not in k)
[228cbd3]449
[17bbadd]450    # If it is a product model P*S, then check the individual forms for special
451    # cases.  Note: despite the structure factor alone not having scale or
452    # background, the product model does, so this is below the test for
453    # models without scale or background.
454    namelist = name.split('*') if '*' in name else [name]
455    for name in namelist:
[256dfe1]456        if name in MODELS_WITHOUT_VOLFRACTION:
457            del oldpars['volfraction']
[f1765a2]458        elif name == 'core_multi_shell':
459            # kill extra shells
460            for k in range(5, 11):
461                oldpars.pop('sld_shell'+str(k), 0)
462                oldpars.pop('thick_shell'+str(k), 0)
463                oldpars.pop('mtheta:sld'+str(k), 0)
464                oldpars.pop('mphi:sld'+str(k), 0)
465                oldpars.pop('M0:sld'+str(k), 0)
466                _remove_pd(oldpars, 'sld_shell'+str(k), 'sld')
467                _remove_pd(oldpars, 'thick_shell'+str(k), 'thickness')
468        elif name == 'core_shell_parallelepiped':
469            _remove_pd(oldpars, 'rimA', name)
[1b49bf8]470            _remove_pd(oldpars, 'rimB', name)
471            _remove_pd(oldpars, 'rimC', name)
472        elif name == 'hollow_cylinder':
473            # now uses radius and thickness
[aea515c]474            thickness = oldpars['core_radius']
475            oldpars['radius'] += thickness
476            oldpars['radius.width'] *= thickness/oldpars['radius']
477        #elif name in ['mono_gauss_coil', 'poly_gauss_coil']:
478        #    del oldpars['i_zero']
[f1765a2]479        elif name == 'onion':
480            oldpars.pop('n_shells', None)
[1b49bf8]481        elif name == 'pearl_necklace':
482            _remove_pd(oldpars, 'num_pearls', name)
483            _remove_pd(oldpars, 'thick_string', name)
484        elif name == 'polymer_micelle':
485            if 'ndensity' in oldpars:
486                oldpars['ndensity'] *= 1e15
[17bbadd]487        elif name == 'rpa':
488            # convert scattering lengths from femtometers to centimeters
[8bd7b77]489            for p in "L1", "L2", "L3", "L4":
[17bbadd]490                if p in oldpars: oldpars[p] *= 1e-13
[51ec7e8]491            if pars['case_num'] < 2:
492                for k in ("a", "b"):
493                    for p in ("L", "N", "Phi", "b", "v"):
494                        oldpars.pop(p+k, None)
495                for k in "Kab,Kac,Kad,Kbc,Kbd".split(','):
496                    oldpars.pop(k, None)
497            elif pars['case_num'] < 5:
498                for k in ("a",):
499                    for p in ("L", "N", "Phi", "b", "v"):
500                        oldpars.pop(p+k, None)
501                for k in "Kab,Kac,Kad".split(','):
502                    oldpars.pop(k, None)
[1b49bf8]503        elif name == 'spherical_sld':
504            oldpars["CONTROL"] -= 1
505            # remove polydispersity from shells
506            for k in range(1, 11):
507                _remove_pd(oldpars, 'thick_flat'+str(k), 'thickness')
508                _remove_pd(oldpars, 'thick_inter'+str(k), 'interface')
509            # remove extra shells
510            for k in range(int(pars['n_shells']), 11):
511                oldpars.pop('sld_flat'+str(k), 0)
512                oldpars.pop('thick_flat'+str(k), 0)
513                oldpars.pop('thick_inter'+str(k), 0)
514                oldpars.pop('func_inter'+str(k), 0)
515                oldpars.pop('nu_inter'+str(k), 0)
516        elif name == 'stacked_disks':
517            _remove_pd(oldpars, 'n_stacking', name)
518        elif name == 'teubner_strey':
519            # basically redoing the entire Teubner-Strey calculations here.
520            volfraction = oldpars.pop('volfraction_a')
521            xi = oldpars.pop('xi')
522            d = oldpars.pop('d')
523            sld_a = oldpars.pop('sld_a')
524            sld_b = oldpars.pop('sld_b')
525            drho = 1e6*(sld_a - sld_b)  # conversion autoscaled these
526            k = 2.0*math.pi*xi/d
527            a2 = (1.0 + k**2)**2
528            c1 = 2.0 * xi**2 * (1.0 - k**2)
529            c2 = xi**4
530            prefactor = 8.0*math.pi*volfraction*(1.0-volfraction)*drho**2*c2/xi
531            scale = 1e-4*prefactor
532            oldpars['scale'] = a2/scale
533            oldpars['c1'] = c1/scale
534            oldpars['c2'] = c2/scale
[17bbadd]535
[54bcd4a]536    #print("convert from",list(sorted(pars)))
537    #print("convert to",list(sorted(oldpars.items())))
[17bbadd]538    return oldpars
539
540def constrain_new_to_old(model_info, pars):
[3964f92]541    """
542    Restrict parameter values to those that will match sasview.
543    """
[6d6508e]544    name = model_info.id
[3964f92]545    # Note: update convert.revert_model to match
[6d6508e]546    if name in MODELS_WITHOUT_SCALE or model_info.structure_factor:
[3964f92]547        pars['scale'] = 1
[6d6508e]548    if name in MODELS_WITHOUT_BACKGROUND or model_info.structure_factor:
[0a4628d]549        pars['background'] = 0
[667a6f2]550    # sasview multiplies background by structure factor
551    if '*' in name:
552        pars['background'] = 0
[9a66e65]553
[f1765a2]554    # Shut off magnetism when comparing non-magnetic sasview models
555    if name not in MAGNETIC_SASVIEW_MODELS:
556        suppress_magnetism = False
557        for key in pars.keys():
558            if key.startswith("M0:"):
[2c74c11]559                suppress_magnetism = suppress_magnetism or (pars[key] != 0)
[f1765a2]560                pars[key] = 0
561        if suppress_magnetism:
562            warnings.warn("suppressing magnetism for comparison with sasview")
563
564    # Shut off theta polydispersity since algorithm has changed
565    if 'theta_pd_n' in pars:
566        if pars['theta_pd_n'] != 0:
567            warnings.warn("suppressing theta polydispersity for comparison with sasview")
568        pars['theta_pd_n'] = 0
569
[17bbadd]570    # If it is a product model P*S, then check the individual forms for special
571    # cases.  Note: despite the structure factor alone not having scale or
572    # background, the product model does, so this is below the test for
573    # models without scale or background.
574    namelist = name.split('*') if '*' in name else [name]
575    for name in namelist:
[256dfe1]576        if name in MODELS_WITHOUT_VOLFRACTION:
577            pars['volfraction'] = 1
[1b49bf8]578        if name == 'core_multi_shell':
579            pars['n'] = min(math.ceil(pars['n']), 4)
580        elif name == 'gel_fit':
581            pars['scale'] = 1
[17bbadd]582        elif name == 'line':
583            pars['scale'] = 1
584            pars['background'] = 0
[228cbd3]585        elif name == 'mono_gauss_coil':
[aea515c]586            pars['scale'] = 1
[d119f34]587        elif name == 'onion':
588            pars['n_shells'] = math.ceil(pars['n_shells'])
[1b49bf8]589        elif name == 'pearl_necklace':
590            pars['string_thickness_pd_n'] = 0
591            pars['number_of_pearls_pd_n'] = 0
592        elif name == 'poly_gauss_coil':
[aea515c]593            pars['scale'] = 1
[1b49bf8]594        elif name == 'rpa':
595            pars['case_num'] = int(pars['case_num'])
[a0494e9]596        elif name == 'spherical_sld':
[0dd4199]597            pars['n_shells'] = math.ceil(pars['n_shells'])
[54bcd4a]598            pars['n_steps'] = math.ceil(pars['n_steps'])
[51241113]599            for k in range(1, 11):
[54bcd4a]600                pars['shape%d'%k] = math.trunc(pars['shape%d'%k]+0.5)
[51241113]601            for k in range(2, 11):
[54bcd4a]602                pars['thickness%d_pd_n'%k] = 0
603                pars['interface%d_pd_n'%k] = 0
[1b49bf8]604        elif name == 'teubner_strey':
605            pars['scale'] = 1
[aea515c]606            if pars['volfraction_a'] > 0.5:
607                pars['volfraction_a'] = 1.0 - pars['volfraction_a']
[ef476e6]608        elif name == 'unified_power_Rg':
[aea515c]609            pars['level'] = int(pars['level'])
610
611def _check_one(name, seed=None):
612    """
613    Generate a random set of parameters for *name*, and check that they can
614    be converted back to SasView 3.x and forward again to sasmodels.  Raises
615    an error if the parameters are changed.
616    """
617    from . import compare
618
619    model_info = load_model_info(name)
620
621    old_name = revert_name(model_info)
622    if old_name is None:
623        return
624
625    pars = compare.get_pars(model_info, use_demo=False)
[32398dc]626    if seed is not None:
627        np.random.seed(seed)
628    pars = compare.randomize_pars(model_info, pars)
[aea515c]629    if name == "teubner_strey":
630        # T-S model is underconstrained, so fix the assumptions.
631        pars['sld_a'], pars['sld_b'] = 1.0, 0.0
632    compare.constrain_pars(model_info, pars)
633    constrain_new_to_old(model_info, pars)
634    old_pars = revert_pars(model_info, pars)
635    new_name, new_pars = convert_model(old_name, old_pars, use_underscore=True)
636    if 1:
637        print("==== %s in ====="%name)
638        print(str(compare.parlist(model_info, pars, True)))
639        print("==== %s ====="%old_name)
640        for k, v in sorted(old_pars.items()):
641            print(k, v)
642        print("==== %s out ====="%new_name)
643        print(str(compare.parlist(model_info, new_pars, True)))
[e65c3ba]644    assert name == new_name, "%r != %r"%(name, new_name)
[aea515c]645    for k, v in new_pars.items():
646        assert k in pars, "%s: %r appeared from conversion"%(name, k)
647        if isinstance(v, float):
[e65c3ba]648            assert abs(v-pars[k]) <= abs(1e-12*v), \
649                "%s: %r  %s != %s"%(name, k, v, pars[k])
[aea515c]650        else:
651            assert v == pars[k], "%s: %r  %s != %s"%(name, k, v, pars[k])
652    for k, v in pars.items():
653        assert k in pars, "%s: %r not converted"%(name, k)
[256dfe1]654
[aea515c]655def test_backward_forward():
[b297ba9]656    """
657    Test conversion of model parameters from 4.x to 3.x and back.
658    """
[aea515c]659    from .core import list_models
[a69d8cd]660    L = lambda name: _check_one(name, seed=1)
[aea515c]661    for name in list_models('all'):
[a69d8cd]662        yield L, name
Note: See TracBrowser for help on using the repository browser.