Changeset 13d86bc in sasmodels

Timestamp:

Aug 26, 2014 7: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:

• compare-new.py (modified) (6 diffs)
• sasmodels/core.py (modified) (4 diffs)

compare-new.py

@@ -6 +6 @@
 import numpy as np
 
-from sasmodels.core import BumpsModel, fake_data2D, set_beam_stop, plot_data, \
-    tic, opencl_model, dll_model
+from sasmodels.core import BumpsModel, plot_data, tic, opencl_model, dll_model
 
 def sasview_model(modelname, **pars):
@@ -27 +26 @@
         elif k.endswith("_pd_nsigma"):
             model.dispersion[k[:-10]]['nsigmas'] = v
+        elif k.endswith("_pd_type"):
+            model.dispersion[k[:-8]]['type'] = v
         else:
             model.setParam(k, v)
     return model
 
+def load_opencl(modelname, dtype='single'):
+    sasmodels = __import__('sasmodels.models.'+modelname)
+    module = getattr(sasmodels.models, modelname, None)
+    kernel = opencl_model(module, dtype=dtype)
+    return kernel
+
+def load_dll(modelname, dtype='single'):
+    sasmodels = __import__('sasmodels.models.'+modelname)
+    module = getattr(sasmodels.models, modelname, None)
+    kernel = dll_model(module, dtype=dtype)
+    return kernel
+
 
 def compare(Ncpu, cpuname, cpupars, Ngpu, gpuname, gpupars):
@@ -36 +49 @@
     #from sasmodels.core import load_data
     #data = load_data('December/DEC07098.DAT')
-    data = fake_data2D(np.linspace(-0.05, 0.05, 128))
-    set_beam_stop(data, 0.004)
+    from sasmodels.core import empty_data1D
+    data = empty_data1D(np.logspace(-4, -1, 128))
+    #from sasmodels.core import empty_2D, set_beam_stop
+    #data = empty_data2D(np.linspace(-0.05, 0.05, 128))
+    #set_beam_stop(data, 0.004)
+    is2D = hasattr(data, 'qx_data')
 
     if Ngpu > 0:
-        gpumodel = opencl_model(gpuname, dtype='single')
+        gpumodel = load_opencl(gpuname, dtype='single')
         model = BumpsModel(data, gpumodel, **gpupars)
         toc = tic()
@@ -52 +69 @@
 
     if 0 and Ncpu > 0: # Hack to compare ctypes vs. opencl
-        dllmodel = dll_model(gpuname)
+        dllmodel = load_dll(gpuname, dtype='double')
         model = BumpsModel(data, dllmodel, **gpupars)
         toc = tic()
@@ -75 +92 @@
         cpumodel = sasview_model(cpuname, **cpupars)
         toc = tic()
-        for i in range(Ncpu):
-            cpu = cpumodel.evalDistribution([data.qx_data, data.qy_data])
+        if is2D:
+            for i in range(Ncpu):
+                cpu = cpumodel.evalDistribution([data.qx_data, data.qy_data])
+        else:
+            for i in range(Ncpu):
+                cpu = cpumodel.evalDistribution(data.x)
         cpu_time = toc()*1000./Ncpu
         print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[model.index]))
@@ -92 +113 @@
     if Ncpu > 0:
         if Ngpu > 0: plt.subplot(131)
-        plot_data(data, cpu)
+        plot_data(data, cpu, scale='log')
         plt.title("omp t=%.1f ms"%cpu_time)
     if Ngpu > 0:
         if Ncpu > 0: plt.subplot(132)
-        plot_data(data, gpu)
+        plot_data(data, gpu, scale='log')
         plt.title("ocl t=%.1f ms"%gpu_time)
     if Ncpu > 0 and Ngpu > 0:
         plt.subplot(133)
-        plot_data(data, 1e8*relerr)
+        plot_data(data, 1e8*relerr, scale='linear')
         plt.title("max rel err = %.3g"%max(abs(relerr)))
-        plt.colorbar()
+        if is2D: plt.colorbar()
     plt.show()
 

sasmodels/core.py

@@ -73 +73 @@
 
 
-def fake_data2D(qx, qy=None):
+def empty_data2D(qx, qy=None):
     from sans.dataloader.data_info import Data2D, Detector
 
@@ -114 +114 @@
 
 
+def empty_data1D(q):
+    from sans.dataloader.data_info import Data1D
+
+    Iq = 100*np.ones_like(q)
+    dIq = np.sqrt(Iq)
+    data = Data1D(q, Iq, dx=0.05*q, dy=dIq)
+    data.filename = "fake data"
+    data.qmin, data.qmax = q.min(), q.max()
+    return data
+
+
 def set_beam_stop(data, radius, outer=None):
     from sans.dataloader.manipulations import Ringcut
@@ -139 +150 @@
 
 
-def plot_data(data, iq, vmin=None, vmax=None):
-    from numpy.ma import masked_array
+def plot_data(data, iq, vmin=None, vmax=None, scale='log'):
+    from numpy.ma import masked_array, masked
     import matplotlib.pyplot as plt
-    img = masked_array(iq, data.mask)
-    xmin, xmax = min(data.qx_data), max(data.qx_data)
-    ymin, ymax = min(data.qy_data), max(data.qy_data)
-    plt.imshow(img.reshape(128,128),
-               interpolation='nearest', aspect=1, origin='upper',
-               extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax)
-
-
-def plot_result2D(data, theory, view='linear'):
+    if hasattr(data, 'qx_data'):
+        img = masked_array(iq, data.mask)
+        if scale == 'log':
+            img[(img <= 0) | ~np.isfinite(img)] = masked
+            img = np.log10(img)
+        xmin, xmax = min(data.qx_data), max(data.qx_data)
+        ymin, ymax = min(data.qy_data), max(data.qy_data)
+        plt.imshow(img.reshape(128,128),
+                   interpolation='nearest', aspect=1, origin='upper',
+                   extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax)
+    else: # 1D data
+        if scale == 'linear':
+            idx = np.isfinite(iq)
+            plt.plot(data.x[idx], iq[idx])
+        else:
+            idx = np.isfinite(iq) & (iq>0)
+            plt.loglog(data.x[idx], iq[idx])
+
+
+def plot_result2D(data, theory, view='log'):
     import matplotlib.pyplot as plt
-    from numpy.ma import masked_array, masked
-    #print "not a number",sum(np.isnan(data.data))
-    #data.data[data.data<0.05] = 0.5
-    mdata = masked_array(data.data, data.mask)
-    mdata[np.isnan(mdata)] = masked
-    if view is 'log':
-        mdata[mdata <= 0] = masked
-        mdata = np.log10(mdata)
-        mtheory = masked_array(np.log10(theory), mdata.mask)
-    else:
-        mtheory = masked_array(theory, mdata.mask)
-    mresid = masked_array((theory-data.data)/data.err_data, data.mask)
-    vmin = min(mdata.min(), mtheory.min())
-    vmax = max(mdata.max(), mtheory.max())
-    print np.exp(np.mean(mtheory)), np.std(mtheory),np.max(mtheory),np.min(mtheory)
-
-    plt.subplot(1, 3, 1)
-    plot_data(data, mdata, vmin=vmin, vmax=vmax)
+    resid = (theory-data.data)/data.err_data
+    plt.subplot(131)
+    plot_data(data, data.data, scale=view)
     plt.colorbar()
-    plt.subplot(1, 3, 2)
-    plot_data(data, mtheory, vmin=vmin, vmax=vmax)
+    plt.subplot(132)
+    plot_data(data, theory, scale=view)
     plt.colorbar()
-    plt.subplot(1, 3, 3)
-    print abs(mresid).max()
-    plot_data(data, mresid)
+    plt.subplot(133)
+    plot_data(data, resid, scale='linear')
     plt.colorbar()
 
 
-def plot_result1D(data, theory, view='linear'):
+def plot_result1D(data, theory, view='log'):
     import matplotlib.pyplot as plt
     from numpy.ma import masked_array, masked
@@ -218 +224 @@
             q_vectors = [data.qx_data, data.qy_data]
         else:
-            self.index = (data.mask==0) & (~np.isnan(data.y))
+            self.index = (data.x>=data.qmin) & (data.x<=data.qmax) & ~np.isnan(data.y)
             self.iq = data.y[self.index]
             self.diq = data.dy[self.index]