Changeset 54bcd4a in sasmodels

Timestamp:

Aug 4, 2016 4:48:37 PM (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:

745b7bb

Parents:

bd49c79

Message:

spherical sld: simplify code so that it works on AMD GPUs

Location:

sasmodels

Files:

• conversion_table.py (modified) (1 diff)
• convert.py (modified) (3 diffs)
• models/spherical_sld.c (modified) (1 diff)
• models/spherical_sld.py (modified) (4 diffs)

sasmodels/conversion_table.py

@@ -688 +688 @@
         }
     ],
-    "magsphere": [
-        "SphereModel",
-        {
-            "sld": "sldSph",
-            "radius": "radius",
-            "sld_solvent": "sldSolv",
-            "sld_M0": "M0_sld_sph",
-            "sphere_theta": "M_theta_sph",
-            "sphere_phi": "M_phi_sph",
-            "sld_solvent_M0": "M0_sld_solv",
-            "solvent_theta": "M_theta_solv",
-            "solvent_phi": "M_phi_solv",
-            "in_spin": "Up_frac_i",
-            "out_spin": "Up_frac_f",
-            "spin_theta": "Up_theta"
-        }
-    ],
-    "_spherepy": [
-        "SphereModel",
-        {
-            "sld": "sldSph",
-            "radius": "radius",
-            "sld_solvent": "sldSolv"
-        }
-    ],
     "spherical_sld": [
         "SphericalSLDModel",
-        {
-            "radius_core": "rad_core0",
-            "sld_core": "sld_core0",
-            "sld_solvent": "sld_solv"
-        }
+        # Be lazy and use a generator expression to define
+        #    sld1: sld_flat0, ...
+        #    thickness1: thick_flat0, ...
+        #    interface1: thick_inter0, ...
+        #    shape1: func_inter0, ...
+        #    nu1: nu_inter0, ...
+        # but override thickness1 => rad_cor0 and sld1 => sld_core0.
+        # Note: explicit key,value pairs given by **{...} override the
+        # keys from the gnerator expression ((k,v) for k,v in seq) when
+        # used as dict((generator), **{...})
+        dict(((field_new+str(index+1),field_old+str(index))
+             for field_new, field_old in [("sld","sld_flat"),
+                                          ("thickness","thick_flat"),
+                                          ("interface","thick_inter"),
+                                          ("shape","func_inter"),
+                                          ("nu","nu_inter"),]
+             for index in range(11)),
+             **{
+            "n_shells": "n_shells",
+            "n_steps": "npts_inter",
+            "sld_solvent": "sld_solv",
+            "thickness1": "rad_core0",
+            "sld1": "sld_core0",
+        })
     ],
     "squarewell": [

sasmodels/convert.py

@@ -229 +229 @@
                 oldpars['ndensity'] *= 1e15
         elif name == 'spherical_sld':
-            for k in range(1, int(pars['n_shells'])+1):
-                _remove_pd(oldpars, 'thick_flat'+str(k), 'thick_flat')
-                _remove_pd(oldpars, 'thick_inter'+str(k), 'thick_inter')
+            oldpars["CONTROL"] -= 1
+            # remove polydispersity from shells
+            for k in range(1, 11):
+                _remove_pd(oldpars, 'thick_flat'+str(k), 'thickness')
+                _remove_pd(oldpars, 'thick_inter'+str(k), 'interface')
+            # remove extra shells
+            for k in range(int(pars['n_shells']), 11):
+                oldpars.pop('sld_flat'+str(k), 0)
+                oldpars.pop('thick_flat'+str(k), 0)
+                oldpars.pop('thick_inter'+str(k), 0)
+                oldpars.pop('func_inter'+str(k), 0)
+                oldpars.pop('nu_inter'+str(k), 0)
         elif name == 'core_multi_shell':
             # kill extra shells
@@ -265 +274 @@
                     oldpars.pop(k, None)
 
+    #print("convert from",list(sorted(pars)))
+    #print("convert to",list(sorted(oldpars.items())))
     return oldpars
 
@@ -323 +334 @@
         elif name == 'spherical_sld':
             pars['n_shells'] = math.ceil(pars['n_shells'])
-            pars['npts_inter'] = math.ceil(pars['npts_inter'])
-            pars['func_inter0'] = math.trunc(pars['func_inter0']+0.5)
-            for k in range(1, 11):
-                pars['func_inter%d'%k] = math.trunc(pars['func_inter%d'%k]+0.5)
-                pars['thick_flat%d_pd_n'%k] = 0
-                pars['thick_inter%d_pd_n'%k] = 0
-
+            pars['n_steps'] = math.ceil(pars['n_steps'])
+            for k in range(1, 12):
+                pars['shape%d'%k] = math.trunc(pars['shape%d'%k]+0.5)
+            for k in range(2, 12):
+                pars['thickness%d_pd_n'%k] = 0
+                pars['interface%d_pd_n'%k] = 0
+

sasmodels/models/spherical_sld.c

@@ -1 +1 @@
 static double form_volume(
     int n_shells,
-    double radius_core,
-    double thick_inter0,
-    double thick_flat[],
-    double thick_inter[])
+    double thickness[],
+    double interface[])
 {
-    int i;
-    double r = radius_core;
-    r += thick_inter0;
-    for (i=0; i < n_shells; i++) {
-        r += thick_inter[i];
-        r += thick_flat[i];
+    double r = 0.0;
+    for (int i=0; i < n_shells; i++) {
+        r += thickness[i] + interface[i];
     }
     return M_4PI_3*cube(r);
 }
 
+static double blend(int shape, double nu, double z)
+{
+    if (shape==0) {
+        const double num = sas_erf(nu * M_SQRT1_2 * (2.0*z - 1.0));
+        const double denom = 2.0 * sas_erf(nu * M_SQRT1_2);
+        return num/denom + 0.5;
+    } else if (shape==1) {
+        return pow(z, nu);
+    } else if (shape==2) {
+        return 1.0 - pow(1. - z, nu);
+    } else if (shape==3) {
+        return expm1(-nu*z)/expm1(-nu);
+    } else if (shape==4) {
+        return expm1(nu*z)/expm1(nu);
+    } else {
+        return NAN;
+    }
+}
 
-static double sphere_sld_kernel(
+static double f_linear(double q, double r, double contrast, double slope)
+{
+    const double qr = q * r;
+    const double qrsq = qr * qr;
+    const double bes = sph_j1c(qr);
+    double sinqr, cosqr;
+    SINCOS(qr, sinqr, cosqr);
+    const double fun = 3.0*r*(2.0*qr*sinqr - (qrsq-2.0)*cosqr)/(qrsq*qrsq);
+    const double vol = M_4PI_3 * cube(r);
+    return vol*(bes*contrast + fun*slope);
+}
+
+static double Iq(
     double q,
     int n_shells,
-    int npts_inter,
-    double radius_core,
-    double sld_core,
     double sld_solvent,
-    double func_inter_core,
-    double thick_inter_core,
-    double nu_inter_core,
-    double sld_flat[],
-    double thick_flat[],
-    double func_inter[],
-    double thick_inter[],
-    double nu_inter[] ) {
+    double sld[],
+    double thickness[],
+    double interface[],
+    double shape[],
+    double nu[],
+    int n_steps)
+{
+    // iteration for # of shells + core + solvent
+    double f=0.0;
+    double r=0.0;
+    for (int shell=0; shell<n_shells; shell++){
+        const double sld_l = sld[shell];
 
-    int i,j,k;
-    int n_s;
+        // uniform shell; r=0 => r^3=0 => f=0, so works for core as well.
+        f -= M_4PI_3 * cube(r) * sld_l * sph_j1c(q*r);
+        r += thickness[shell];
+        f += M_4PI_3 * cube(r) * sld_l * sph_j1c(q*r);
 
-    double sld_i,sld_f,dz,bes,fun,f,vol,qr,r,contr,f2;
-    double sign,slope=0.0;
-    double pi;
+        // iterate over sub_shells in the interface
+        const double dr = interface[shell]/n_steps;
+        const double delta = (shell==n_shells-1 ? sld_solvent : sld[shell+1]) - sld_l;
+        const double nu_shell = fmax(fabs(nu[shell]), 1.e-14);
+        const int shape_shell = (int)(shape[shell]);
 
-    double total_thick=0.0;
+        // if there is no interface the equations don't work
+        if (dr == 0.) continue;
 
-    //TODO: This part can be further simplified
-    int fun_type[12];
-    double sld[12];
-    double thick_internal[12];
-    double thick[12];
-    double fun_coef[12];
+        double sld_in = sld_l;
+        for (int step=1; step <= n_steps; step++) {
+            // find sld_i at the outer boundary of sub-shell step
+            //const double z = (double)step/(double)n_steps;
+            const double z = (double)step/(double)n_steps;
+            const double fraction = blend(shape_shell, nu_shell, z);
+            const double sld_out = fraction*delta + sld_l;
+            // calculate slope
+            const double slope = (sld_out - sld_in)/dr;
+            const double contrast = sld_in - slope*r;
 
-    fun_type[0] = func_inter_core;
-    fun_coef[0] = fabs(nu_inter_core);
-    sld[0] = sld_core;
-    thick[0] = radius_core;
-    thick_internal[0] = thick_inter_core;
-
-    for (i =1; i<=n_shells; i++){
-        sld[i] = sld_flat[i-1];
-        thick_internal[i]= thick_inter[i-1];
-        thick[i] = thick_flat[i-1];
-        fun_type[i] = func_inter[i-1];
-        fun_coef[i] = fabs(nu_inter[i-1]);
-        total_thick += thick[i];
-        total_thick += thick_internal[i]; //doesn't account for core layer
-    }
-
-    sld[n_shells+1] = sld_solvent;
-    thick[n_shells+1] = total_thick/5.0;
-    thick_internal[n_shells+1] = 0.0;
-    fun_coef[n_shells+1] = 0.0;
-    fun_type[n_shells+1] = 0;
-
-    pi = 4.0*atan(1.0);
-    f = 0.0;
-    r = 0.0;
-    vol = 0.0;
-    sld_f = sld_core;
-
-    dz = 0.0;
-    // iteration for # of shells + core + solvent
-    for (i=0;i<=n_shells+1; i++){
-        // iteration for flat and interface
-        for (j=0;j<2;j++){
-            // iteration for sub_shells in the interface
-            // starts from #1 sub-layer
-            for (n_s=1;n_s<=npts_inter; n_s++){
-                // for solvent, it doesn't have an interface
-                if (i==n_shells+1 && j==1)
-                    break;
-                // for flat layers
-                if (j==0){
-                    dz = thick[i];
-                    sld_i = sld[i];
-                    slope = 0.0;
-                }
-                // for interfacial sub_shells
-                else{
-                    dz = thick_internal[i]/npts_inter;
-                    // find sld_i at the outer boundary of sub-layer #n_s
-                    sld_i = intersldfunc(fun_type[i], npts_inter, n_s,
-                            fun_coef[i], sld[i], sld[i+1]);
-                    // calculate slope
-                    slope= (sld_i -sld_f)/dz;
-                }
-                contr = sld_f-slope*r;
-                // iteration for the left and right boundary of the shells
-                for (k=0; k<2; k++){
-                    // At r=0, the contribution to I is zero, so skip it.
-                    if ( i == 0 && j == 0 && k == 0){
-                        continue;
-                    }
-                    // On the top of slovent there is no interface; skip it.
-                    if (i == n_shells+1 && k == 1){
-                        continue;
-                    }
-                    // At the right side (outer) boundary
-                    if ( k == 1){
-                        sign = 1.0;
-                        r += dz;
-                    }
-                    // At the left side (inner) boundary
-                    else{
-                        sign = -1.0;
-                    }
-
-                    qr = q * r;
-                    fun = 0.0;
-
-                    if(qr == 0.0){
-                        bes = sign * 1.0;
-                    }
-                    else{
-                        // for flat sub-layer
-                        //TODO: Single precision calculation fails here
-                        bes = sign *  sph_j1c(qr);
-                        if (fabs(slope) > 0.0 ){
-                            const double qrsq = qr*qr;
-                            double sinqr, cosqr;
-                            SINCOS(qr, sinqr, cosqr);
-                            fun = sign * 3.0 * r *
-                            (2.0*qr*sinqr - (qrsq-2.0)*cosqr)/(qrsq * qrsq);
-                            // In the onion model Jae-He's formula is rederived
-                            // and gives following:
-                            //fun = 3.0 * sign * r *
-                            //(2.0*cosqr + qr*sinqr)/(qrsq*qrsq);
-                            //But this seems not to be working in this case...
-                        }
-                    }
-
-                    // update total volume
-                    vol = M_4PI_3 * cube(r);
-                    f += vol * (bes * contr + fun * slope);
-                }
-                sld_f = sld_i;
-                // no sub-layer iteration (n_s loop) for the flat layer
-                if (j==0)
-                    break;
-            }
+            // iteration for the left and right boundary of the shells
+            f -= f_linear(q, r, contrast, slope);
+            r += dr;
+            f += f_linear(q, r, contrast, slope);
+            sld_in = sld_out;
         }
     }
+    // add in solvent effect
+    f -= M_4PI_3 * cube(r) * sld_solvent * sph_j1c(q*r);
 
-    f2 = f * f * 1.0e-4;
-    return (f2);
+    const double f2 = f * f * 1.0e-4;
+    return f2;
 }
 
-
-/**
- * Function to evaluate 1D SphereSLD function
- * @param q: q-value
- * @return: function value
- */
-static double Iq(double q,
-    int n_shells,
-    int npts_inter,
-    double radius_core,
-    double sld_core,
-    double sld_solvent,
-    double func_inter0,
-    double thick_inter0,
-    double nu_inter0,
-    double sld_flat[],
-    double thick_flat[],
-    double func_inter[],
-    double thick_inter[],
-    double nu_inter[] ) {
-
-    double intensity = sphere_sld_kernel(q, n_shells, npts_inter, radius_core,
-                sld_core, sld_solvent, func_inter0, thick_inter0, nu_inter0,
-                sld_flat, thick_flat, func_inter, thick_inter, nu_inter);
-
-    return intensity;
-}
-
-/**
- * Function to evaluate 2D SphereSLD function
- * @param q_x: value of Q along x
- * @param q_y: value of Q along y
- * @return: function value
- */
-
-/*static double Iqxy(double qx, double qy,
-    int n_shells,
-    int npts_inter,
-    double radius_core
-    double sld_core,
-    double sld_solvent,
-    double sld_flat[],
-    double thick_flat[],
-    double func_inter[],
-    double thick_inter[],
-    double nu_inter[],
-    ) {
-
-    double q = sqrt(qx*qx + qy*qy);
-    return Iq(q, n_shells, npts_inter, radius_core, sld_core, sld_solvent,
-    sld_flat[], thick_flat[], func_inter[], thick_inter[], nu_inter[])
-}*/
-

sasmodels/models/spherical_sld.py

@@ -190 +190 @@
 
 import numpy as np
-from numpy import inf
+from numpy import inf, expm1, sqrt
+from scipy.special import erf
 
 name = "spherical_sld"
@@ -200 +201 @@
 category = "shape:sphere"
 
+SHAPES = ["erf(|nu|*z)", "Rpow(z^|nu|)", "Lpow(z^|nu|)",
+          "Rexp(-|nu|z)", "Lexp(-|nu|z)"],
+
 # pylint: disable=bad-whitespace, line-too-long
 #            ["name", "units", default, [lower, upper], "type", "description"],
-parameters = [["n_shells",          "",               1,      [0, 10],        "volume", "number of shells"],
-              ["npts_inter",        "",               35,     [0, inf],       "", "number of points in each sublayer Must be odd number"],
-              ["radius_core",       "Ang",            50.0,   [0, inf],       "volume", "intern layer thickness"],
-              ["sld_core",          "1e-6/Ang^2",     2.07,   [-inf, inf],    "sld", "sld function flat"],
-              ["sld_solvent",       "1e-6/Ang^2",     1.0,    [-inf, inf],    "sld", "sld function solvent"],
-              ["func_inter0",       "",               0,      [0, 5],         "", "Erf:0, RPower:1, LPower:2, RExp:3, LExp:4, Linear:5"],
-              ["thick_inter0",      "Ang",            50.0,   [0, inf],       "volume", "intern layer thickness for core layer"],
-              ["nu_inter0",         "",               2.5,    [-inf, inf],    "", "steepness parameter for core layer"],
-              ["sld_flat[n_shells]",      "1e-6/Ang^2",     4.06,   [-inf, inf],    "sld", "sld function flat"],
-              ["thick_flat[n_shells]",    "Ang",            100.0,  [0, inf],       "volume", "flat layer_thickness"],
-              ["func_inter[n_shells]",    "",               0,      [0, 5],         "", "Erf:0, RPower:1, LPower:2, RExp:3, LExp:4, Linear:5"],
-              ["thick_inter[n_shells]",   "Ang",            50.0,   [0, inf],       "volume", "intern layer thickness"],
-              ["nu_inter[n_shells]",      "",               2.5,    [-inf, inf],    "", "steepness parameter"],
+parameters = [["n_shells",             "",           1,      [1, 11],        "volume", "number of shells"],
+              ["sld_solvent",          "1e-6/Ang^2", 1.0,    [-inf, inf],    "sld", "solvent sld"],
+              ["sld[n_shells]",        "1e-6/Ang^2", 4.06,   [-inf, inf],    "sld", "sld of the shell"],
+              ["thickness[n_shells]",  "Ang",        100.0,  [0, inf],       "volume", "thickness shell"],
+              ["interface[n_shells]",  "Ang",        50.0,   [0, inf],       "volume", "thickness of the interface"],
+              ["shape[n_shells]",      "",           0,      SHAPES,         "", "interface shape"],
+              ["nu[n_shells]",         "",           2.5,    [0, inf],       "", "interface shape exponent"],
+              ["n_steps",              "",           35,     [0, inf],       "", "number of steps in each interface (must be an odd integer)"],
               ]
 # pylint: enable=bad-whitespace, line-too-long
-source = ["lib/polevl.c", "lib/sas_erf.c", "lib/librefl.c",  "lib/sph_j1c.c", "spherical_sld.c"]
+source = ["lib/polevl.c", "lib/sas_erf.c", "lib/sph_j1c.c", "spherical_sld.c"]
 single = False  # TODO: fix low q behaviour
 
 profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)']
-def profile(n_shells, radius_core,  sld_core,  sld_solvent, sld_flat,
-            thick_flat, func_inter, thick_inter, nu_inter, npts_inter):
+
+SHAPE_FUNCTIONS = [
+    lambda z, nu: erf(nu/sqrt(2)*(2*z-1))/(2*erf(nu/sqrt(2))) + 0.5,  # erf
+    lambda z, nu: z**nu,                    # Rpow
+    lambda z, nu: 1 - (1-z)**nu,            # Lpow
+    lambda z, nu: expm1(-nu*z)/expm1(-nu),  # Rexp
+    lambda z, nu: expm1(nu*z)/expm1(nu),    # Lexp
+]
+
+def profile(n_shells, sld_solvent, sld, thickness,
+            interface, shape, nu, n_steps):
     """
     Returns shape profile with x=radius, y=SLD.
@@ -228 +236 @@
 
     z = []
-    beta = []
+    rho = []
     z0 = 0
     # two sld points for core
     z.append(0)
-    beta.append(sld_core)
-    z.append(radius_core)
-    beta.append(sld_core)
-    z0 += radius_core
-
-    for i in range(1, n_shells+2):
-        dz = thick_inter[i-1]/npts_inter
-        # j=0 for interface, j=1 for flat layer
-        for j in range(0, 2):
-            # interation for sub-layers
-            for n_s in range(0, npts_inter+1):
-                if j == 1:
-                    if i == n_shells+1:
-                        break
-                    # shift half sub thickness for the first point
-                    z0 -= dz#/2.0
-                    z.append(z0)
-                    #z0 -= dz/2.0
-                    z0 += thick_flat[i]
-                    sld_i = sld_flat[i]
-                    beta.append(sld_flat[i])
-                    dz = 0
-                else:
-                    nu = nu_inter[i-1]
-                    # decide which sld is which, sld_r or sld_l
-                    if i == 1:
-                        sld_l = sld_core
-                    else:
-                        sld_l = sld_flat[i-1]
-                    if i == n_shells+1:
-                        sld_r = sld_solvent
-                    else:
-                        sld_r = sld_flat[i]
-                    # get function type
-                    func_idx = func_inter[i-1]
-                    # calculate the sld
-                    sld_i = intersldfunc(func_idx, npts_inter, n_s, nu,
-                                            sld_l, sld_r)
-                # append to the list
-                z.append(z0)
-                beta.append(sld_i)
-                z0 += dz
-                if j == 1:
-                    break
-    z.append(z0)
-    beta.append(sld_solvent)
-    z_ext = z0/5.0
-    z.append(z0+z_ext)
-    beta.append(sld_solvent)
+    rho.append(sld[0])
+
+    for i in range(0, n_shells):
+        z0 += thickness[i]
+        z.append(z0)
+        rho.append(sld[i])
+        dz = interface[i]/n_steps
+        sld_l = sld[i]
+        sld_r = sld[i+1] if i < n_shells-1 else sld_solvent
+        interface = SHAPE_FUNCTIONS[int(np.clip(shape[i], 0, len(SHAPES)-1))]
+        for step in range(1, n_steps+1):
+            portion = interface(float(step)/n_steps, max(abs(nu[i]), 1e-14))
+            z0 += dz
+            z.append(z0)
+            rho.append((sld_r - sld_l)*portion + sld_l)
+    z.append(z0*1.2)
+    rho.append(sld_solvent)
     # return sld profile (r, beta)
-    return np.asarray(z), np.asarray(beta)*1e-6
-
-def ER(n_shells, radius_core, thick_inter0, thick_inter, thick_flat):
+    return np.asarray(z), np.asarray(rho)*1e-6
+
+
+def ER(n_shells, thickness, interface):
     n_shells = int(n_shells)
-    total_thickness = thick_inter0
-    total_thickness += np.sum(thick_inter[:n_shells], axis=0)
-    total_thickness += np.sum(thick_flat[:n_shells], axis=0)
-    return total_thickness + radius_core
+    total = (np.sum(thickness[:n_shells], axis=1)
+             + np.sum(interface[:n_shells], axis=1))
+    return total
 
 
 demo = {
-    "n_shells": 4,
-    "npts_inter": 35.0,
-    "radius_core": 50.0,
-    "sld_core": 2.07,
+    "n_shells": 5,
+    "n_steps": 35.0,
     "sld_solvent": 1.0,
-    "thick_inter0": 50.0,
-    "func_inter0": 0,
-    "nu_inter0": 2.5,
-    "sld_flat":[4.0,3.5,4.0,3.5],
-    "thick_flat":[100.0,100.0,100.0,100.0],
-    "func_inter":[0,0,0,0],
-    "thick_inter":[50.0,50.0,50.0,50.0],
-    "nu_inter":[2.5,2.5,2.5,2.5],
+    "sld":[2.07,4.0,3.5,4.0,3.5],
+    "thickness":[50.0,100.0,100.0,100.0,100.0],
+    "interface":[50.0,50.0,50.0,50.0],
+    "shape": [0,0,0,0,0],
+    "nu":[2.5,2.5,2.5,2.5,2.5],
     }
 
@@ -312 +283 @@
 tests = [
     # Accuracy tests based on content in test/utest_extra_models.py
-    [{"n_shells":4,
-        'npts_inter':35,
-        "radius_core":50.0,
-        "sld_core":2.07,
+    [{"n_shells": 5,
+        "n_steps": 35,
         "sld_solvent": 1.0,
-        "sld_flat":[4.0,3.5,4.0,3.5],
-        "thick_flat":[100.0,100.0,100.0,100.0],
-        "func_inter":[0,0,0,0],
-        "thick_inter":[50.0,50.0,50.0,50.0],
-        "nu_inter":[2.5,2.5,2.5,2.5]
+        "sld": [2.07, 4.0, 3.5, 4.0, 3.5],
+        "thickness": [50.0, 100.0, 100.0, 100.0, 100.0],
+        "interface": [50]*5,
+        "shape": [0]*5,
+        "nu": [2.5]*5,
     }, 0.001, 0.001],
 ]