Changes in / [17a8c94:9a7ef88] in sasmodels


Ignore:
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • doc/guide/pd/polydispersity.rst

    rf4ae8c4 r29afc50  
    2828sigmas $N_\sigma$ to include from the tails of the distribution, and the 
    2929number of points used to compute the average. The center of the distribution 
    30 is set by the value of the model parameter. 
    31  
    32 Volume parameters have polydispersity *PD* (not to be confused with a 
    33 molecular weight distributions in polymer science), but orientation parameters 
    34 use angular distributions of width $\sigma$. 
     30is set by the value of the model parameter. The meaning of a polydispersity  
     31parameter *PD* (not to be confused with a molecular weight distributions  
     32in polymer science) in a model depends on the type of parameter it is being  
     33applied too. 
     34 
     35The distribution width applied to *volume* (ie, shape-describing) parameters  
     36is relative to the center value such that $\sigma = \mathrm{PD} \cdot \bar x$.  
     37However, the distribution width applied to *orientation* (ie, angle-describing)  
     38parameters is just $\sigma = \mathrm{PD}$. 
    3539 
    3640$N_\sigma$ determines how far into the tails to evaluate the distribution, 
    … …  
    6973or angular orientations, use the Gaussian or Boltzmann distributions. 
    7074 
     75If applying polydispersion to parameters describing angles, use the Uniform  
     76distribution. Beware of using distributions that are always positive (eg, the  
     77Lognormal) because angles can be negative! 
     78 
    7179The array distribution allows a user-defined distribution to be applied. 
    7280 
    … …  
    215223The polydispersity in sasmodels is given by 
    216224 
    217 .. math:: \text{PD} = p = \sigma / x_\text{med} 
    218  
    219 The mean value of the distribution is given by $\bar x = \exp(\mu+ p^2/2)$ 
    220 and the peak value by $\max x = \exp(\mu - p^2)$. 
     225.. math:: \text{PD} = \sigma = p / x_\text{med} 
     226 
     227The mean value of the distribution is given by $\bar x = \exp(\mu+ \sigma^2/2)$ 
     228and the peak value by $\max x = \exp(\mu - \sigma^2)$. 
    221229 
    222230The variance (the square of the standard deviation) of the *lognormal* 
    … …  
    232240.. figure:: pd_lognormal.jpg 
    233241 
    234     Lognormal distribution. 
     242    Lognormal distribution for PD=0.1. 
    235243 
    236244For further information on the Lognormal distribution see: 
    … …  
    334342| 2017-05-08 Paul Kienzle 
    335343| 2018-03-20 Steve King 
     344| 2018-04-04 Steve King 

  • explore/asymint.py

    ra1c32c2 rc462169  
    3030import numpy as np 
    3131import mpmath as mp 
    32 from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10 
     32from numpy import pi, sin, cos, sqrt, exp, expm1, degrees, log10, arccos 
    3333from numpy.polynomial.legendre import leggauss 
    3434from scipy.integrate import dblquad, simps, romb, romberg 
    … …  
    4141shape = 'parallelepiped' if len(sys.argv) < 2 else sys.argv[1] 
    4242Qstr = '0.005' if len(sys.argv) < 3 else sys.argv[2] 
     43 
     44DTYPE = 'd' 
    4345 
    4446class MPenv: 
    … …  
    7274    sas_sinx_x = staticmethod(sp.sas_sinx_x) 
    7375    pi = np.pi 
    74     mpf = staticmethod(float) 
     76    #mpf = staticmethod(float) 
     77    mpf = staticmethod(lambda x: np.array(x, DTYPE)) 
    7578 
    7679SLD = 3 
    … …  
    220223    #A, B, C = 4450, 14000, 47 
    221224    #A, B, C = 445, 140, 47  # integer for the sake of mpf 
    222     A, B, C = 6800, 114, 1380 
    223     DA, DB, DC = 2300, 21, 58 
     225    A, B, C = 114, 1380, 6800 
     226    DA, DB, DC = 21, 58, 2300 
    224227    SLDA, SLDB, SLDC = "5", "-0.3", "11.5" 
     228    ## default parameters from sasmodels 
     229    #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 400,75,35,10,10,10,2,4,2 
     230    ## swap A-B-C to C-B-A 
     231    #A, B, C, DA, DB, DC, SLDA, SLDB, SLDC = C, B, A, DC, DB, DA, SLDC, SLDB, SLDA 
    225232    #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 10,20,30,100,200,300,1,2,3 
    226233    #SLD_SOLVENT,CONTRAST = 0, 4 
    … …  
    233240        DA, DC = DC, DA 
    234241        SLDA, SLDC = SLDC, SLDA 
     242    #NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 
    235243    NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 
    236244    NORM_MP, KERNEL_MP = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC, env=MPenv) 
    … …  
    348356    return n**2, Iq 
    349357 
     358def gauss_quad_usub(q, n=150, dtype=DTYPE): 
     359    """ 
     360    Compute the integral using gaussian quadrature for n = 20, 76 or 150. 
     361 
     362    Use *u = sin theta* substitution, and restrict integration over a single 
     363    quadrant for shapes that are mirror symmetric about AB, AC and BC planes. 
     364 
     365    Note that this doesn't work for fcc/bcc paracrystals, which instead step 
     366    over the entire 4 pi surface uniformly in theta-phi. 
     367    """ 
     368    z, w = leggauss(n) 
     369    cos_theta = 0.5 * (z + 1) 
     370    theta = arccos(cos_theta) 
     371    phi = pi/2*(0.5 * (z + 1)) 
     372    Atheta, Aphi = np.meshgrid(theta, phi) 
     373    Aw = w[None, :] * w[:, None] 
     374    q, Atheta, Aphi, Aw = [np.asarray(v, dtype=dtype) for v in (q, Atheta, Aphi, Aw)] 
     375    Zq = kernel_2d(q=q, theta=Atheta, phi=Aphi) 
     376    Iq = np.sum(Zq*Aw)*0.25 
     377    return n**2, Iq 
     378 
    350379def gridded_2d(q, n=300): 
    351380    """ 
    … …  
    395424    print("gauss-1025", *gauss_quad_2d(Q, n=1025)) 
    396425    print("gauss-2049", *gauss_quad_2d(Q, n=2049)) 
     426    print("gauss-20 usub", *gauss_quad_usub(Q, n=20)) 
     427    print("gauss-76 usub", *gauss_quad_usub(Q, n=76)) 
     428    print("gauss-150 usub", *gauss_quad_usub(Q, n=150)) 
    397429    #gridded_2d(Q, n=2**8+1) 
    398430    gridded_2d(Q, n=2**10+1) 

  • sasmodels/compare.py

    rd86f0fc r5770493  
    752752    """ 
    753753    for k in range(opts['sets']): 
    754         if k > 1: 
     754        if k > 0: 
    755755            # print a separate seed for each dataset for better reproducibility 
    756756            new_seed = np.random.randint(1000000) 
    757             print("Set %d uses -random=%i"%(k+1, new_seed)) 
     757            print("=== Set %d uses -random=%i ==="%(k+1, new_seed)) 
    758758            np.random.seed(new_seed) 
    759759        opts['pars'] = parse_pars(opts, maxdim=maxdim) 
    … …  
    762762        result = run_models(opts, verbose=True) 
    763763        if opts['plot']: 
     764            if opts['is2d'] and k > 0: 
     765                import matplotlib.pyplot as plt 
     766                plt.figure() 
    764767            limits = plot_models(opts, result, limits=limits, setnum=k) 
    765768        if opts['show_weights']: 
    … …  
    13291332 
    13301333    # Evaluate preset parameter expressions 
     1334    # Note: need to replace ':' with '_' in parameter names and expressions 
     1335    # in order to support math on magnetic parameters. 
    13311336    context = MATH.copy() 
    13321337    context['np'] = np 
    1333     context.update(pars) 
     1338    context.update((k.replace(':', '_'), v) for k, v in pars.items()) 
    13341339    context.update((k, v) for k, v in presets.items() if isinstance(v, float)) 
     1340    #for k,v in sorted(context.items()): print(k, v) 
    13351341    for k, v in presets.items(): 
    13361342        if not isinstance(v, float) and not k.endswith('_type'): 
    1337             presets[k] = eval(v, context) 
     1343            presets[k] = eval(v.replace(':', '_'), context) 
    13381344    context.update(presets) 
    1339     context.update((k, v) for k, v in presets2.items() if isinstance(v, float)) 
     1345    context.update((k.replace(':', '_'), v) for k, v in presets2.items() if isinstance(v, float)) 
    13401346    for k, v in presets2.items(): 
    13411347        if not isinstance(v, float) and not k.endswith('_type'): 
    1342             presets2[k] = eval(v, context) 
     1348            presets2[k] = eval(v.replace(':', '_'), context) 
    13431349 
    13441350    # update parameters with presets 

  • sasmodels/details.py

    r108e70e r885753a  
    243243    offset = np.cumsum(np.hstack((0, length))) 
    244244    call_details = make_details(kernel.info, length, offset[:-1], offset[-1]) 
    245     # Pad value array to a 32 value boundaryd 
     245    # Pad value array to a 32 value boundary 
    246246    data_len = nvalues + 2*sum(len(v) for v in dispersity) 
    247247    extra = (32 - data_len%32)%32 
    … …  
    250250    is_magnetic = convert_magnetism(kernel.info.parameters, data) 
    251251    #call_details.show() 
     252    #print("data", data) 
    252253    return call_details, data, is_magnetic 
    253254 
    … …  
    296297    mag = values[parameters.nvalues-3*parameters.nmagnetic:parameters.nvalues] 
    297298    mag = mag.reshape(-1, 3) 
    298     scale = mag[:, 0] 
    299     if np.any(scale): 
     299    if np.any(mag[:, 0] != 0.0): 
     300        M0 = mag[:, 0].copy() 
    300301        theta, phi = radians(mag[:, 1]), radians(mag[:, 2]) 
    301         cos_theta = cos(theta) 
    302         mag[:, 0] = scale*cos_theta*cos(phi)  # mx 
    303         mag[:, 1] = scale*sin(theta)  # my 
    304         mag[:, 2] = -scale*cos_theta*sin(phi)  # mz 
     302        mag[:, 0] = +M0*cos(theta)*cos(phi)  # mx 
     303        mag[:, 1] = +M0*sin(theta) # my 
     304        mag[:, 2] = -M0*cos(theta)*sin(phi)  # mz 
    305305        return True 
    306306    else: 

  • sasmodels/kernel_iq.c

    rd86f0fc rdc6f601  
    7979//     du * (m_sigma_y + 1j*m_sigma_z); 
    8080// weights for spin crosssections: dd du real, ud real, uu, du imag, ud imag 
    81 static void set_spin_weights(double in_spin, double out_spin, double spins[4]) 
     81static void set_spin_weights(double in_spin, double out_spin, double weight[6]) 
    8282{ 
    8383  in_spin = clip(in_spin, 0.0, 1.0); 
    8484  out_spin = clip(out_spin, 0.0, 1.0); 
    85   spins[0] = sqrt(sqrt((1.0-in_spin) * (1.0-out_spin))); // dd 
    86   spins[1] = sqrt(sqrt((1.0-in_spin) * out_spin));       // du real 
    87   spins[2] = sqrt(sqrt(in_spin * (1.0-out_spin)));       // ud real 
    88   spins[3] = sqrt(sqrt(in_spin * out_spin));             // uu 
    89   spins[4] = spins[1]; // du imag 
    90   spins[5] = spins[2]; // ud imag 
     85  // Note: sasview 3.1 scaled all slds by sqrt(weight) and assumed that 
     86  //     w*I(q, rho1, rho2, ...) = I(q, sqrt(w)*rho1, sqrt(w)*rho2, ...) 
     87  // which is likely to be the case for simple models. 
     88  weight[0] = sqrt((1.0-in_spin) * (1.0-out_spin)); // dd 
     89  weight[1] = sqrt((1.0-in_spin) * out_spin);       // du.real 
     90  weight[2] = sqrt(in_spin * (1.0-out_spin));       // ud.real 
     91  weight[3] = sqrt(in_spin * out_spin);             // uu 
     92  weight[4] = weight[1]; // du.imag 
     93  weight[5] = weight[2]; // ud.imag 
    9194} 
    9295 
    9396// Compute the magnetic sld 
    9497static double mag_sld( 
    95   const unsigned int xs, // 0=dd, 1=du real, 2=ud real, 3=uu, 4=du imag, 5=up imag 
     98  const unsigned int xs, // 0=dd, 1=du.real, 2=ud.real, 3=uu, 4=du.imag, 5=ud.imag 
    9699  const double qx, const double qy, 
    97100  const double px, const double py, 
    … …  
    106109      case 0: // uu => sld - D M_perpx 
    107110          return sld - px*perp; 
    108       case 1: // ud real => -D M_perpy 
     111      case 1: // ud.real => -D M_perpy 
    109112          return py*perp; 
    110       case 2: // du real => -D M_perpy 
     113      case 2: // du.real => -D M_perpy 
    111114          return py*perp; 
    112       case 3: // dd real => sld + D M_perpx 
     115      case 3: // dd => sld + D M_perpx 
    113116          return sld + px*perp; 
    114117    } 
    115118  } else { 
    116119    if (xs== 4) { 
    117       return -mz;  // ud imag => -D M_perpz 
     120      return -mz;  // du.imag => +D M_perpz 
    118121    } else { // index == 5 
    119       return mz;   // du imag => D M_perpz 
     122      return +mz;  // ud.imag => -D M_perpz 
    120123    } 
    121124  } 
    … …  
    312315  //     up_angle = values[NUM_PARS+4]; 
    313316  // TODO: could precompute more magnetism parameters before calling the kernel. 
    314   double spins[8];  // uu, ud real, du real, dd, ud imag, du imag, fill, fill 
     317  double xs_weights[8];  // uu, ud real, du real, dd, ud imag, du imag, fill, fill 
    315318  double cos_mspin, sin_mspin; 
    316   set_spin_weights(values[NUM_PARS+2], values[NUM_PARS+3], spins); 
     319  set_spin_weights(values[NUM_PARS+2], values[NUM_PARS+3], xs_weights); 
    317320  SINCOS(-values[NUM_PARS+4]*M_PI_180, sin_mspin, cos_mspin); 
    318321#endif // MAGNETIC 
    … …  
    661664            // loop over uu, ud real, du real, dd, ud imag, du imag 
    662665            for (unsigned int xs=0; xs<6; xs++) { 
    663               const double xs_weight = spins[xs]; 
     666              const double xs_weight = xs_weights[xs]; 
    664667              if (xs_weight > 1.e-8) { 
    665668                // Since the cross section weight is significant, set the slds 
    … …  
    674677                  local_values.vector[sld_index] = 
    675678                    mag_sld(xs, qx, qy, px, py, values[sld_index+2], mx, my, mz); 
     679//if (q_index==0) printf("%d: (qx,qy)=(%g,%g) xs=%d sld%d=%g p=(%g,%g) m=(%g,%g,%g)\n", 
     680//  q_index, qx, qy, xs, sk, local_values.vector[sld_index], px, py, mx, my, mz); 
    676681                } 
    677682                scattering += xs_weight * CALL_KERNEL(); 

  • sasmodels/sasview_model.py

    r3221de0 r7c3fb15  
    593593            # Check whether we have a list of ndarrays [qx,qy] 
    594594            qx, qy = qdist 
    595             if not self._model_info.parameters.has_2d: 
    596                 return self.calculate_Iq(np.sqrt(qx ** 2 + qy ** 2)) 
    597             else: 
    598                 return self.calculate_Iq(qx, qy) 
     595            return self.calculate_Iq(qx, qy) 
    599596 
    600597        elif isinstance(qdist, np.ndarray): 
    … …  
    677674        call_details, values, is_magnetic = make_kernel_args(calculator, pairs) 
    678675        #call_details.show() 
    679         #print("pairs", pairs) 
     676        #print("================ parameters ==================") 
     677        #for p, v in zip(parameters.call_parameters, pairs): print(p.name, v[0]) 
    680678        #for k, p in enumerate(self._model_info.parameters.call_parameters): 
    681679        #    print(k, p.name, *pairs[k]) 
    … …  
    721719 
    722720    def set_dispersion(self, parameter, dispersion): 
    723         # type: (str, weights.Dispersion) -> Dict[str, Any] 
     721        # type: (str, weights.Dispersion) -> None 
    724722        """ 
    725723        Set the dispersion object for a model parameter 
    … …  
    744742            self.dispersion[parameter] = dispersion.get_pars() 
    745743        else: 
    746             raise ValueError("%r is not a dispersity or orientation parameter") 
     744            raise ValueError("%r is not a dispersity or orientation parameter" 
     745                             % parameter) 
    747746 
    748747    def _dispersion_mesh(self): 
    … …  
    871870    CylinderModel().evalDistribution([0.1, 0.1]) 
    872871 
     872def magnetic_demo(): 
     873    Model = _make_standard_model('sphere') 
     874    model = Model() 
     875    model.setParam('M0:sld', 8) 
     876    q = np.linspace(-0.35, 0.35, 500) 
     877    qx, qy = np.meshgrid(q, q) 
     878    result = model.calculate_Iq(qx.flatten(), qy.flatten()) 
     879    result = result.reshape(qx.shape) 
     880 
     881    import pylab 
     882    pylab.imshow(np.log(result + 0.001)) 
     883    pylab.show() 
     884 
    873885if __name__ == "__main__": 
    874886    print("cylinder(0.1,0.1)=%g"%test_cylinder()) 
     887    #magnetic_demo() 
    875888    #test_product() 
    876889    #test_structure_factor() 
Note: See TracChangeset for help on using the changeset viewer.