Changeset e3571cb in sasmodels

Timestamp:

Oct 21, 2017 10:57:43 PM (7 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

2bccb5a

Parents:

6773b02

Message:

allow comparison of 1D with integrated 2D

Files:

• explore/check1d.py (added)
• sasmodels/compare.py (modified) (16 diffs)
• sasmodels/data.py (modified) (3 diffs)
• sasmodels/modelinfo.py (modified) (1 diff)

sasmodels/compare.py

@@ -83 +83 @@
     -pars/-nopars* prints the parameter set or not
     -default/-demo* use demo vs default parameters
+    -sphere[=150] set up spherical integration over theta/phi using n points
 
     === calculation options ===
-    -mono*/-poly force monodisperse or allow polydisperse demo parameters
+    -mono*/-poly force monodisperse or allow polydisperse random parameters
     -cutoff=1e-5* cutoff value for including a point in polydispersity
     -magnetic/-nonmagnetic* suppress magnetism
@@ -279 +280 @@
         # orientation in [-180,180], orientation pd in [0,45]
         if p.endswith('_pd'):
-            return 0., 45.
+            return 0., 180.
         else:
             return -180., 180.
@@ -434 +435 @@
 
 
+def limit_dimensions(model_info, pars, maxdim):
+    # type: (ModelInfo, ParameterSet, float) -> None
+    """
+    Limit parameters of units of Ang to maxdim.
+    """
+    for p in model_info.parameters.call_parameters:
+        value = pars[p.name]
+        if p.units == 'Ang' and value > maxdim:
+            pars[p.name] = maxdim*10**np.random.uniform(-3,0)
+
 def constrain_pars(model_info, pars):
     # type: (ModelInfo, ParameterSet) -> None
@@ -496 +507 @@
     Format the parameter list for printing.
     """
+    is2d = True
     lines = []
     parameters = model_info.parameters
@@ -817 +829 @@
 
 
-def compare(opts, limits=None):
+def compare(opts, limits=None, maxdim=np.inf):
     # type: (Dict[str, Any], Optional[Tuple[float, float]]) -> Tuple[float, float]
     """
@@ -828 +840 @@
     as the values are adjusted, making it easier to see the effects of the
     parameters.
+
+    *maxdim* is the maximum value for any parameter with units of Angstrom.
     """
     for k in range(opts['sets']):
-        opts['pars'] = parse_pars(opts)
+        if k > 1:
+            # print a separate seed for each dataset for better reproducibility
+            new_seed = np.random.randint(1000000)
+            print("Set %d uses -random=%i"%(k+1,new_seed))
+            np.random.seed(new_seed)
+        opts['pars'] = parse_pars(opts, maxdim=maxdim)
         if opts['pars'] is None:
             return
@@ -971 +990 @@
         #err,errstr = base/comp,"ratio"
         plot_theory(data, None, resid=err, view=errview, use_data=use_data)
-        if view == 'linear':
-            plt.xscale('linear')
+        plt.xscale('log' if view == 'log' else linear)
+        plt.legend(['P%d'%(k+1) for k in range(setnum+1)], loc='best')
         plt.title("max %s = %.3g"%(errstr, abs(err).max()))
         #cbar_title = errstr if errview=="linear" else "log "+errstr
@@ -1019 +1038 @@
     'demo', 'default',  # TODO: remove demo/default
     'nopars', 'pars',
+    'sphere', 'sphere=', # integrate over a sphere in 2d with n points
 
     # Calculation options
@@ -1153 +1173 @@
         'sets'      : 0,
         'engine'    : 'default',
-        'evals'     : '1',
+        'count'     : '1',
         'show_weights' : False,
+        'sphere'    : 0,
     }
     for arg in flags:
@@ -1179 +1200 @@
         elif arg.startswith('-data='):     opts['datafile'] = arg[6:]
         elif arg.startswith('-engine='):   opts['engine'] = arg[8:]
-        elif arg.startswith('-neval='):    opts['evals'] = arg[7:]
+        elif arg.startswith('-neval='):    opts['count'] = arg[7:]
+        elif arg.startswith('-sphere'):
+            opts['sphere'] = int(arg[8:]) if len(arg) > 7 else 150
+            opts['is2d'] = True
         elif arg == '-random':  opts['seed'] = np.random.randint(1000000)
         elif arg == '-preset':  opts['seed'] = -1
@@ -1238 +1262 @@
 
     if PAR_SPLIT in opts['engine']:
-        engine_types = opts['engine'].split(PAR_SPLIT, 2)
+        opts['engine'] = opts['engine'].split(PAR_SPLIT, 2)
         comparison = True
     else:
-        engine_types = [opts['engine']]*2
-
-    if PAR_SPLIT in opts['evals']:
-        evals = [int(k) for k in opts['evals'].split(PAR_SPLIT, 2)]
+        opts['engine'] = [opts['engine']]*2
+
+    if PAR_SPLIT in opts['count']:
+        opts['count'] = [int(k) for k in opts['count'].split(PAR_SPLIT, 2)]
         comparison = True
     else:
-        evals = [int(opts['evals'])]*2
+        opts['count'] = [int(opts['count'])]*2
 
     if PAR_SPLIT in opts['cutoff']:
-        cutoff = [float(k) for k in opts['cutoff'].split(PAR_SPLIT, 2)]
+        opts['cutoff'] = [float(k) for k in opts['cutoff'].split(PAR_SPLIT, 2)]
         comparison = True
     else:
-        cutoff = [float(opts['cutoff'])]*2
-
-    base = make_engine(model_info[0], data, engine_types[0], cutoff[0])
+        opts['cutoff'] = [float(opts['cutoff'])]*2
+
+    base = make_engine(model_info[0], data, opts['engine'][0], opts['cutoff'][0])
     if comparison:
-        comp = make_engine(model_info[1], data, engine_types[1], cutoff[1])
+        comp = make_engine(model_info[1], data, opts['engine'][1], opts['cutoff'][1])
     else:
         comp = None
@@ -1266 +1290 @@
         'data'      : data,
         'name'      : names,
-        'def'       : model_info,
-        'count'     : evals,
+        'info'      : model_info,
         'engines'   : [base, comp],
         'values'    : values,
@@ -1273 +1296 @@
     # pylint: enable=bad-whitespace
 
+    # Set the integration parameters to the half sphere
+    if opts['sphere'] > 0:
+        set_spherical_integration_parameters(opts, opts['sphere'])
+
     return opts
 
-def parse_pars(opts):
-    model_info, model_info2 = opts['def']
+def set_spherical_integration_parameters(opts, steps):
+    """
+    Set integration parameters for spherical integration over the entire
+    surface in theta-phi coordinates.
+    """
+    # Set the integration parameters to the half sphere
+    opts['values'].extend([
+        'theta=90',
+        'theta_pd=%g'%(90/np.sqrt(3)),
+        'theta_pd_n=%d'%steps,
+        'theta_pd_type=rectangle',
+        'phi=0',
+        'phi_pd=%g'%(180/np.sqrt(3)),
+        'phi_pd_n=%d'%(2*steps),
+        'phi_pd_type=rectangle',
+        #'background=0',
+    ])
+    if 'psi' in opts['info'][0].parameters:
+        opts['values'].append('psi=0')
+
+def parse_pars(opts, maxdim=np.inf):
+    model_info, model_info2 = opts['info']
 
     # Get demo parameters from model definition, or use default parameters
@@ -1287 +1334 @@
     if opts['seed'] > -1:
         pars = randomize_pars(model_info, pars)
+        limit_dimensions(model_info, pars, maxdim)
         if model_info != model_info2:
             pars2 = randomize_pars(model_info2, pars2)
+            limit_dimensions(model_info, pars2, maxdim)
             # Share values for parameters with the same name
             for k, v in pars.items():
@@ -1365 +1414 @@
     from . import rst2html
 
-    info = opts['def'][0]
+    info = opts['info'][0]
     html = make_html(info)
     path = os.path.dirname(info.filename)
@@ -1416 +1465 @@
         opts['pars'] = list(opts['pars'])
         p1, p2 = opts['pars']
-        m1, m2 = opts['def']
+        m1, m2 = opts['info']
         self.fix_p2 = m1 != m2 or p1 != p2
         model_info = m1

sasmodels/data.py

@@ -363 +363 @@
     if hasattr(data, 'isSesans') and data.isSesans:
         _plot_result_sesans(data, None, None, use_data=True, limits=limits)
-    elif hasattr(data, 'qx_data'):
+    elif hasattr(data, 'qx_data') and not getattr(data, 'radial', False):
         _plot_result2D(data, None, None, view, use_data=True, limits=limits)
     else:
@@ -391 +391 @@
     if hasattr(data, 'isSesans') and data.isSesans:
         _plot_result_sesans(data, theory, resid, use_data=True, limits=limits)
-    elif hasattr(data, 'qx_data'):
+    elif hasattr(data, 'qx_data') and not getattr(data, 'radial', False):
         _plot_result2D(data, theory, resid, view, use_data, limits=limits)
     else:
@@ -425 +425 @@
     import matplotlib.pyplot as plt  # type: ignore
     from numpy.ma import masked_array, masked  # type: ignore
+
+    if getattr(data, 'radial', False):
+        radial_data.x = radial_data.q_data
+        radial_data.y = radial_data.data
 
     use_data = use_data and data.y is not None

sasmodels/modelinfo.py

@@ -502 +502 @@
             raise KeyError("unknown parameter %r"%key)
         return par
+
+    def __contains__(self, key):
+        for par in self.call_parameters:
+            if par.name == key:
+                return True
+        else:
+            return False
 
     def _set_vector_lengths(self):