Changeset 216a9e1 in sasmodels for compare.py

Timestamp:

Dec 3, 2014 11:38:49 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:

8fc26f8

Parents:

31819c5

Message:

batch generate random old/new comparisons across all models

File:

• compare.py (modified) (13 diffs)

compare.py

@@ -75 +75 @@
         return 200*np.random.rand()
 
+def randomize_model(name, pars, seed=None):
+    if seed is None:
+        seed = np.random.randint(1e9)
+    np.random.seed(seed)
+    # Note: the sort guarantees order of calls to random number generator
+    pars = dict((p,randomize(p,v)) for p,v in sorted(pars.items()))
+    # The capped cylinder model has a constraint on its parameters
+    if name == 'capped_cylinder' and pars['cap_radius'] < pars['radius']:
+        pars['radius'],pars['cap_radius'] = pars['cap_radius'],pars['radius']
+    return pars, seed
+
 def parlist(pars):
     return "\n".join("%s: %s"%(p,v) for p,v in sorted(pars.items()))
@@ -88 +99 @@
         if p.endswith("_pd"): pars[p] = 0
 
+def eval_sasview(name, pars, data, Nevals=1):
+    model = sasview_model(name, **pars)
+    toc = tic()
+    for _ in range(Nevals):
+        if hasattr(data, 'qx_data'):
+            value = model.evalDistribution([data.qx_data, data.qy_data])
+        else:
+            value = model.evalDistribution(data.x)
+    average_time = toc()*1000./Nevals
+    return value, average_time
+
+def eval_opencl(name, pars, data, dtype='single', Nevals=1, cutoff=0):
+    try:
+        model = load_opencl(name, dtype=dtype)
+    except Exception,exc:
+        print exc
+        print "... trying again with single precision"
+        model = load_opencl(name, dtype='single')
+    problem = BumpsModel(data, model, cutoff=cutoff, **pars)
+    toc = tic()
+    for _ in range(Nevals):
+        #pars['scale'] = np.random.rand()
+        problem.update()
+        value = problem.theory()
+    average_time = toc()*1000./Nevals
+    return value, average_time
+
+def eval_ctypes(name, pars, data, dtype='double', Nevals=1, cutoff=0):
+    model = load_ctypes(name, dtype=dtype)
+    problem = BumpsModel(data, model, cutoff=cutoff, **pars)
+    toc = tic()
+    for _ in range(Nevals):
+        problem.update()
+        value = problem.theory()
+    average_time = toc()*1000./Nevals
+    return value, average_time
+
+def make_data(qmax, is2D, Nq=128):
+    if is2D:
+        from sasmodels.bumps_model import empty_data2D, set_beam_stop
+        data = empty_data2D(np.linspace(-qmax, qmax, Nq))
+        set_beam_stop(data, 0.004)
+        index = ~data.mask
+    else:
+        from sasmodels.bumps_model import empty_data1D
+        qmax = math.log10(qmax)
+        data = empty_data1D(np.logspace(qmax-3, qmax, Nq))
+        index = slice(None, None)
+    return data, index
+
 def compare(name, pars, Ncpu, Ngpu, opts, set_pars):
+    opt_values = dict(split
+                      for s in opts for split in ((s.split('='),))
+                      if len(split) == 2)
     # Sort out data
     qmax = 1.0 if '-highq' in opts else (0.2 if '-midq' in opts else 0.05)
-    if "-1d" in opts:
-        from sasmodels.bumps_model import empty_data1D
-        qmax = math.log10(qmax)
-        data = empty_data1D(np.logspace(qmax-3, qmax, 128))
-        index = slice(None, None)
-    else:
-        from sasmodels.bumps_model import empty_data2D, set_beam_stop
-        data = empty_data2D(np.linspace(-qmax, qmax, 128))
-        set_beam_stop(data, 0.004)
-        index = ~data.mask
-    is2D = hasattr(data, 'qx_data')
+    Nq = int(opt_values.get('-Nq', '128'))
+    is2D = not "-1d" in opts
+    data, index = make_data(qmax, is2D, Nq)
+
 
     # modelling accuracy is determined by dtype and cutoff
     dtype = 'double' if '-double' in opts else 'single'
-    cutoff_opts = [s for s in opts if s.startswith('-cutoff')]
-    cutoff = float(cutoff_opts[0].split('=')[1]) if cutoff_opts else 1e-5
+    cutoff = float(opt_values.get('-cutoff','1e-5'))
 
     # randomize parameters
-    random_opts = [s for s in opts if s.startswith('-random')]
-    if random_opts:
-        if '=' in random_opts[0]:
-            seed = int(random_opts[0].split('=')[1])
-        else:
-            seed = int(np.random.rand()*10000)
-        np.random.seed(seed)
+    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
-        # Note: the sort guarantees order of calls to random number generator
-        pars = dict((p,randomize(p,v)) for p,v in sorted(pars.items()))
-        # The capped cylinder model has a constraint on its parameters
-        if name == 'capped_cylinder' and pars['cap_radius'] < pars['radius']:
-            pars['radius'],pars['cap_radius'] = pars['cap_radius'],pars['radius']
     pars.update(set_pars)
 
@@ -132 +179 @@
     # OpenCl calculation
     if Ngpu > 0:
-        try:
-            model = load_opencl(name, dtype=dtype)
-        except Exception,exc:
-            print exc
-            print "... trying again with single precision"
-            model = load_opencl(name, dtype='single')
-        problem = BumpsModel(data, model, cutoff=cutoff, **pars)
-        toc = tic()
-        for _ in range(Ngpu):
-            #pars['scale'] = np.random.rand()
-            problem.update()
-            gpu = problem.theory()
-        gpu_time = toc()*1000./Ngpu
+        gpu, gpu_time = eval_opencl(name, pars, data, dtype, Ngpu)
         print "opencl t=%.1f ms, intensity=%.0f"%(gpu_time, sum(gpu[index]))
         #print max(gpu), min(gpu)
@@ -150 +185 @@
     # ctypes/sasview calculation
     if Ncpu > 0 and "-ctypes" in opts:
-        model = load_ctypes(name, dtype=dtype)
-        problem = BumpsModel(data, model, cutoff=cutoff, **pars)
-        toc = tic()
-        for _ in range(Ncpu):
-            problem.update()
-            cpu = problem.theory()
-        cpu_time = toc()*1000./Ncpu
+        cpu, cpu_time = eval_ctypes(name, pars, data, dtype=dtype, cutoff=cutoff, Nevals=Ncpu)
         comp = "ctypes"
         print "ctypes t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[index]))
     elif Ncpu > 0:
-        model = sasview_model(name, **pars)
-        toc = tic()
-        for _ in range(Ncpu):
-            if is2D:
-                cpu = model.evalDistribution([data.qx_data, data.qy_data])
-            else:
-                cpu = model.evalDistribution(data.x)
-        cpu_time = toc()*1000./Ncpu
+        cpu, cpu_time = eval_sasview(name, pars, data, Ncpu)
         comp = "sasview"
         print "sasview t=%.1f ms, intensity=%.0f"%(cpu_time, sum(cpu[index]))
@@ -235 +257 @@
     -single*/-double uses double precision for comparison
     -lowq*/-midq/-highq use q values up to 0.05, 0.2 or 1.0
-    -1d/-2d* uses 1d or 2d random data
+    -Nq=128 sets the number of Q points in the data set
+    -1d/-2d* computes 1d or 2d data
     -preset*/-random[=seed] preset or random parameters
     -mono/-poly* force monodisperse/polydisperse
@@ -252 +275 @@
 """
 
-VALID_OPTIONS = [
+NAME_OPTIONS = set([
     'plot','noplot',
     'single','double',
@@ -263 +286 @@
     'rel','abs',
     'hist','nohist',
+    ])
+VALUE_OPTIONS = [
+    # Note: random is both a name option and a value option
+    'cutoff', 'random', 'Nq',
     ]
 
@@ -277 +304 @@
         sys.exit(1)
 
-    valid_opts = set(VALID_OPTIONS)
     invalid = [o[1:] for o in opts
-               if o[1:] not in valid_opts
-                  and not o.startswith('-cutoff=')
-                  and not o.startswith('-random=')]
+               if o[1:] not in NAME_OPTIONS
+                  and not any(o.startswith('-%s='%t) for t in VALUE_OPTIONS)]
     if invalid:
         print "Invalid options: %s"%(", ".join(invalid))
@@ -328 +353 @@
         radius=260, length=290,
         theta=30, phi=0,
-        radius_pd=.2, radius_pd_n=1,
-        length_pd=.2,length_pd_n=1,
+        radius_pd=.2, radius_pd_n=9,
+        length_pd=.2,length_pd_n=10,
         theta_pd=15, theta_pd_n=45,
         phi_pd=15, phi_pd_n=1,
@@ -344 +369 @@
         radius_pd=.2, radius_pd_n=1,
         cap_radius_pd=.2, cap_radius_pd_n=1,
-        length_pd=.2, length_pd_n=1,
+        length_pd=.2, length_pd_n=10,
         theta_pd=15, theta_pd_n=45,
         phi_pd=15, phi_pd_n=1,
@@ -359 +384 @@
         theta=30, phi=15,
         radius_pd=.2, radius_pd_n=1,
-        length_pd=.2, length_pd_n=1,
-        thickness_pd=.2, thickness_pd_n=1,
+        length_pd=.2, length_pd_n=10,
+        thickness_pd=.2, thickness_pd_n=10,
         theta_pd=15, theta_pd_n=45,
         phi_pd=15, phi_pd_n=1,
@@ -374 +399 @@
         rpolar=50, requatorial=30,
         theta=30, phi=15,
-        rpolar_pd=.2, rpolar_pd_n=1,
-        requatorial_pd=.2, requatorial_pd_n=1,
+        rpolar_pd=.2, rpolar_pd_n=15,
+        requatorial_pd=.2, requatorial_pd_n=15,
         theta_pd=15, theta_pd_n=45,
         phi_pd=15, phi_pd_n=1,
@@ -391 +416 @@
         req_minor_pd=0, req_minor_pd_n=1,
         req_major_pd=0, req_major_pd_n=1,
-        rpolar_pd=.2, rpolar_pd_n=1,
+        rpolar_pd=.2, rpolar_pd_n=30,
         theta_pd=15, theta_pd_n=45,
         phi_pd=15, phi_pd_n=1,