Changeset e78edc4 in sasmodels

Timestamp:

Apr 14, 2016 8:41:17 AM (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:

7abcc59

Parents:

62cf915

git-author:

Paul Kienzle <pkienzle@…> (04/14/16 08:37:03)

git-committer:

Paul Kienzle <pkienzle@…> (04/14/16 08:41:17)

Message:

add -zero option to test q=0

Location:

sasmodels

Files:

• compare.py (modified) (11 diffs)
• data.py (modified) (4 diffs)
• model_test.py (modified) (2 diffs)

sasmodels/compare.py

@@ -58 +58 @@
     -lowq*/-midq/-highq/-exq use q values up to 0.05, 0.2, 1.0, 10.0
     -nq=128 sets the number of Q points in the data set
+    -zero indicates that q=0 should be included
     -1d*/-2d computes 1d or 2d data
     -preset*/-random[=seed] preset or random parameters
@@ -476 +477 @@
         index = ~data.mask
     else:
-        if opts['view'] == 'log':
+        if opts['view'] == 'log' and not opts['zero']:
             qmax = math.log10(qmax)
             q = np.logspace(qmax-3, qmax, nq)
         else:
             q = np.linspace(0.001*qmax, qmax, nq)
+        if opts['zero']:
+            q = np.hstack((0, q))
         data = empty_data1D(q, resolution=res)
         index = slice(None, None)
@@ -498 +501 @@
     else:
         return eval_opencl(model_info, data, dtype=dtype, cutoff=cutoff)
+
+def _show_invalid(data, theory):
+    if not theory.mask.any():
+        return
+
+    if hasattr(data, 'x'):
+        bad = zip(data.x[theory.mask], theory[theory.mask])
+        print("   *** ", ", ".join("I(%g)=%g"%(x, y) for x,y in bad))
+
 
 def compare(opts, limits=None):
@@ -519 +531 @@
         try:
             base_value, base_time = time_calculation(base, pars, Nbase)
+            base_value = np.ma.masked_invalid(base_value)
             print("%s t=%.2f ms, intensity=%.0f"
-                  % (base.engine, base_time, sum(base_value)))
+                  % (base.engine, base_time, base_value.sum()))
+            _show_invalid(data, base_value)
         except ImportError:
             traceback.print_exc()
@@ -530 +544 @@
         try:
             comp_value, comp_time = time_calculation(comp, pars, Ncomp)
+            comp_value = np.ma.masked_invalid(comp_value)
             print("%s t=%.2f ms, intensity=%.0f"
-                  % (comp.engine, comp_time, sum(comp_value)))
+                  % (comp.engine, comp_time, comp_value.sum()))
+            _show_invalid(data, comp_value)
         except ImportError:
             traceback.print_exc()
@@ -554 +570 @@
         vmin, vmax = np.Inf, -np.Inf
         if Nbase > 0:
-            vmin = min(vmin, min(base_value))
-            vmax = max(vmax, max(base_value))
+            vmin = min(vmin, base_value.min())
+            vmax = max(vmax, base_value.max())
         if Ncomp > 0:
-            vmin = min(vmin, min(comp_value))
-            vmax = max(vmax, max(comp_value))
+            vmin = min(vmin, comp_value.min())
+            vmax = max(vmax, comp_value.max())
         limits = vmin, vmax
 
@@ -581 +597 @@
         if view == 'linear':
             plt.xscale('linear')
-        plt.title("max %s = %.3g"%(errstr, max(abs(err))))
+        plt.title("max %s = %.3g"%(errstr, abs(err).max()))
         #cbar_title = errstr if errview=="linear" else "log "+errstr
     #if is2D:
@@ -603 +619 @@
 
 def _print_stats(label, err):
-    sorted_err = np.sort(abs(err))
+    sorted_err = np.sort(abs(err.compressed()))
     p50 = int((len(err)-1)*0.50)
     p98 = int((len(err)-1)*0.98)
@@ -623 +639 @@
     'half', 'fast', 'single', 'double',
     'single!', 'double!', 'quad!', 'sasview',
-    'lowq', 'midq', 'highq', 'exq',
+    'lowq', 'midq', 'highq', 'exq', 'zero',
     '2d', '1d',
     'preset', 'random',
@@ -700 +716 @@
     try:
         model_info = core.load_model_info(name)
-    except ImportError, exc:
+    except ImportError as exc:
         print(str(exc))
         print("Could not find model; use one of:\n    " + models)
@@ -745 +761 @@
         elif arg == '-midq':    opts['qmax'] = 0.2
         elif arg == '-lowq':    opts['qmax'] = 0.05
+        elif arg == '-zero':    opts['zero'] = True
         elif arg.startswith('-nq='):       opts['nq'] = int(arg[4:])
         elif arg.startswith('-res='):      opts['res'] = float(arg[5:])

sasmodels/data.py

@@ -380 +380 @@
     use_calc = use_theory and Iq_calc is not None
     num_plots = (use_data or use_theory) + use_calc + use_resid
-
+    non_positive_x = (data.x<=0.0).any()
 
     scale = data.x**4 if view == 'q4' else 1.0
@@ -402 +402 @@
 
         if use_theory:
+            # Note: masks merge, so any masked theory points will stay masked,
+            # and the data mask will be added to it.
             mtheory = masked_array(theory, data.mask.copy())
             mtheory[~np.isfinite(mtheory)] = masked
@@ -413 +415 @@
             plt.ylim(*limits)
 
-        plt.xscale('linear' if not some_present else view)
+        plt.xscale('linear' if not some_present or non_positive_x  else view)
         plt.yscale('linear'
                    if view == 'q4' or not some_present or not all_positive
@@ -441 +443 @@
         plt.xlabel("$q$/A$^{-1}$")
         plt.ylabel('residuals')
-        plt.xscale('linear' if not some_present else view)
+        plt.xscale('linear' if not some_present or non_positive_x else view)
 
 

sasmodels/model_test.py

@@ -145 +145 @@
         def _runTest(self):
             smoke_tests = [
+                # test validity at reasonable values
                 [{}, 0.1, None],
                 [{}, (0.1, 0.1), None],
+                # test validity at q = 0
+                #[{}, 0.0, None],
+                #[{}, (0.0, 0.0), None],
+                # test that ER/VR will run if they exist
                 [{}, 'ER', None],
                 [{}, 'VR', None],
@@ -194 +199 @@
                 actual = call_kernel(kernel, pars)
 
-            self.assertGreater(len(actual), 0)
+            self.assertTrue(len(actual) > 0)
             self.assertEqual(len(y), len(actual))