Changeset 87985ca in sasmodels


Ignore:
Timestamp:
Sep 2, 2014 12:35:23 PM (10 years ago)
Author:
Paul Kienzle <pkienzle@…>
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:
19dcb933
Parents:
f4cf580
Message:

clean up source tree

Files:
36 deleted
5 edited

Legend:

Unmodified
Added
Removed

  • compare.py

    • Property mode changed from 100644 to 100755
    ra953943 r87985ca  
    22# -*- coding: utf-8 -*- 
    33 
    4 from sasmodel import SasModel, load_data, set_beam_stop, plot_data, tic, toc 
     4import sys 
     5import math 
     6 
    57import numpy as np 
    68 
     9from sasmodels.bumps_model import BumpsModel, plot_data, tic 
     10from sasmodels import gpu, dll 
     11from sasmodels.convert import revert_model 
     12 
    713 
    814def sasview_model(modelname, **pars): 
    9     modelname = modelname+"Model" 
     15    """ 
     16    Load a sasview model given the model name. 
     17    """ 
     18    # convert model parameters from sasmodel form to sasview form 
     19    #print "old",sorted(pars.items()) 
     20    modelname, pars = revert_model(modelname, pars) 
     21    #print "new",sorted(pars.items()) 
    1022    sans = __import__('sans.models.'+modelname) 
    1123    ModelClass = getattr(getattr(sans.models,modelname,None),modelname,None) 
    … …  
    2133        elif k.endswith("_pd_nsigma"): 
    2234            model.dispersion[k[:-10]]['nsigmas'] = v 
     35        elif k.endswith("_pd_type"): 
     36            model.dispersion[k[:-8]]['type'] = v 
    2337        else: 
    2438            model.setParam(k, v) 
    2539    return model 
    2640 
    27  
    28 def sasview_eval(model, data): 
    29     theory = model.evalDistribution([data.qx_data, data.qy_data]) 
    30     return theory 
    31  
    32 def plot(data, cpumodel, gpumodel, N, pars): 
    33  
    34     model = SasModel(data, gpumodel, dtype='f', **pars) 
    35     tic() 
    36     for i in range(N): 
    37         #pars['scale'] = np.random.rand() 
    38         model.update() 
    39         gpu = model.theory() 
    40     gpu_time = toc()*1000./N 
    41     print "ocl t=%.1f ms"%gpu_time 
    42     tic() 
    43     for i in range(N): 
    44         cpu = sasview_eval(cpumodel, data) 
    45     cpu_time = toc()*1000./N 
    46     print "omp t=%.1f ms"%cpu_time 
    47  
     41def load_opencl(modelname, dtype='single'): 
     42    sasmodels = __import__('sasmodels.models.'+modelname) 
     43    module = getattr(sasmodels.models, modelname, None) 
     44    kernel = gpu.load_model(module, dtype=dtype) 
     45    return kernel 
     46 
     47def load_ctypes(modelname, dtype='single'): 
     48    sasmodels = __import__('sasmodels.models.'+modelname) 
     49    module = getattr(sasmodels.models, modelname, None) 
     50    kernel = dll.load_model(module, dtype=dtype) 
     51    return kernel 
     52 
     53def randomize(p, v): 
     54    """ 
     55    Randomizing parameter. 
     56 
     57    Guess the parameter type from name. 
     58    """ 
     59    if any(p.endswith(s) for s in ('_pd_n','_pd_nsigma','_pd_type')): 
     60        return v 
     61    elif any(s in p for s in ('theta','phi','psi')): 
     62        # orientation in [-180,180], orientation pd in [0,45] 
     63        if p.endswith('_pd'): 
     64            return 45*np.random.rand() 
     65        else: 
     66            return 360*np.random.rand() - 180 
     67    elif 'sld' in p: 
     68        # sld in in [-0.5,10] 
     69        return 10.5*np.random.rand() - 0.5 
     70    elif p.endswith('_pd'): 
     71        # length pd in [0,1] 
     72        return np.random.rand() 
     73    else: 
     74        # length, scale, background in [0,200] 
     75        return 200*np.random.rand() 
     76 
     77def parlist(pars): 
     78    return "\n".join("%s: %s"%(p,v) for p,v in sorted(pars.items())) 
     79 
     80def suppress_pd(pars): 
     81    """ 
     82    Suppress theta_pd for now until the normalization is resolved. 
     83 
     84    May also suppress complete polydispersity of the model to test 
     85    models more quickly. 
     86    """ 
     87    for p in pars: 
     88        if p.endswith("_pd"): pars[p] = 0 
     89 
     90def compare(name, pars, Ncpu, Ngpu, opts, set_pars): 
     91 
     92    # Sort out data 
     93    qmax = 1.0 if '-highq' in opts else (0.2 if '-midq' in opts else 0.05) 
     94    if "-1d" in opts: 
     95        from sasmodels.bumps_model import empty_data1D 
     96        qmax = math.log10(qmax) 
     97        data = empty_data1D(np.logspace(qmax-3, qmax, 128)) 
     98        index = slice(None, None) 
     99    else: 
     100        from sasmodels.bumps_model import empty_data2D, set_beam_stop 
     101        data = empty_data2D(np.linspace(-qmax, qmax, 128)) 
     102        set_beam_stop(data, 0.004) 
     103        index = ~data.mask 
     104    is2D = hasattr(data, 'qx_data') 
     105 
     106    # modelling accuracy is determined by dtype and cutoff 
     107    dtype = 'double' if '-double' in opts else 'single' 
     108    cutoff_opts = [s for s in opts if s.startswith('-cutoff')] 
     109    cutoff = float(cutoff_opts[0].split('=')[1]) if cutoff_opts else 1e-5 
     110 
     111    # randomize parameters 
     112    random_opts = [s for s in opts if s.startswith('-random')] 
     113    if random_opts: 
     114        if '=' in random_opts[0]: 
     115            seed = int(random_opts[0].split('=')[1]) 
     116        else: 
     117            seed = int(np.random.rand()*10000) 
     118        np.random.seed(seed) 
     119        print "Randomize using -random=%i"%seed 
     120        # Note: the sort guarantees order of calls to random number generator 
     121        pars = dict((p,randomize(p,v)) for p,v in sorted(pars.items())) 
     122        # The capped cylinder model has a constraint on its parameters 
     123        if name == 'capped_cylinder' and pars['cap_radius'] < pars['radius']: 
     124            pars['radius'],pars['cap_radius'] = pars['cap_radius'],pars['radius'] 
     125    pars.update(set_pars) 
     126 
     127    # parameter selection 
     128    if '-mono' in opts: 
     129        suppress_pd(pars) 
     130    if '-pars' in opts: 
     131        print "pars",parlist(pars) 
     132 
     133    # OpenCl calculation 
     134    if Ngpu > 0: 
     135        try: 
     136            model = load_opencl(name, dtype=dtype) 
     137        except Exception,exc: 
     138            print exc 
     139            print "... trying again with single precision" 
     140            model = load_opencl(name, dtype='single') 
     141        problem = BumpsModel(data, model, cutoff=cutoff, **pars) 
     142        toc = tic() 
     143        for _ in range(Ngpu): 
     144            #pars['scale'] = np.random.rand() 
     145            problem.update() 
     146            gpu = problem.theory() 
     147        gpu_time = toc()*1000./Ngpu 
     148        print "opencl t=%.1f ms, intensity=%.0f"%(gpu_time, sum(gpu[index])) 
     149        #print max(gpu), min(gpu) 
     150 
     151    # ctypes/sasview calculation 
     152    if Ncpu > 0 and "-ctypes" in opts: 
     153        model = load_ctypes(name, dtype=dtype) 
     154        problem = BumpsModel(data, model, cutoff=cutoff, **pars) 
     155        toc = tic() 
     156        for _ in range(Ncpu): 
     157            problem.update() 
     158            cpu = problem.theory() 
     159        cpu_time = toc()*1000./Ncpu 
     160        comp = "ctypes" 
     161        print "ctypes t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[index])) 
     162    elif Ncpu > 0: 
     163        model = sasview_model(name, **pars) 
     164        toc = tic() 
     165        for _ in range(Ncpu): 
     166            if is2D: 
     167                cpu = model.evalDistribution([data.qx_data, data.qy_data]) 
     168            else: 
     169                cpu = model.evalDistribution(data.x) 
     170        cpu_time = toc()*1000./Ncpu 
     171        comp = "sasview" 
     172        print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[index])) 
     173 
     174    # Compare, but only if computing both forms 
     175    if Ngpu > 0 and Ncpu > 0: 
     176        #print "speedup %.2g"%(cpu_time/gpu_time) 
     177        #print "max |gpu/cpu|", max(abs(gpu/cpu)), "%.15g"%max(abs(gpu)), "%.15g"%max(abs(cpu)) 
     178        #cpu *= max(gpu/cpu) 
     179        resid, relerr = np.zeros_like(gpu), np.zeros_like(gpu) 
     180        resid[index] = (gpu - cpu)[index] 
     181        relerr[index] = resid[index]/cpu[index] 
     182        print "max(|ocl-%s|)"%comp, max(abs(resid[index])) 
     183        print "max(|(ocl-%s)/ocl|)"%comp, max(abs(relerr[index])) 
     184 
     185    # Plot if requested 
     186    if '-noplot' in opts: return 
    48187    import matplotlib.pyplot as plt 
    49  
    50     print "gpu/cpu", max(gpu/cpu) 
    51     cpu *= max(gpu/cpu) 
    52     relerr = (gpu - cpu)/cpu 
    53     print "max(|(ocl-omp)/ocl|)", max(abs(relerr)) 
    54  
    55  
    56     plt.subplot(131); plot_data(data, cpu); plt.title("omp t=%.1f ms"%cpu_time) 
    57     plt.subplot(132); plot_data(data, gpu); plt.title("ocl t=%.1f ms"%gpu_time) 
    58     plt.subplot(133); plot_data(data, relerr); plt.title("relerr x 10^8"); plt.colorbar() 
     188    if Ncpu > 0: 
     189        if Ngpu > 0: plt.subplot(131) 
     190        plot_data(data, cpu, scale='log') 
     191        plt.title("%s t=%.1f ms"%(comp,cpu_time)) 
     192    if Ngpu > 0: 
     193        if Ncpu > 0: plt.subplot(132) 
     194        plot_data(data, gpu, scale='log') 
     195        plt.title("opencl t=%.1f ms"%gpu_time) 
     196    if Ncpu > 0 and Ngpu > 0: 
     197        plt.subplot(133) 
     198        err = resid if '-abs' in opts else relerr 
     199        errstr = "abs err" if '-abs' in opts else "rel err" 
     200        #err,errstr = gpu/cpu,"ratio" 
     201        plot_data(data, err, scale='linear') 
     202        plt.title("max %s = %.3g"%(errstr, max(abs(err[index])))) 
     203        if is2D: plt.colorbar() 
    59204    plt.show() 
    60205 
    61 def newlen(N): 
    62     import sys 
    63  
    64     if len(sys.argv) > 1: 
    65         N = int(sys.argv[1]) 
    66     data = load_data('December/DEC07098.DAT') 
    67     set_beam_stop(data, 0.004) 
    68  
    69     return data, N 
    70  
    71 def cyl(N=1): 
    72  
    73     dtype = "float" 
    74     pars = dict( 
    75         scale=.003, radius=64.1, length=66.96, sldCyl=.291e-6, sldSolv=5.77e-6, background=.1, 
    76         cyl_theta=90, cyl_phi=0, 
    77         radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3, 
    78         length_pd=0.1,length_pd_n=5, length_pd_nsigma=3, 
    79         cyl_theta_pd=0.1, cyl_theta_pd_n=5, cyl_theta_pd_nsigma=3, 
    80         cyl_phi_pd=0.1, cyl_phi_pd_n=10, cyl_phi_pd_nsigma=3) 
    81  
    82     from Models.code_cylinder_f import GpuCylinder as gpumodel 
    83     model = sasview_model('Cylinder', **pars) 
    84     data, N = newlen(N) 
    85  
    86     plot(data, model, gpumodel, N, pars) 
    87  
    88  
    89 def ellipse(N=1): 
    90  
    91     pars = dict(scale=.027, radius_a=60, radius_b=180, sldEll=.297e-6, sldSolv=5.773e-6, background=4.9, 
    92                 axis_theta=0, axis_phi=90, 
    93                 radius_a_pd=0.1, radius_a_pd_n=10, radius_a_pd_nsigma=3, 
    94                 radius_b_pd=0.1, radius_b_pd_n=10, radius_b_pd_nsigma=3, 
    95                 axis_theta_pd=0.1, axis_theta_pd_n=6, axis_theta_pd_nsigma=3, 
    96                 axis_phi_pd=0.1, axis_phi_pd_n=6, axis_phi_pd_nsigma=3,) 
    97  
    98     from Models.code_ellipse import GpuEllipse as gpumodel 
    99     model = sasview_model('Ellipsoid', **pars) 
    100  
    101     data, N = newlen(N) 
    102     plot(data, model, gpumodel, N, pars) 
    103  
    104 def coreshell(N=1): 
    105  
    106     data, N = newlen(N) 
    107  
    108     pars = dict(scale= 1.77881e-06, radius=325, thickness=25, length=34.2709, 
    109                  core_sld=1e-6, shell_sld=.291e-6, solvent_sld=7.105e-6, 
    110                  background=223.827, axis_theta=90, axis_phi=0, 
    111                  axis_theta_pd=15.8, 
    112                  radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3, 
    113                  length_pd=0.1, length_pd_n=10, length_pd_nsigma=3, 
    114                  thickness_pd=0.1, thickness_pd_n=5, thickness_pd_nsigma=3, 
    115                  axis_theta_pd_n=20, axis_theta_pd_nsigma=5, 
    116                  axis_phi_pd=0.0008748, axis_phi_pd_n=5, axis_phi_pd_nsigma=3,) 
    117  
    118     model = sasview_model('CoreShellCylinder', **pars) 
    119     from Models.code_coreshellcyl_f import GpuCoreShellCylinder as gpumodel 
    120  
    121     plot(data, model, gpumodel, N, pars) 
    122  
    123 def trellipse(N=1): 
    124  
    125     data, N = newlen(N) 
    126  
    127     pars = dict(scale=0.08, semi_axisA=15, semi_axisB=20, semi_axisC=500, 
    128                 sldEll=7.105e-6, sldSolv=.291e-6, 
    129                 background=5, axis_theta=0, axis_phi=0, axis_psi=0, 
    130                 axis_theta_pd=20, axis_theta_pd_n=10, axis_theta_pd_nsigma=3, 
    131                 axis_phi_pd=.1, axis_phi_pd_n=10, axis_phi_pd_nsigma=3, 
    132                 axis_psi_pd=30, axis_psi_pd_n=5, axis_psi_pd_nsigma=3, 
    133                 semi_axisA_pd=.1, semi_axisA_pd_n=5, semi_axisA_pd_nsigma=3, 
    134                 semi_axisB_pd=.1, semi_axisB_pd_n=5, semi_axisB_pd_nsigma=3, 
    135                 semi_axisC_pd=.1, semi_axisC_pd_n=5, semi_axisC_pd_nsigma=3,) 
    136  
    137     model = sasview_model('TriaxialEllipsoid', **pars) 
    138     from Models.code_triaxialellipse import GpuTriEllipse as gpumodel 
    139  
    140     plot(data, model, gpumodel, N, pars) 
    141  
    142 def lamellar(N=1): 
    143  
    144     data, N = newlen(N) 
    145  
    146     pars = dict(scale=0.08, 
    147                 bi_thick=19.2946, 
    148                 sld_bi=5.38e-6,sld_sol=7.105e-6, 
    149                 background=0.003, 
    150                 bi_thick_pd= 0.37765, bi_thick_pd_n=40, bi_thick_pd_nsigma=3) 
    151  
    152     model = sasview_model('Lamellar', **pars) 
    153     from Models.code_lamellar import GpuLamellar as gpumodel 
    154  
    155     plot(data, model, gpumodel, N, pars) 
    156  
    157 def cap(N=1): 
    158  
    159     data, N = newlen(N) 
    160  
    161     pars = dict(scale=.08, rad_cyl=20, rad_cap=40, len_cyl=400, 
    162                 sld_capcyl=1e-6, sld_solv=6.3e-6, 
    163                 background=0, theta=0, phi=0, 
    164                 rad_cyl_pd=.1, rad_cyl_pd_n=10, rad_cyl_pd_nsigma=3, 
    165                 rad_cap_pd=.1, rad_cap_pd_n=10, rad_cap_pd_nsigma=3, 
    166                 len_cyl_pd=.1, len_cyl_pd_n=3, len_cyl_pd_nsigma=3, 
    167                 theta_pd=.1, theta_pd_n=3, theta_pd_nsigma=3, 
    168                 phi_pd=.1, phi_pd_n=3, phi_pd_nsigma=3) 
    169  
    170  
    171     model = sasview_model('CappedCylinder', **pars) 
    172     from Models.code_capcyl import GpuCapCylinder as gpumodel 
    173  
    174     plot(data, model, gpumodel, N, pars) 
    175  
     206# =========================================================================== 
     207# 
     208USAGE=""" 
     209usage: compare.py model [Nopencl] [Nsasview] [options...] [key=val] 
     210 
     211Compare the speed and value for a model between the SasView original and the 
     212OpenCL rewrite. 
     213 
     214model is the name of the model to compare (see below). 
     215Nopencl is the number of times to run the OpenCL model (default=5) 
     216Nsasview is the number of times to run the Sasview model (default=1) 
     217 
     218Options (* for default): 
     219 
     220    -plot*/-noplot plots or suppress the plot of the model 
     221    -single/-double uses double precision for comparison 
     222    -lowq/-midq/-highq use q values up to 0.05, 0.2 or 1.0 
     223    -2d/-1d uses 1d or 2d random data 
     224    -preset/-random[=seed] preset or random parameters 
     225    -poly/-mono force monodisperse/polydisperse 
     226    -sasview/-ctypes whether cpu is tested using sasview or ctypes 
     227    -cutoff=1e-5/value cutoff for including a point in polydispersity 
     228    -nopars/-pars prints the parameter set or not 
     229    -rel/-abs plot relative or absolute error 
     230 
     231Key=value pairs allow you to set specific values to any of the model 
     232parameters. 
     233 
     234Available models: 
     235 
     236    %s 
     237""" 
     238 
     239VALID_OPTIONS = [ 
     240    'plot','noplot', 
     241    'single','double', 
     242    'lowq','midq','highq', 
     243    '2d','1d', 
     244    'preset','random', 
     245    'poly','mono', 
     246    'sasview','ctypes', 
     247    'nopars','pars', 
     248    'rel','abs', 
     249    ] 
     250 
     251def main(): 
     252    opts = [arg for arg in sys.argv[1:] if arg.startswith('-')] 
     253    args = [arg for arg in sys.argv[1:] if not arg.startswith('-')] 
     254    models = "\n    ".join("%-7s: %s"%(k,v.__name__.replace('_',' ')) 
     255                           for k,v in sorted(MODELS.items())) 
     256    if len(args) == 0: 
     257        print(USAGE%models) 
     258        sys.exit(1) 
     259    if args[0] not in MODELS: 
     260        print "Model %r not available. Use one of:\n    %s"%(args[0],models) 
     261        sys.exit(1) 
     262 
     263    valid_opts = set(VALID_OPTIONS) 
     264    invalid = [o[1:] for o in opts 
     265               if o[1:] not in valid_opts 
     266                  and not o.startswith('-cutoff=') 
     267                  and not o.startswith('-random=')] 
     268    if invalid: 
     269        print "Invalid options: %s"%(", ".join(invalid)) 
     270        sys.exit(1) 
     271 
     272    name, pars = MODELS[args[0]]() 
     273    Nopencl = int(args[1]) if len(args) > 1 else 5 
     274    Nsasview = int(args[2]) if len(args) > 3 else 1 
     275 
     276    # Fill in default polydispersity parameters 
     277    pds = set(p.split('_pd')[0] for p in pars if p.endswith('_pd')) 
     278    for p in pds: 
     279        if p+"_pd_nsigma" not in pars: pars[p+"_pd_nsigma"] = 3 
     280        if p+"_pd_type" not in pars: pars[p+"_pd_type"] = "gaussian" 
     281 
     282    # Fill in parameters given on the command line 
     283    set_pars = {} 
     284    for arg in args[3:]: 
     285        k,v = arg.split('=') 
     286        if k not in pars: 
     287            # extract base name without distribution 
     288            s = set(p.split('_pd')[0] for p in pars) 
     289            print "%r invalid; parameters are: %s"%(k,", ".join(sorted(s))) 
     290            sys.exit(1) 
     291        set_pars[k] = float(v) if not v.endswith('type') else v 
     292 
     293    compare(name, pars, Nsasview, Nopencl, opts, set_pars) 
     294 
     295# =========================================================================== 
     296# 
     297 
     298MODELS = {} 
     299def model(name): 
     300    def gather_function(fn): 
     301        MODELS[name] = fn 
     302        return fn 
     303    return gather_function 
     304 
     305 
     306@model('cyl') 
     307def cylinder(): 
     308    pars = dict( 
     309        scale=1, background=0, 
     310        sld=6, solvent_sld=1, 
     311        #radius=5, length=20, 
     312        radius=260, length=290, 
     313        theta=30, phi=0, 
     314        radius_pd=.2, radius_pd_n=1, 
     315        length_pd=.2,length_pd_n=1, 
     316        theta_pd=15, theta_pd_n=45, 
     317        phi_pd=15, phi_pd_n=1, 
     318        ) 
     319    return 'cylinder', pars 
     320 
     321@model('capcyl') 
     322def capped_cylinder(): 
     323    pars = dict( 
     324        scale=1, background=0, 
     325        sld=6, solvent_sld=1, 
     326        radius=260, cap_radius=290, length=290, 
     327        theta=30, phi=15, 
     328        radius_pd=.2, radius_pd_n=1, 
     329        cap_radius_pd=.2, cap_radius_pd_n=1, 
     330        length_pd=.2, length_pd_n=1, 
     331        theta_pd=15, theta_pd_n=45, 
     332        phi_pd=15, phi_pd_n=1, 
     333        ) 
     334    return 'capped_cylinder', pars 
     335 
     336 
     337@model('cscyl') 
     338def core_shell_cylinder(): 
     339    pars = dict( 
     340        scale=1, background=0, 
     341        core_sld=6, shell_sld=8, solvent_sld=1, 
     342        radius=45, thickness=25, length=340, 
     343        theta=30, phi=15, 
     344        radius_pd=.2, radius_pd_n=1, 
     345        length_pd=.2, length_pd_n=1, 
     346        thickness_pd=.2, thickness_pd_n=1, 
     347        theta_pd=15, theta_pd_n=45, 
     348        phi_pd=15, phi_pd_n=1, 
     349        ) 
     350    return 'core_shell_cylinder', pars 
     351 
     352 
     353@model('ell') 
     354def ellipsoid(): 
     355    pars = dict( 
     356        scale=1, background=0, 
     357        sld=6, solvent_sld=1, 
     358        rpolar=50, requatorial=30, 
     359        theta=30, phi=15, 
     360        rpolar_pd=.2, rpolar_pd_n=1, 
     361        requatorial_pd=.2, requatorial_pd_n=1, 
     362        theta_pd=15, theta_pd_n=45, 
     363        phi_pd=15, phi_pd_n=1, 
     364        ) 
     365    return 'ellipsoid', pars 
     366 
     367 
     368@model('ell3') 
     369def triaxial_ellipsoid(): 
     370    pars = dict( 
     371        scale=1, background=0, 
     372        sld=6, solvent_sld=1, 
     373        theta=30, phi=15, psi=5, 
     374        req_minor=25, req_major=36, rpolar=50, 
     375        req_minor_pd=0, req_minor_pd_n=1, 
     376        req_major_pd=0, req_major_pd_n=1, 
     377        rpolar_pd=.2, rpolar_pd_n=1, 
     378        theta_pd=15, theta_pd_n=45, 
     379        phi_pd=15, phi_pd_n=1, 
     380        psi_pd=15, psi_pd_n=1, 
     381        ) 
     382    return 'triaxial_ellipsoid', pars 
     383 
     384@model('sph') 
     385def sphere(): 
     386    pars = dict( 
     387        scale=1, background=0, 
     388        sld=6, solvent_sld=1, 
     389        radius=120, 
     390        radius_pd=.2, radius_pd_n=45, 
     391        ) 
     392    return 'sphere', pars 
     393 
     394@model('lam') 
     395def lamellar(): 
     396    pars = dict( 
     397        scale=1, background=0, 
     398        sld=6, solvent_sld=1, 
     399        thickness=40, 
     400        thickness_pd= 0.2, thickness_pd_n=40, 
     401        ) 
     402    return 'lamellar', pars 
    176403 
    177404if __name__ == "__main__": 
    178     coreshell(1) 
    179  
    180  
    181  
    182  
    183  
    184  
    185  
    186  
    187  
    188  
    189  
    190  
    191  
    192  
    193  
    194  
    195  
    196  
    197  
    198  
    199  
    200  
    201  
    202  
    203  
    204  
    205  
    206  
    207  
     405    main() 

  • fit.py

    r4001d6e r87985ca  
    22# -*- coding: utf-8 -*- 
    33 
     4import sys 
    45from bumps.names import * 
    5 from sasmodel import SasModel, load_data, set_beam_stop, set_half, set_top 
    6 from Models.code_capcyl import GpuCapCylinder 
    7 from Models.code_coreshellcyl_f import GpuCoreShellCylinder 
    8 from Models.code_cylinder_f import GpuCylinder 
    9 #from Models.code_cylinder import GpuCylinder, OneDGpuCylinder 
    10 from Models.code_ellipse_f import GpuEllipse 
    11 from Models.code_lamellar import GpuLamellar 
    12 from Models.code_triaxialellipse import GpuTriEllipse 
     6from sasmodels import bumps_model as sas 
    137 
    148""" IMPORT THE DATA USED """ 
    15 #data = load_data('December/Tangential/Sector0/DEC07133.ABS') 
    16 data = load_data('December/DEC07098.DAT') 
     9radial_data = sas.load_data('December/DEC07267.DAT') 
     10sas.set_beam_stop(radial_data, 0.00669, outer=0.025) 
     11sas.set_top(radial_data, -.0185) 
    1712 
    18 """ SET INNER BEAM STOP, OUTER RING, AND MASK HALF OF THE DATA """ 
    19 set_beam_stop(data, 0.004)#, outer=0.025) 
    20 set_top(data, -.018) 
    21 #set_half(data, 'right') 
     13tan_data = sas.load_data('December/DEC07266.DAT') 
     14sas.set_beam_stop(tan_data, 0.00669, outer=0.025) 
     15sas.set_top(tan_data, -.0185) 
     16#sas.set_half(tan_data, 'right') 
    2217 
     18name = "ellipsoid" if len(sys.argv) < 2 else sys.argv[1] 
     19section = "radial" if len(sys.argv) < 3 else sys.argv[2] 
     20data = radial_data if section is not "tangent" else tan_data 
     21kernel = sas.load_model(name, dtype="single") 
    2322 
    24 if 0: 
    25     model = SasModel(data, OneDGpuCylinder, 
    26     scale=0.0013, radius=105, length=1000, 
    27     background=21, sldCyl=.291e-6,sldSolv=7.105e-6, 
    28     radius_pd=0.1,radius_pd_n=10,radius_pd_nsigma=0, 
    29     length_pd=0.1,length_pd_n=5,length_pd_nsigma=0, 
    30     bolim=0.0, uplim=90) #bottom limit, upper limit of angle integral 
    31  
    32  
    33 if 0: 
    34     model = SasModel(data, GpuEllipse, 
     23if name == "ellipsoid": 
     24    model = sas.BumpsModel(data, kernel, 
    3525    scale=0.08, 
    36     radius_a=15, radius_b=800, 
    37     sldEll=.291e-6, sldSolv=7.105e-6, 
     26    rpolar=15, requatorial=800, 
     27    sld=.291, solvent_sld=7.105, 
    3828    background=0, 
    39     axis_theta=90, axis_phi=0, 
    40     axis_theta_pd=15, axis_theta_pd_n=40, axis_theta_pd_nsigma=3, 
    41     radius_a_pd=0.222296, radius_a_pd_n=1, radius_a_pd_nsigma=0, 
    42     radius_b_pd=.000128, radius_b_pd_n=1, radius_b_pd_nsigma=0, 
    43     axis_phi_pd=2.63698e-05, axis_phi_pd_n=20, axis_phi_pd_nsigma=3, 
    44     dtype='float32') 
     29    theta=90, phi=0, 
     30    theta_pd=15, theta_pd_n=40, theta_pd_nsigma=3, 
     31    rpolar_pd=0.222296, rpolar_pd_n=1, rpolar_pd_nsigma=0, 
     32    requatorial_pd=.000128, requatorial_pd_n=1, requatorial_pd_nsigma=0, 
     33    phi_pd=0, phi_pd_n=20, phi_pd_nsigma=3, 
     34    ) 
    4535 
    4636 
    4737    # SET THE FITTING PARAMETERS 
    48     #model.radius_a.range(15, 1000) 
    49     #model.radius_b.range(15, 1000) 
    50     #model.axis_theta_pd.range(0, 360) 
     38    #model.rpolar.range(15, 1000) 
     39    #model.requatorial.range(15, 1000) 
     40    #model.theta_pd.range(0, 360) 
    5141    #model.background.range(0,1000) 
    5242    model.scale.range(0, 1) 
    5343 
    5444 
    55 if 0: 
    56     model = SasModel(data, GpuLamellar, 
     45elif name == "lamellar": 
     46    model = sas.BumpsModel(data, kernel, 
    5747    scale=0.08, 
    58     bi_thick=19.2946, 
    59     sld_bi=5.38e-6,sld_sol=7.105e-6, 
     48    thickness=19.2946, 
     49    sld=5.38,sld_sol=7.105, 
    6050    background=0.003, 
    61     bi_thick_pd= 0.37765, bi_thick_pd_n=10, bi_thick_pd_nsigma=3, 
    62     dtype='float32') 
     51    thickness_pd= 0.37765, thickness_pd_n=10, thickness_pd_nsigma=3, 
     52    ) 
    6353 
    6454    # SET THE FITTING PARAMETERS 
    65     #model.bi_thick.range(0, 1000) 
     55    #model.thickness.range(0, 1000) 
    6656    #model.scale.range(0, 1) 
    67     #model.bi_thick_pd.range(0, 1000) 
     57    #model.thickness_pd.range(0, 1000) 
    6858    #model.background.range(0, 1000) 
    69     model.sld_bi.range(0, 1) 
     59    model.sld.range(0, 1) 
    7060 
    7161 
    72 if 0: 
     62elif name == "cylinder": 
    7363    """ 
    7464    pars = dict(scale=0.0023, radius=92.5, length=798.3, 
    75         sldCyl=.29e-6, sldSolv=7.105e-6, background=5, 
    76         cyl_theta=0, cyl_phi=0, 
    77         cyl_theta_pd=22.11, cyl_theta_pd_n=5, cyl_theta_pd_nsigma=3, 
     65        sld=.29, solvent_sld=7.105, background=5, 
     66        theta=0, phi=0, 
     67        theta_pd=22.11, theta_pd_n=5, theta_pd_nsigma=3, 
    7868        radius_pd=.0084, radius_pd_n=10, radius_pd_nsigma=3, 
    7969        length_pd=0.493, length_pd_n=10, length_pd_nsigma=3, 
    80         cyl_phi_pd=0, cyl_phi_pd_n=5, cyl_phi_pd_nsigma=3,) 
     70        phi_pd=0, phi_pd_n=5, phi_pd_nsigma=3,) 
    8171        """ 
    8272    pars = dict( 
    83         scale=.01, radius=64.1, length=66.96, sldCyl=.291e-6, sldSolv=5.77e-6, background=.1, 
    84         cyl_theta=90, cyl_phi=0, 
     73        scale=.01, radius=64.1, length=66.96, sld=.291, solvent_sld=5.77, background=.1, 
     74        theta=90, phi=0, 
    8575        radius_pd=0.1, radius_pd_n=10, radius_pd_nsigma=3, 
    8676        length_pd=0.1,length_pd_n=5, length_pd_nsigma=3, 
    87         cyl_theta_pd=0.1, cyl_theta_pd_n=5, cyl_theta_pd_nsigma=3, 
    88         cyl_phi_pd=0.1, cyl_phi_pd_n=10, cyl_phi_pd_nsigma=3) 
    89     model = SasModel(data, GpuCylinder, dtype="float32", **pars) 
     77        theta_pd=0.1, theta_pd_n=5, theta_pd_nsigma=3, 
     78        phi_pd=0.1, phi_pd_n=10, phi_pd_nsigma=3) 
     79    model = sas.BumpsModel(data, kernel, **pars) 
    9080 
    9181 
    … …  
    9484    model.radius.range(1, 500) 
    9585    model.length.range(1, 5000) 
    96     #model.cyl_theta.range(-90,100) 
    97     #model.cyl_theta_pd.range(0, 90) 
    98     #model.cyl_theta_pd_n = model.cyl_theta_pd + 5 
     86    #model.theta.range(-90,100) 
     87    #model.theta_pd.range(0, 90) 
     88    #model.theta_pd_n = model.theta_pd + 5 
    9989    #model.radius_pd.range(0, 90) 
    10090    #model.length_pd.range(0, 90) 
    10191    model.scale.range(0, 1) 
    10292    #model.background.range(0, 100) 
    103     #model.sldCyl.range(0, 1) 
     93    #model.sld.range(0, 1) 
    10494 
    10595 
    106 if 1: 
    107     model = SasModel(data, GpuCoreShellCylinder, 
     96elif name == "core_shell_cylinder": 
     97    model = sas.BumpsModel(data, kernel, 
    10898    scale= .031, radius=19.5, thickness=30, length=22, 
    109     core_sld=7.105e-6, shell_sld=.291e-6, solvent_sld=7.105e-6, 
    110     background=0, axis_theta=0, axis_phi=0, 
     99    core_sld=7.105, shell_sld=.291, solvent_sld=7.105, 
     100    background=0, theta=0, phi=0, 
    111101 
    112102    radius_pd=0.26, radius_pd_n=10, radius_pd_nsigma=3, 
    113103    length_pd=0.26, length_pd_n=10, length_pd_nsigma=3, 
    114104    thickness_pd=1, thickness_pd_n=1, thickness_pd_nsigma=1, 
    115     axis_theta_pd=1, axis_theta_pd_n=10, axis_theta_pd_nsigma=3, 
    116     axis_phi_pd=0.1, axis_phi_pd_n=1, axis_phi_pd_nsigma=1, 
    117     dtype='float32') 
     105    theta_pd=1, theta_pd_n=10, theta_pd_nsigma=3, 
     106    phi_pd=0.1, phi_pd_n=1, phi_pd_nsigma=1, 
     107    ) 
    118108 
    119109    # SET THE FITTING PARAMETERS 
    … …  
    122112    #model.thickness.range(18, 38) 
    123113    #model.thickness_pd.range(0, 1) 
    124     #model.axis_phi.range(0, 90) 
     114    #model.phi.range(0, 90) 
    125115    #model.radius_pd.range(0, 1) 
    126116    #model.length_pd.range(0, 1) 
    127     #model.axis_theta_pd.range(0, 360) 
     117    #model.theta_pd.range(0, 360) 
    128118    #model.background.range(0,5) 
    129119    model.scale.range(0, 1) 
    … …  
    131121 
    132122 
    133 if 0: 
    134     model = SasModel(data, GpuCapCylinder, 
    135     scale=.08, rad_cyl=20, rad_cap=40, len_cyl=400, 
    136     sld_capcyl=1e-6, sld_solv=6.3e-6, 
     123elif name == "capped_cylinder": 
     124    model = sas.BumpsModel(data, kernel, 
     125    scale=.08, radius=20, cap_radius=40, length=400, 
     126    sld_capcyl=1, sld_solv=6.3, 
    137127    background=0, theta=0, phi=0, 
    138     rad_cyl_pd=.1, rad_cyl_pd_n=5, rad_cyl_pd_nsigma=3, 
    139     rad_cap_pd=.1, rad_cap_pd_n=5, rad_cap_pd_nsigma=3, 
    140     len_cyl_pd=.1, len_cyl_pd_n=1, len_cyl_pd_nsigma=0, 
     128    radius_pd=.1, radius_pd_n=5, radius_pd_nsigma=3, 
     129    cap_radius_pd=.1, cap_radius_pd_n=5, cap_radius_pd_nsigma=3, 
     130    length_pd=.1, length_pd_n=1, length_pd_nsigma=0, 
    141131    theta_pd=.1, theta_pd_n=1, theta_pd_nsigma=0, 
    142132    phi_pd=.1, phi_pd_n=1, phi_pd_nsigma=0, 
    143     dtype='float32') 
     133    dtype=dtype) 
    144134 
    145135    model.scale.range(0, 1) 
    146136 
    147137 
    148 if 0: 
    149     model = SasModel(data, GpuTriEllipse, 
    150     scale=0.08, axisA=15, axisB=20, axisC=500, 
    151     sldEll=7.105e-6, sldSolv=.291e-6, 
     138elif name == "triaxial_ellipsoid": 
     139    model = sas.BumpsModel(data, kernel, 
     140    scale=0.08, req_minor=15, req_major=20, rpolar=500, 
     141    sldEll=7.105, solvent_sld=.291, 
    152142    background=5, theta=0, phi=0, psi=0, 
    153143    theta_pd=20, theta_pd_n=40, theta_pd_nsigma=3, 
    154144    phi_pd=.1, phi_pd_n=1, phi_pd_nsigma=0, 
    155145    psi_pd=30, psi_pd_n=1, psi_pd_nsigma=0, 
    156     axisA_pd=.1, axisA_pd_n=1, axisA_pd_nsigma=0, 
    157     axisB_pd=.1, axisB_pd_n=1, axisB_pd_nsigma=0, 
    158     axisC_pd=.1, axisC_pd_n=1, axisC_pd_nsigma=0, dtype='float32') 
     146    req_minor_pd=.1, req_minor_pd_n=1, req_minor_pd_nsigma=0, 
     147    req_major_pd=.1, req_major_pd_n=1, req_major_pd_nsigma=0, 
     148    rpolar_pd=.1, rpolar_pd_n=1, rpolar_pd_nsigma=0, dtype=dtype) 
    159149 
    160150    # SET THE FITTING PARAMETERS 
    161     model.axisA.range(15, 1000) 
    162     model.axisB.range(15, 1000) 
    163     #model.axisC.range(15, 1000) 
     151    model.req_minor.range(15, 1000) 
     152    model.req_major.range(15, 1000) 
     153    #model.rpolar.range(15, 1000) 
    164154    #model.background.range(0,1000) 
    165155    #model.theta_pd.range(0, 360) 
    … …  
    167157    #model.psi_pd.range(0, 360) 
    168158 
     159else: 
     160    print "No parameters for %s"%name 
     161    sys.exit(1) 
    169162 
     163if section is "both": 
     164   tan_model = sas.BumpsModel(tan_data, model.model, model.parameters()) 
     165   tan_model.phi = model.phi - 90 
     166   problem = FitProblem([model, tan_model]) 
     167else: 
     168   problem = FitProblem(model) 
    170169 
    171  
    172 problem = FitProblem(model) 
    173  

  • sasmodels/bumps_model.py

    r5d4777d r87985ca  
    66 
    77import numpy as np 
     8 
     9try: 
     10    from .gpu import load_model as _loader 
     11except ImportError,exc: 
     12    warnings.warn(str(exc)) 
     13    warnings.warn("OpenCL not available --- using ctypes instead") 
     14    from .dll import load_model as _loader 
     15 
     16def load_model(modelname, dtype='single'): 
     17    """ 
     18    Load model by name. 
     19    """ 
     20    sasmodels = __import__('sasmodels.models.'+modelname) 
     21    module = getattr(sasmodels.models, modelname, None) 
     22    model = _loader(module, dtype=dtype) 
     23    return model 
    824 
    925 
    … …  
    235251        partype = model.info['partype'] 
    236252 
     253        # remember inputs so we can inspect from outside 
     254        self.data = data 
     255        self.model = model 
     256 
    237257        # interpret data 
    238         self.data = data 
    239258        if hasattr(data, 'qx_data'): 
    240259            self.index = (data.mask==0) & (~np.isnan(data.data)) 

  • sasmodels/sasview_model.py

    r5d4777d r87985ca  
     1""" 
     2Sasview model constructor. 
     3 
     4Given a module defining an OpenCL kernel such as sasmodels.models.cylinder, 
     5create a sasview model class to run that kernel as follows:: 
     6 
     7    from sasmodels.sasview_model import make_class 
     8    from sasmodels.models import cylinder 
     9    CylinderModel = make_class(cylinder, dtype='single') 
     10 
     11The model parameters for sasmodels are different from those in sasview. 
     12When reloading previously saved models, the parameters should be converted 
     13using :func:`sasmodels.convert.convert`. 
     14""" 
     15 
     16# TODO: add a sasview=>sasmodels parameter translation layer 
     17# this will allow us to use the new sasmodels as drop in replacements, and 
     18# delay renaming parameters until all models have been converted. 
     19 
    120import math 
    221from copy import deepcopy 
    … …  
    1635    """ 
    1736    Load the sasview model defined in *kernel_module*. 
    18     :param kernel_module: 
    19     :param dtype: 
    20     :return: 
     37 
     38    Returns a class that can be used directly as a sasview model. 
    2139    """ 
    2240    model =  load_model(kernel_module, dtype=dtype) 

  • setup_py2exe.py

    r14de349 r87985ca  
    4040import sys 
    4141 
    42 sys.dont_write_bytecode = True 
     42#sys.dont_write_bytecode = True 
    4343 
    4444# Force build before continuing 
    … …  
    228228# Specify required packages to bundle in the executable image. 
    229229packages = ['numpy', 'scipy', 'matplotlib', 'pytz', 'pyparsing', 
    230             'periodictable', 'bumps.names', 'dream' 
     230            'periodictable', 'bumps', 'sasmodels', 'pyopencl', 
    231231            ] 
    232232 
Note: See TracChangeset for help on using the changeset viewer.