Changeset 1a580fb in sasmodels


Ignore:
Timestamp:
Feb 14, 2017 3:32:05 PM (8 years ago)
Author:
krzywon
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
e1aa129, 23054a1
Parents:
94d13f1 (diff), 9ed53e7 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' into sesans41

Files:
1 added
10 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/conversion_table.py

    rfa6d6fc r790b16bb  
    55names for each parameter in sasmodels.  This is used by :mod:`convert` to 
    66determine the equivalent parameter set when comparing a sasmodels model to 
    7 the models defined in SasView 3.1. 
     7the models defined in previous versions of SasView and sasmodels. This is now 
     8versioned based on the version number of SasView. 
     9 
     10When any sasmodels parameter or model name is changed, this must be modified to 
     11account for that. 
     12 
     13Usage: 
     14<old_Sasview_version> : { 
     15    <new_model_name> : [ 
     16        <old_model_name> , 
     17        { 
     18            <new_param_name_1> : <old_param_name_1>, 
     19            ... 
     20            <new_param_name_n> : <old_param_name_n> 
     21        } 
     22    ] 
     23} 
     24 
     25Any future parameter and model name changes can and should be given in this 
     26table for future compatibility. 
    827""" 
    928 
    10  
    1129CONVERSION_TABLE = { 
     30    (3,1,2) : { 
    1231    "adsorbed_layer": [ 
    1332        "Core2ndMomentModel", 
    … …  
    211230    ], 
    212231    "correlation_length": [ 
    213         "CorrLengthModel", 
     232        "CorrLength", 
    214233        { 
    215234            "porod_scale": "scale_p", 
    … …  
    218237            "lorentz_exp": "exponent_l", 
    219238            "cor_length": "length_l" 
    220         } 
     239        }, 
     240        "CorrLengthModel" 
    221241    ], 
    222242    "cylinder": [ 
    … …  
    299319    ], 
    300320    "fractal_core_shell": [ 
    301         "FractalCoreShellModel", 
     321        "FractalCoreShell", 
    302322        { 
    303323            "sld_core": "core_sld", 
    … …  
    309329            "cor_length": "cor_length", 
    310330            "volfraction": "volfraction", 
    311         } 
     331        }, 
     332        "FractalCoreShellModel" 
    312333    ], 
    313334    "fuzzy_sphere": [ 
    … …  
    321342    ], 
    322343    "gauss_lorentz_gel": [ 
    323         "GaussLorentzGelModel", 
     344        "GaussLorentzGel", 
    324345        { 
    325346            "gauss_scale": "scale_g", 
    … …  
    328349            "background": "background", 
    329350            "lorentz_scale": "scale_l" 
    330         } 
     351        }, 
     352        "GaussLorentzGelModel" 
    331353    ], 
    332354    "gaussian_peak": [ 
    333         "PeakGaussModel", 
     355        "Peak Gauss Model", 
    334356        { 
    335357            "peak_pos": "q0", 
    336358            "sigma": "B", 
    337         } 
     359        }, 
     360        "PeakGaussModel", 
    338361    ], 
    339362    "gel_fit": [ 
    … …  
    343366            "lorentz_scale": "lScale", 
    344367            "guinier_scale": "gScale", 
    345             "fractal_dim": "scale", 
     368            "fractal_dim": "FractalExp", 
    346369            "cor_length": "zeta", 
    347370        } 
    348371    ], 
    349372    "guinier": [ 
     373        "Guinier", 
     374        { 
     375            "rg": "rg" 
     376        }, 
    350377        "GuinierModel", 
    351         { 
    352             "rg": "rg" 
    353         } 
    354378    ], 
    355379    "guinier_porod": [ 
    356         "GuinierPorodModel", 
     380        "GuinierPorod", 
    357381        { 
    358382            "s": "dim", 
    … …  
    361385            "scale": "scale", 
    362386            "background": "background" 
    363         } 
     387        }, 
     388        "GuinierPorodModel", 
    364389    ], 
    365390    "hardsphere": [ 
    … …  
    454479            "d_spacing": "spacing", 
    455480            "Caille_parameter": "caille", 
    456             "Nlayers": "N_plates", 
     481            "Nlayers": "n_plates", 
    457482        } 
    458483    ], 
    … …  
    486511    ], 
    487512    "lorentz": [ 
     513        "Lorentz", 
     514        { 
     515            "cor_length": "length" 
     516        }, 
    488517        "LorentzModel", 
    489         { 
    490             "cor_length": "length" 
    491         } 
    492518    ], 
    493519    "mass_fractal": [ 
    … …  
    510536    ], 
    511537    "mono_gauss_coil": [ 
    512         "DebyeModel", 
     538        "Debye", 
    513539        { 
    514540            "rg": "rg", 
    515541            "i_zero": "scale", 
    516542            "background": "background", 
    517         } 
     543        }, 
     544        "DebyeModel", 
    518545    ], 
    519546    "multilayer_vesicle": [ 
    … …  
    564591    ], 
    565592    "peak_lorentz": [ 
    566         "PeakLorentzModel", 
     593        "Peak Lorentz Model", 
    567594        { 
    568595            "peak_pos": "q0", 
    569596            "peak_hwhm": "B" 
    570         } 
     597        }, 
     598        "PeakLorentzModel", 
    571599    ], 
    572600    "pearl_necklace": [ 
    … …  
    802830    ], 
    803831    "two_lorentzian": [ 
    804         "TwoLorentzianModel", 
     832        "TwoLorentzian", 
    805833        { 
    806834            "lorentz_scale_1": "scale_1", 
    … …  
    811839            "lorentz_length_1": "length_1", 
    812840            "background": "background" 
    813         } 
     841        }, 
     842        "TwoLorentzianModel", 
    814843    ], 
    815844    "two_power_law": [ 
    816         "TwoPowerLawModel", 
     845        "TwoPowerLaw", 
    817846        { 
    818847            "coefficent_1": "coef_A", 
    … …  
    821850            "background": "background", 
    822851            "crossover": "qc" 
    823         } 
     852        }, 
     853        "TwoPowerLawModel", 
    824854    ], 
    825855    "unified_power_Rg": [ 
     856        "UnifiedPowerRg", 
     857        dict(((field_new+str(index), field_old+str(index)) 
     858              for field_new, field_old in [("rg", "Rg"), 
     859                                           ("power", "power"), 
     860                                           ("G", "G"), 
     861                                           ("B", "B"),] 
     862              for index in range(11)), 
     863             **{ 
     864                   "background": "background", 
     865                   "scale": "scale", 
     866               }), 
    826867        "UnifiedPowerRgModel", 
    827         { 
    828         } 
    829868    ], 
    830869    "vesicle": [ 
    … …  
    835874        } 
    836875    ] 
     876    } 
    837877} 

  • sasmodels/convert.py

    rf80f334 r07c8d46  
    5353    ("_pd_nsigma", ".nsigmas"), 
    5454    ("_pd_type", ".type"), 
     55    (".lower", ".lower"), 
     56    (".upper", ".upper"), 
     57    (".fittable", ".fittable"), 
     58    (".std", ".std"), 
     59    (".units", ".units"), 
     60    ("", "") 
    5561    ] 
    5662 
    … …  
    6470    if id.startswith('M0:'): 
    6571        return True 
    66     if id.startswith('volfraction') or id.startswith('radius_effective'): 
    67         return False 
    6872    if '_pd' in id or '.' in id: 
    6973        return False 
    … …  
    7579        if p.id == id: 
    7680            return p.type == 'sld' 
    77     raise ValueError("unknown parameter %r in conversion"%id) 
     81    return False 
    7882 
    7983def _rescale_sld(model_info, pars, scale): 
    … …  
    8892 
    8993 
    90 def _get_translation_table(model_info): 
    91     _, translation = CONVERSION_TABLE.get(model_info.id, [None, {}]) 
    92     translation = translation.copy() 
     94def _get_translation_table(model_info, version=(3,1,2)): 
     95    conv_param = CONVERSION_TABLE.get(version, {}).get(model_info.id, [None, {}]) 
     96    translation = conv_param[1].copy() 
    9397    for p in model_info.parameters.kernel_parameters: 
    9498        if p.length > 1: 
    … …  
    119123def _pd_to_underscores(pars): 
    120124    return dict((_dot_pd_to_underscore_pd(k), v) for k, v in pars.items()) 
    121  
    122 def _convert_name(conv_dict, pars): 
    123     """ 
    124     Renames parameter values (upper, lower, etc) to v4.0 names 
    125     :param conv_dict: conversion dictionary mapping new name : old name 
    126     :param pars: parameters to convert 
    127     :return: 
    128     """ 
    129     new_pars = {} 
    130     i = 0 
    131     j = 0 
    132     for key_par, value_par in pars.iteritems(): 
    133         j += 1 
    134         for key_conv, value_conv in conv_dict.iteritems(): 
    135             if re.search(value_conv, key_par): 
    136                 new_pars[key_par.replace(value_conv, key_conv)] = value_par 
    137                 i += 1 
    138                 break 
    139             elif re.search("background", key_par) or re.search("scale", key_par): 
    140                 new_pars[key_par] = value_par 
    141                 i += 1 
    142                 break 
    143         if i != j: 
    144             new_pars[key_par] = value_par 
    145             i += 1 
    146     return new_pars 
    147125 
    148126def _convert_pars(pars, mapping): 
    … …  
    167145    return newpars 
    168146 
    169  
    170 def _conversion_target(model_name): 
     147def _conversion_target(model_name, version=(3,1,2)): 
    171148    """ 
    172149    Find the sasmodel name which translates into the sasview name. 
    … …  
    176153    two variants in sasview. 
    177154    """ 
    178     for sasmodels_name, [sasview_name, _] in CONVERSION_TABLE.items(): 
    179         if sasview_name == model_name: 
     155    for sasmodels_name, sasview_dict in \ 
     156            CONVERSION_TABLE.get(version, {}).items(): 
     157        if sasview_dict[0] == model_name: 
    180158            return sasmodels_name 
    181159    return None 
    182160 
    183  
    184 def _hand_convert(name, oldpars): 
     161def _hand_convert(name, oldpars, version=(3,1,2)): 
     162    if version == (3,1,2): 
     163        oldpars = _hand_convert_3_1_2_to_4_1(name, oldpars) 
     164    return oldpars 
     165 
     166def _hand_convert_3_1_2_to_4_1(name, oldpars): 
    185167    if name == 'core_shell_parallelepiped': 
    186168        # Make sure pd on rim parameters defaults to zero 
    … …  
    215197            pd = oldpars['radius.width']*oldpars['radius']/thickness 
    216198            oldpars['radius.width'] = pd 
     199    elif name == 'multilayer_vesicle': 
     200        if 'scale' in oldpars: 
     201            oldpars['volfraction'] = oldpars['scale'] 
     202            oldpars['scale'] = 1.0 
     203        if 'scale.lower' in oldpars: 
     204            oldpars['volfraction.lower'] = oldpars['scale.lower'] 
     205        if 'scale.upper' in oldpars: 
     206            oldpars['volfraction.upper'] = oldpars['scale.upper'] 
     207        if 'scale.fittable' in oldpars: 
     208            oldpars['volfraction.fittable'] = oldpars['scale.fittable'] 
     209        if 'scale.std' in oldpars: 
     210            oldpars['volfraction.std'] = oldpars['scale.std'] 
     211        if 'scale.units' in oldpars: 
     212            oldpars['volfraction.units'] = oldpars['scale.units'] 
    217213    elif name == 'pearl_necklace': 
    218214        pass 
    … …  
    236232                oldpars[p + ".upper"] /= 1e-13 
    237233    elif name == 'spherical_sld': 
    238         oldpars["CONTROL"] += 1 
     234        j = 0 
     235        while "func_inter" + str(j) in oldpars: 
     236            name = "func_inter" + str(j) 
     237            new_name = "shape" + str(j + 1) 
     238            if oldpars[name] == 'Erf(|nu|*z)': 
     239                oldpars[new_name] = int(0) 
     240            elif oldpars[name] == 'RPower(z^|nu|)': 
     241                oldpars[new_name] = int(1) 
     242            elif oldpars[name] == 'LPower(z^|nu|)': 
     243                oldpars[new_name] = int(2) 
     244            elif oldpars[name] == 'RExp(-|nu|*z)': 
     245                oldpars[new_name] = int(3) 
     246            elif oldpars[name] == 'LExp(-|nu|*z)': 
     247                oldpars[new_name] = int(4) 
     248            else: 
     249                oldpars[new_name] = int(0) 
     250            oldpars.pop(name) 
     251            oldpars['n_shells'] = str(j + 1) 
     252            j += 1 
    239253    elif name == 'teubner_strey': 
    240254        # basically undoing the entire Teubner-Strey calculations here. 
    … …  
    281295    return oldpars 
    282296 
    283 def convert_model(name, pars, use_underscore=False): 
     297def convert_model(name, pars, use_underscore=False, model_version=(3,1,2)): 
    284298    """ 
    285299    Convert model from old style parameter names to new style. 
    286300    """ 
    287     newname = _conversion_target(name) 
    288     if newname is None: 
    289         return name, pars 
    290     if ':' in newname:   # core_shell_ellipsoid:1 
    291         model_info = load_model_info(newname[:-2]) 
    292         # Know that the table exists and isn't multiplicity so grab it directly 
    293         # Can't use _get_translation_table since that will return the 'bare' 
    294         # version. 
    295         translation = CONVERSION_TABLE[newname][1] 
    296     else: 
    297         model_info = load_model_info(newname) 
    298         translation = _get_translation_table(model_info) 
    299     newpars = _hand_convert(newname, pars.copy()) 
    300     newpars = _convert_name(translation, newpars) 
    301     newpars = _convert_pars(newpars, translation) 
    302     if not model_info.structure_factor: 
    303         newpars = _rescale_sld(model_info, newpars, 1e6) 
    304     newpars.setdefault('scale', 1.0) 
    305     newpars.setdefault('background', 0.0) 
    306     if use_underscore: 
    307         newpars = _pd_to_underscores(newpars) 
     301    newpars = pars 
     302    keys = sorted(CONVERSION_TABLE.keys()) 
     303    for i, version in enumerate(keys): 
     304        # Don't allow indices outside list 
     305        next_i = i + 1 
     306        if next_i == len(keys): 
     307            next_i = i 
     308        # If the save state is from a later version, skip the check 
     309        if model_version <= keys[next_i]: 
     310            newname = _conversion_target(name, version) 
     311        else: 
     312            newname = None 
     313        # If no conversion is found, move on 
     314        if newname is None: 
     315            newname = name 
     316            continue 
     317        if ':' in newname:   # core_shell_ellipsoid:1 
     318            model_info = load_model_info(newname[:-2]) 
     319            # Know the table exists and isn't multiplicity so grab it directly 
     320            # Can't use _get_translation_table since that will return the 'bare' 
     321            # version. 
     322            translation = CONVERSION_TABLE.get(version, {})[newname][1] 
     323        else: 
     324            model_info = load_model_info(newname) 
     325            translation = _get_translation_table(model_info, version) 
     326        newpars = _hand_convert(newname, newpars, version) 
     327        newpars = _convert_pars(newpars, translation) 
     328        # TODO: Still not convinced this is the best check 
     329        if not model_info.structure_factor and version == (3,1,2): 
     330            newpars = _rescale_sld(model_info, newpars, 1e6) 
     331        newpars.setdefault('scale', 1.0) 
     332        newpars.setdefault('background', 0.0) 
     333        if use_underscore: 
     334            newpars = _pd_to_underscores(newpars) 
     335        name = newname 
    308336    return newname, newpars 
    309  
    310337 
    311338# ========= BACKWARD CONVERSION sasmodels => sasview 3.x =========== 

  • sasmodels/kernelcl.py

    r7fcdc9f rd2b939c  
    465465        # architectures tested so far. 
    466466        if self.is_2d: 
    467             # Note: 17 rather than 15 because results is 2 elements 
    468             # longer than input. 
    469             width = ((self.nq+17)//16)*16 
     467            # Note: 16 rather than 15 because result is 1 longer than input. 
     468            width = ((self.nq+16)//16)*16 
    470469            self.q = np.empty((width, 2), dtype=dtype) 
    471470            self.q[:self.nq, 0] = q_vectors[0] 
    472471            self.q[:self.nq, 1] = q_vectors[1] 
    473472        else: 
    474             # Note: 33 rather than 31 because results is 2 elements 
    475             # longer than input. 
    476             width = ((self.nq+33)//32)*32 
     473            # Note: 32 rather than 31 because result is 1 longer than input. 
     474            width = ((self.nq+32)//32)*32 
    477475            self.q = np.empty(width, dtype=dtype) 
    478476            self.q[:self.nq] = q_vectors[0] 
    … …  
    524522        self.dim = '2d' if q_input.is_2d else '1d' 
    525523        # plus three for the normalization values 
    526         self.result = np.empty(q_input.nq+3, dtype) 
     524        self.result = np.empty(q_input.nq+1, dtype) 
    527525 
    528526        # Inputs and outputs for each kernel call 

  • sasmodels/models/core_shell_bicelle.py

    r4b541ac r3b9a526  
    156156            phi=0) 
    157157 
    158 qx, qy = 0.4 * cos(90), 0.5 * sin(0) 
     158#qx, qy = 0.4 * cos(pi/2.0), 0.5 * sin(0) 

  • sasmodels/models/core_shell_bicelle_elliptical.py

    r8e68ea0 r3b9a526  
    8383 
    8484    Definition of the angles for the oriented core_shell_bicelle_elliptical model. 
     85    Note that *theta* and *phi* are currently defined differently to those for the core_shell_bicelle model. 
    8586 
    8687 
    … …  
    149150            psi=0) 
    150151 
    151 qx, qy = 0.4 * cos(90), 0.5 * sin(0) 
     152#qx, qy = 0.4 * cos(pi/2.0), 0.5 * sin(0) 
    152153 
    153154tests = [ 

  • sasmodels/models/ellipsoid.c

    r925ad6e r130d4c7  
    44    double radius_polar, double radius_equatorial, double theta, double phi); 
    55 
    6 double _ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double sin_alpha); 
    7 double _ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double sin_alpha) 
     6static double 
     7_ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double cos_alpha) 
    88{ 
    99    double ratio = radius_polar/radius_equatorial; 
    10     // Given the following under the radical: 
    11     //     1 + sin^2(T) (v^2 - 1) 
    12     // we can expand to match the form given in Guinier (1955) 
    13     //     = (1 - sin^2(T)) + v^2 sin^2(T) = cos^2(T) + sin^2(T) 
     10    // Using ratio v = Rp/Re, we can expand the following to match the 
     11    // form given in Guinier (1955) 
     12    //     r = Re * sqrt(1 + cos^2(T) (v^2 - 1)) 
     13    //       = Re * sqrt( (1 - cos^2(T)) + v^2 cos^2(T) ) 
     14    //       = Re * sqrt( sin^2(T) + v^2 cos^2(T) ) 
     15    //       = sqrt( Re^2 sin^2(T) + Rp^2 cos^2(T) ) 
     16    // 
    1417    // Instead of using pythagoras we could pass in sin and cos; this may be 
    1518    // slightly better for 2D which has already computed it, but it introduces 
    … …  
    1720    // leave it as is. 
    1821    const double r = radius_equatorial 
    19                      * sqrt(1.0 + sin_alpha*sin_alpha*(ratio*ratio - 1.0)); 
     22                     * sqrt(1.0 + cos_alpha*cos_alpha*(ratio*ratio - 1.0)); 
    2023    const double f = sas_3j1x_x(q*r); 
    2124 
    … …  
    3942    double total = 0.0; 
    4043    for (int i=0;i<76;i++) { 
    41         //const double sin_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2; 
    42         const double sin_alpha = Gauss76Z[i]*zm + zb; 
    43         total += Gauss76Wt[i] * _ellipsoid_kernel(q, radius_polar, radius_equatorial, sin_alpha); 
     44        //const double cos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2; 
     45        const double cos_alpha = Gauss76Z[i]*zm + zb; 
     46        total += Gauss76Wt[i] * _ellipsoid_kernel(q, radius_polar, radius_equatorial, cos_alpha); 
    4447    } 
    4548    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    5962    double q, sin_alpha, cos_alpha; 
    6063    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    61     const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, sin_alpha); 
     64    const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, cos_alpha); 
    6265    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    6366 

  • sasmodels/models/guinier_porod.py

    ra807206 rcdcebf1  
    44and dimensionality of scattering objects, including asymmetric objects 
    55such as rods or platelets, and shapes intermediate between spheres 
    6 and rods or between rods and platelets. 
     6and rods or between rods and platelets, and overcomes some of the  
     7deficiencies of the (Beaucage) Unified_Power_Rg model (see Hammouda, 2010). 
    78 
    89Definition 
    … …  
    5960--------- 
    6061 
    61 A Guinier, G Fournet, *Small-Angle Scattering of X-Rays*, 
    62 John Wiley and Sons, New York, (1955) 
     62B Hammouda, *A new Guinier-Porod model, J. Appl. Cryst.*, (2010), 43, 716-719 
    6363 
    64 O Glatter, O Kratky, *Small-Angle X-Ray Scattering*, Academic Press (1982) 
    65 Check out Chapter 4 on Data Treatment, pages 155-156. 
     64B Hammouda, *Analysis of the Beaucage model, J. Appl. Cryst.*, (2010), 43, 1474-1478 
     65 
    6666""" 
    6767 

  • sasmodels/models/unified_power_Rg.py

    rb3f2a24 r66ca2a6  
    33---------- 
    44 
    5 The Beaucage model employs the empirical multiple level unified 
    6 Exponential/Power-law fit method developed by G. Beaucage. Four functions 
    7 are included so that 1, 2, 3, or 4 levels can be used. In addition a 0 level 
    8 has been added which simply calculates 
     5This model employs the empirical multiple level unified Exponential/Power-law  
     6fit method developed by Beaucage. Four functions are included so that 1, 2, 3,  
     7or 4 levels can be used. In addition a 0 level has been added which simply  
     8calculates 
    99 
    1010.. math:: 
    … …  
    1515many different types of particles, including fractal clusters, random coils 
    1616(Debye equation), ellipsoidal particles, etc. 
     17 
     18The model works best for mass fractal systems characterized by Porod exponents  
     19between 5/3 and 3. It should not be used for surface fractal systems. Hammouda  
     20(2010) has pointed out a deficiency in the way this model handles the  
     21transitioning between the Guinier and Porod regimes and which can create  
     22artefacts that appear as kinks in the fitted model function. 
     23 
     24Also see the Guinier_Porod model. 
    1725 
    1826The empirical fit function is: 
    … …  
    3038.. math:: 
    3139 
    32     q_i^* = \frac{q}{\operatorname{erf}^3(q R_{gi}/\sqrt{6}} 
     40    q_i^* = q \left[\operatorname{erf} 
     41            \left(\frac{q R_{gi}}{\sqrt{6}}\right) 
     42        \right]^{-3} 
    3343 
    3444 
    … …  
    5666 
    5767G Beaucage, *J. Appl. Cryst.*, 29 (1996) 134-146 
     68 
     69B Hammouda, *Analysis of the Beaucage model, J. Appl. Cryst.*, (2010), 43, 1474-1478 
    5870 
    5971""" 

  • sasmodels/sasview_model.py

    r64614ad ra38b065  
    1616import logging 
    1717from os.path import basename, splitext 
     18import thread 
    1819 
    1920import numpy as np  # type: ignore 
    … …  
    3940    pass 
    4041 
     42calculation_lock = thread.allocate_lock() 
     43 
    4144SUPPORT_OLD_STYLE_PLUGINS = True 
    4245 
    … …  
    5760    import sas.models 
    5861    from sasmodels.conversion_table import CONVERSION_TABLE 
    59     for new_name, conversion in CONVERSION_TABLE.items(): 
     62    for new_name, conversion in CONVERSION_TABLE.get((3,1,2), {}).items(): 
    6063        # CoreShellEllipsoidModel => core_shell_ellipsoid:1 
    6164        new_name = new_name.split(':')[0] 
    62         old_name = conversion[0] 
     65        old_name = conversion[0] if len(conversion) < 3 else conversion[2] 
    6366        module_attrs = {old_name: find_model(new_name)} 
    6467        ConstructedModule = type(old_name, (), module_attrs) 
    … …  
    605608        to the card for each evaluation. 
    606609        """ 
     610        ## uncomment the following when trying to debug the uncoordinated calls 
     611        ## to calculate_Iq 
     612        #if calculation_lock.locked(): 
     613        #    logging.info("calculation waiting for another thread to complete") 
     614        #    logging.info("\n".join(traceback.format_stack())) 
     615 
     616        with calculation_lock: 
     617            return self._calculate_Iq(qx, qy) 
     618 
     619    def _calculate_Iq(self, qx, qy=None): 
    607620        #core.HAVE_OPENCL = False 
    608621        if self._model is None: 

  • sasmodels/sesans.py

    rb397165 r94d13f1  
    1414import numpy as np  # type: ignore 
    1515from numpy import pi, exp  # type: ignore 
    16 from scipy.special import jv as besselj 
    17 #import direct_model.DataMixin as model 
    18          
    19 def make_q(q_max, Rmax): 
    20     r""" 
    21     Return a $q$ vector suitable for SESANS covering from $2\pi/ (10 R_{\max})$ 
    22     to $q_max$. This is the integration range of the Hankel transform; bigger range and  
    23     more points makes a better numerical integration. 
    24     Smaller q_min will increase reliable spin echo length range.  
    25     Rmax is the "radius" of the largest expected object and can be set elsewhere. 
    26     q_max is determined by the acceptance angle of the SESANS instrument. 
     16from scipy.special import j0 
     17 
     18class SesansTransform(object): 
    2719    """ 
    28     from sas.sascalc.data_util.nxsunit import Converter 
     20    Spin-Echo SANS transform calculator.  Similar to a resolution function, 
     21    the SesansTransform object takes I(q) for the set of *q_calc* values and 
     22    produces a transformed dataset 
    2923 
    30     q_min = dq = 0.1 * 2*pi / Rmax 
    31     return np.arange(q_min, 
    32                      Converter(q_max[1])(q_max[0], 
    33                                          units="1/A"), 
    34                      dq) 
    35      
    36 def make_all_q(data): 
     24    *SElength* (A) is the set of spin-echo lengths in the measured data. 
     25 
     26    *zaccept* (1/A) is the maximum acceptance of scattering vector in the spin 
     27    echo encoding dimension (for ToF: Q of min(R) and max(lam)). 
     28 
     29    *Rmax* (A) is the maximum size sensitivity; larger radius requires more 
     30    computation time. 
    3731    """ 
    38     Return a $q$ vector suitable for calculating the total scattering cross section for 
    39     calculating the effect of finite acceptance angles on Time of Flight SESANS instruments. 
    40     If no acceptance is given, or unwanted (set "unwanted" flag in paramfile), no all_q vector is needed. 
    41     If the instrument has a rectangular acceptance, 2 all_q vectors are needed. 
    42     If the instrument has a circular acceptance, 1 all_q vector is needed 
    43      
    44     """ 
    45     if not data.has_no_finite_acceptance: 
    46         return [] 
    47     elif data.has_yz_acceptance(data): 
    48         # compute qx, qy 
    49         Qx, Qy = np.meshgrid(qx, qy) 
    50         return [Qx, Qy] 
    51     else: 
    52         # else only need q 
    53         # data.has_z_acceptance 
    54         return [q] 
     32    #: SElength from the data in the original data units; not used by transform 
     33    #: but the GUI uses it, so make sure that it is present. 
     34    q = None  # type: np.ndarray 
    5535 
    56 def transform(data, q_calc, Iq_calc, qmono, Iq_mono): 
    57     """ 
    58     Decides which transform type is to be used, based on the experiment data file contents (header) 
    59     (2016-03-19: currently controlled from parameters script) 
    60     nqmono is the number of q vectors to be used for the detector integration 
    61     """ 
    62     nqmono = len(qmono) 
    63     if nqmono == 0: 
    64         result = call_hankel(data, q_calc, Iq_calc) 
    65     elif nqmono == 1: 
    66         q = qmono[0] 
    67         result = call_HankelAccept(data, q_calc, Iq_calc, q, Iq_mono) 
    68     else: 
    69         Qx, Qy = [qmono[0], qmono[1]] 
    70         Qx = np.reshape(Qx, nqx, nqy) 
    71         Qy = np.reshape(Qy, nqx, nqy) 
    72         Iq_mono = np.reshape(Iq_mono, nqx, nqy) 
    73         qx = Qx[0, :] 
    74         qy = Qy[:, 0] 
    75         result = call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono) 
     36    #: q values to calculate when computing transform 
     37    q_calc = None  # type: np.ndarray 
    7638 
    77     return result 
     39    # transform arrays 
     40    _H = None  # type: np.ndarray 
     41    _H0 = None # type: np.ndarray 
    7842 
    79 def call_hankel(data, q_calc, Iq_calc): 
    80     return hankel((data.x, data.x_unit), 
    81                   (data.lam, data.lam_unit), 
    82                   (data.sample.thickness, 
    83                    data.sample.thickness_unit), 
    84                   q_calc, Iq_calc) 
    85    
    86 def call_HankelAccept(data, q_calc, Iq_calc, q_mono, Iq_mono): 
    87     return hankel(data.x, data.lam * 1e-9, 
    88                   data.sample.thickness / 10, 
    89                   q_calc, Iq_calc) 
    90                    
    91 def call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono): 
    92     return hankel(data.x, data.y, data.lam * 1e-9, 
    93                   data.sample.thickness / 10, 
    94                   q_calc, Iq_calc) 
    95                          
    96 def TotalScatter(model, parameters):  #Work in progress!! 
    97 #    Calls a model with existing model parameters already in place, then integrate the product of q and I(q) from 0 to (4*pi/lambda) 
    98     allq = np.linspace(0,4*pi/wavelength,1000) 
    99     allIq = 1 
    100     integral = allq*allIq 
    101      
     43    def __init__(self, z, SElength, zaccept, Rmax): 
     44        # type: (np.ndarray, float, float) -> None 
     45        #import logging; logging.info("creating SESANS transform") 
     46        self.q = z 
     47        self._set_hankel(SElength, zaccept, Rmax) 
    10248 
     49    def apply(self, Iq): 
     50        # tye: (np.ndarray) -> np.ndarray 
     51        G0 = np.dot(self._H0, Iq) 
     52        G = np.dot(self._H.T, Iq) 
     53        P = G - G0 
     54        return P 
    10355 
    104 def Cosine2D(wavelength, magfield, thickness, qy, qz, Iqy, Iqz, modelname): #Work in progress!! Needs to call model still 
    105 #============================================================================== 
    106 #     2D Cosine Transform if "wavelength" is a vector 
    107 #============================================================================== 
    108 #allq is the q-space needed to create the total scattering cross-section 
     56    def _set_hankel(self, SElength, zaccept, Rmax): 
     57        # type: (np.ndarray, float, float) -> None 
     58        # Force float32 arrays, otherwise run into memory problems on some machines 
     59        SElength = np.asarray(SElength, dtype='float32') 
    10960 
    110     Gprime = np.zeros_like(wavelength, 'd') 
    111     s = np.zeros_like(wavelength, 'd') 
    112     sd = np.zeros_like(wavelength, 'd') 
    113     Gprime = np.zeros_like(wavelength, 'd') 
    114     f = np.zeros_like(wavelength, 'd') 
    115     for i, wavelength_i in enumerate(wavelength): 
    116         z = magfield*wavelength_i 
    117         allq=np.linspace() #for calculating the Q-range of the  scattering power integral 
    118         allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow 
    119         alldq = (allq[1]-allq[0])*1e10 
    120         sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) 
    121         s[i]=1-exp(-sigma) 
    122         for j, Iqy_j, qy_j in enumerate(qy): 
    123             for k, Iqz_k, qz_k in enumerate(qz): 
    124                 Iq = np.sqrt(Iqy_j^2+Iqz_k^2) 
    125                 q = np.sqrt(qy_j^2 + qz_k^2) 
    126                 Gintegral = Iq*cos(z*Qz_k) 
    127                 Gprime[i] += Gintegral 
    128 #                sigma = wavelength^2*thickness/2/pi* allq[i]*allIq[i] 
    129 #                s[i] += 1-exp(Totalscatter(modelname)*thickness) 
    130 #                For now, work with standard 2-phase scatter 
     61        #Rmax = #value in text box somewhere in FitPage? 
     62        q_max = 2*pi / (SElength[1] - SElength[0]) 
     63        q_min = 0.1 * 2*pi / (np.size(SElength) * SElength[-1]) 
     64        q = np.arange(q_min, q_max, q_min, dtype='float32') 
     65        dq = q_min 
    13166 
     67        H0 = np.float32(dq/(2*pi)) * q 
    13268 
    133                 sd[i] += Iq 
    134         f[i] = 1-s[i]+sd[i] 
    135         P[i] = (1-sd[i]/f[i])+1/f[i]*Gprime[i] 
     69        repq = np.tile(q, (SElength.size, 1)).T 
     70        repSE = np.tile(SElength, (q.size, 1)) 
     71        H = np.float32(dq/(2*pi)) * j0(repSE*repq) * repq 
    13672 
    137  
    138  
    139  
    140 def HankelAccept(wavelength, magfield, thickness, q, Iq, theta, modelname): 
    141 #============================================================================== 
    142 #     HankelTransform with fixed circular acceptance angle (circular aperture) for Time of Flight SESANS 
    143 #============================================================================== 
    144 #acceptq is the q-space needed to create limited acceptance effect 
    145     SElength= wavelength*magfield 
    146     G = np.zeros_like(SElength, 'd') 
    147     threshold=2*pi*theta/wavelength 
    148     for i, SElength_i in enumerate(SElength): 
    149         allq=np.linspace() #for calculating the Q-range of the  scattering power integral 
    150         allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow 
    151         alldq = (allq[1]-allq[0])*1e10 
    152         sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) 
    153         s[i]=1-exp(-sigma) 
    154  
    155         dq = (q[1]-q[0])*1e10 
    156         a = (x<threshold) 
    157         acceptq = a*q 
    158         acceptIq = a*Iq 
    159  
    160         G[i] = np.sum(besselj(0, acceptq*SElength_i)*acceptIq*acceptq*dq) 
    161  
    162 #        G[i]=np.sum(integral) 
    163  
    164     G *= dq*1e10*2*pi 
    165  
    166     P = exp(thickness*wavelength**2/(4*pi**2)*(G-G[0])) 
    167      
    168 def hankel(SElength, wavelength, thickness, q, Iq): 
    169     r""" 
    170     Compute the expected SESANS polarization for a given SANS pattern. 
    171  
    172     Uses the hankel transform followed by the exponential.  The values for *zz* 
    173     (or spin echo length, or delta), wavelength and sample thickness should 
    174     come from the dataset.  $q$ should be chosen such that the oscillations 
    175     in $I(q)$ are well sampled (e.g., $5 \cdot 2 \pi/d_{\max}$). 
    176  
    177     *SElength* [A] is the set of $z$ points at which to compute the 
    178     Hankel transform 
    179  
    180     *wavelength* [m]  is the wavelength of each individual point *zz* 
    181  
    182     *thickness* [cm] is the sample thickness. 
    183  
    184     *q* [A$^{-1}$] is the set of $q$ points at which the model has been 
    185     computed. These should be equally spaced. 
    186  
    187     *I* [cm$^{-1}$] is the value of the SANS model at *q* 
    188     """ 
    189  
    190     from sas.sascalc.data_util.nxsunit import Converter 
    191     wavelength = Converter(wavelength[1])(wavelength[0],"A") 
    192     thickness = Converter(thickness[1])(thickness[0],"A") 
    193     Iq = Converter("1/cm")(Iq,"1/A") # All models default to inverse centimeters 
    194     SElength = Converter(SElength[1])(SElength[0],"A") 
    195  
    196     G = np.zeros_like(SElength, 'd') 
    197 #============================================================================== 
    198 #     Hankel Transform method if "wavelength" is a scalar; mono-chromatic SESANS 
    199 #============================================================================== 
    200     for i, SElength_i in enumerate(SElength): 
    201         integral = besselj(0, q*SElength_i)*Iq*q 
    202         G[i] = np.sum(integral) 
    203     G0 = np.sum(Iq*q) 
    204  
    205     # [m^-1] step size in q, needed for integration 
    206     dq = (q[1]-q[0]) 
    207  
    208     # integration step, convert q into [m**-1] and 2 pi circle integration 
    209     G *= dq*2*pi 
    210     G0 = np.sum(Iq*q)*dq*2*np.pi 
    211  
    212     P = exp(thickness*wavelength**2/(4*pi**2)*(G-G0)) 
    213  
    214     return P 
     73        self.q_calc = q 
     74        self._H, self._H0 = H, H0 
Note: See TracChangeset for help on using the changeset viewer.