Changeset 7cf2cfd in sasmodels for compare.py

Timestamp:

Nov 22, 2015 9:37:15 PM (8 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:

3b4243d

Parents:

677ccf1

Message:

refactor compare.py so that bumps/sasview not required for simple tests

File:

• compare.py (modified) (13 diffs)

compare.py

@@ -6 +6 @@
 from os.path import basename, dirname, join as joinpath
 import glob
+import datetime
 
 import numpy as np
@@ -13 +14 @@
 
 
-from sasmodels.bumps_model import Model, Experiment, plot_theory, tic
 from sasmodels import core
 from sasmodels import kerneldll
+from sasmodels.data import plot_theory, empty_data1D, empty_data2D
+from sasmodels.direct_model import DirectModel
 from sasmodels.convert import revert_model
 kerneldll.ALLOW_SINGLE_PRECISION_DLLS = True
@@ -22 +24 @@
 MODELS = [basename(f)[:-3]
           for f in sorted(glob.glob(joinpath(ROOT,"sasmodels","models","[a-zA-Z]*.py")))]
+
+# CRUFT python 2.6
+if not hasattr(datetime.timedelta, 'total_seconds'):
+    def delay(dt):
+        """Return number date-time delta as number seconds"""
+        return dt.days * 86400 + dt.seconds + 1e-6 * dt.microseconds
+else:
+    def delay(dt):
+        """Return number date-time delta as number seconds"""
+        return dt.total_seconds()
+
+
+def tic():
+    """
+    Timer function.
+
+    Use "toc=tic()" to start the clock and "toc()" to measure
+    a time interval.
+    """
+    then = datetime.datetime.now()
+    return lambda: delay(datetime.datetime.now() - then)
+
+
+def set_beam_stop(data, radius, outer=None):
+    """
+    Add a beam stop of the given *radius*.  If *outer*, make an annulus.
+
+    Note: this function does not use the sasview package
+    """
+    if hasattr(data, 'qx_data'):
+        q = np.sqrt(data.qx_data**2 + data.qy_data**2)
+        data.mask = (q < radius)
+        if outer is not None:
+            data.mask |= (q >= outer)
+    else:
+        data.mask = (data.x < radius)
+        if outer is not None:
+            data.mask |= (data.x >= outer)
 
 
@@ -99 +139 @@
         if p.endswith("_pd"): pars[p] = 0
 
-def eval_sasview(name, pars, data, Nevals=1):
+def eval_sasview(model_definition, pars, data, Nevals=1):
     from sas.models.qsmearing import smear_selection
-    model = sasview_model(name, **pars)
+    model = sasview_model(model_definition, **pars)
     smearer = smear_selection(data, model=model)
     value = None  # silence the linter
@@ -131 +171 @@
         print "... trying again with single precision"
         model = core.load_model(model_definition, dtype='single', platform="ocl")
-    problem = Experiment(data, Model(model, **pars), cutoff=cutoff)
+    calculator = DirectModel(data, model, cutoff=cutoff)
     value = None  # silence the linter
     toc = tic()
     for _ in range(max(Nevals, 1)):  # force at least one eval
-        #pars['scale'] = np.random.rand()
-        problem.update()
-        value = problem.theory()
+        value = calculator(**pars)
     average_time = toc()*1000./Nevals
     return value, average_time
 
+
 def eval_ctypes(model_definition, pars, data, dtype='double', Nevals=1, cutoff=0.):
     model = core.load_model(model_definition, dtype=dtype, platform="dll")
-    problem = Experiment(data, Model(model, **pars), cutoff=cutoff)
+    calculator = DirectModel(data, model, cutoff=cutoff)
     value = None  # silence the linter
     toc = tic()
     for _ in range(max(Nevals, 1)):  # force at least one eval
-        problem.update()
-        value = problem.theory()
+        value = calculator(**pars)
     average_time = toc()*1000./Nevals
     return value, average_time
 
+
 def make_data(qmax, is2D, Nq=128, resolution=0.0, accuracy='Low', view='log'):
     if is2D:
-        from sasmodels.bumps_model import empty_data2D, set_beam_stop
         data = empty_data2D(np.linspace(-qmax, qmax, Nq), resolution=resolution)
         data.accuracy = accuracy
@@ -160 +198 @@
         index = ~data.mask
     else:
-        from sasmodels.bumps_model import empty_data1D
         if view == 'log':
             qmax = math.log10(qmax)
@@ -190 +227 @@
 
     # randomize parameters
-    pars.update(set_pars)
+    #pars.update(set_pars)  # set value before random to control range
     if '-random' in opts or '-random' in opt_values:
         seed = int(opt_values['-random']) if '-random' in opt_values else None
         pars, seed = randomize_model(name, pars, seed=seed)
         print "Randomize using -random=%i"%seed
+    pars.update(set_pars)  # set value after random to control value
 
     # parameter selection
@@ -217 +255 @@
         print "ctypes t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu))
     elif Ncpu > 0:
-        cpu, cpu_time = eval_sasview(model_definition, pars, data, Ncpu)
-        comp = "sasview"
-        print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu))
+        try:
+            cpu, cpu_time = eval_sasview(model_definition, pars, data, Ncpu)
+            comp = "sasview"
+            #print "ocl/sasview", (ocl-pars['background'])/(cpu-pars['background'])
+            print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu))
+        except ImportError:
+            Ncpu = 0
 
     # Compare, but only if computing both forms
@@ -238 +280 @@
     if Ncpu > 0:
         if Nocl > 0: plt.subplot(131)
-        plot_theory(data, cpu, view=view)
+        plot_theory(data, cpu, view=view, plot_data=False)
         plt.title("%s t=%.1f ms"%(comp,cpu_time))
-        cbar_title = "log I"
+        #cbar_title = "log I"
     if Nocl > 0:
         if Ncpu > 0: plt.subplot(132)
-        plot_theory(data, ocl, view=view)
+        plot_theory(data, ocl, view=view, plot_data=False)
         plt.title("opencl t=%.1f ms"%ocl_time)
-        cbar_title = "log I"
+        #cbar_title = "log I"
     if Ncpu > 0 and Nocl > 0:
         plt.subplot(133)
@@ -253 +295 @@
             err,errstr,errview = abs(relerr), "rel err", "log"
         #err,errstr = ocl/cpu,"ratio"
-        plot_theory(data, err, view=errview)
+        plot_theory(data, None, resid=err, view=errview, plot_data=False)
         plt.title("max %s = %.3g"%(errstr, max(abs(err))))
-        cbar_title = errstr if errview=="linear" else "log "+errstr
-    if is2D:
-        h = plt.colorbar()
-        h.ax.set_title(cbar_title)
+        #cbar_title = errstr if errview=="linear" else "log "+errstr
+    #if is2D:
+    #    h = plt.colorbar()
+    #    h.ax.set_title(cbar_title)
 
     if Ncpu > 0 and Nocl > 0 and '-hist' in opts:
@@ -320 +362 @@
 
 Available models:
-
-    %s
 """
+
 
 NAME_OPTIONS = set([
@@ -342 +383 @@
     ]
 
+def columnize(L, indent="", width=79):
+    column_width = max(len(w) for w in L) + 1
+    num_columns = (width - len(indent)) // column_width
+    num_rows = len(L) // num_columns
+    L = L + [""] * (num_rows*num_columns - len(L))
+    columns = [L[k*num_rows:(k+1)*num_rows] for k in range(num_columns)]
+    lines = [" ".join("%-*s"%(column_width, entry) for entry in row)
+             for row in zip(*columns)]
+    output = indent + ("\n"+indent).join(lines)
+    return output
+
+
 def get_demo_pars(name):
     import sasmodels.models
@@ -355 +408 @@
     models = "\n    ".join("%-15s"%v for v in MODELS)
     if len(args) == 0:
-        print(USAGE%models)
+        print(USAGE)
+        print(columnize(MODELS, indent="  "))
         sys.exit(1)
     if args[0] not in MODELS: