Changes in / [0c2da4b:d3e3f756] in sasmodels


Ignore:
Files:
1 deleted
9 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/conversion_table.py

    rd3e3f756 rd3e3f756  
    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 previous versions of SasView and sasmodels. This is now 
    8 versioned based on the version number of SasView. 
    9  
    10 When any sasmodels parameter or model name is changed, this must be modified to 
    11 account for that. 
    12  
    13 Usage: 
    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  
    25 Any future parameter and model name changes can and should be given in this 
    26 table for future compatibility. 
     7the models defined in SasView 3.1. 
    278""" 
    289 
     10 
    2911CONVERSION_TABLE = { 
    30     (3,1,2) : { 
    3112    "adsorbed_layer": [ 
    3213        "Core2ndMomentModel", 
    … …  
    230211    ], 
    231212    "correlation_length": [ 
    232         "CorrLength", 
     213        "CorrLengthModel", 
    233214        { 
    234215            "porod_scale": "scale_p", 
    … …  
    237218            "lorentz_exp": "exponent_l", 
    238219            "cor_length": "length_l" 
    239         }, 
    240         "CorrLengthModel" 
     220        } 
    241221    ], 
    242222    "cylinder": [ 
    … …  
    319299    ], 
    320300    "fractal_core_shell": [ 
    321         "FractalCoreShell", 
     301        "FractalCoreShellModel", 
    322302        { 
    323303            "sld_core": "core_sld", 
    … …  
    329309            "cor_length": "cor_length", 
    330310            "volfraction": "volfraction", 
    331         }, 
    332         "FractalCoreShellModel" 
     311        } 
    333312    ], 
    334313    "fuzzy_sphere": [ 
    … …  
    342321    ], 
    343322    "gauss_lorentz_gel": [ 
    344         "GaussLorentzGel", 
     323        "GaussLorentzGelModel", 
    345324        { 
    346325            "gauss_scale": "scale_g", 
    … …  
    349328            "background": "background", 
    350329            "lorentz_scale": "scale_l" 
    351         }, 
    352         "GaussLorentzGelModel" 
     330        } 
    353331    ], 
    354332    "gaussian_peak": [ 
    355         "Peak Gauss Model", 
     333        "PeakGaussModel", 
    356334        { 
    357335            "peak_pos": "q0", 
    358336            "sigma": "B", 
    359         }, 
    360         "PeakGaussModel", 
     337        } 
    361338    ], 
    362339    "gel_fit": [ 
    … …  
    366343            "lorentz_scale": "lScale", 
    367344            "guinier_scale": "gScale", 
    368             "fractal_dim": "FractalExp", 
     345            "fractal_dim": "scale", 
    369346            "cor_length": "zeta", 
    370347        } 
    371348    ], 
    372349    "guinier": [ 
    373         "Guinier", 
     350        "GuinierModel", 
    374351        { 
    375352            "rg": "rg" 
    376         }, 
    377         "GuinierModel", 
     353        } 
    378354    ], 
    379355    "guinier_porod": [ 
    380         "GuinierPorod", 
     356        "GuinierPorodModel", 
    381357        { 
    382358            "s": "dim", 
    … …  
    385361            "scale": "scale", 
    386362            "background": "background" 
    387         }, 
    388         "GuinierPorodModel", 
     363        } 
    389364    ], 
    390365    "hardsphere": [ 
    … …  
    479454            "d_spacing": "spacing", 
    480455            "Caille_parameter": "caille", 
    481             "Nlayers": "n_plates", 
     456            "Nlayers": "N_plates", 
    482457        } 
    483458    ], 
    … …  
    511486    ], 
    512487    "lorentz": [ 
    513         "Lorentz", 
     488        "LorentzModel", 
    514489        { 
    515490            "cor_length": "length" 
    516         }, 
    517         "LorentzModel", 
     491        } 
    518492    ], 
    519493    "mass_fractal": [ 
    … …  
    536510    ], 
    537511    "mono_gauss_coil": [ 
    538         "Debye", 
     512        "DebyeModel", 
    539513        { 
    540514            "rg": "rg", 
    541515            "i_zero": "scale", 
    542516            "background": "background", 
    543         }, 
    544         "DebyeModel", 
     517        } 
    545518    ], 
    546519    "multilayer_vesicle": [ 
    … …  
    591564    ], 
    592565    "peak_lorentz": [ 
    593         "Peak Lorentz Model", 
     566        "PeakLorentzModel", 
    594567        { 
    595568            "peak_pos": "q0", 
    596569            "peak_hwhm": "B" 
    597         }, 
    598         "PeakLorentzModel", 
     570        } 
    599571    ], 
    600572    "pearl_necklace": [ 
    … …  
    830802    ], 
    831803    "two_lorentzian": [ 
    832         "TwoLorentzian", 
     804        "TwoLorentzianModel", 
    833805        { 
    834806            "lorentz_scale_1": "scale_1", 
    … …  
    839811            "lorentz_length_1": "length_1", 
    840812            "background": "background" 
    841         }, 
    842         "TwoLorentzianModel", 
     813        } 
    843814    ], 
    844815    "two_power_law": [ 
    845         "TwoPowerLaw", 
     816        "TwoPowerLawModel", 
    846817        { 
    847818            "coefficent_1": "coef_A", 
    … …  
    850821            "background": "background", 
    851822            "crossover": "qc" 
    852         }, 
    853         "TwoPowerLawModel", 
     823        } 
    854824    ], 
    855825    "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                }), 
    867826        "UnifiedPowerRgModel", 
     827        { 
     828        } 
    868829    ], 
    869830    "vesicle": [ 
    … …  
    874835        } 
    875836    ] 
    876     } 
    877837} 

  • sasmodels/convert.py

    r07c8d46 rf80f334  
    5353    ("_pd_nsigma", ".nsigmas"), 
    5454    ("_pd_type", ".type"), 
    55     (".lower", ".lower"), 
    56     (".upper", ".upper"), 
    57     (".fittable", ".fittable"), 
    58     (".std", ".std"), 
    59     (".units", ".units"), 
    60     ("", "") 
    6155    ] 
    6256 
    … …  
    7064    if id.startswith('M0:'): 
    7165        return True 
     66    if id.startswith('volfraction') or id.startswith('radius_effective'): 
     67        return False 
    7268    if '_pd' in id or '.' in id: 
    7369        return False 
    … …  
    7975        if p.id == id: 
    8076            return p.type == 'sld' 
    81     return False 
     77    raise ValueError("unknown parameter %r in conversion"%id) 
    8278 
    8379def _rescale_sld(model_info, pars, scale): 
    … …  
    9288 
    9389 
    94 def _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() 
     90def _get_translation_table(model_info): 
     91    _, translation = CONVERSION_TABLE.get(model_info.id, [None, {}]) 
     92    translation = translation.copy() 
    9793    for p in model_info.parameters.kernel_parameters: 
    9894        if p.length > 1: 
    … …  
    123119def _pd_to_underscores(pars): 
    124120    return dict((_dot_pd_to_underscore_pd(k), v) for k, v in pars.items()) 
     121 
     122def _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 
    125147 
    126148def _convert_pars(pars, mapping): 
    … …  
    145167    return newpars 
    146168 
    147 def _conversion_target(model_name, version=(3,1,2)): 
     169 
     170def _conversion_target(model_name): 
    148171    """ 
    149172    Find the sasmodel name which translates into the sasview name. 
    … …  
    153176    two variants in sasview. 
    154177    """ 
    155     for sasmodels_name, sasview_dict in \ 
    156             CONVERSION_TABLE.get(version, {}).items(): 
    157         if sasview_dict[0] == model_name: 
     178    for sasmodels_name, [sasview_name, _] in CONVERSION_TABLE.items(): 
     179        if sasview_name == model_name: 
    158180            return sasmodels_name 
    159181    return None 
    160182 
    161 def _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  
    166 def _hand_convert_3_1_2_to_4_1(name, oldpars): 
     183 
     184def _hand_convert(name, oldpars): 
    167185    if name == 'core_shell_parallelepiped': 
    168186        # Make sure pd on rim parameters defaults to zero 
    … …  
    197215            pd = oldpars['radius.width']*oldpars['radius']/thickness 
    198216            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'] 
    213217    elif name == 'pearl_necklace': 
    214218        pass 
    … …  
    232236                oldpars[p + ".upper"] /= 1e-13 
    233237    elif name == 'spherical_sld': 
    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 
     238        oldpars["CONTROL"] += 1 
    253239    elif name == 'teubner_strey': 
    254240        # basically undoing the entire Teubner-Strey calculations here. 
    … …  
    295281    return oldpars 
    296282 
    297 def convert_model(name, pars, use_underscore=False, model_version=(3,1,2)): 
     283def convert_model(name, pars, use_underscore=False): 
    298284    """ 
    299285    Convert model from old style parameter names to new style. 
    300286    """ 
    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 
     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) 
    336308    return newname, newpars 
     309 
    337310 
    338311# ========= BACKWARD CONVERSION sasmodels => sasview 3.x =========== 

  • sasmodels/kernelcl.py

    rb8ddf2e r7fcdc9f  
    465465        # architectures tested so far. 
    466466        if self.is_2d: 
    467             # Note: 16 rather than 15 because result is 1 longer than input. 
    468             width = ((self.nq+16)//16)*16 
     467            # Note: 17 rather than 15 because results is 2 elements 
     468            # longer than input. 
     469            width = ((self.nq+17)//16)*16 
    469470            self.q = np.empty((width, 2), dtype=dtype) 
    470471            self.q[:self.nq, 0] = q_vectors[0] 
    471472            self.q[:self.nq, 1] = q_vectors[1] 
    472473        else: 
    473             # Note: 32 rather than 31 because result is 1 longer than input. 
    474             width = ((self.nq+32)//32)*32 
     474            # Note: 33 rather than 31 because results is 2 elements 
     475            # longer than input. 
     476            width = ((self.nq+33)//32)*32 
    475477            self.q = np.empty(width, dtype=dtype) 
    476478            self.q[:self.nq] = q_vectors[0] 
    … …  
    522524        self.dim = '2d' if q_input.is_2d else '1d' 
    523525        # plus three for the normalization values 
    524         self.result = np.empty(q_input.nq+1, dtype) 
     526        self.result = np.empty(q_input.nq+3, dtype) 
    525527 
    526528        # Inputs and outputs for each kernel call 

  • sasmodels/models/ellipsoid.c

    r130d4c7 r925ad6e  
    44    double radius_polar, double radius_equatorial, double theta, double phi); 
    55 
    6 static double 
    7 _ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double cos_alpha) 
     6double _ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double sin_alpha); 
     7double _ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double sin_alpha) 
    88{ 
    99    double ratio = radius_polar/radius_equatorial; 
    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     // 
     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) 
    1714    // Instead of using pythagoras we could pass in sin and cos; this may be 
    1815    // slightly better for 2D which has already computed it, but it introduces 
    … …  
    2017    // leave it as is. 
    2118    const double r = radius_equatorial 
    22                      * sqrt(1.0 + cos_alpha*cos_alpha*(ratio*ratio - 1.0)); 
     19                     * sqrt(1.0 + sin_alpha*sin_alpha*(ratio*ratio - 1.0)); 
    2320    const double f = sas_3j1x_x(q*r); 
    2421 
    … …  
    4239    double total = 0.0; 
    4340    for (int i=0;i<76;i++) { 
    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); 
     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); 
    4744    } 
    4845    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    6259    double q, sin_alpha, cos_alpha; 
    6360    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 
    64     const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, cos_alpha); 
     61    const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, sin_alpha); 
    6562    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
    6663 

  • sasmodels/models/hayter_msa.c

    r3fc5d27 r4962519  
    7070                SofQ=sqhcal(Qdiam, gMSAWave); 
    7171        }else{ 
    72         SofQ=NAN; 
     72        //SofQ=NaN; 
     73                SofQ=-1.0; 
    7374                //      print "Error Level = ",ierr 
    7475                //      print "Please report HPMSA problem with above error code" 

  • sasmodels/models/rpa.c

    racfb094 r6351bfa  
    3939  double S31,S32,S33,S34,S41,S42,S43,S44; 
    4040  double Lad,Lbd,Lcd,Nav,Intg; 
    41    
    42   // Set values for non existent parameters (eg. no A or B in case 0 and 1 etc) 
     41 
    4342  //icase was shifted to N-1 from the original code 
    4443  if (icase <= 1){ 
    … …  
    5857    Kab = Kac = Kad = -0.0004; 
    5958  } 
    60   
    61   // Set volume fraction of component D based on constraint that sum of vol frac =1 
    62   Phi[3]=1.0-Phi[0]-Phi[1]-Phi[2]; 
    63  
    64   //set up values for cross terms in case of block copolymers (1,3,4,6,7,8,9) 
     59 
    6560  Nab=sqrt(N[0]*N[1]); 
    6661  Nac=sqrt(N[0]*N[2]); 
    … …  
    8479  Phicd=sqrt(Phi[2]*Phi[3]); 
    8580 
    86   // Calculate Q^2 * Rg^2 for each homopolymer assuming random walk  
    8781  Xa=q*q*b[0]*b[0]*N[0]/6.0; 
    8882  Xb=q*q*b[1]*b[1]*N[1]/6.0; 
    … …  
    9084  Xd=q*q*b[3]*b[3]*N[3]/6.0; 
    9185 
    92   //calculate all partial structure factors Pij and normalize n^2 
    93   Paa=2.0*(exp(-Xa)-1.0+Xa)/(Xa*Xa); // free A chain form factor 
    94   S0aa=N[0]*Phi[0]*v[0]*Paa; // Phi * Vp * P(Q)= I(Q0)/delRho^2 
    95   Pab=((1.0-exp(-Xa))/Xa)*((1.0-exp(-Xb))/Xb); //AB diblock (anchored Paa * anchored Pbb) partial form factor 
     86  Paa=2.0*(exp(-Xa)-1.0+Xa)/(Xa*Xa); 
     87  S0aa=N[0]*Phi[0]*v[0]*Paa; 
     88  Pab=((1.0-exp(-Xa))/Xa)*((1.0-exp(-Xb))/Xb); 
    9689  S0ab=(Phiab*vab*Nab)*Pab; 
    97   Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc); //ABC triblock AC partial form factor 
     90  Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc); 
    9891  S0ac=(Phiac*vac*Nac)*Pac; 
    99   Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd); //ABCD four block 
     92  Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd); 
    10093  S0ad=(Phiad*vad*Nad)*Pad; 
    10194 
    10295  S0ba=S0ab; 
    103   Pbb=2.0*(exp(-Xb)-1.0+Xb)/(Xb*Xb); // free B chain 
     96  Pbb=2.0*(exp(-Xb)-1.0+Xb)/(Xb*Xb); 
    10497  S0bb=N[1]*Phi[1]*v[1]*Pbb; 
    105   Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc); // BC diblock 
     98  Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc); 
    10699  S0bc=(Phibc*vbc*Nbc)*Pbc; 
    107   Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd); // BCD triblock 
     100  Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd); 
    108101  S0bd=(Phibd*vbd*Nbd)*Pbd; 
    109102 
    110103  S0ca=S0ac; 
    111104  S0cb=S0bc; 
    112   Pcc=2.0*(exp(-Xc)-1.0+Xc)/(Xc*Xc); // Free C chain 
     105  Pcc=2.0*(exp(-Xc)-1.0+Xc)/(Xc*Xc); 
    113106  S0cc=N[2]*Phi[2]*v[2]*Pcc; 
    114   Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd); // CD diblock 
     107  Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd); 
    115108  S0cd=(Phicd*vcd*Ncd)*Pcd; 
    116109 
    … …  
    118111  S0db=S0bd; 
    119112  S0dc=S0cd; 
    120   Pdd=2.0*(exp(-Xd)-1.0+Xd)/(Xd*Xd); // free D chain 
     113  Pdd=2.0*(exp(-Xd)-1.0+Xd)/(Xd*Xd); 
    121114  S0dd=N[3]*Phi[3]*v[3]*Pdd; 
    122  
    123   // Reset all unused partial structure factors to 0 (depends on case) 
    124115  //icase was shifted to N-1 from the original code 
    125116  switch(icase){ 
    … …  
    202193  S0dc=S0cd; 
    203194 
    204   // self chi parameter is 0 ... of course 
    205195  Kaa=0.0; 
    206196  Kbb=0.0; 
    … …  
    253243  ZZ=S0ad*(T11*S0ad+T12*S0bd+T13*S0cd)+S0bd*(T21*S0ad+T22*S0bd+T23*S0cd)+S0cd*(T31*S0ad+T32*S0bd+T33*S0cd); 
    254244 
    255   // D is considered the matrix or background component so enters here 
    256245  m=1.0/(S0dd-ZZ); 
    257246 
    … …  
    308297  S44=S11+S22+S33+2.0*S12+2.0*S13+2.0*S23; 
    309298 
    310   //calculate contrast where L[i] is the scattering length of i and D is the matrix 
    311   //need to verify why the sqrt of Nav rather than just Nav (assuming v is molar volume) 
    312299  Nav=6.022045e+23; 
    313300  Lad=(L[0]/v[0]-L[3]/v[3])*sqrt(Nav); 
    … …  
    316303 
    317304  Intg=Lad*Lad*S11+Lbd*Lbd*S22+Lcd*Lcd*S33+2.0*Lad*Lbd*S12+2.0*Lbd*Lcd*S23+2.0*Lad*Lcd*S13; 
    318  
    319   //rescale for units of Lij^2 (in 10e-12 m^2 to m^2 ?) 
    320   Intg *= 1.0e-24;     
    321305 
    322306  return Intg; 

  • sasmodels/models/rpa.py

    rbb73096 r40a87fa  
    107107 
    108108control = "case_num" 
    109 HIDE_ALL = set("Phi4".split()) 
    110 HIDE_A = set("N1 Phi1 v1 L1 b1 K12 K13 K14".split()).union(HIDE_ALL) 
     109HIDE_NONE = set() 
     110HIDE_A = set("N1 Phi1 v1 L1 b1 K12 K13 K14".split()) 
    111111HIDE_AB = set("N2 Phi2 v2 L2 b2 K23 K24".split()).union(HIDE_A) 
    112112def hidden(case_num): 
    … …  
    119119        return HIDE_A 
    120120    else: 
    121         return HIDE_ALL 
     121        return HIDE_NONE 
    122122 

  • sasmodels/sasview_model.py

    r749a7d4 r64614ad  
    1616import logging 
    1717from os.path import basename, splitext 
    18 import thread 
    1918 
    2019import numpy as np  # type: ignore 
    … …  
    4039    pass 
    4140 
    42 calculation_lock = thread.allocate_lock() 
    43  
    4441SUPPORT_OLD_STYLE_PLUGINS = True 
    4542 
    … …  
    6057    import sas.models 
    6158    from sasmodels.conversion_table import CONVERSION_TABLE 
    62     for new_name, conversion in CONVERSION_TABLE.get((3,1,2), {}).items(): 
     59    for new_name, conversion in CONVERSION_TABLE.items(): 
    6360        # CoreShellEllipsoidModel => core_shell_ellipsoid:1 
    6461        new_name = new_name.split(':')[0] 
    65         old_name = conversion[0] if len(conversion) < 3 else conversion[2] 
     62        old_name = conversion[0] 
    6663        module_attrs = {old_name: find_model(new_name)} 
    6764        ConstructedModule = type(old_name, (), module_attrs) 
    … …  
    608605        to the card for each evaluation. 
    609606        """ 
    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): 
    620607        #core.HAVE_OPENCL = False 
    621608        if self._model is None: 
    … …  
    734721            return [self.params[par.name]], [1.0] 
    735722 
    736 def test_cylinder(): 
     723def test_model(): 
    737724    # type: () -> float 
    738725    """ 
    739     Test that the cylinder model runs, returning the value at [0.1,0.1]. 
     726    Test that a sasview model (cylinder) can be run. 
    740727    """ 
    741728    Cylinder = _make_standard_model('cylinder') 
    … …  
    746733    # type: () -> float 
    747734    """ 
    748     Test that 2-D hardsphere model runs and doesn't produce NaN. 
     735    Test that a sasview model (cylinder) can be run. 
    749736    """ 
    750737    Model = _make_standard_model('hardsphere') 
    … …  
    757744    # type: () -> float 
    758745    """ 
    759     Test that the 2-D RPA model runs 
     746    Test that a sasview model (cylinder) can be run. 
    760747    """ 
    761748    RPA = _make_standard_model('rpa') 
    … …  
    763750    return rpa.evalDistribution([0.1, 0.1]) 
    764751 
    765 def test_empty_distribution(): 
    766     # type: () -> None 
    767     """ 
    768     Make sure that sasmodels returns NaN when there are no polydispersity points 
    769     """ 
    770     Cylinder = _make_standard_model('cylinder') 
    771     cylinder = Cylinder() 
    772     cylinder.setParam('radius', -1.0) 
    773     cylinder.setParam('background', 0.) 
    774     Iq = cylinder.evalDistribution(np.asarray([0.1])) 
    775     assert np.isnan(Iq[0]), "empty distribution fails" 
    776752 
    777753def test_model_list(): 
    778754    # type: () -> None 
    779755    """ 
    780     Make sure that all models build as sasview models 
     756    Make sure that all models build as sasview models. 
    781757    """ 
    782758    from .exception import annotate_exception 
    … …  
    805781 
    806782if __name__ == "__main__": 
    807     print("cylinder(0.1,0.1)=%g"%test_cylinder()) 
    808     #test_empty_distribution() 
     783    print("cylinder(0.1,0.1)=%g"%test_model()) 

  • sasmodels/sesans.py

    r94d13f1 rb397165  
    1414import numpy as np  # type: ignore 
    1515from numpy import pi, exp  # type: ignore 
    16 from scipy.special import j0 
    17  
    18 class SesansTransform(object): 
    19     """ 
    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 
    23  
    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. 
    31     """ 
    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 
    35  
    36     #: q values to calculate when computing transform 
    37     q_calc = None  # type: np.ndarray 
    38  
    39     # transform arrays 
    40     _H = None  # type: np.ndarray 
    41     _H0 = None # type: np.ndarray 
    42  
    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) 
    48  
    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 
    55  
    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') 
    60  
    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 
    66  
    67         H0 = np.float32(dq/(2*pi)) * q 
    68  
    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 
    72  
    73         self.q_calc = q 
    74         self._H, self._H0 = H, H0 
     16from scipy.special import jv as besselj 
     17#import direct_model.DataMixin as model 
     18         
     19def 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. 
     27    """ 
     28    from sas.sascalc.data_util.nxsunit import Converter 
     29 
     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     
     36def make_all_q(data): 
     37    """ 
     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] 
     55 
     56def 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) 
     76 
     77    return result 
     78 
     79def 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   
     86def 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                   
     91def 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                         
     96def 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     
     102 
     103 
     104def 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 
     109 
     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 
     131 
     132 
     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] 
     136 
     137 
     138 
     139 
     140def 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     
     168def 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 
Note: See TracChangeset for help on using the changeset viewer.