Changeset f4cf580 in sasmodels for compare-new.py


Ignore:
Timestamp:
Sep 2, 2014 11:15:40 AM (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:
87985ca
Parents:
5d4777d
Message:

resolve remaining differences between sasview and sasmodels

File:
1 edited

Legend:

Unmodified
Added
Removed

  • compare-new.py

    r5d4777d rf4cf580  
    33 
    44import sys 
     5import math 
    56 
    67import numpy as np 
    … …  
    1516    Load a sasview model given the model name. 
    1617    """ 
    17     modelname = modelname 
     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()) 
    1822    sans = __import__('sans.models.'+modelname) 
    1923    ModelClass = getattr(getattr(sans.models,modelname,None),modelname,None) 
    … …  
    4145    return kernel 
    4246 
    43 def load_dll(modelname, dtype='single'): 
     47def load_ctypes(modelname, dtype='single'): 
    4448    sasmodels = __import__('sasmodels.models.'+modelname) 
    4549    module = getattr(sasmodels.models, modelname, None) 
    … …  
    4953def randomize(p, v): 
    5054    """ 
    51     Randomizing pars using sasview names. 
    52     Guessing at angles and slds. 
    53     """ 
    54     if any(p.endswith(s) for s in ('_pd_n','_pd_nsigmas','_pd_type')): 
     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')): 
    5560        return v 
    56     if any(s in p for s in ('theta','phi','psi')): 
    57         return np.random.randint(0,90) 
    58     elif any(s in p for s in ('sld','Sld','SLD')): 
    59         return np.random.rand()*1e-5 
     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() 
    6073    else: 
    61         return np.random.rand() 
     74        # length, scale, background in [0,200] 
     75        return 200*np.random.rand() 
    6276 
    6377def parlist(pars): 
    … …  
    7488        if p.endswith("_pd"): pars[p] = 0 
    7589 
    76 def compare(gpuname, gpupars, Ncpu, Ngpu, opts): 
    77     #from sasmodels.core import load_data 
    78     #data = load_data('December/DEC07098.DAT') 
     90def compare(name, pars, Ncpu, Ngpu, opts, set_pars): 
     91 
     92    # Sort out data 
    7993    qmax = 1.0 if '-highq' in opts else (0.2 if '-midq' in opts else 0.05) 
    8094    if "-1d" in opts: 
    8195        from sasmodels.bumps_model import empty_data1D 
    82         qmax = np.log10(qmax) 
     96        qmax = math.log10(qmax) 
    8397        data = empty_data1D(np.logspace(qmax-3, qmax, 128)) 
     98        index = slice(None, None) 
    8499    else: 
    85100        from sasmodels.bumps_model import empty_data2D, set_beam_stop 
    86101        data = empty_data2D(np.linspace(-qmax, qmax, 128)) 
    87102        set_beam_stop(data, 0.004) 
     103        index = ~data.mask 
    88104    is2D = hasattr(data, 'qx_data') 
     105 
     106    # modelling accuracy is determined by dtype and cutoff 
    89107    dtype = 'double' if '-double' in opts else 'single' 
    90108    cutoff_opts = [s for s in opts if s.startswith('-cutoff')] 
    91109    cutoff = float(cutoff_opts[0].split('=')[1]) if cutoff_opts else 1e-5 
    92110 
    93  
    94     if '-random' in opts: 
    95         gpupars = dict((p,randomize(p,v)) for p,v in gpupars.items()) 
    96     if '-pars' in opts: print "pars",parlist(gpupars) 
    97     if '-mono' in opts: suppress_pd(gpupars) 
    98  
    99     cpuname, cpupars = revert_model(gpuname, gpupars) 
    100  
    101     try: 
    102         gpumodel = load_opencl(gpuname, dtype=dtype) 
    103     except Exception,exc: 
    104         print exc 
    105         print "... trying again with single precision" 
    106         gpumodel = load_opencl(gpuname, dtype='single') 
    107     model = BumpsModel(data, gpumodel, cutoff=cutoff, **gpupars) 
     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 
    108134    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) 
    109142        toc = tic() 
    110         for i in range(Ngpu): 
     143        for _ in range(Ngpu): 
    111144            #pars['scale'] = np.random.rand() 
    112             model.update() 
    113             gpu = model.theory() 
     145            problem.update() 
     146            gpu = problem.theory() 
    114147        gpu_time = toc()*1000./Ngpu 
    115         print "ocl t=%.1f ms, intensity=%f"%(gpu_time, sum(gpu[~np.isnan(gpu)])) 
     148        print "opencl t=%.1f ms, intensity=%.0f"%(gpu_time, sum(gpu[index])) 
    116149        #print max(gpu), min(gpu) 
    117150 
    118     comp = None 
    119     if Ncpu > 0 and "-dll" in opts: 
    120         dllmodel = load_dll(gpuname, dtype=dtype) 
    121         model = BumpsModel(data, dllmodel, cutoff=cutoff, **gpupars) 
     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) 
    122155        toc = tic() 
    123         for i in range(Ncpu): 
    124             model.update() 
    125             cpu = model.theory() 
     156        for _ in range(Ncpu): 
     157            problem.update() 
     158            cpu = problem.theory() 
    126159        cpu_time = toc()*1000./Ncpu 
    127         comp = "dll" 
    128  
    129     elif 0: # Hack to check new vs old for GpuCylinder 
    130         from Models.code_cylinder_f import GpuCylinder as oldgpu 
    131         from sasmodel import SasModel 
    132         oldmodel = SasModel(data, oldgpu, dtype=dtype, **cpupars) 
     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) 
    133164        toc = tic() 
    134         for i in range(Ngpu): 
    135             oldmodel.update() 
    136             cpu = oldmodel.theory() 
    137         cpu_time = toc()*1000./Ngpu 
    138         comp = "old" 
    139  
    140     elif Ncpu > 0: 
    141         cpumodel = sasview_model(cpuname, **cpupars) 
    142         toc = tic() 
    143         if is2D: 
    144             for i in range(Ncpu): 
    145                 cpu = cpumodel.evalDistribution([data.qx_data, data.qy_data]) 
    146         else: 
    147             for i in range(Ncpu): 
    148                 cpu = cpumodel.evalDistribution(data.x) 
     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) 
    149170        cpu_time = toc()*1000./Ncpu 
    150         comp = 'sasview' 
    151         #print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[model.index])) 
    152  
    153     if comp: 
    154         print "%s t=%.1f ms, intensity=%f"%(comp, cpu_time, sum(cpu[model.index])) 
     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 
    155175    if Ngpu > 0 and Ncpu > 0: 
    156176        #print "speedup %.2g"%(cpu_time/gpu_time) 
    … …  
    158178        #cpu *= max(gpu/cpu) 
    159179        resid, relerr = np.zeros_like(gpu), np.zeros_like(gpu) 
    160         resid[model.index] = (gpu - cpu)[model.index] 
    161         relerr[model.index] = resid[model.index]/cpu[model.index] 
    162         print "max(|ocl-%s|)"%comp, max(abs(resid[model.index])) 
    163         print "max(|(ocl-%s)/ocl|)"%comp, max(abs(relerr[model.index])) 
    164  
     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 
    165186    if '-noplot' in opts: return 
    166  
    167187    import matplotlib.pyplot as plt 
    168188    if Ncpu > 0: 
    … …  
    173193        if Ncpu > 0: plt.subplot(132) 
    174194        plot_data(data, gpu, scale='log') 
    175         plt.title("ocl t=%.1f ms"%gpu_time) 
     195        plt.title("opencl t=%.1f ms"%gpu_time) 
    176196    if Ncpu > 0 and Ngpu > 0: 
    177197        plt.subplot(133) 
    178198        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" 
    179201        plot_data(data, err, scale='linear') 
    180         plt.title("max rel err = %.3g"%max(abs(err[model.index]))) 
     202        plt.title("max %s = %.3g"%(errstr, max(abs(err[index])))) 
    181203        if is2D: plt.colorbar() 
    182204    plt.show() 
    … …  
    200222    -lowq/-midq/-highq use q values up to 0.05, 0.2 or 1.0 
    201223    -2d/-1d uses 1d or 2d random data 
    202     -preset/-random randomizes the parameters 
     224    -preset/-random[=seed] preset or random parameters 
    203225    -poly/-mono force monodisperse/polydisperse 
    204     -sasview/-dll whether cpu is tested using sasview or dll 
     226    -sasview/-ctypes whether cpu is tested using sasview or ctypes 
    205227    -cutoff=1e-5/value cutoff for including a point in polydispersity 
    206228    -nopars/-pars prints the parameter set or not 
    … …  
    214236    %s 
    215237""" 
     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    ] 
    216250 
    217251def main(): 
    … …  
    227261        sys.exit(1) 
    228262 
     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 
    229272    name, pars = MODELS[args[0]]() 
    230273    Nopencl = int(args[1]) if len(args) > 1 else 5 
    … …  
    238281 
    239282    # Fill in parameters given on the command line 
     283    set_pars = {} 
    240284    for arg in args[3:]: 
    241285        k,v = arg.split('=') 
    242286        if k not in pars: 
    243287            # extract base name without distribution 
    244             # style may be either a.d or a_pd_d 
    245             s = set((p.split('_pd')[0]).split('.')[0] for p in pars) 
     288            s = set(p.split('_pd')[0] for p in pars) 
    246289            print "%r invalid; parameters are: %s"%(k,", ".join(sorted(s))) 
    247290            sys.exit(1) 
    248         pars[k] = float(v) if not v.endswith('type') else v 
    249  
    250     compare(name, pars, Nsasview, Nopencl, opts) 
     291        set_pars[k] = float(v) if not v.endswith('type') else v 
     292 
     293    compare(name, pars, Nsasview, Nopencl, opts, set_pars) 
    251294 
    252295# =========================================================================== 
    … …  
    265308    pars = dict( 
    266309        scale=1, background=0, 
    267         sld=6e-6, solvent_sld=1e-6, 
     310        sld=6, solvent_sld=1, 
    268311        #radius=5, length=20, 
    269312        radius=260, length=290, 
    270         theta=30, phi=15, 
     313        theta=30, phi=0, 
    271314        radius_pd=.2, radius_pd_n=1, 
    272315        length_pd=.2,length_pd_n=1, 
    273         theta_pd=15, theta_pd_n=25, 
     316        theta_pd=15, theta_pd_n=45, 
    274317        phi_pd=15, phi_pd_n=1, 
    275318        ) 
    … …  
    280323    pars = dict( 
    281324        scale=1, background=0, 
    282         sld=6e-6, solvent_sld=1e-6, 
    283         radius=260, cap_radius=80000, length=290, 
     325        sld=6, solvent_sld=1, 
     326        radius=260, cap_radius=290, length=290, 
    284327        theta=30, phi=15, 
    285328        radius_pd=.2, radius_pd_n=1, 
    286329        cap_radius_pd=.2, cap_radius_pd_n=1, 
    287330        length_pd=.2, length_pd_n=1, 
    288         theta_pd=15, theta_pd_n=25, 
     331        theta_pd=15, theta_pd_n=45, 
    289332        phi_pd=15, phi_pd_n=1, 
    290333        ) 
    … …  
    296339    pars = dict( 
    297340        scale=1, background=0, 
    298         core_sld=6e-6, shell_sld=8e-6, solvent_sld=1e-6, 
    299         radius=325, thickness=25, length=34.2709, 
    300         theta=90, phi=0, 
    301         radius_pd=0.1, radius_pd_n=10, 
    302         length_pd=0.1, length_pd_n=5, 
    303         thickness_pd=0.1, thickness_pd_n=5, 
    304         theta_pd=15.8, theta_pd_n=5, 
    305         phi_pd=0.0008748, phi_pd_n=5, 
     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, 
    306349        ) 
    307350    return 'core_shell_cylinder', pars 
    … …  
    312355    pars = dict( 
    313356        scale=1, background=0, 
    314         sld=6e-6, solvent_sld=1e-6, 
     357        sld=6, solvent_sld=1, 
    315358        rpolar=50, requatorial=30, 
    316         theta=0, phi=0, 
    317         rpolar_pd=0.3, rpolar_pd_n=10, 
    318         requatorial_pd=0, requatorial_pd_n=10, 
    319         theta_pd=0, theta_pd_n=45, 
    320         phi_pd=0, phi_pd_n=45, 
     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, 
    321364        ) 
    322365    return 'ellipsoid', pars 
    … …  
    327370    pars = dict( 
    328371        scale=1, background=0, 
    329         sld=6e-6, solvent_sld=1e-6, 
    330         theta=0, phi=0, psi=0, 
     372        sld=6, solvent_sld=1, 
     373        theta=30, phi=15, psi=5, 
    331374        req_minor=25, req_major=36, rpolar=50, 
    332         theta_pd=0, theta_pd_n=5, 
    333         phi_pd=0, phi_pd_n=5, 
    334         psi_pd=0, psi_pd_n=5, 
    335         req_minor_pd=0, req_minor_pd_n=5, 
    336         req_major_pd=0, req_major_pd_n=5, 
    337         rpolar_pd=.3, rpolar_pd_n=25, 
     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, 
    338381        ) 
    339382    return 'triaxial_ellipsoid', pars 
    … …  
    343386    pars = dict( 
    344387        scale=1, background=0, 
    345         sld=6e-6, solvent_sld=1e-6, 
     388        sld=6, solvent_sld=1, 
    346389        radius=120, 
    347         radius_pd=.3, radius_pd_n=5, 
     390        radius_pd=.2, radius_pd_n=45, 
    348391        ) 
    349392    return 'sphere', pars 
    … …  
    353396    pars = dict( 
    354397        scale=1, background=0, 
    355         sld=6e-6, solvent_sld=1e-6, 
     398        sld=6, solvent_sld=1, 
    356399        thickness=40, 
    357         thickness_pd= 0.3, thickness_pd_n=40, 
     400        thickness_pd= 0.2, thickness_pd_n=40, 
    358401        ) 
    359402    return 'lamellar', pars 
Note: See TracChangeset for help on using the changeset viewer.