# Changeset 9a7ef88 in sasmodels

Ignore:
Timestamp:
Apr 17, 2018 6:16:45 AM (6 years ago)
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
Parents:
17a8c94 (diff), 7c3fb15 (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 ticket-896

Files:
10 edited

Unmodified
Removed
• ## doc/guide/pd/polydispersity.rst

 rf4ae8c4 sigmas $N_\sigma$ to include from the tails of the distribution, and the number of points used to compute the average. The center of the distribution is set by the value of the model parameter. Volume parameters have polydispersity *PD* (not to be confused with a molecular weight distributions in polymer science), but orientation parameters use angular distributions of width $\sigma$. is set by the value of the model parameter. The meaning of a polydispersity parameter *PD* (not to be confused with a molecular weight distributions in polymer science) in a model depends on the type of parameter it is being applied too. The distribution width applied to *volume* (ie, shape-describing) parameters is relative to the center value such that $\sigma = \mathrm{PD} \cdot \bar x$. However, the distribution width applied to *orientation* (ie, angle-describing) parameters is just $\sigma = \mathrm{PD}$. $N_\sigma$ determines how far into the tails to evaluate the distribution, or angular orientations, use the Gaussian or Boltzmann distributions. If applying polydispersion to parameters describing angles, use the Uniform distribution. Beware of using distributions that are always positive (eg, the Lognormal) because angles can be negative! The array distribution allows a user-defined distribution to be applied. The polydispersity in sasmodels is given by .. math:: \text{PD} = p = \sigma / x_\text{med} The mean value of the distribution is given by $\bar x = \exp(\mu+ p^2/2)$ and the peak value by $\max x = \exp(\mu - p^2)$. .. math:: \text{PD} = \sigma = p / x_\text{med} The mean value of the distribution is given by $\bar x = \exp(\mu+ \sigma^2/2)$ and the peak value by $\max x = \exp(\mu - \sigma^2)$. The variance (the square of the standard deviation) of the *lognormal* .. figure:: pd_lognormal.jpg Lognormal distribution. Lognormal distribution for PD=0.1. For further information on the Lognormal distribution see: | 2017-05-08 Paul Kienzle | 2018-03-20 Steve King | 2018-04-04 Steve King

• ## sasmodels/compare.py

 rd86f0fc """ for k in range(opts['sets']): if k > 1: if k > 0: # print a separate seed for each dataset for better reproducibility new_seed = np.random.randint(1000000) print("Set %d uses -random=%i"%(k+1, new_seed)) print("=== Set %d uses -random=%i ==="%(k+1, new_seed)) np.random.seed(new_seed) opts['pars'] = parse_pars(opts, maxdim=maxdim) result = run_models(opts, verbose=True) if opts['plot']: if opts['is2d'] and k > 0: import matplotlib.pyplot as plt plt.figure() limits = plot_models(opts, result, limits=limits, setnum=k) if opts['show_weights']: # Evaluate preset parameter expressions # Note: need to replace ':' with '_' in parameter names and expressions # in order to support math on magnetic parameters. context = MATH.copy() context['np'] = np context.update(pars) context.update((k.replace(':', '_'), v) for k, v in pars.items()) context.update((k, v) for k, v in presets.items() if isinstance(v, float)) #for k,v in sorted(context.items()): print(k, v) for k, v in presets.items(): if not isinstance(v, float) and not k.endswith('_type'): presets[k] = eval(v, context) presets[k] = eval(v.replace(':', '_'), context) context.update(presets) context.update((k, v) for k, v in presets2.items() if isinstance(v, float)) context.update((k.replace(':', '_'), v) for k, v in presets2.items() if isinstance(v, float)) for k, v in presets2.items(): if not isinstance(v, float) and not k.endswith('_type'): presets2[k] = eval(v, context) presets2[k] = eval(v.replace(':', '_'), context) # update parameters with presets
• ## sasmodels/details.py

 r108e70e offset = np.cumsum(np.hstack((0, length))) call_details = make_details(kernel.info, length, offset[:-1], offset[-1]) # Pad value array to a 32 value boundaryd # Pad value array to a 32 value boundary data_len = nvalues + 2*sum(len(v) for v in dispersity) extra = (32 - data_len%32)%32 is_magnetic = convert_magnetism(kernel.info.parameters, data) #call_details.show() #print("data", data) return call_details, data, is_magnetic mag = values[parameters.nvalues-3*parameters.nmagnetic:parameters.nvalues] mag = mag.reshape(-1, 3) scale = mag[:, 0] if np.any(scale): if np.any(mag[:, 0] != 0.0): M0 = mag[:, 0].copy() theta, phi = radians(mag[:, 1]), radians(mag[:, 2]) cos_theta = cos(theta) mag[:, 0] = scale*cos_theta*cos(phi)  # mx mag[:, 1] = scale*sin(theta)  # my mag[:, 2] = -scale*cos_theta*sin(phi)  # mz mag[:, 0] = +M0*cos(theta)*cos(phi)  # mx mag[:, 1] = +M0*sin(theta) # my mag[:, 2] = -M0*cos(theta)*sin(phi)  # mz return True else:
• ## sasmodels/kernel_iq.c

 rd86f0fc //     du * (m_sigma_y + 1j*m_sigma_z); // weights for spin crosssections: dd du real, ud real, uu, du imag, ud imag static void set_spin_weights(double in_spin, double out_spin, double spins[4]) static void set_spin_weights(double in_spin, double out_spin, double weight[6]) { in_spin = clip(in_spin, 0.0, 1.0); out_spin = clip(out_spin, 0.0, 1.0); spins[0] = sqrt(sqrt((1.0-in_spin) * (1.0-out_spin))); // dd spins[1] = sqrt(sqrt((1.0-in_spin) * out_spin));       // du real spins[2] = sqrt(sqrt(in_spin * (1.0-out_spin)));       // ud real spins[3] = sqrt(sqrt(in_spin * out_spin));             // uu spins[4] = spins[1]; // du imag spins[5] = spins[2]; // ud imag // Note: sasview 3.1 scaled all slds by sqrt(weight) and assumed that //     w*I(q, rho1, rho2, ...) = I(q, sqrt(w)*rho1, sqrt(w)*rho2, ...) // which is likely to be the case for simple models. weight[0] = sqrt((1.0-in_spin) * (1.0-out_spin)); // dd weight[1] = sqrt((1.0-in_spin) * out_spin);       // du.real weight[2] = sqrt(in_spin * (1.0-out_spin));       // ud.real weight[3] = sqrt(in_spin * out_spin);             // uu weight[4] = weight[1]; // du.imag weight[5] = weight[2]; // ud.imag } // Compute the magnetic sld static double mag_sld( const unsigned int xs, // 0=dd, 1=du real, 2=ud real, 3=uu, 4=du imag, 5=up imag const unsigned int xs, // 0=dd, 1=du.real, 2=ud.real, 3=uu, 4=du.imag, 5=ud.imag const double qx, const double qy, const double px, const double py, case 0: // uu => sld - D M_perpx return sld - px*perp; case 1: // ud real => -D M_perpy case 1: // ud.real => -D M_perpy return py*perp; case 2: // du real => -D M_perpy case 2: // du.real => -D M_perpy return py*perp; case 3: // dd real => sld + D M_perpx case 3: // dd => sld + D M_perpx return sld + px*perp; } } else { if (xs== 4) { return -mz;  // ud imag => -D M_perpz return -mz;  // du.imag => +D M_perpz } else { // index == 5 return mz;   // du imag => D M_perpz return +mz;  // ud.imag => -D M_perpz } } //     up_angle = values[NUM_PARS+4]; // TODO: could precompute more magnetism parameters before calling the kernel. double spins[8];  // uu, ud real, du real, dd, ud imag, du imag, fill, fill double xs_weights[8];  // uu, ud real, du real, dd, ud imag, du imag, fill, fill double cos_mspin, sin_mspin; set_spin_weights(values[NUM_PARS+2], values[NUM_PARS+3], spins); set_spin_weights(values[NUM_PARS+2], values[NUM_PARS+3], xs_weights); SINCOS(-values[NUM_PARS+4]*M_PI_180, sin_mspin, cos_mspin); #endif // MAGNETIC // loop over uu, ud real, du real, dd, ud imag, du imag for (unsigned int xs=0; xs<6; xs++) { const double xs_weight = spins[xs]; const double xs_weight = xs_weights[xs]; if (xs_weight > 1.e-8) { // Since the cross section weight is significant, set the slds local_values.vector[sld_index] = mag_sld(xs, qx, qy, px, py, values[sld_index+2], mx, my, mz); //if (q_index==0) printf("%d: (qx,qy)=(%g,%g) xs=%d sld%d=%g p=(%g,%g) m=(%g,%g,%g)\n", //  q_index, qx, qy, xs, sk, local_values.vector[sld_index], px, py, mx, my, mz); } scattering += xs_weight * CALL_KERNEL();
• ## sasmodels/sasview_model.py

 r3221de0 # Check whether we have a list of ndarrays [qx,qy] qx, qy = qdist if not self._model_info.parameters.has_2d: return self.calculate_Iq(np.sqrt(qx ** 2 + qy ** 2)) else: return self.calculate_Iq(qx, qy) return self.calculate_Iq(qx, qy) elif isinstance(qdist, np.ndarray): call_details, values, is_magnetic = make_kernel_args(calculator, pairs) #call_details.show() #print("pairs", pairs) #print("================ parameters ==================") #for p, v in zip(parameters.call_parameters, pairs): print(p.name, v[0]) #for k, p in enumerate(self._model_info.parameters.call_parameters): #    print(k, p.name, *pairs[k]) def set_dispersion(self, parameter, dispersion): # type: (str, weights.Dispersion) -> Dict[str, Any] # type: (str, weights.Dispersion) -> None """ Set the dispersion object for a model parameter self.dispersion[parameter] = dispersion.get_pars() else: raise ValueError("%r is not a dispersity or orientation parameter") raise ValueError("%r is not a dispersity or orientation parameter" % parameter) def _dispersion_mesh(self): CylinderModel().evalDistribution([0.1, 0.1]) def magnetic_demo(): Model = _make_standard_model('sphere') model = Model() model.setParam('M0:sld', 8) q = np.linspace(-0.35, 0.35, 500) qx, qy = np.meshgrid(q, q) result = model.calculate_Iq(qx.flatten(), qy.flatten()) result = result.reshape(qx.shape) import pylab pylab.imshow(np.log(result + 0.001)) pylab.show() if __name__ == "__main__": print("cylinder(0.1,0.1)=%g"%test_cylinder()) #magnetic_demo() #test_product() #test_structure_factor()
• ## sasmodels/models/core_shell_parallelepiped.c

 re077231 // outer integral (with gauss points), integration limits = 0, 1 // substitute d_cos_alpha for sin_alpha d_alpha double outer_sum = 0; //initialize integral for( int i=0; i
• ## sasmodels/models/core_shell_parallelepiped.py

 r97be877 .. math:: I(q) = \text{scale}\frac{\langle f^2 \rangle}{V} + \text{background} I(q) = \frac{\text{scale}}{V} \langle P(q,\alpha,\beta) \rangle + \text{background} where $\langle \ldots \rangle$ is an average over all possible orientations of the rectangular solid. The function calculated is the form factor of the rectangular solid below. of the rectangular solid, and the usual $\Delta \rho^2 \ V^2$ term cannot be pulled out of the form factor term due to the multiple slds in the model. The core of the solid is defined by the dimensions $A$, $B$, $C$ such that $A < B < C$. .. image:: img/core_shell_parallelepiped_geometry.jpg .. figure:: img/parallelepiped_geometry.jpg Core of the core shell parallelepiped with the corresponding definition of sides. There are rectangular "slabs" of thickness $t_A$ that add to the $A$ dimension (on the $BC$ faces). There are similar slabs on the $AC$ $(=t_B)$ and $AB$ $(=t_C)$ faces. The projection in the $AB$ plane is then .. image:: img/core_shell_parallelepiped_projection.jpg The volume of the solid is $(=t_C)$ faces. The projection in the $AB$ plane is .. figure:: img/core_shell_parallelepiped_projection.jpg AB cut through the core-shell parallelipiped showing the cross secion of four of the six shell slabs. As can be seen, this model leaves **"gaps"** at the corners of the solid. The total volume of the solid is thus given as .. math:: V = ABC + 2t_ABC + 2t_BAC + 2t_CAB **meaning that there are "gaps" at the corners of the solid.** The intensity calculated follows the :ref:parallelepiped model, with the core-shell intensity being calculated as the square of the sum of the amplitudes of the core and the slabs on the edges. the scattering amplitude is computed for a particular orientation of the core-shell parallelepiped with respect to the scattering vector and then averaged over all possible orientations, where $\alpha$ is the angle between the $z$ axis and the $C$ axis of the parallelepiped, $\beta$ is the angle between projection of the particle in the $xy$ detector plane and the $y$ axis. .. math:: F(Q) amplitudes of the core and the slabs on the edges. The scattering amplitude is computed for a particular orientation of the core-shell parallelepiped with respect to the scattering vector and then averaged over all possible orientations, where $\alpha$ is the angle between the $z$ axis and the $C$ axis of the parallelepiped, and $\beta$ is the angle between the projection of the particle in the $xy$ detector plane and the $y$ axis. .. math:: P(q)=\frac {\int_{0}^{\pi/2}\int_{0}^{\pi/2}F^2(q,\alpha,\beta) \ sin\alpha \ d\alpha \ d\beta} {\int_{0}^{\pi/2} \ sin\alpha \ d\alpha \ d\beta} and .. math:: F(q,\alpha,\beta) &= (\rho_\text{core}-\rho_\text{solvent}) S(Q_A, A) S(Q_B, B) S(Q_C, C) \\ &+ (\rho_\text{A}-\rho_\text{solvent}) \left[S(Q_A, A+2t_A) - S(Q_A, Q)\right] S(Q_B, B) S(Q_C, C) \\ \left[S(Q_A, A+2t_A) - S(Q_A, A)\right] S(Q_B, B) S(Q_C, C) \\ &+ (\rho_\text{B}-\rho_\text{solvent}) S(Q_A, A) \left[S(Q_B, B+2t_B) - S(Q_B, B)\right] S(Q_C, C) \\ .. math:: S(Q, L) = L \frac{\sin \tfrac{1}{2} Q L}{\tfrac{1}{2} Q L} S(Q_X, L) = L \frac{\sin (\tfrac{1}{2} Q_X L)}{\tfrac{1}{2} Q_X L} and .. math:: Q_A &= \sin\alpha \sin\beta \\ Q_B &= \sin\alpha \cos\beta \\ Q_C &= \cos\alpha Q_A &= q \sin\alpha \sin\beta \\ Q_B &= q \sin\alpha \cos\beta \\ Q_C &= q \cos\alpha where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$ are the scattering length of the parallelepiped core, and the rectangular are the scattering lengths of the parallelepiped core, and the rectangular slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$ is the scattering length of the solvent. .. note:: the code actually implements two substitutions: $d(cos\alpha)$ is substituted for -$sin\alpha \ d\alpha$ (note that in the :ref:parallelepiped code this is explicitly implemented with $\sigma = cos\alpha$), and $\beta$ is set to $\beta = u \pi/2$ so that $du = \pi/2 \ d\beta$.  Thus both integrals go from 0 to 1 rather than 0 to $\pi/2$. FITTING NOTES based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ and length $(C+2t_C)$ values, after appropriately sorting the three dimensions to give an oblate or prolate particle, to give an effective radius, for $S(Q)$ when $P(Q) * S(Q)$ is applied. to give an oblate or prolate particle, to give an effective radius for $S(q)$ when $P(q) * S(q)$ is applied. For 2d data the orientation of the particle is required, described using angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below. For further details of the calculation and angular dispersions see :ref:orientation. The angle $\Psi$ is the rotational angle around the *long_c* axis. For example, $\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector. For 2d, constraints must be applied during fitting to ensure that the inequality $A < B < C$ is not violated, and hence the correct definition of angles is preserved. The calculation will not report an error, but the results may be not correct. .. note:: For 2d, constraints must be applied during fitting to ensure that the inequality $A < B < C$ is not violated, and hence the correct definition of angles is preserved. The calculation will not report an error, but the results may be not correct. .. figure:: img/parallelepiped_angle_definition.png Note that rotation $\theta$, initially in the $xz$ plane, is carried out first, then rotation $\phi$ about the $z$ axis, finally rotation $\Psi$ is now around the axis of the cylinder. The neutron or X-ray $\Psi$ is now around the axis of the particle. The neutron or X-ray beam is along the $z$ axis. Examples of the angles for oriented core-shell parallelepipeds against the detector plane. Validation ---------- Cross-checked against hollow rectangular prism and rectangular prism for equal thickness overlapping sides, and by Monte Carlo sampling of points within the shape for non-uniform, non-overlapping sides. References * **Author:** NIST IGOR/DANSE **Date:** pre 2010 * **Converted to sasmodels by:** Miguel Gonzales **Date:** February 26, 2016 * **Converted to sasmodels by:** Miguel Gonzalez **Date:** February 26, 2016 * **Last Modified by:** Paul Kienzle **Date:** October 17, 2017 * Cross-checked against hollow rectangular prism and rectangular prism for equal thickness overlapping sides, and by Monte Carlo sampling of points within the shape for non-uniform, non-overlapping sides. """
• ## sasmodels/models/parallelepiped.c

 r108e70e inner_total += GAUSS_W[j] * square(si1 * si2); } // now complete change of inner integration variable (1-0)/(1-(-1))= 0.5 inner_total *= 0.5; outer_total += GAUSS_W[i] * inner_total * si * si; } // now complete change of outer integration variable (1-0)/(1-(-1))= 0.5 outer_total *= 0.5;
• ## sasmodels/models/parallelepiped.py

 ref07e95 # Note: model title and parameter table are inserted automatically r""" The form factor is normalized by the particle volume. For information about polarised and magnetic scattering, see the :ref:magnetism documentation. Definition ---------- This model calculates the scattering from a rectangular parallelepiped (\:numref:parallelepiped-image\). If you need to apply polydispersity, see also :ref:rectangular-prism. (:numref:parallelepiped-image). If you need to apply polydispersity, see also :ref:rectangular-prism. For information about polarised and magnetic scattering, see the :ref:magnetism documentation. .. _parallelepiped-image: error, or fixing of some dimensions at expected values, may help. The 1D scattering intensity $I(q)$ is calculated as: The form factor is normalized by the particle volume and the 1D scattering intensity $I(q)$ is then calculated as: .. Comment by Miguel Gonzalez: I(q) = \frac{\text{scale}}{V} (\Delta\rho \cdot V)^2 \left< P(q, \alpha) \right> + \text{background} \left< P(q, \alpha, \beta) \right> + \text{background} where the volume $V = A B C$, the contrast is defined as $\Delta\rho = \rho_\text{p} - \rho_\text{solvent}$, $P(q, \alpha)$ is the form factor corresponding to a parallelepiped oriented at an angle $\alpha$ (angle between the long axis C and $\vec q$), and the averaging $\left<\ldots\right>$ is applied over all orientations. $\Delta\rho = \rho_\text{p} - \rho_\text{solvent}$, $P(q, \alpha, \beta)$ is the form factor corresponding to a parallelepiped oriented at an angle $\alpha$ (angle between the long axis C and $\vec q$), and $\beta$ ( the angle between the projection of the particle in the $xy$ detector plane and the $y$ axis) and the averaging $\left<\ldots\right>$ is applied over all orientations. Assuming $a = A/B < 1$, $b = B /B = 1$, and $c = C/B > 1$, the form factor is given by (Mittelbach and Porod, 1961) form factor is given by (Mittelbach and Porod, 1961 [#Mittelbach]_) .. math:: \mu &= qB The scattering intensity per unit volume is returned in units of |cm^-1|. where substitution of $\sigma = cos\alpha$ and $\beta = \pi/2 \ u$ have been applied. NB: The 2nd virial coefficient of the parallelepiped is calculated based on .. math:: P(q_x, q_y) = \left[\frac{\sin(\tfrac{1}{2}qA\cos\alpha)}{(\tfrac{1}{2}qA\cos\alpha)}\right]^2 \left[\frac{\sin(\tfrac{1}{2}qB\cos\beta)}{(\tfrac{1}{2}qB\cos\beta)}\right]^2 \left[\frac{\sin(\tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1}{2}qC\cos\gamma)}\right]^2 P(q_x, q_y) = \left[\frac{\sin(\tfrac{1}{2}qA\cos\alpha)}{(\tfrac{1} {2}qA\cos\alpha)}\right]^2 \left[\frac{\sin(\tfrac{1}{2}qB\cos\beta)}{(\tfrac{1} {2}qB\cos\beta)}\right]^2 \left[\frac{\sin(\tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1} {2}qC\cos\gamma)}\right]^2 with ---------- P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211 R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 .. [#Mittelbach] P Mittelbach and G Porod, *Acta Physica Austriaca*, 14 (1961) 185-211 .. [#] R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 Authorship and Verification
Note: See TracChangeset for help on using the changeset viewer.