Changeset aea515c in sasmodels


Ignore:
Timestamp:
Oct 24, 2016 6:17:19 PM (7 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
999c87f
Parents:
97f9b46
Message:

convert model parameters from SasView? 3.1 to sasmodels

Location:
sasmodels
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/conversion_table.py

    r1b49bf8 raea515c  
    143143        } 
    144144    ], 
    145     "core_shell_ellipsoid(nonXT)": [ 
     145    "core_shell_ellipsoid:1": [ 
    146146        "CoreShellEllipsoidModel", 
    147147        { 
     
    150150            "sld_solvent": "sld_solvent", 
    151151            "equat_core": "equat_core", 
    152             "equat_shell": "equat_shell", 
    153             "polar_core": "polar_core", 
    154             "polar_shell": "polar_shell", 
     152            "thick_shell": "equat_shell", 
     153            "x_core": "polar_core", 
     154            "x_polar_shell": "polar_shell", 
    155155            "theta": "axis_theta", 
    156156            "phi": "axis_phi", 
     
    399399        "RectangularHollowPrismModel", 
    400400        { 
     401            "sld": "sldPipe", 
     402            "sld_solvent": "sldSolv", 
     403            "length_a": "short_side", 
    401404            "b2a_ratio": "b2a_ratio", 
    402             "length_a": "short_side", 
    403             "sld": "sldPipe", 
    404             "length_c": "c2a_ratio", 
    405             "sld_solvent": "sldSolv", 
    406             "thickness": "thickness" 
     405            "c2a_ratio": "c2a_ratio", 
     406            "thickness": "thickness", 
    407407        } 
    408408    ], 
     
    411411        { 
    412412            "sld": "sldPipe", 
     413            "sld_solvent": "sldSolv", 
     414            "length_a": "short_side", 
    413415            "b2a_ratio": "b2a_ratio", 
    414             "length_a": "short_side", 
    415             "length_c": "c2a_ratio", 
    416             "sld_solvent": "sldSolv" 
     416            "c2a_ratio": "c2a_ratio", 
    417417        } 
    418418    ], 
     
    515515        { 
    516516            "rg": "rg", 
    517             "scale": "scale", 
    518             "background": "background" 
     517            "i_zero": "scale", 
     518            "background": "background", 
    519519        } 
    520520    ], 
     
    576576        { 
    577577            "scale": "scale", 
    578             "string_thickness": "thick_string", 
     578            "thick_string": "thick_string", 
    579579            "sld_string": "sld_string", 
    580580            "sld_solvent": "sld_solv", 
    581581            "edge_sep": "edge_separation", 
    582             "number_of_pearls": "num_pearls", 
     582            "num_pearls": "num_pearls", 
    583583            "radius": "radius", 
    584584            "background": "background", 
     
    591591            "rg": "rg", 
    592592            "polydispersity": "poly_m", 
    593             "scale": "scale", 
    594             "background": "background" 
     593            "i_zero": "scale", 
     594            "background": "background", 
    595595        } 
    596596    ], 
     
    638638        { 
    639639            "scale": "scale", 
    640             "solvent_sld": "sld_solvent", 
     640            "sld_solvent": "sld_solvent", 
    641641            "thickness": "thickness", 
    642642            "beta": "beta", 
     
    665665        { 
    666666            "sld": "sldPipe", 
     667            "length_a": "short_side", 
    667668            "b2a_ratio": "b2a_ratio", 
    668             "length_a": "short_side", 
    669             "length_c": "c2a_ratio", 
     669            "c2a_ratio": "c2a_ratio", 
    670670            "sld_solvent": "sldSolv" 
    671671        } 
  • sasmodels/convert.py

    r1b49bf8 raea515c  
    88 
    99from .conversion_table import CONVERSION_TABLE 
     10from .core import load_model_info 
    1011 
    1112# List of models which SasView versions don't contain the explicit 'scale' argument. 
     
    5455    ] 
    5556 
    56 def _convert_pars(pars, mapping): 
    57     """ 
    58     Rename the parameters and any associated polydispersity attributes. 
    59     """ 
    60     newpars = pars.copy() 
    61     for new, old in mapping.items(): 
    62         if old == new: continue 
    63         for underscore, dot in PD_DOT: 
    64             if old+dot in newpars: 
    65                 if new is not None: 
    66                     newpars[new+underscore] = pars[old+dot] 
    67                 del newpars[old+dot] 
    68     return newpars 
    69  
    70 def convert_model(name, pars): 
    71     """ 
    72     Convert model from old style parameter names to new style. 
    73     """ 
    74     _, _ = name, pars # lint 
    75     raise NotImplementedError 
    76     # need to load all new models in order to determine old=>new 
    77     # model name mapping 
    78  
    79 def _unscale(par, scale): 
     57def _rescale(par, scale): 
    8058    return [pk*scale for pk in par] if isinstance(par, list) else par*scale 
    8159 
    82 def _is_sld(modelinfo, id): 
     60def _is_sld(model_info, id): 
     61    """ 
     62    Return True if parameter is a magnetic magnitude or SLD parameter. 
     63    """ 
    8364    if id.startswith('M0:'): 
    8465        return True 
    85     if (id.endswith('_pd') or id.endswith('_pd_n') or id.endswith('_pd_nsigma') 
    86             or id.endswith('_pd_width') or id.endswith('_pd_type')): 
     66    if '_pd' in id or '.' in id: 
    8767        return False 
    88     for p in modelinfo.parameters.call_parameters: 
     68    for p in model_info.parameters.call_parameters: 
    8969        if p.id == id: 
    9070            return p.type == 'sld' 
    9171    # check through kernel parameters in case it is a named as a vector 
    92     for p in modelinfo.parameters.kernel_parameters: 
     72    for p in model_info.parameters.kernel_parameters: 
    9373        if p.id == id: 
    9474            return p.type == 'sld' 
    9575    raise ValueError("unknown parameter %r in conversion"%id) 
    9676 
    97 def _unscale_sld(modelinfo, pars): 
    98     """ 
    99     rescale all sld parameters in the new model definition by 1e6 so the 
     77def _rescale_sld(model_info, pars, scale): 
     78    """ 
     79    rescale all sld parameters in the new model definition by *scale* so the 
    10080    numbers are nicer.  Relies on the fact that all sld parameters in the 
    101     new model definition end with sld. 
    102     """ 
    103     return dict((id, (_unscale(v, 1e-6) if _is_sld(modelinfo, id) else v)) 
     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)) 
    10485                for id, v in pars.items()) 
    10586 
    106 def _remove_pd(pars, key, name): 
    107     """ 
    108     Remove polydispersity from the parameter list. 
    109  
    110     Note: operates in place 
    111     """ 
    112     # Bumps style parameter names 
    113     width = pars.pop(key+".width", 0.0) 
    114     n_points = pars.pop(key+".npts", 0) 
    115     if width != 0.0 and n_points != 0: 
    116         warnings.warn("parameter %s not polydisperse in sasview %s"%(key, name)) 
    117     pars.pop(key+".nsigmas", None) 
    118     pars.pop(key+".type", None) 
    119     return pars 
    120  
    121 def _revert_pars(pars, mapping): 
    122     """ 
    123     Rename the parameters and any associated polydispersity attributes. 
    124     """ 
    125     newpars = pars.copy() 
    126  
    127     for new, old in mapping.items(): 
    128         for underscore, dot in PD_DOT: 
    129             if old and old+underscore == new+dot: 
    130                 continue 
    131             if new+underscore in newpars: 
    132                 if old is not None: 
    133                     newpars[old+dot] = pars[new+underscore] 
    134                 del newpars[new+underscore] 
    135     for k in list(newpars.keys()): 
    136         for underscore, dot in PD_DOT[1:]:  # skip "" => "" 
    137             if k.endswith(underscore): 
    138                 newpars[k[:-len(underscore)]+dot] = newpars[k] 
    139                 del newpars[k] 
    140     return newpars 
    141  
    142 def revert_name(model_info): 
    143     oldname, _ = CONVERSION_TABLE.get(model_info.id, [None, {}]) 
    144     return oldname 
    14587 
    14688def _get_translation_table(model_info): 
     
    160102    return translation 
    161103 
     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 
    162289def _trim_vectors(model_info, pars, oldpars): 
    163290    _, translation = CONVERSION_TABLE.get(model_info.id, [None, {}]) 
     
    185312    else: 
    186313        translation = _get_translation_table(model_info) 
    187     oldpars = _revert_pars(_unscale_sld(model_info, pars), translation) 
     314    oldpars = _revert_pars(_rescale_sld(model_info, pars, 1e-6), translation) 
    188315    oldpars = _trim_vectors(model_info, pars, oldpars) 
    189316 
     
    229356        elif name == 'hollow_cylinder': 
    230357            # now uses radius and thickness 
    231             oldpars['radius'] += oldpars['core_radius'] 
    232         elif name in ['mono_gauss_coil', 'poly_gauss_coil']: 
    233             del oldpars['i_zero'] 
     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'] 
    234363        elif name == 'onion': 
    235364            oldpars.pop('n_shells', None) 
     
    339468            pars['background'] = 0 
    340469        elif name == 'mono_gauss_coil': 
    341             pars['i_zero'] = 1 
     470            pars['scale'] = 1 
    342471        elif name == 'onion': 
    343472            pars['n_shells'] = math.ceil(pars['n_shells']) 
     
    346475            pars['number_of_pearls_pd_n'] = 0 
    347476        elif name == 'poly_gauss_coil': 
    348             pars['i_zero'] = 1 
     477            pars['scale'] = 1 
    349478        elif name == 'rpa': 
    350479            pars['case_num'] = int(pars['case_num']) 
     
    359488        elif name == 'teubner_strey': 
    360489            pars['scale'] = 1 
    361  
     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 TracChangeset for help on using the changeset viewer.