Changeset a98958b in sasmodels


Ignore:
Timestamp:
Mar 16, 2016 1:01:34 PM (9 years ago)
Author:
krzywon
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:
2f0c07d, 092933d
Parents:
5041682
Message:

Good progress on SESANS scripting.

Location:
example
Files:
6 added
1 edited
1 moved

Legend:

Unmodified
Added
Removed
  • example/sesansfit.py

    r346bc88 ra98958b  
     1#TODO: Convert units properly (nm -> A) 
     2#TODO: Implement constraints 
     3 
    14from bumps.names import * 
    2  
    35from sasmodels import core, bumps_model 
    46 
    5 if True: # fix when data loader exists 
    6 #    from sas.dataloader.readers\ 
    7     from sas.dataloader.loader import Loader 
    8     loader = Loader() 
    9     filename = 'testsasview1.ses' 
    10     data = loader.load(filename) 
    11     if data is None: raise IOError("Could not load file %r"%(filename,)) 
    12     data.x /= 10 
    13 #    print data 
    14 #    data = load_sesans('mydatfile.pz') 
    15 #    sans_data = load_sans('mysansfile.xml') 
     7HAS_CONVERTER = True 
     8try: 
     9    from sas.sascalc.data_util.nxsunit import Converter 
     10except ImportError: 
     11    HAS_CONVERTER = False 
    1612 
    17 else: 
    18     SElength = np.linspace(0, 2400, 61) # [A] 
    19     data = np.ones_like(SElength) 
    20     err_data = np.ones_like(SElength)*0.03 
     13def sesans_fit(file, model_name, initial_vals={}, custom_params={}, param_range=[]): 
     14    """ 
    2115 
    22     class Sample: 
    23         zacceptance = 0.1 # [A^-1] 
    24         thickness = 0.2 # [cm] 
    25          
    26     class SESANSData1D: 
    27         #q_zmax = 0.23 # [A^-1] 
    28         lam = 0.2 # [nm] 
    29         x = SElength 
    30         y = data 
    31         dy = err_data 
    32         sample = Sample() 
    33     data = SESANSData1D() 
     16    @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 initial_vals: dictionary of {param_name : initial_value} 
     19    @param custom_params: dictionary of {custom_parameter_name : Parameter() object} 
     20    @param param_range: dictionary of {parameter_name : [minimum, maximum]} 
     21    @return: FitProblem for Bumps usage 
     22    """ 
     23    try: 
     24        from sas.sascalc.dataloader.loader import Loader 
     25        loader = Loader() 
     26        data = loader.load(file) 
     27        if data is None: raise IOError("Could not load file %r"%(file)) 
     28        if HAS_CONVERTER == True: 
     29            default_unit = "A" 
     30            data_conv_q = Converter(data._xunit) 
     31            data.x = data_conv_q(data.x, units=default_unit) 
     32            data._xunit = default_unit 
    3433 
    35 radius = 1000 
    36 data.Rmax = 3*radius # [A] 
     34    except: 
     35        # If no loadable data file, generate random data 
     36        SElength = np.linspace(0, 2400, 61) # [A] 
     37        data = np.ones_like(SElength) 
     38        err_data = np.ones_like(SElength)*0.03 
    3739 
    38 ##  Sphere parameters 
     40        class Sample: 
     41            zacceptance = 0.1 # [A^-1] 
     42            thickness = 0.2 # [cm] 
    3943 
    40 kernel = core.load_model("sphere", dtype='single') 
    41 phi = Parameter(0.1, name="phi") 
    42 model = bumps_model.Model(kernel, 
    43     scale=phi*(1-phi), sld=7.0, solvent_sld=1.0, radius=radius, 
    44     ) 
    45 phi.range(0.001,0.5) 
    46 #model.radius.pmp(40) 
    47 model.radius.range(1,10000) 
    48 #model.sld.pm(5) 
    49 #model.background 
    50 #model.radius_pd=0 
    51 #model.radius_pd_n=0 
     44        class SESANSData1D: 
     45            #q_zmax = 0.23 # [A^-1] 
     46            lam = 0.2 # [nm] 
     47            x = SElength 
     48            y = data 
     49            dy = err_data 
     50            sample = Sample() 
     51        data = SESANSData1D() 
    5252 
    53 ### Tri-Axial Ellipsoid 
    54 # 
    55 #kernel = core.load_model("triaxial_ellipsoid", dtype='single') 
    56 #phi = Parameter(0.1, name='phi') 
    57 #model = bumps_model.Model(kernel, 
    58 #    scale=phi*(1-phi), sld=7.0, solvent_sld=1.0, radius=radius, 
    59 #    ) 
    60 #phi.range(0.001,0.90) 
    61 ##model.radius.pmp(40) 
    62 #model.radius.range(100,10000) 
    63 ##model.sld.pmp(5) 
    64 ##model.background 
    65 ##model.radius_pd = 0 
    66 ##model.radius_pd_n = 0 
     53    radius = 1000 
     54    data.Rmax = 3*radius # [A] 
    6755 
    68 if False: # have sans data 
    69     M_sesans = bumps_model.Experiment(data=data, model=model) 
    70     M_sans = bumps_model.Experiment(data=sans_data, model=model) 
    71     problem = FitProblem([M_sesans, M_sans]) 
    72 else: 
    73     M_sesans = bumps_model.Experiment(data=data, model=model) 
    74     problem = FitProblem(M_sesans) 
     56    kernel = core.load_model(model_name) 
     57    model = bumps_model.Model(kernel) 
    7558 
     59    # Load custom parameters, initial values and parameter constraints 
     60    for k, v in custom_params.items(): 
     61        setattr(model, k, v) 
     62        model._parameter_names.append(k) 
     63    for k, v in initial_vals.items(): 
     64        param = model.parameters().get(k) 
     65        setattr(param, "value", v) 
     66    for k, v in param_range.items(): 
     67        param = model.parameters().get(k) 
     68        if param is not None: 
     69            setattr(param.bounds, "limits", v) 
     70 
     71    if False: # have sans data 
     72        M_sesans = bumps_model.Experiment(data=data, model=model) 
     73        M_sans = bumps_model.Experiment(data=sans_data, model=model) 
     74        problem = FitProblem([M_sesans, M_sans]) 
     75    else: 
     76        M_sesans = bumps_model.Experiment(data=data, model=model) 
     77        problem = FitProblem(M_sesans) 
     78    return problem 
Note: See TracChangeset for help on using the changeset viewer.