Changeset 735507b in sasmodels


Ignore:
Timestamp:
Mar 18, 2016 4:12:48 AM (9 years ago)
Author:
gonzalezm
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
2622b3f
Parents:
03e8a6e (diff), c094758 (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' of ​https://github.com/SasView/sasmodels

Conflicts:

doc/genmodel.py

Files:
8 added
2 deleted
14 edited

Legend:

Unmodified
Added
Removed

  • doc/genmodel.py

    r03e8a6e r735507b  
    1010# Convert ../sasmodels/models/name.py to name 
    1111model_name = os.path.basename(sys.argv[1])[:-3] 
     12model_info = core.load_model_info(model_name) 
     13model = core.build_model(model_info) 
    1214 
    1315# Load the doc string from the module definition file and store it in rst 
    14 docstr = generate.make_doc(core.load_model_info(model_name)) 
    15      
    16 # Generate automatically plot of the model and add it to rst documentation 
     16docstr = generate.make_doc(model_info) 
    1717 
    18 info = core.load_model_info(model_name) 
    1918 
    2019# Calculate 1D curve for default parameters 
    21 pars = dict((p[0], p[2]) for p in info['parameters']) 
     20pars = dict((p.name, p.default) for p in model_info['parameters']) 
    2221 
    2322# Plotting ranges and options 
    2423opts = { 
    25         'xscale'    : 'log', 
    26         'yscale'    : 'log' if not info['structure_factor'] else 'linear', 
    27         'xlabel'    : '$Q \/(\AA^{-1})$', 
    28         'ylabel'    : '$I(Q) \/(\mathrm{cm}^{-1})$' if not info['structure_factor'] else '$S(Q) \/(\mathrm{cm}^{-1})$', 
    29         'qmin'      : 0.001, 
    30         'qmax'      : 1.0, 
    31         'nq'        : 1000, 
    32         'nq2d'      : 100, 
     24    'xscale'    : 'log', 
     25    'yscale'    : 'log' if not model_info['structure_factor'] else 'linear', 
     26    'zscale'    : 'log' if not model_info['structure_factor'] else 'linear', 
     27    'q_min'     : 0.001, 
     28    'q_max'     : 1.0, 
     29    'nq'        : 1000, 
     30    'nq2d'      : 400, 
     31    'vmin'      : 1e-3,  # floor for the 2D data results 
     32    'qx_max'    : 0.5, 
    3333} 
    3434 
    35 qmin, qmax, nq = opts['qmin'], opts['qmax'], opts['nq'] 
    36 qmin = math.log10(qmin) 
    37 qmax = math.log10(qmax) 
    38 q = np.logspace(qmin, qmax, nq) 
    39 data = empty_data1D(q) 
    40 model = core.load_model(model_name) 
    41 calculator = DirectModel(data, model) 
    42 Iq1D = calculator() 
    4335 
    44 # TO DO: Generation of 2D plots  
    45 # Problem in sasmodels.direct_model._calc_theory 
    46 # There self._kernel.q_input.nq  gets a value of 0 in the 2D case 
    47 # and returns a 0 numpy array (it does not call the C code) 
     36def plot_1d(model, opts, ax): 
     37    q_min, q_max, nq = opts['q_min'], opts['q_max'], opts['nq'] 
     38    q_min = math.log10(q_min) 
     39    q_max = math.log10(q_max) 
     40    q = np.logspace(q_min, q_max, nq) 
     41    data = empty_data1D(q) 
     42    calculator = DirectModel(data, model) 
     43    Iq1D = calculator() 
    4844 
    49 # If 2D model, compute 2D image 
    50 #if info['has_2d'] != []: 
    51 #    qmax, nq2d = opts['qmax'], opts['nq2d'] 
    52 #    data2d = empty_data2D(np.linspace(-qmax, qmax, nq2d), resolution=0.0)  
    53 #    #model = core.load_model(model_name) 
    54 #    calculator = DirectModel(data2d, model) 
    55 #    Iq2D = calculator() 
     45    ax.plot(q, Iq1D, color='blue', lw=2, label=model_info['name']) 
     46    ax.set_xlabel(r'$Q \/(\AA^{-1})$') 
     47    ax.set_ylabel(r'$I(Q) \/(\mathrm{cm}^{-1})$') 
     48    ax.set_xscale(opts['xscale']) 
     49    ax.set_yscale(opts['yscale']) 
     50    #ax.legend(loc='best') 
    5651 
    57 # Generate image (comment IF for 1D/2D for the moment) and generate only 1D 
    58 #if info['has_2d'] == []: 
    59 #    fig = plt.figure() 
    60 #    ax = fig.add_subplot(1,1,1) 
    61 #    ax.plot(q, Iq1D, color='blue', lw=2, label=model_name) 
    62 #    ax.set_xlabel(r'$Q \/(\AA^{-1})$') 
    63 #    ax.set_xscale(opts['xscale'])    
    64 #    ax.set_ylabel(r'$I(Q) \/(\mathrm{cm}^{-1})$') 
    65 #    ax.set_yscale(opts['yscale'])   
    66 #    ax.legend() 
    67 #else: 
    68 #    # need figure with 1D + 2D 
    69 #    pass 
    70 fig = plt.figure() 
    71 ax = fig.add_subplot(1,1,1) 
    72 ax.plot(q, Iq1D, color='blue', lw=2, label=info['name']) 
    73 ax.set_xlabel(opts['xlabel']) 
    74 ax.set_xscale(opts['xscale'])    
    75 ax.set_ylabel(opts['ylabel']) 
    76 ax.set_yscale(opts['yscale'])   
    77 ax.legend() 
    78   
     52def plot_2d(model, opts, ax): 
     53    qx_max, nq2d = opts['qx_max'], opts['nq2d'] 
     54    q = np.linspace(-qx_max, qx_max, nq2d) 
     55    data2d = empty_data2D(q, resolution=0.0) 
     56    calculator = DirectModel(data2d, model) 
     57    Iq2D = calculator() #background=0) 
     58    Iq2D = Iq2D.reshape(nq2d, nq2d) 
     59    if opts['zscale'] == 'log': 
     60        Iq2D = np.log(np.clip(Iq2D, opts['vmin'], np.inf)) 
     61    h = ax.imshow(Iq2D, interpolation='nearest', aspect=1, origin='upper', 
     62           extent=[-qx_max, qx_max, -qx_max, qx_max], cmap=ice_cm()) 
     63    # , vmin=vmin, vmax=vmax) 
     64    ax.set_xlabel(r'$Q_x \/(\AA^{-1})$') 
     65    ax.set_ylabel(r'$Q_y \/(\AA^{-1})$') 
     66 
     67def ice_cm(): 
     68    from matplotlib._cm import _Blues_data 
     69    from matplotlib import colors 
     70    from matplotlib import rcParams 
     71    def from_white(segments): 
     72        scale = 1.0/segments[0][1] 
     73        return [(k, v*scale, w*scale) for k, v, w in segments] 
     74    ice_data = dict((k,from_white(v)) for k,v in _Blues_data.items()) 
     75    ice = colors.LinearSegmentedColormap("ice", ice_data, rcParams['image.lut']) 
     76    return ice 
     77 
     78 
     79# Generate image  
     80fig_height = 3.0 # in 
     81fig_left = 0.6 # in 
     82fig_right = 0.5 # in 
     83fig_top = 0.6*0.25 # in 
     84fig_bottom = 0.6*0.75 
     85if model_info['has_2d']: 
     86    plot_height = fig_height - (fig_top+fig_bottom) 
     87    plot_width = plot_height 
     88    fig_width = 2*(plot_width + fig_left + fig_right) 
     89    aspect = (fig_width, fig_height) 
     90    ratio = aspect[0]/aspect[1] 
     91    ax_left = fig_left/fig_width 
     92    ax_bottom = fig_bottom/fig_height 
     93    ax_height = plot_height/fig_height 
     94    ax_width = ax_height/ratio # square axes 
     95    fig = plt.figure(figsize=aspect) 
     96    ax2d = fig.add_axes([0.5+ax_left, ax_bottom, ax_width, ax_height]) 
     97    plot_2d(model, opts, ax2d) 
     98    ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height]) 
     99    #ax.set_aspect('square') 
     100else: 
     101    plot_height = fig_height - (fig_top+fig_bottom) 
     102    plot_width = (1+np.sqrt(5))/2*fig_height 
     103    fig_width = plot_width + fig_left + fig_right 
     104    ax_left = fig_left/fig_width 
     105    ax_bottom = fig_bottom/fig_height 
     106    ax_width = plot_width/fig_width 
     107    ax_height = plot_height/fig_height 
     108    aspect = (fig_width, fig_height) 
     109    fig = plt.figure(figsize=aspect) 
     110    ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height]) 
     111    plot_1d(model, opts, ax1d) 
    79112 
    80113# Save image in model/img 
    81114figname = model_name + '_autogenfig.png' 
    82115filename = os.path.join('model', 'img', figname) 
    83 plt.savefig(filename) 
     116plt.savefig(filename, bbox_inches='tight') 
     117#print "figure saved in",filename 
    84118 
    85119# Auto caption for figure 
    86120captionstr = '\n' 
    87 captionstr += '.. figure:: img/' + model_name + '_autogenfig.png\n' 
     121captionstr += '.. figure:: img/' + model_info['id'] + '_autogenfig.png\n' 
    88122captionstr += '\n' 
    89 #if info['has_2d'] == []: 
    90 #    captionstr += '    1D plot corresponding to the default parameters of the model.\n' 
    91 #else: 
    92 #    captionstr += '    1D and 2D plots corresponding to the default parameters of the model.\n' 
    93123captionstr += '    1D plot corresponding to the default parameters of the model.\n' 
    94124captionstr += '\n' 
    … …  
    104134else: 
    105135    print '------------------------------------------------------------------' 
    106     print 'References NOT FOUND for model: ', model_name 
     136    print 'References NOT FOUND for model: ', model_info['id'] 
    107137    print '------------------------------------------------------------------' 
    108138    docstr = docstr + captionstr 

  • example/sesans_parameters_css-hs.py

    ra98958b r84db7a5  
    88# Enter the model name to use 
    99model_name = "core_shell_sphere*hardsphere" 
     10 
     11# DO NOT MODIFY THIS LINE 
     12model = sesansfit.get_bumps_model(model_name) 
    1013 
    1114# Enter any custom parameters 
    … …  
    1922# Initial parameter values (if other than defaults) 
    2023initial_vals = { 
    21     "scale" : 0.09, 
    2224    "core_sld" : 1.0592, 
    2325    "solvent_sld" : 2.88, 
    24     "shell_sld" : 2.88, 
    2526    "radius" : 890, 
    26     "thickness" : 130, 
    27     "volfraction" : 0.45 
     27    "thickness" : 130 
    2828} 
    2929 
    … …  
    3636} 
    3737 
     38# Constraints 
     39# model.param_name = f(other params) 
     40# EXAMPLE: model.scale = model.radius*model.radius*(1 - phi) - where radius and scale are model functions and phi is 
     41# a custom parameter 
     42model.scale = phi*(1-phi) 
     43model.volfraction = phi 
     44model.shell_sld = pen*2.88 
     45 
    3846# Send to the fitting engine 
    3947problem = sesansfit.sesans_fit(sesans_file, model_name, initial_vals, custom_params, param_range) 

  • example/sesans_parameters_sphere.py

    ra98958b r84db7a5  
    99model_name = "sphere" 
    1010 
     11# DO NOT MODIFY THIS LINE 
     12model = sesansfit.get_bumps_model(model_name) 
     13 
    1114# Enter any custom parameters 
     15# name = Parameter(initial_value, name='name') 
    1216phi = Parameter(0.10, name='phi') 
     17# Add the parameters to this list that should be displayed in the fitting window 
    1318custom_params = {"phi" : phi} 
    1419 
    15 # SESANS data file 
     20# SESANS data file name 
    1621sesans_file = "sphere.ses" 
    1722 
    1823# Initial parameter values (if other than defaults) 
     24# "model_parameter_name" : value 
    1925initial_vals = { 
    20     "scale" : phi*(1 - phi), 
    2126    "sld" : 7.0, 
     27    "radius" : 1000, 
    2228    "solvent_sld" : 1.0, 
    23     "radius" : 1000, 
    2429} 
    2530 
    2631# Ranges for parameters if other than default 
     32# "model_parameter_name" : [min, max] 
    2733param_range = { 
    2834    "phi" : [0.001, 0.5], 
    … …  
    3036} 
    3137 
     38# Constraints 
     39# model.param_name = f(other params) 
     40# EXAMPLE: model.scale = model.radius*model.radius*(1 - phi) - where radius and scale are model functions and phi is 
     41# a custom parameter 
     42model.scale = phi*(1-phi) 
     43 
    3244# Send to the fitting engine 
    33 problem = sesansfit.sesans_fit(sesans_file, model_name, initial_vals, custom_params, param_range) 
     45# DO NOT MODIFY THIS LINE 
     46problem = sesansfit.sesans_fit(sesans_file, model, initial_vals, custom_params, param_range) 
     47 

  • example/sesansfit.py

    ra98958b r84db7a5  
    1 #TODO: Convert units properly (nm -> A) 
    2 #TODO: Implement constraints 
    3  
    41from bumps.names import * 
    5 from sasmodels import core, bumps_model 
     2from sasmodels import core, bumps_model, sesans 
    63 
    74HAS_CONVERTER = True 
    … …  
    118    HAS_CONVERTER = False 
    129 
    13 def sesans_fit(file, model_name, initial_vals={}, custom_params={}, param_range=[]): 
     10def get_bumps_model(model_name): 
     11    kernel = core.load_model(model_name) 
     12    model = bumps_model.Model(kernel) 
     13    return model 
     14 
     15def sesans_fit(file, model, initial_vals={}, custom_params={}, param_range=[]): 
    1416    """ 
    15  
    1617    @param file: SESANS file location 
    17     @param model_name: model name string - can be model, model_1 * model_2, and/or model_1 + model_2 
     18    @param model: Bumps model object or model name - can be model, model_1 * model_2, and/or model_1 + model_2 
    1819    @param initial_vals: dictionary of {param_name : initial_value} 
    1920    @param custom_params: dictionary of {custom_parameter_name : Parameter() object} 
    2021    @param param_range: dictionary of {parameter_name : [minimum, maximum]} 
     22    @param constraints: dictionary of {parameter_name : constraint} 
    2123    @return: FitProblem for Bumps usage 
    2224    """ 
    … …  
    2931            default_unit = "A" 
    3032            data_conv_q = Converter(data._xunit) 
     33            for x in data.x: 
     34                print x 
    3135            data.x = data_conv_q(data.x, units=default_unit) 
     36            for x in data.x: 
     37                print x 
    3238            data._xunit = default_unit 
    3339 
    … …  
    5157        data = SESANSData1D() 
    5258 
    53     radius = 1000 
     59    if "radius" in initial_vals: 
     60        radius = initial_vals.get("radius") 
     61    else: 
     62        radius = 1000 
    5463    data.Rmax = 3*radius # [A] 
    5564 
    56     kernel = core.load_model(model_name) 
    57     model = bumps_model.Model(kernel) 
     65    if isinstance(model, basestring): 
     66        model = get_bumps_model(model) 
    5867 
    59     # Load custom parameters, initial values and parameter constraints 
     68    # Load custom parameters, initial values and parameter ranges 
    6069    for k, v in custom_params.items(): 
    6170        setattr(model, k, v) 
    … …  
    6978            setattr(param.bounds, "limits", v) 
    7079 
    71     if False: # have sans data 
     80    if False: # for future implementation 
    7281        M_sesans = bumps_model.Experiment(data=data, model=model) 
    7382        M_sans = bumps_model.Experiment(data=sans_data, model=model) 

  • example/sesansfit.sh

    r15bd6e7 rbdb653c  
    88set -x 
    99 
    10 pythonw -m bumps.cli $* 
     10python -m bumps.cli $* 

  • explore/J1c.py

    r0a6da3c rcbd37a7  
    99 
    1010 
    11 SHOW_DIFF = True  # True if show diff rather than function value 
     11SHOW_DIFF = True # True if show diff rather than function value 
     12#SHOW_DIFF = False # True if show diff rather than function value 
    1213LINEAR_X = False  # True if q is linearly spaced instead of log spaced 
     14#LINEAR_X = True # True if q is linearly spaced instead of log spaced 
     15FUNCTION = "2*J1(x)/x" 
    1316 
    14 def mp_J1c(vec, bits=500): 
     17def mp_fn(vec, bits=500): 
    1518    """ 
    1619    Direct calculation using sympy multiprecision library. 
    1720    """ 
    1821    with mp.workprec(bits): 
    19         return [_mp_J1c(mp.mpf(x)) for x in vec] 
     22        return [_mp_fn(mp.mpf(x)) for x in vec] 
    2023 
    21 def _mp_J1c(x): 
     24def _mp_fn(x): 
    2225    """ 
    23     Helper funciton for mp_j1c 
     26    Actual function that gets evaluated.  The caller just vectorizes. 
    2427    """ 
    2528    return mp.mpf(2)*mp.j1(x)/x 
    2629 
    27 def np_J1c(x, dtype): 
     30def np_fn(x, dtype): 
    2831    """ 
    2932    Direct calculation using scipy. 
    … …  
    3336    return np.asarray(2, dtype)*J1(x)/x 
    3437 
    35 def cephes_J1c(x, dtype, n): 
     38def sasmodels_fn(x, dtype, platform='ocl'): 
    3639    """ 
    3740    Calculation using pade approximant. 
    3841    """ 
    39     f = np.float64 if np.dtype(dtype) == np.float64 else np.float32 
    40     x = np.asarray(x, dtype) 
    41     ans = np.empty_like(x) 
    42     ax = abs(x) 
    43  
    44     # Branch a 
    45     a_idx = ax < f(8.0) 
    46     a_xsq = x[a_idx]**2 
    47     a_coeff1 = list(reversed((72362614232.0, -7895059235.0, 242396853.1, -2972611.439, 15704.48260, -30.16036606))) 
    48     a_coeff2 = list(reversed((144725228442.0, 2300535178.0, 18583304.74, 99447.43394, 376.9991397, 1.0))) 
    49     a_ans1 = np.polyval(np.asarray(a_coeff1[n:], dtype), a_xsq) 
    50     a_ans2 = np.polyval(np.asarray(a_coeff2[n:], dtype), a_xsq) 
    51     ans[a_idx] = f(2.0)*a_ans1/a_ans2 
    52  
    53     # Branch b 
    54     b_idx = ~a_idx 
    55     b_ax = ax[b_idx] 
    56     b_x = x[b_idx] 
    57  
    58     b_y = f(64.0)/(b_ax**2) 
    59     b_xx = b_ax - f(2.356194491) 
    60     b_coeff1 = list(reversed((1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6))) 
    61     b_coeff2 = list(reversed((0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6))) 
    62     b_ans1 = np.polyval(np.asarray(b_coeff1[n:], dtype),b_y) 
    63     b_ans2 = np.polyval(np.asarray(b_coeff2[n:], dtype), b_y) 
    64     b_sn, b_cn = np.sin(b_xx), np.cos(b_xx) 
    65     ans[b_idx] = np.sign(b_x)*np.sqrt(f(0.636619772)/b_ax) * (b_cn*b_ans1 - (f(8.0)/b_ax)*b_sn*b_ans2)*f(2.0)/b_x 
    66  
    67     return ans 
    68  
    69 def div_J1c(x, dtype): 
    70     f = np.float64 if np.dtype(dtype) == np.float64 else np.float32 
    71     x = np.asarray(x, dtype) 
    72     return f(2.0)*np.asarray([_J1(xi, f)/xi for xi in x], dtype) 
    73  
    74 def _J1(x, f): 
    75     ax = abs(x) 
    76     if ax < f(8.0): 
    77         y = x*x 
    78         ans1 = x*(f(72362614232.0) 
    79                   + y*(f(-7895059235.0) 
    80                   + y*(f(242396853.1) 
    81                   + y*(f(-2972611.439) 
    82                   + y*(f(15704.48260) 
    83                   + y*(f(-30.16036606))))))) 
    84         ans2 = (f(144725228442.0) 
    85                   + y*(f(2300535178.0) 
    86                   + y*(f(18583304.74) 
    87                   + y*(f(99447.43394) 
    88                   + y*(f(376.9991397) 
    89                   + y))))) 
    90         return ans1/ans2 
    91     else: 
    92         y = f(64.0)/(ax*ax) 
    93         xx = ax - f(2.356194491) 
    94         ans1 = (f(1.0) 
    95                   + y*(f(0.183105e-2) 
    96                   + y*(f(-0.3516396496e-4) 
    97                   + y*(f(0.2457520174e-5) 
    98                   + y*f(-0.240337019e-6))))) 
    99         ans2 = (f(0.04687499995) 
    100                   + y*(f(-0.2002690873e-3) 
    101                   + y*(f(0.8449199096e-5) 
    102                   + y*(f(-0.88228987e-6) 
    103                   + y*f(0.105787412e-6))))) 
    104         sn, cn = np.sin(xx), np.cos(xx) 
    105         ans = np.sqrt(f(0.636619772)/ax) * (cn*ans1 - (f(8.0)/ax)*sn*ans2) 
    106         return -ans if (x < f(0.0)) else ans 
     42    from sasmodels import core, data, direct_model 
     43    model = core.load_model('bessel', dtype=dtype) 
     44    calculator = direct_model.DirectModel(data.empty_data1D(x), model) 
     45    return calculator(background=0) 
    10746 
    10847def plotdiff(x, target, actual, label): 
    … …  
    11352    """ 
    11453    if SHOW_DIFF: 
    115         err = np.clip(abs((target-actual)/target), 0, 1) 
     54        err = abs((target-actual)/target) 
     55        #err = np.clip(err, 0, 1) 
    11656        pylab.loglog(x, err, '-', label=label) 
    11757    else: 
    … …  
    11959        pylab.semilogx(x, np.clip(actual,*limits),  '-', label=label) 
    12060 
    121 def compare(x, precision): 
     61def compare(x, precision, target): 
    12262    r""" 
    12363    Compare the different computation methods using the given precision. 
    12464    """ 
    125     target = np.asarray(mp_J1c(x), 'double') 
    126     #plotdiff(x, target, mp_J1c(x, 11), 'mp 11 bits') 
    127     plotdiff(x, target, np_J1c(x, precision), 'direct '+precision) 
    128     plotdiff(x, target, cephes_J1c(x, precision, 0), 'cephes '+precision) 
    129     #plotdiff(x, target, cephes_J1c(x, precision, 1), 'cephes '+precision) 
    130     #plotdiff(x, target, div_J1c(x, precision), 'cephes 2 J1(x)/x '+precision) 
     65    #plotdiff(x, target, mp_fn(x, 11), 'mp 11 bits') 
     66    plotdiff(x, target, np_fn(x, precision), 'numpy '+precision) 
     67    plotdiff(x, target, sasmodels_fn(x, precision, 0), 'sasmodels '+precision) 
    13168    pylab.xlabel("qr (1/Ang)") 
    13269    if SHOW_DIFF: 
    13370        pylab.ylabel("relative error") 
    13471    else: 
    135         pylab.ylabel("2 J1(x)/x") 
     72        pylab.ylabel(FUNCTION) 
    13673        pylab.semilogx(x, target,  '-', label="true value") 
    13774    if LINEAR_X: 
    … …  
    14784    else: 
    14885        qr = np.logspace(-3,5,400) 
     86    target = np.asarray(mp_fn(qr), 'double') 
    14987    pylab.subplot(121) 
    150     compare(qr, 'single') 
     88    compare(qr, 'single', target) 
    15189    pylab.legend(loc='best') 
    15290    pylab.subplot(122) 
    153     compare(qr, 'double') 
     91    compare(qr, 'double', target) 
    15492    pylab.legend(loc='best') 
    155     pylab.suptitle('2 J1(x)/x') 
     93    pylab.suptitle(FUNCTION) 
    15694 
    15795if __name__ == "__main__": 

  • sasmodels/convert.py

    r228cbd3 r78d3341  
    7272    new model definition end with sld. 
    7373    """ 
    74     return dict((p, (v*1e-6 if p.endswith('sld') 
     74    return dict((p, (v*1e-6 if p.startswith('sld') or p.endswith('sld') 
    7575                     else v*1e15 if 'ndensity' in p 
    7676                     else v)) 

  • sasmodels/data.py

    r7824276 rc094758  
    178178        self.qy_data, self.dqy_data = y, dy 
    179179        self.data, self.err_data = z, dz 
    180         self.mask = (~np.isnan(z) if z is not None 
    181                      else np.ones_like(x) if x is not None 
     180        self.mask = (np.isnan(z) if z is not None 
     181                     else np.zeros_like(x, dtype='bool') if x is not None 
    182182                     else None) 
    183183        self.q_data = np.sqrt(x**2 + y**2) 
    … …  
    440440 
    441441    if use_data or use_theory: 
     442        is_tof = np.any(data.lam!=data.lam[0]) 
    442443        if num_plots > 1: 
    443444            plt.subplot(1, num_plots, 1) 
    444445        if use_data: 
    445             plt.errorbar(data.x, data.y, yerr=data.dy) 
     446            if is_tof: 
     447                plt.errorbar(data.x, np.log(data.y)/(data.lam*data.lam), yerr=data.dy/data.y/(data.lam*data.lam)) 
     448            else: 
     449                plt.errorbar(data.x, data.y, yerr=data.dy) 
    446450        if theory is not None: 
    447             plt.plot(data.x, theory, '-', hold=True) 
     451            if is_tof: 
     452                plt.plot(data.x, np.log(theory)/(data.lam*data.lam), '-', hold=True) 
     453            else: 
     454                plt.plot(data.x, theory, '-', hold=True) 
    448455        if limits is not None: 
    449456            plt.ylim(*limits) 
    450         plt.xlabel('spin echo length (nm)') 
    451         plt.ylabel('polarization (P/P0)') 
     457 
     458        plt.xlabel('spin echo length ({})'.format(data._xunit)) 
     459        if is_tof: 
     460            plt.ylabel('(Log (P/P$_0$))/$\lambda^2$') 
     461        else: 
     462            plt.ylabel('polarization (P/P0)') 
     463 
    452464 
    453465    if resid is not None: 
    … …  
    455467            plt.subplot(1, num_plots, (use_data or use_theory) + 1) 
    456468        plt.plot(data.x, resid, 'x') 
    457         plt.xlabel('spin echo length (nm)') 
     469        plt.xlabel('spin echo length ({})'.format(data._xunit)) 
    458470        plt.ylabel('residuals (P/P0)') 
    459471 

  • sasmodels/generate.py

    rfcd7bbd r78d3341  
    9191    in *Iqxy* and *Imagnetic*.  "magnetic* parameters will be used in 
    9292    *Imagnetic* only.  If *type* is the empty string, the parameter will 
    93     be used in all of *Iq*, *Iqxy* and *Imagnetic*. 
     93    be used in all of *Iq*, *Iqxy* and *Imagnetic*.  "sld" parameters 
     94    can automatically be promoted to magnetic parameters, each of which 
     95    will have a magnitude and a direction, which may be different from 
     96    other sld parameters. 
    9497 
    9598    *description* is a short description of the parameter.  This will 
    … …  
    216219import re 
    217220import string 
     221import warnings 
    218222from collections import namedtuple 
    219223 
    … …  
    604608    """ 
    605609    partype = { 
    606         'volume': [], 'orientation': [], 'magnetic': [], '': [], 
     610        'volume': [], 'orientation': [], 'magnetic': [], 'sld': [], '': [], 
    607611        'fixed-1d': [], 'fixed-2d': [], 'pd-1d': [], 'pd-2d': [], 
    608612        'pd-rel': set(), 
    … …  
    618622        elif p.type == 'orientation': 
    619623            partype['pd-2d'].append(p.name) 
    620         elif p.type == '': 
     624        elif p.type in ('', 'sld'): 
    621625            partype['fixed-1d'].append(p.name) 
    622626            partype['fixed-2d'].append(p.name) 
    … …  
    632636    """ 
    633637    # convert parameters into named tuples 
     638    for p in model_info['parameters']: 
     639        if p[4] == '' and (p[0].startswith('sld') or p[0].endswith('sld')): 
     640            p[4] = 'sld' 
     641            # TODO: make sure all models explicitly label their sld parameters 
     642            #raise ValueError("%s.%s needs to be explicitly set to type 'sld'" %(model_info['id'], p[0])) 
     643 
    634644    pars = [Parameter(*p) for p in model_info['parameters']] 
    635645    # Fill in the derived attributes 

  • sasmodels/kernelcl.py

    r17bbadd rc094758  
    367367        self.q_vectors = [_stretch_input(q, self.dtype, 32) for q in q_vectors] 
    368368        context = env.get_context(self.dtype) 
     369        self.global_size = [self.q_vectors[0].size] 
     370        #print("creating inputs of size", self.global_size) 
    369371        self.q_buffers = [ 
    370372            cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q) 
    371373            for q in self.q_vectors 
    372374        ] 
    373         self.global_size = [self.q_vectors[0].size] 
    374375 
    375376    def release(self): 

  • sasmodels/models/adsorbed_layer.py

    r54954e1 rf10bc52  
    8989               background = 'background') 
    9090 
     91# these unit test values taken from SasView 3.1.2 
    9192tests =  [ 
    9293    [{'scale': 1.0, 'second_moment': 23.0, 'adsorbed_amount': 1.9,  

  • sasmodels/models/mono_gauss_coil.py

    r15bd6e7 rf10bc52  
    8686               background = 'background') 
    8787 
     88# these unit test values taken from SasView 3.1.2 
    8889tests =  [ 
    89     [{'scale': 70.0, 'radius_gyration': 75.0, 'background': 0.0}, 
     90    [{'scale': 1.0, 'i_zero': 70.0, 'radius_gyration': 75.0, 'background': 0.0}, 
    9091     [0.0106939, 0.469418], [57.1241, 0.112859]], 
    9192    ] 

  • sasmodels/models/poly_gauss_coil.py

    r15bd6e7 r09b84ed  
    5050""" 
    5151 
    52 from numpy import inf, sqrt, power 
     52from numpy import inf, sqrt, exp, power 
    5353 
    5454name =  "poly_gauss_coil" 
    … …  
    6969def Iq(q, i_zero, radius_gyration, polydispersity): 
    7070    # pylint: disable = missing-docstring 
    71     # need to trap the case of the polydispersity being 1 (ie, monodispersity) 
    7271    u = polydispersity - 1.0 
    73     if polydispersity == 1: 
    74        minusoneonu = -1.0 / u 
    75     else: 
    76        minusoneonu = -1.0 / u 
    7772    z = ((q * radius_gyration) * (q * radius_gyration)) / (1.0 + 2.0 * u) 
    7873    if (q == 0).any(): 
    79        inten = i_zero 
     74        inten = i_zero 
    8075    else: 
    81        inten = i_zero * 2.0 * (power((1.0 + u * z),minusoneonu) + z - 1.0 ) / ((1.0 + u) * (z * z)) 
     76    # need to trap the case of the polydispersity being 1 (ie, monodispersity!) 
     77        if polydispersity == 1: 
     78            inten = i_zero * 2.0 * (exp(-z) + z - 1.0 ) / (z * z) 
     79        else: 
     80            minusoneonu = -1.0 / u 
     81            inten = i_zero * 2.0 * (power((1.0 + u * z),minusoneonu) + z - 1.0 ) / ((1.0 + u) * (z * z)) 
    8282    return inten 
    83 Iq.vectorized =  True  # Iq accepts an array of q values 
     83#Iq.vectorized =  True  # Iq accepts an array of q values 
    8484 
    8585def Iqxy(qx, qy, *args): 
    … …  
    100100               background = 'background') 
    101101 
     102# these unit test values taken from SasView 3.1.2 
    102103tests =  [ 
    103     [{'scale': 70.0, 'radius_gyration': 75.0, 'polydispersity': 2.0, 'background': 0.0}, 
     104    [{'scale': 1.0, 'i_zero': 70.0, 'radius_gyration': 75.0, 'polydispersity': 2.0, 'background': 0.0}, 
    104105     [0.0106939, 0.469418], [57.6405, 0.169016]], 
    105106    ] 

  • sasmodels/product.py

    rfcd7bbd r3bcd03d  
    1818SCALE=0 
    1919BACKGROUND=1 
    20 EFFECT_RADIUS=2 
     20RADIUS_EFFECTIVE=2 
    2121VOLFRACTION=3 
    2222 
    … …  
    3131    assert s_pars[BACKGROUND].name == 'background' 
    3232    # We require structure factors to start with effect radius and volfraction 
    33     assert s_pars[EFFECT_RADIUS].name == 'effect_radius' 
     33    assert s_pars[RADIUS_EFFECTIVE].name == 'radius_effective' 
    3434    assert s_pars[VOLFRACTION].name == 'volfraction' 
    3535    # Combine the parameter sets.  We are skipping the first three 
Note: See TracChangeset for help on using the changeset viewer.