Changeset 26a6608 in sasmodels for sasmodels

Timestamp:

Nov 20, 2017 11:49:03 AM (6 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:

a66b004

Parents:

a70959a (diff), fa70e04 (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' into boltzmann

Location:

sasmodels

Files:

• models/core_shell_parallelepiped.c (modified) (11 diffs)
• models/core_shell_parallelepiped.py (modified) (1 diff)
• product.py (modified) (2 diffs)
• sasview_model.py (modified) (3 diffs)
• weights.py (modified) (5 diffs)

sasmodels/models/core_shell_parallelepiped.c

@@ -1 @@
-double form_volume(double length_a, double length_b, double length_c, 
+double form_volume(double length_a, double length_b, double length_c,
                    double thick_rim_a, double thick_rim_b, double thick_rim_c);
 double Iq(double q, double core_sld, double arim_sld, double brim_sld, double crim_sld,
@@ -9 +9 @@
             double thick_rim_c, double theta, double phi, double psi);
 
-double form_volume(double length_a, double length_b, double length_c, 
+double form_volume(double length_a, double length_b, double length_c,
                    double thick_rim_a, double thick_rim_b, double thick_rim_c)
 {
     //return length_a * length_b * length_c;
-    return length_a * length_b * length_c + 
-           2.0 * thick_rim_a * length_b * length_c + 
+    return length_a * length_b * length_c +
+           2.0 * thick_rim_a * length_b * length_c +
            2.0 * thick_rim_b * length_a * length_c +
            2.0 * thick_rim_c * length_a * length_b;
@@ -34 +34 @@
     // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled
     // Did not understand the code completely, it should be rechecked (Miguel Gonzalez)
-    
+    //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker)
+
     const double mu = 0.5 * q * length_b;
-    
+
     //calculate volume before rescaling (in original code, but not used)
-    //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);        
-    //double vol = length_a * length_b * length_c + 
-    //       2.0 * thick_rim_a * length_b * length_c + 
+    //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);
+    //double vol = length_a * length_b * length_c +
+    //       2.0 * thick_rim_a * length_b * length_c +
     //       2.0 * thick_rim_b * length_a * length_c +
     //       2.0 * thick_rim_c * length_a * length_b;
-    
+
     // Scale sides by B
     const double a_scaled = length_a / length_b;
     const double c_scaled = length_c / length_b;
 
-    // ta and tb correspond to the definitions in CSPPKernel, but they don't make sense to me (MG)
-    // the a_scaled in the definition of tb was present in CSPPKernel in libCylinder.c,
-    // while in cspkernel in csparallelepiped.cpp (used for the 2D), all the definitions
-    // for ta, tb, tc use also A + 2*rim_thickness (but not scaled by B!!!)
-    double ta = (a_scaled + 2.0*thick_rim_a)/length_b;
-    double tb = (a_scaled + 2.0*thick_rim_b)/length_b;
+    double ta = a_scaled + 2.0*thick_rim_a/length_b; // incorrect ta = (a_scaled + 2.0*thick_rim_a)/length_b;
+    double tb = 1+ 2.0*thick_rim_b/length_b; // incorrect tb = (a_scaled + 2.0*thick_rim_b)/length_b;
+    double tc = c_scaled + 2.0*thick_rim_c/length_b; //not present
 
     double Vin = length_a * length_b * length_c;
@@ -62 +60 @@
     double V1 = (2.0 * thick_rim_a * length_b * length_c);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc)
     double V2 = (2.0 * length_a * thick_rim_b * length_c);    // incorrect V2(aa*bb*cc+2*aa*tb*cc)
+    double V3 = (2.0 * length_a * length_b * thick_rim_c);    //not present
+    double Vot = Vin + V1 + V2 + V3;
 
     // Scale factors (note that drC is not used later)
@@ -67 +67 @@
     const double drhoA = (arim_sld-solvent_sld);
     const double drhoB = (brim_sld-solvent_sld);
-    //const double drC_Vot = (crim_sld-solvent_sld)*Vot;
+    const double drhoC = (crim_sld-solvent_sld);  // incorrect const double drC_Vot = (crim_sld-solvent_sld)*Vot;
+
 
     // Precompute scale factors for combining cross terms from the shape
     const double scale23 = drhoA*V1;
     const double scale14 = drhoB*V2;
-    const double scale12 = drho0*Vin - scale23 - scale14;
+    const double scale24 = drhoC*V3;
+    const double scale11 = drho0*Vin;
+    const double scale12 = drho0*Vin - scale23 - scale14 - scale24;
 
     // outer integral (with gauss points), integration limits = 0, 1
@@ -83 +86 @@
         // inner integral (with gauss points), integration limits = 0, 1
         double inner_total = 0.0;
+        double inner_total_crim = 0.0;
         for(int j=0; j<76; j++) {
             const double uu = 0.5 * ( Gauss76Z[j] + 1.0 );
@@ -88 +92 @@
             SINCOS(M_PI_2*uu, sin_uu, cos_uu);
             const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled);
-            const double si2 = sas_sinx_x(mu_proj * cos_uu);
+            const double si2 = sas_sinx_x(mu_proj * cos_uu );
             const double si3 = sas_sinx_x(mu_proj * sin_uu * ta);
             const double si4 = sas_sinx_x(mu_proj * cos_uu * tb);
@@ -94 +98 @@
             // Expression in libCylinder.c (neither drC nor Vot are used)
             const double form = scale12*si1*si2 + scale23*si2*si3 + scale14*si1*si4;
+            const double form_crim = scale11*si1*si2;
 
-            // To note also that in csparallelepiped.cpp, there is a function called
-            // cspkernel, which seems to make more sense and has the following comment:
-            //   This expression is different from NIST/IGOR package (I strongly believe the IGOR is wrong!!!). 10/15/2010.
-            //   tmp =( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3)*
-            //   ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3);   //  correct FF : square of sum of phase factors
-            // This is the function called by csparallelepiped_analytical_2D_scaled,
-            // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c        
-            
+
             //  correct FF : sum of square of phase factors
             inner_total += Gauss76Wt[j] * form * form;
+            inner_total_crim += Gauss76Wt[j] * form_crim * form_crim;
         }
         inner_total *= 0.5;
-
+        inner_total_crim *= 0.5;
         // now sum up the outer integral
         const double si = sas_sinx_x(mu * c_scaled * sigma);
-        outer_total += Gauss76Wt[i] * inner_total * si * si;
+        const double si_crim = sas_sinx_x(mu * tc * sigma);
+        outer_total += Gauss76Wt[i] * (inner_total * si * si + inner_total_crim * si_crim * si_crim);
     }
     outer_total *= 0.5;
@@ -154 +154 @@
 
     // The definitions of ta, tb, tc are not the same as in the 1D case because there is no
-    // the scaling by B. The use of length_a for the 3 of them seems clearly a mistake to me,
-    // but for the moment I let it like this until understanding better the code.
+    // the scaling by B.
     double ta = length_a + 2.0*thick_rim_a;
-    double tb = length_a + 2.0*thick_rim_b;
-    double tc = length_a + 2.0*thick_rim_c;
+    double tb = length_b + 2.0*thick_rim_b;
+    double tc = length_c + 2.0*thick_rim_c;
     //handle arg=0 separately, as sin(t)/t -> 1 as t->0
     double siA = sas_sinx_x(0.5*q*length_a*xhat);
@@ -166 +165 @@
     double siBt = sas_sinx_x(0.5*q*tb*yhat);
     double siCt = sas_sinx_x(0.5*q*tc*zhat);
-    
+
 
     // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed
@@ -173 +172 @@
                + drA*(siAt-siA)*siB*siC*V1
                + drB*siA*(siBt-siB)*siC*V2
-               + drC*siA*siB*(siCt*siCt-siC)*V3);
-   
+               + drC*siA*siB*(siCt-siC)*V3);
+
     return 1.0e-4 * f * f;
 }

sasmodels/models/core_shell_parallelepiped.py

@@ -211 +211 @@
 
 # rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models
-qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.)
-tests = [[{}, 0.2, 0.533149288477],
-         [{}, [0.2], [0.533149288477]],
-         [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0853299803222],
-         [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0853299803222]],
-        ]
-del qx, qy  # not necessary to delete, but cleaner
+if 0:  # pak: model rewrite; need to update tests
+    qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.)
+    tests = [[{}, 0.2, 0.533149288477],
+            [{}, [0.2], [0.533149288477]],
+            [{'theta':10.0, 'phi':20.0}, (qx, qy), 0.0853299803222],
+            [{'theta':10.0, 'phi':20.0}, [(qx, qy)], [0.0853299803222]],
+            ]
+    del qx, qy  # not necessary to delete, but cleaner

sasmodels/product.py

@@ -100 +100 @@
     # Remember the component info blocks so we can build the model
     model_info.composition = ('product', [p_info, s_info])
+    model_info.control = p_info.control
+    model_info.hidden = p_info.hidden
+    if getattr(p_info, 'profile', None) is not None:
+        profile_pars = set(p.id for p in p_info.parameters.kernel_parameters)
+        def profile(**kwargs):
+            # extract the profile args
+            kwargs = dict((k, v) for k, v in kwargs.items() if k in profile_pars)
+            return p_info.profile(**kwargs)
+    else:
+        profile = None
+    model_info.profile = profile
+    model_info.profile_axes = p_info.profile_axes
+
     # TODO: delegate random to p_info, s_info
     #model_info.random = lambda: {}
@@ -129 +142 @@
     def __init__(self, model_info, P, S):
         # type: (ModelInfo, KernelModel, KernelModel) -> None
+        #: Combined info plock for the product model
         self.info = model_info
+        #: Form factor modelling individual particles.
         self.P = P
+        #: Structure factor modelling interaction between particles.
         self.S = S
-        self.dtype = P.dtype
+        #: Model precision. This is not really relevant, since it is the
+        #: individual P and S models that control the effective dtype,
+        #: converting the q-vectors to the correct type when the kernels
+        #: for each are created. Ideally this should be set to the more
+        #: precise type to avoid loss of precision, but precision in q is
+        #: not critical (single is good enough for our purposes), so it just
+        #: uses the precision of the form factor.
+        self.dtype = P.dtype  # type: np.dtype
 
     def make_kernel(self, q_vectors):

sasmodels/sasview_model.py

@@ -205 +205 @@
                                            structure_factor._model_info)
     ConstructedModel = make_model_from_info(model_info)
-    return ConstructedModel()
+    return ConstructedModel(form_factor.multiplicity)
 
 
@@ -323 +323 @@
     #: True if model has multiplicity
     is_multiplicity_model = False
-    #: Mulitplicity information
+    #: Multiplicity information
     multiplicity_info = None # type: MultiplicityInfoType
 
@@ -354 +354 @@
         # and lines to plot.
 
-        # Get the list of hidden parameters given the mulitplicity
+        # Get the list of hidden parameters given the multiplicity
         # Don't include multiplicity in the list of parameters
         self.multiplicity = multiplicity

sasmodels/weights.py

@@ -55 +55 @@
         """
         sigma = self.width * center if relative else self.width
+        if not relative:
+            # For orientation, the jitter is relative to 0 not the angle
+            center = 0
+            pass
         if sigma == 0 or self.npts < 2:
             if lb <= center <= ub:
@@ -93 +97 @@
         return x, px
 
+class UniformDispersion(Dispersion):
+    r"""
+    Uniform dispersion, with width $\sigma$.
+
+    .. math::
+
+        w = 1
+    """
+    type = "uniform"
+    default = dict(npts=35, width=0, nsigmas=1)
+    def _weights(self, center, sigma, lb, ub):
+        x = self._linspace(center, sigma, lb, ub)
+        x = x[np.fabs(x-center) <= np.fabs(sigma)]
+        return x, np.ones_like(x)
 
 class RectangleDispersion(Dispersion):
@@ -186 +204 @@
         return x, px
 
+class BoltzmannDispersion(Dispersion):
+    r"""
+    Boltzmann dispersion, with $\sigma=k T/E$.
+
+    .. math::
+
+        w = \exp\left( -|x - c|/\sigma\right)
+    """
+    type = "boltzmann"
+    default = dict(npts=35, width=0, nsigmas=3)
+    def _weights(self, center, sigma, lb, ub):
+        x = self._linspace(center, sigma, lb, ub)
+        px = np.exp(-np.abs(x-center) / np.abs(sigma))
+        return x, px
 
 # dispersion name -> disperser lookup table.
@@ -192 +224 @@
 MODELS = OrderedDict((d.type, d) for d in (
     RectangleDispersion,
+    UniformDispersion,
     ArrayDispersion,
     LogNormalDispersion,
     GaussianDispersion,
     SchulzDispersion,
+    BoltzmannDispersion
 ))
 
@@ -225 +259 @@
     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
+    return v, w/np.sum(w)
+
+
+def plot_weights(model_info, mesh):
+    # type: (ModelInfo, List[Tuple[float, np.ndarray, np.ndarray]]) -> None
     """
     Plot the weights returned by :func:`get_weights`.
 
-    *model_info* is
-    :param model_info:
-    :param pairs:
-    :return:
+    *model_info* defines model parameters, etc.
+
+    *mesh* is a list of tuples containing (*value*, *dispersity*, *weights*)
+    for each parameter, where (*dispersity*, *weights*) pairs are the
+    distributions to be plotted.
     """
     import pylab
 
-    if any(len(values)>1 for values, weights in pairs):
+    if any(len(dispersity)>1 for value, dispersity, weights in mesh):
         labels = [p.name for p in model_info.parameters.call_parameters]
-        pylab.interactive(True)
+        #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)
+        for (v,x,w), s in zip(mesh, labels):
+            if len(x) > 1:
+                pylab.plot(x, w, '-o', label=s)
         pylab.grid(True)
         pylab.legend()