Changeset e589e9a in sasmodels for sasmodels

Timestamp:

Mar 5, 2019 10:36:34 AM (5 years ago)

Author:

GitHub <noreply@…>

Branches:

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

Children:

3b6567f, bd91e8f, c11d09f, 9150036

Parents:

a8a1d48 (diff), 8064d5e (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

git-author:

Paul Butler <butlerpd@…> (03/05/19 10:36:34)

git-committer:

GitHub <noreply@…> (03/05/19 10:36:34)

Message:

Merge pull request #86 from SasView?/py3

fix test warnings from py3 and numpy

as agreed at today's call

Location:

sasmodels

Files:

• jitter.py (modified) (2 diffs)
• kernelcl.py (modified) (3 diffs)
• resolution.py (modified) (2 diffs)
• weights.py (modified) (1 diff)
• compare.py (modified) (1 diff)
• direct_model.py (modified) (1 diff)
• generate.py (modified) (3 diffs)
• kernel.py (modified) (1 diff)
• modelinfo.py (modified) (6 diffs)
• models/hardsphere.py (modified) (1 diff)
• models/pearl_necklace.c (modified) (2 diffs)
• sasview_model.py (modified) (7 diffs)

sasmodels/jitter.py

@@ -15 +15 @@
     pass
 
+import matplotlib as mpl
 import matplotlib.pyplot as plt
 from matplotlib.widgets import Slider
@@ -746 +747 @@
         pass
 
-    axcolor = 'lightgoldenrodyellow'
+    # CRUFT: use axisbg instead of facecolor for matplotlib<2
+    facecolor_prop = 'facecolor' if mpl.__version__ > '2' else 'axisbg'
+    props = {facecolor_prop: 'lightgoldenrodyellow'}
 
     ## add control widgets to plot
-    axes_theta = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor)
-    axes_phi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor)
-    axes_psi = plt.axes([0.1, 0.05, 0.45, 0.04], axisbg=axcolor)
+    axes_theta = plt.axes([0.1, 0.15, 0.45, 0.04], **props)
+    axes_phi = plt.axes([0.1, 0.1, 0.45, 0.04], **props)
+    axes_psi = plt.axes([0.1, 0.05, 0.45, 0.04], **props)
     stheta = Slider(axes_theta, 'Theta', -90, 90, valinit=theta)
     sphi = Slider(axes_phi, 'Phi', -180, 180, valinit=phi)
     spsi = Slider(axes_psi, 'Psi', -180, 180, valinit=psi)
 
-    axes_dtheta = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor)
-    axes_dphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor)
-    axes_dpsi = plt.axes([0.75, 0.05, 0.15, 0.04], axisbg=axcolor)
+    axes_dtheta = plt.axes([0.75, 0.15, 0.15, 0.04], **props)
+    axes_dphi = plt.axes([0.75, 0.1, 0.15, 0.04], **props)
+    axes_dpsi = plt.axes([0.75, 0.05, 0.15, 0.04], **props)
     # Note: using ridiculous definition of rectangle distribution, whose width
     # in sasmodels is sqrt(3) times the given width.  Divide by sqrt(3) to keep

sasmodels/kernelcl.py

@@ -58 +58 @@
 import time
 
+try:
+    from time import perf_counter as clock
+except ImportError: # CRUFT: python < 3.3
+    import sys
+    if sys.platform.count("darwin") > 0:
+        from time import time as clock
+    else:
+        from time import clock
+
 import numpy as np  # type: ignore
-
 
 # Attempt to setup opencl. This may fail if the pyopencl package is not
@@ -611 +619 @@
         #Call kernel and retrieve results
         wait_for = None
-        last_nap = time.clock()
+        last_nap = clock()
         step = 1000000//self.q_input.nq + 1
         for start in range(0, call_details.num_eval, step):
@@ -622 +630 @@
                 # Allow other processes to run
                 wait_for[0].wait()
-                current_time = time.clock()
+                current_time = clock()
                 if current_time - last_nap > 0.5:
                     time.sleep(0.001)

sasmodels/resolution.py

@@ -445 +445 @@
     q = np.sort(q)
     if q_min + 2*MINIMUM_RESOLUTION < q[0]:
-        n_low = np.ceil((q[0]-q_min) / (q[1]-q[0])) if q[1] > q[0] else 15
+        n_low = int(np.ceil((q[0]-q_min) / (q[1]-q[0]))) if q[1] > q[0] else 15
         q_low = np.linspace(q_min, q[0], n_low+1)[:-1]
     else:
         q_low = []
     if q_max - 2*MINIMUM_RESOLUTION > q[-1]:
-        n_high = np.ceil((q_max-q[-1]) / (q[-1]-q[-2])) if q[-1] > q[-2] else 15
+        n_high = int(np.ceil((q_max-q[-1]) / (q[-1]-q[-2]))) if q[-1] > q[-2] else 15
         q_high = np.linspace(q[-1], q_max, n_high+1)[1:]
     else:
@@ -499 +499 @@
             q_min = q[0]*MINIMUM_ABSOLUTE_Q
         n_low = log_delta_q * (log(q[0])-log(q_min))
-        q_low = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:-1]
+        q_low = np.logspace(log10(q_min), log10(q[0]), int(np.ceil(n_low))+1)[:-1]
     else:
         q_low = []
     if q_max > q[-1]:
         n_high = log_delta_q * (log(q_max)-log(q[-1]))
-        q_high = np.logspace(log10(q[-1]), log10(q_max), np.ceil(n_high)+1)[1:]
+        q_high = np.logspace(log10(q[-1]), log10(q_max), int(np.ceil(n_high))+1)[1:]
     else:
         q_high = []

sasmodels/weights.py

@@ -23 +23 @@
     default = dict(npts=35, width=0, nsigmas=3)
     def __init__(self, npts=None, width=None, nsigmas=None):
-        self.npts = self.default['npts'] if npts is None else npts
+        self.npts = self.default['npts'] if npts is None else int(npts)
         self.width = self.default['width'] if width is None else width
         self.nsigmas = self.default['nsigmas'] if nsigmas is None else nsigmas

sasmodels/compare.py

@@ -1152 +1152 @@
         'rel_err'   : True,
         'explore'   : False,
-        'use_demo'  : True,
+        'use_demo'  : False,
         'zero'      : False,
         'html'      : False,

sasmodels/direct_model.py

@@ -332 +332 @@
 
         # Need to pull background out of resolution for multiple scattering
-        background = pars.get('background', DEFAULT_BACKGROUND)
+        default_background = self._model.info.parameters.common_parameters[1].default
+        background = pars.get('background', default_background)
         pars = pars.copy()
         pars['background'] = 0.

sasmodels/generate.py

@@ -703 +703 @@
     """
     for code in source:
-        m = _FQ_PATTERN.search(code)
-        if m is not None:
+        if _FQ_PATTERN.search(code) is not None:
             return True
     return False
@@ -712 +711 @@
     # type: (List[str]) -> bool
     """
-    Return True if C source defines "void Fq(".
+    Return True if C source defines "double shell_volume(".
     """
     for code in source:
-        m = _SHELL_VOLUME_PATTERN.search(code)
-        if m is not None:
+        if _SHELL_VOLUME_PATTERN.search(code) is not None:
             return True
     return False
@@ -1008 +1006 @@
         pars = model_info.parameters.kernel_parameters
     else:
-        pars = model_info.parameters.COMMON + model_info.parameters.kernel_parameters
+        pars = (model_info.parameters.common_parameters
+                + model_info.parameters.kernel_parameters)
     partable = make_partable(pars)
     subst = dict(id=model_info.id.replace('_', '-'),

sasmodels/kernel.py

@@ -133 +133 @@
         nout = 2 if self.info.have_Fq and self.dim == '1d' else 1
         total_weight = self.result[nout*self.q_input.nq + 0]
+        # Note: total_weight = sum(weight > cutoff), with cutoff >= 0, so it
+        # is okay to test directly against zero.  If weight is zero then I(q),
+        # etc. must also be zero.
         if total_weight == 0.:
             total_weight = 1.

sasmodels/modelinfo.py

@@ -404 +404 @@
       parameters counted as n individual parameters p1, p2, ...
 
+    * *common_parameters* is the list of common parameters, with a unique
+      copy for each model so that structure factors can have a default
+      background of 0.0.
+
     * *call_parameters* is the complete list of parameters to the kernel,
       including scale and background, with vector parameters recorded as
@@ -416 +420 @@
     parameters don't use vector notation, and instead use p1, p2, ...
     """
-    # scale and background are implicit parameters
-    COMMON = [Parameter(*p) for p in COMMON_PARAMETERS]
-
     def __init__(self, parameters):
         # type: (List[Parameter]) -> None
+
+        # scale and background are implicit parameters
+        # Need them to be unique to each model in case they have different
+        # properties, such as default=0.0 for structure factor backgrounds.
+        self.common_parameters = [Parameter(*p) for p in COMMON_PARAMETERS]
+
         self.kernel_parameters = parameters
         self._set_vector_lengths()
@@ -468 +475 @@
                          if p.polydisperse and p.type not in ('orientation', 'magnetic'))
         self.pd_2d = set(p.name for p in self.call_parameters if p.polydisperse)
+
+    def set_zero_background(self):
+        """
+        Set the default background to zero for this model.  This is done for
+        structure factor models.
+        """
+        # type: () -> None
+        # Make sure background is the second common parameter.
+        assert self.common_parameters[1].id == "background"
+        self.common_parameters[1].default = 0.0
+        self.defaults = self._get_defaults()
 
     def check_angles(self):
@@ -567 +585 @@
     def _get_call_parameters(self):
         # type: () -> List[Parameter]
-        full_list = self.COMMON[:]
+        full_list = self.common_parameters[:]
         for p in self.kernel_parameters:
             if p.length == 1:
@@ -670 +688 @@
 
         # Gather the user parameters in order
-        result = control + self.COMMON
+        result = control + self.common_parameters
         for p in self.kernel_parameters:
             if not is2d and p.type in ('orientation', 'magnetic'):
@@ -770 +788 @@
 
     info = ModelInfo()
+
+    # Build the parameter table
     #print("make parameter table", kernel_module.parameters)
     parameters = make_parameter_table(getattr(kernel_module, 'parameters', []))
+
+    # background defaults to zero for structure factor models
+    structure_factor = getattr(kernel_module, 'structure_factor', False)
+    if structure_factor:
+        parameters.set_zero_background()
+
+    # TODO: remove demo parameters
+    # The plots in the docs are generated from the model default values.
+    # Sascomp set parameters from the command line, and so doesn't need
+    # demo values for testing.
     demo = expand_pars(parameters, getattr(kernel_module, 'demo', None))
+
     filename = abspath(kernel_module.__file__).replace('.pyc', '.py')
     kernel_id = splitext(basename(filename))[0]

sasmodels/models/hardsphere.py

@@ -162 +162 @@
     return pars
 
-demo = dict(radius_effective=200, volfraction=0.2,
-            radius_effective_pd=0.1, radius_effective_pd_n=40)
 # Q=0.001 is in the Taylor series, low Q part, so add Q=0.1,
 # assuming double precision sasview is correct

sasmodels/models/pearl_necklace.c

@@ -40 +40 @@
     const double si = sas_sinx_x(q*A_s);
     const double omsi = 1.0 - si;
-    const double pow_si = pow(si, num_pearls);
+    const double pow_si = pown(si, num_pearls);
 
     // form factor for num_pearls
@@ -81 +81 @@
 radius_from_volume(double radius, double edge_sep, double thick_string, double fp_num_pearls)
 {
-    const int num_pearls = (int) fp_num_pearls +0.5;
     const double vol_tot = form_volume(radius, edge_sep, thick_string, fp_num_pearls);
     return cbrt(vol_tot/M_4PI_3);

sasmodels/sasview_model.py

@@ -382 +382 @@
             hidden.add('scale')
             hidden.add('background')
-            self._model_info.parameters.defaults['background'] = 0.
 
         # Update the parameter lists to exclude any hidden parameters
@@ -695 +694 @@
             return self._calculate_Iq(qx, qy)
 
-    def _calculate_Iq(self, qx, qy=None, Fq=False, effective_radius_type=1):
+    def _calculate_Iq(self, qx, qy=None):
         if self._model is None:
             self._model = core.build_model(self._model_info)
@@ -715 +714 @@
         #print("values", values)
         #print("is_mag", is_magnetic)
-        if Fq:
-            result = calculator.Fq(call_details, values, cutoff=self.cutoff,
-                                   magnetic=is_magnetic,
-                                   effective_radius_type=effective_radius_type)
         result = calculator(call_details, values, cutoff=self.cutoff,
                             magnetic=is_magnetic)
@@ -736 +731 @@
         Calculate the effective radius for P(q)*S(q)
 
+        *mode* is the R_eff type, which defaults to 1 to match the ER
+        calculation for sasview models from version 3.x.
+
         :return: the value of the effective radius
         """
-        Fq = self._calculate_Iq([0.1], True, mode)
-        return Fq[2]
+        # ER and VR are only needed for old multiplication models, based on
+        # sas.sascalc.fit.MultiplicationModel.  Fail for now.  If we want to
+        # continue supporting them then add some test cases so that the code
+        # is exercised.  We can access ER/VR using the kernel Fq function by
+        # extending _calculate_Iq so that it calls:
+        #    if er_mode > 0:
+        #        res = calculator.Fq(call_details, values, cutoff=self.cutoff,
+        #                            magnetic=False, effective_radius_type=mode)
+        #        R_eff, form_shell_ratio = res[2], res[4]
+        #        return R_eff, form_shell_ratio
+        # Then use the following in calculate_ER:
+        #    ER, VR = self._calculate_Iq(q=[0.1], er_mode=mode)
+        #    return ER
+        # Similarly, for calculate_VR:
+        #    ER, VR = self._calculate_Iq(q=[0.1], er_mode=1)
+        #    return VR
+        # Obviously a combined calculate_ER_VR method would be better, but
+        # we only need them to support very old models, so ignore the 2x
+        # performance hit.
+        raise NotImplementedError("ER function is no longer available.")
 
     def calculate_VR(self):
@@ -748 +764 @@
         :return: the value of the form:shell volume ratio
         """
-        Fq = self._calculate_Iq([0.1], True, mode)
-        return Fq[4]
+        # See comments in calculate_ER.
+        raise NotImplementedError("VR function is no longer available.")
 
     def set_dispersion(self, parameter, dispersion):
@@ -914 +930 @@
     CylinderModel().evalDistribution([0.1, 0.1])
 
+def test_structure_factor_background():
+    # type: () -> None
+    """
+    Check that sasview model and direct model match, with background=0.
+    """
+    from .data import empty_data1D
+    from .core import load_model_info, build_model
+    from .direct_model import DirectModel
+
+    model_name = "hardsphere"
+    q = [0.0]
+
+    sasview_model = _make_standard_model(model_name)()
+    sasview_value = sasview_model.evalDistribution(np.array(q))[0]
+
+    data = empty_data1D(q)
+    model_info = load_model_info(model_name)
+    model = build_model(model_info)
+    direct_model = DirectModel(data, model)
+    direct_value_zero_background = direct_model(background=0.0)
+
+    assert sasview_value == direct_value_zero_background
+
+    # Additionally check that direct value background defaults to zero
+    direct_value_default = direct_model()
+    assert sasview_value == direct_value_default
+
+
 def magnetic_demo():
     Model = _make_standard_model('sphere')
@@ -934 +978 @@
     #print("rpa:", test_rpa())
     #test_empty_distribution()
+    #test_structure_factor_background()