Changeset 9acade6 in sasmodels

Timestamp:

Sep 27, 2016 9:52:14 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:

bb46723

Parents:

6cefbc9 (diff), 5fd684d (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.

Message:

Merge branch 'master' of github.com:sasview/sasmodels

Location:

sasmodels

Files:

• models/pringle.c (modified) (3 diffs)
• models/pringle.py (modified) (6 diffs)
• sasview_model.py (modified) (3 diffs)
• weights.py (modified) (7 diffs)

sasmodels/models/pringle.c

@@ -42 +42 @@
     }
 
-    #if 1
-    // TODO: should the normalization be to R^2 or the last value, r^2
-    // the gaussian window does not go all the way from 0 to 1.
-    //radius = Gauss76Z[75] * zm + zb;
-    *Sn = zm*sumS / (r*r);
-    *Cn = zm*sumC / (r*r);
-    #else
     *Sn = zm*sumS / (radius*radius);
     *Cn = zm*sumC / (radius*radius);
-    #endif
 }
 
@@ -109 +101 @@
 double form_volume(double radius, double thickness, double alpha, double beta)
 {
-    // TODO: Normalize by form volume
-    //return M_PI*radius*radius*thickness;
-    return 1.0;
+    return M_PI*radius*radius*thickness;
 }
 
@@ -126 +116 @@
     double contrast = sld_pringle - sld_solvent;
     double volume = M_PI*radius*radius*thickness;
-    // TODO: If normalize by form volume, need an extra volume here
-    //return 1.0e-4*form * square(contrast * volume);
-    return 1.0e-4*form * square(contrast) * volume;
+    return 1.0e-4*form * square(contrast * volume);
 }

sasmodels/models/pringle.py

@@ -17 +17 @@
 .. math::
 
-    C_n = \int^{R}_{0} r dr\cos(qr^2\alpha \cos{\psi})
+    C_n = \frac{1}{r^2}\int^{R}_{0} r dr\cos(qr^2\alpha \cos{\psi})
     J_n\left( qr^2\beta \cos{\psi}\right)
     J_{2n}\left( qr \sin{\psi}\right)
@@ -23 +23 @@
 .. math::
 
-    S_n = \int^{R}_{0} r dr\sin(qr^2\alpha \cos{\psi})
+    S_n = \frac{1}{r^2}\int^{R}_{0} r dr\sin(qr^2\alpha \cos{\psi})
     J_n\left( qr^2\beta \cos{\psi}\right)
     J_{2n}\left( qr \sin{\psi}\right)
@@ -47 +47 @@
 **Last Modified by:** Wojciech Wpotrzebowski **on:** March 20, 2016
 
-**Last Reviewed by:** Paul Butler **on:** March 21, 2016
+**Last Reviewed by:** Andrew Jackson **on:** September 26, 2016
 
 """
@@ -102 +102 @@
       'sld_pringle': 1.0,
       'sld_solvent': 6.3,
-      'background': 6.3,
-     }, 0.1, 16.185532],
+      'background': 0.001,
+     }, 0.1, 9.87676],
 
     [{'scale' : 1.0,
@@ -112 +112 @@
       'sld_pringle': 1.0,
       'sld_solvent': 6.3,
-      'background': 6.3,
-     }, 0.01, 297.153496],
+      'background': 0.001,
+     }, 0.01, 290.56723],
 
     [{'scale' : 1.0,
@@ -122 +122 @@
       'sld_pringle': 1.0,
       'sld_solvent': 6.3,
-      'background': 6.3,
-     }, 0.001, 324.021256415],
-
-    [{'scale' : 1.0,
-      'radius': 60.0,
-      'thickness': 10.0,
-      'alpha': 0.001,
-      'beta': 0.02,
-      'sld_pringle': 1.0,
-      'sld_solvent': 6.3,
-      'background': 6.3,
-     }, (0.001, 90.0), 6.30000026876],
+      'background': 0.001,
+     }, 0.001, 317.40847],
 ]

sasmodels/sasview_model.py

@@ -565 +565 @@
         parameters = self._model_info.parameters
         pairs = [self._get_weights(p) for p in parameters.call_parameters]
+        #weights.plot_weights(self._model_info, pairs)
         call_details, values, is_magnetic = make_kernel_args(calculator, pairs)
         #call_details.show()
@@ -618 +619 @@
             # remember them is kind of funky.
             # Note: can't seem to get disperser parameters from sasview
-            # (1) Could create a sasview model that has not yet # been
+            # (1) Could create a sasview model that has not yet been
             # converted, assign the disperser to one of its polydisperse
             # parameters, then retrieve the disperser parameters from the
-            # sasview model.  (2) Could write a disperser parameter retriever
-            # in sasview.  (3) Could modify sasview to use sasmodels.weights
-            # dispersers.
+            # sasview model.
+            # (2) Could write a disperser parameter retriever in sasview.
+            # (3) Could modify sasview to use sasmodels.weights dispersers.
             # For now, rely on the fact that the sasview only ever uses
             # new dispersers in the set_dispersion call and create a new
             # one instead of trying to assign parameters.
-            dispersion = weights.MODELS[dispersion.type]()
             self.dispersion[parameter] = dispersion.get_pars()
         else:
@@ -658 +658 @@
         elif par.polydisperse:
             dis = self.dispersion[par.name]
-            value, weight = weights.get_weights(
-                dis['type'], dis['npts'], dis['width'], dis['nsigmas'],
-                self.params[par.name], par.limits, par.relative_pd)
+            if dis['type'] == 'array':
+                value, weight = dis['values'], dis['weights']
+            else:
+                value, weight = weights.get_weights(
+                    dis['type'], dis['npts'], dis['width'], dis['nsigmas'],
+                    self.params[par.name], par.limits, par.relative_pd)
             return value, weight / np.sum(weight)
         else:

sasmodels/weights.py

@@ -3 +3 @@
 """
 # TODO: include dispersion docs with the disperser models
+from __future__ import division, print_function
+
 from math import sqrt  # type: ignore
+from collections import OrderedDict
+
 import numpy as np  # type: ignore
 from scipy.special import gammaln  # type: ignore
@@ -54 +58 @@
             else:
                 return np.array([], 'd'), np.array([], 'd')
-        return self._weights(center, sigma, lb, ub)
+        x, px = self._weights(center, sigma, lb, ub)
+        return x, px
 
     def _weights(self, center, sigma, lb, ub):
@@ -79 +84 @@
     default = dict(npts=35, width=0, nsigmas=3)
     def _weights(self, center, sigma, lb, ub):
+        # TODO: sample high probability regions more densely
+        # i.e., step uniformly in cumulative density rather than x value
+        # so weight = 1/Npts for all weights, but values are unevenly spaced
         x = self._linspace(center, sigma, lb, ub)
         px = np.exp((x-center)**2 / (-2.0 * sigma * sigma))
@@ -165 +173 @@
 
     def _weights(self, center, sigma, lb, ub):
-        # TODO: interpolate into the array dispersion using npts
-        x = center + self.values*sigma
+        # TODO: rebin the array dispersion using npts
+        # TODO: use a distribution that can be recentered and scaled
+        x = self.values
+        #x = center + self.values*sigma
         idx = (x >= lb) & (x <= ub)
         x = x[idx]
@@ -174 +184 @@
 
 # dispersion name -> disperser lookup table.
-MODELS = dict((d.type, d) for d in (
-    GaussianDispersion, RectangleDispersion,
-    ArrayDispersion, SchulzDispersion, LogNormalDispersion
+# Maintain order since this is used by sasview GUI to order the options in
+# the dispersion type combobox.
+MODELS = OrderedDict((d.type, d) for d in (
+    RectangleDispersion,
+    ArrayDispersion,
+    LogNormalDispersion,
+    GaussianDispersion,
+    SchulzDispersion,
 ))
 
@@ -194 +209 @@
     *value* is the value of the parameter in the model.
 
-    *limits* is [lb, ub], the lower and upper bound of the weight value.
+    *limits* is [lb, ub], the lower and upper bound on the possible values.
 
     *relative* is true if *width* is defined in proportion to the value
@@ -201 +216 @@
     Returns *(value, weight)*, where *value* and *weight* are vectors.
     """
+    if disperser == "array":
+        raise NotImplementedError("Don't handle arrays through get_weights; use values and weights directly")
     cls = MODELS[disperser]
     obj = cls(n, width, nsigmas)
     v, w = obj.get_weights(value, limits[0], limits[1], relative)
     return v, w
+
+
+def plot_weights(model_info, pairs):
+    # type: (ModelInfo, List[Tuple[np.ndarray, np.ndarray]]) -> None
+    """
+    Plot the weights returned by :func:`get_weights`.
+
+    *model_info* is
+    :param model_info:
+    :param pairs:
+    :return:
+    """
+    import pylab
+
+    if any(len(values)>1 for values, weights in pairs):
+        labels = [p.name for p in model_info.parameters.call_parameters]
+        pylab.interactive(True)
+        pylab.figure()
+        for (v,w), s in zip(pairs, labels):
+            if len(v) > 1:
+                #print("weights for", s, v, w)
+                pylab.plot(v, w, '-o', label=s)
+        pylab.grid(True)
+        pylab.legend()