Changeset 13d86bc in sasmodels


Ignore:
Timestamp:
Aug 26, 2014 9:06:58 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:
a7684e5
Parents:
32c160a
Message:

1D model comparison with sasview

Files:
2 edited

Legend:

Unmodified
Added
Removed
  • compare-new.py

    rce27e21 r13d86bc  
    66import numpy as np 
    77 
    8 from sasmodels.core import BumpsModel, fake_data2D, set_beam_stop, plot_data, \ 
    9     tic, opencl_model, dll_model 
     8from sasmodels.core import BumpsModel, plot_data, tic, opencl_model, dll_model 
    109 
    1110def sasview_model(modelname, **pars): 
     
    2726        elif k.endswith("_pd_nsigma"): 
    2827            model.dispersion[k[:-10]]['nsigmas'] = v 
     28        elif k.endswith("_pd_type"): 
     29            model.dispersion[k[:-8]]['type'] = v 
    2930        else: 
    3031            model.setParam(k, v) 
    3132    return model 
    3233 
     34def load_opencl(modelname, dtype='single'): 
     35    sasmodels = __import__('sasmodels.models.'+modelname) 
     36    module = getattr(sasmodels.models, modelname, None) 
     37    kernel = opencl_model(module, dtype=dtype) 
     38    return kernel 
     39 
     40def load_dll(modelname, dtype='single'): 
     41    sasmodels = __import__('sasmodels.models.'+modelname) 
     42    module = getattr(sasmodels.models, modelname, None) 
     43    kernel = dll_model(module, dtype=dtype) 
     44    return kernel 
     45 
    3346 
    3447def compare(Ncpu, cpuname, cpupars, Ngpu, gpuname, gpupars): 
     
    3649    #from sasmodels.core import load_data 
    3750    #data = load_data('December/DEC07098.DAT') 
    38     data = fake_data2D(np.linspace(-0.05, 0.05, 128)) 
    39     set_beam_stop(data, 0.004) 
     51    from sasmodels.core import empty_data1D 
     52    data = empty_data1D(np.logspace(-4, -1, 128)) 
     53    #from sasmodels.core import empty_2D, set_beam_stop 
     54    #data = empty_data2D(np.linspace(-0.05, 0.05, 128)) 
     55    #set_beam_stop(data, 0.004) 
     56    is2D = hasattr(data, 'qx_data') 
    4057 
    4158    if Ngpu > 0: 
    42         gpumodel = opencl_model(gpuname, dtype='single') 
     59        gpumodel = load_opencl(gpuname, dtype='single') 
    4360        model = BumpsModel(data, gpumodel, **gpupars) 
    4461        toc = tic() 
     
    5269 
    5370    if 0 and Ncpu > 0: # Hack to compare ctypes vs. opencl 
    54         dllmodel = dll_model(gpuname) 
     71        dllmodel = load_dll(gpuname, dtype='double') 
    5572        model = BumpsModel(data, dllmodel, **gpupars) 
    5673        toc = tic() 
     
    7592        cpumodel = sasview_model(cpuname, **cpupars) 
    7693        toc = tic() 
    77         for i in range(Ncpu): 
    78             cpu = cpumodel.evalDistribution([data.qx_data, data.qy_data]) 
     94        if is2D: 
     95            for i in range(Ncpu): 
     96                cpu = cpumodel.evalDistribution([data.qx_data, data.qy_data]) 
     97        else: 
     98            for i in range(Ncpu): 
     99                cpu = cpumodel.evalDistribution(data.x) 
    79100        cpu_time = toc()*1000./Ncpu 
    80101        print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[model.index])) 
     
    92113    if Ncpu > 0: 
    93114        if Ngpu > 0: plt.subplot(131) 
    94         plot_data(data, cpu) 
     115        plot_data(data, cpu, scale='log') 
    95116        plt.title("omp t=%.1f ms"%cpu_time) 
    96117    if Ngpu > 0: 
    97118        if Ncpu > 0: plt.subplot(132) 
    98         plot_data(data, gpu) 
     119        plot_data(data, gpu, scale='log') 
    99120        plt.title("ocl t=%.1f ms"%gpu_time) 
    100121    if Ncpu > 0 and Ngpu > 0: 
    101122        plt.subplot(133) 
    102         plot_data(data, 1e8*relerr) 
     123        plot_data(data, 1e8*relerr, scale='linear') 
    103124        plt.title("max rel err = %.3g"%max(abs(relerr))) 
    104         plt.colorbar() 
     125        if is2D: plt.colorbar() 
    105126    plt.show() 
    106127 
  • sasmodels/core.py

    r32c160a r13d86bc  
    7373 
    7474 
    75 def fake_data2D(qx, qy=None): 
     75def empty_data2D(qx, qy=None): 
    7676    from sans.dataloader.data_info import Data2D, Detector 
    7777 
     
    114114 
    115115 
     116def empty_data1D(q): 
     117    from sans.dataloader.data_info import Data1D 
     118 
     119    Iq = 100*np.ones_like(q) 
     120    dIq = np.sqrt(Iq) 
     121    data = Data1D(q, Iq, dx=0.05*q, dy=dIq) 
     122    data.filename = "fake data" 
     123    data.qmin, data.qmax = q.min(), q.max() 
     124    return data 
     125 
     126 
    116127def set_beam_stop(data, radius, outer=None): 
    117128    from sans.dataloader.manipulations import Ringcut 
     
    139150 
    140151 
    141 def plot_data(data, iq, vmin=None, vmax=None): 
    142     from numpy.ma import masked_array 
     152def plot_data(data, iq, vmin=None, vmax=None, scale='log'): 
     153    from numpy.ma import masked_array, masked 
    143154    import matplotlib.pyplot as plt 
    144     img = masked_array(iq, data.mask) 
    145     xmin, xmax = min(data.qx_data), max(data.qx_data) 
    146     ymin, ymax = min(data.qy_data), max(data.qy_data) 
    147     plt.imshow(img.reshape(128,128), 
    148                interpolation='nearest', aspect=1, origin='upper', 
    149                extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax) 
    150  
    151  
    152 def plot_result2D(data, theory, view='linear'): 
     155    if hasattr(data, 'qx_data'): 
     156        img = masked_array(iq, data.mask) 
     157        if scale == 'log': 
     158            img[(img <= 0) | ~np.isfinite(img)] = masked 
     159            img = np.log10(img) 
     160        xmin, xmax = min(data.qx_data), max(data.qx_data) 
     161        ymin, ymax = min(data.qy_data), max(data.qy_data) 
     162        plt.imshow(img.reshape(128,128), 
     163                   interpolation='nearest', aspect=1, origin='upper', 
     164                   extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax) 
     165    else: # 1D data 
     166        if scale == 'linear': 
     167            idx = np.isfinite(iq) 
     168            plt.plot(data.x[idx], iq[idx]) 
     169        else: 
     170            idx = np.isfinite(iq) & (iq>0) 
     171            plt.loglog(data.x[idx], iq[idx]) 
     172 
     173 
     174def plot_result2D(data, theory, view='log'): 
    153175    import matplotlib.pyplot as plt 
    154     from numpy.ma import masked_array, masked 
    155     #print "not a number",sum(np.isnan(data.data)) 
    156     #data.data[data.data<0.05] = 0.5 
    157     mdata = masked_array(data.data, data.mask) 
    158     mdata[np.isnan(mdata)] = masked 
    159     if view is 'log': 
    160         mdata[mdata <= 0] = masked 
    161         mdata = np.log10(mdata) 
    162         mtheory = masked_array(np.log10(theory), mdata.mask) 
    163     else: 
    164         mtheory = masked_array(theory, mdata.mask) 
    165     mresid = masked_array((theory-data.data)/data.err_data, data.mask) 
    166     vmin = min(mdata.min(), mtheory.min()) 
    167     vmax = max(mdata.max(), mtheory.max()) 
    168     print np.exp(np.mean(mtheory)), np.std(mtheory),np.max(mtheory),np.min(mtheory) 
    169  
    170     plt.subplot(1, 3, 1) 
    171     plot_data(data, mdata, vmin=vmin, vmax=vmax) 
     176    resid = (theory-data.data)/data.err_data 
     177    plt.subplot(131) 
     178    plot_data(data, data.data, scale=view) 
    172179    plt.colorbar() 
    173     plt.subplot(1, 3, 2) 
    174     plot_data(data, mtheory, vmin=vmin, vmax=vmax) 
     180    plt.subplot(132) 
     181    plot_data(data, theory, scale=view) 
    175182    plt.colorbar() 
    176     plt.subplot(1, 3, 3) 
    177     print abs(mresid).max() 
    178     plot_data(data, mresid) 
     183    plt.subplot(133) 
     184    plot_data(data, resid, scale='linear') 
    179185    plt.colorbar() 
    180186 
    181187 
    182 def plot_result1D(data, theory, view='linear'): 
     188def plot_result1D(data, theory, view='log'): 
    183189    import matplotlib.pyplot as plt 
    184190    from numpy.ma import masked_array, masked 
     
    218224            q_vectors = [data.qx_data, data.qy_data] 
    219225        else: 
    220             self.index = (data.mask==0) & (~np.isnan(data.y)) 
     226            self.index = (data.x>=data.qmin) & (data.x<=data.qmax) & ~np.isnan(data.y) 
    221227            self.iq = data.y[self.index] 
    222228            self.diq = data.dy[self.index] 
Note: See TracChangeset for help on using the changeset viewer.