Changeset ba51e00 in sasmodels
- Timestamp:
- Sep 8, 2018 10:30:57 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:
- c44b611
- Parents:
- 2773c66 (diff), fbaef04 (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. - Files:
-
- 17 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/models/core_multi_shell.c
rd277229 ra94046f 32 32 static double 33 33 effective_radius(int mode, double core_radius, double fp_n, double thickness[]) 34 // this seems regardless to always give the result for outer radius for n=1 shells; why?? 35 // printf shows fp_n is always 1, not 0,1,2 34 36 { 35 if (mode == 1) { 37 // printf("fp_n =%g \n",fp_n); 38 if (mode == 1) { 36 39 double r = core_radius; 37 40 int n = (int)(fp_n+0.5); 38 for (int i=0; i < n; i++) { 39 r += thickness[i]; 41 if ( n > 0) { 42 for (int i=0; i < n; i++) { 43 r += thickness[i]; 44 } 40 45 } 41 46 return r; -
sasmodels/models/core_shell_cylinder.c
rd277229 ra94046f 24 24 { 25 25 const double radius_outer = radius + thickness; 26 const double length_outer = length + thickness;26 const double length_outer = length + 2.0*thickness; 27 27 return sqrt(radius_outer*radius_outer + 0.25*length_outer*length_outer); 28 28 } … … 30 30 static double 31 31 effective_radius(int mode, double radius, double thickness, double length) 32 //effective_radius_type = ["equivalent sphere","outer radius","half outer length","half min outer dimension", 33 // "half max outer dimension","half outer diagonal"] 32 34 { 33 35 if (mode == 1) { -
sasmodels/models/core_shell_parallelepiped.c
rd277229 ra94046f 45 45 effective_radius(int mode, double length_a, double length_b, double length_c, 46 46 double thick_rim_a, double thick_rim_b, double thick_rim_c) 47 //effective_radius_type = ["equivalent sphere","half outer length_a", "half outer length_b", "half outer length_c", 48 // "equivalent circular cross-section","half outer ab diagonal","half outer diagonal"] 49 // note the core box is A*B*C with slabs ta, tb & tc on each face but missing the corners, though that fact is ignored here 50 // in the equvalent sphere option 47 51 { 48 52 if (mode == 1) { 49 53 return radius_from_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 50 54 } else if (mode == 2) { 51 return 0.5 * (length_a + thick_rim_a);55 return 0.5 * length_a + thick_rim_a; 52 56 } else if (mode == 3) { 53 return 0.5 * (length_b + thick_rim_b);57 return 0.5 * length_b + thick_rim_b; 54 58 } else if (mode == 4) { 55 return 0.5 * (length_c + thick_rim_c);59 return 0.5 * length_c + thick_rim_c; 56 60 } else if (mode == 5) { 57 61 return radius_from_crosssection(length_a, length_b, thick_rim_a, thick_rim_b); 58 62 } else if (mode == 6) { 59 return 0.5*sqrt(square(length_a+ thick_rim_a) + square(length_b+thick_rim_b));63 return 0.5*sqrt(square(length_a+ 2.0*thick_rim_a) + square(length_b+ 2.0*thick_rim_b)); 60 64 } else { 61 return 0.5*sqrt(square(length_a+ thick_rim_a) + square(length_b+thick_rim_b) + square(length_c+thick_rim_c));65 return 0.5*sqrt(square(length_a+ 2.0*thick_rim_a) + square(length_b+ 2.0*thick_rim_b) + square(length_c+ 2.0*thick_rim_c)); 62 66 } 63 67 } -
sasmodels/models/elliptical_cylinder.c
rd277229 rfbaef04 13 13 14 14 static double 15 radius_from_min_dimension(double radius_minor, double r_ratio, double length)15 radius_from_min_dimension(double radius_minor, double r_ratio, double hlength) 16 16 { 17 17 const double rad_min = (r_ratio > 1.0 ? radius_minor : r_ratio*radius_minor); 18 return (rad_min < length ? rad_min :length);18 return (rad_min < hlength ? rad_min : hlength); 19 19 } 20 20 21 21 static double 22 radius_from_max_dimension(double radius_minor, double r_ratio, double length)22 radius_from_max_dimension(double radius_minor, double r_ratio, double hlength) 23 23 { 24 24 const double rad_max = (r_ratio < 1.0 ? radius_minor : r_ratio*radius_minor); 25 return (rad_max > length ? rad_max :length);25 return (rad_max > hlength ? rad_max : hlength); 26 26 } 27 27 … … 35 35 static double 36 36 effective_radius(int mode, double radius_minor, double r_ratio, double length) 37 //effective_radius_type = ["equivalent sphere","average radius","min radius","max radius", 38 // "equivalent circular cross-section","half length","half min dimension","half max dimension","half diagonal"] 37 39 { 38 40 if (mode == 1) { … … 49 51 return 0.5*length; 50 52 } else if (mode == 7) { 51 return radius_from_min_dimension(radius_minor,r_ratio, length);53 return radius_from_min_dimension(radius_minor,r_ratio,0.5*length); 52 54 } else if (mode == 8) { 53 return radius_from_max_dimension(radius_minor,r_ratio, length);55 return radius_from_max_dimension(radius_minor,r_ratio,0.5*length); 54 56 } else { 55 57 return radius_from_diagonal(radius_minor,r_ratio,length); -
sasmodels/models/hollow_rectangular_prism.c
rd277229 ra94046f 15 15 static double 16 16 effective_radius(int mode, double length_a, double b2a_ratio, double c2a_ratio, double thickness) 17 //effective_radius_type = ["equivalent sphere","half length_a", "half length_b", "half length_c", 18 // "equivalent outer circular cross-section","half ab diagonal","half diagonal"] 19 // NOTE length_a is external dimension! 17 20 { 18 21 if (mode == 1) { -
sasmodels/models/hollow_rectangular_prism.py
rd277229 ra94046f 126 126 "Solvent scattering length density"], 127 127 ["length_a", "Ang", 35, [0, inf], "volume", 128 "Shorte r side of the parallelepiped"],128 "Shortest, external, size of the parallelepiped"], 129 129 ["b2a_ratio", "Ang", 1, [0, inf], "volume", 130 130 "Ratio sides b/a"], -
doc/guide/scripting.rst
r4aa5dce rbd7630d 10 10 The key functions are :func:`sasmodels.core.load_model` for loading the 11 11 model definition and compiling the kernel and 12 :func:`sasmodels.data.load_data` for calling sasview to load the data. Need 13 the data because that defines the resolution function and the q values to 14 evaluate. If there is no data, then use :func:`sasmodels.data.empty_data1D` 15 or :func:`sasmodels.data.empty_data2D` to create some data with a given $q$. 16 17 Using sasmodels through bumps 18 ============================= 19 20 With the data and the model, you can wrap it in a *bumps* model with 12 :func:`sasmodels.data.load_data` for calling sasview to load the data. 13 14 Preparing data 15 ============== 16 17 Usually you will load data via the sasview loader, with the 18 :func:`sasmodels.data.load_data` function. For example:: 19 20 from sasmodels.data import load_data 21 data = load_data("sasmodels/example/093191_201.dat") 22 23 You may want to apply a data mask, such a beam stop, and trim high $q$:: 24 25 from sasmodels.data import set_beam_stop 26 set_beam_stop(data, qmin, qmax) 27 28 The :func:`sasmodels.data.set_beam_stop` method simply sets the *mask* 29 attribute for the data. 30 31 The data defines the resolution function and the q values to evaluate, so 32 even if you simulating experiments prior to making measurements, you still 33 need a data object for reference. Use :func:`sasmodels.data.empty_data1D` 34 or :func:`sasmodels.data.empty_data2D` to create a container with a 35 given $q$ and $\Delta q/q$. For example:: 36 37 import numpy as np 38 from sasmodels.data import empty_data1D 39 40 # 120 points logarithmically spaced from 0.005 to 0.2, with dq/q = 5% 41 q = np.logspace(np.log10(5e-3), np.log10(2e-1), 120) 42 data = empty_data1D(q, resolution=0.05) 43 44 To use a more realistic model of resolution, or to load data from a file 45 format not understood by SasView, you can use :class:`sasmodels.data.Data1D` 46 or :class:`sasmodels.data.Data2D` directly. The 1D data uses 47 *x*, *y*, *dx* and *dy* for $x = q$ and $y = I(q)$, and 2D data uses 48 *x*, *y*, *z*, *dx*, *dy*, *dz* for $x, y = qx, qy$ and $z = I(qx, qy)$. 49 [Note: internally, the Data2D object uses SasView conventions, 50 *qx_data*, *qy_data*, *data*, *dqx_data*, *dqy_data*, and *err_data*.] 51 52 For USANS data, use 1D data, but set *dxl* and *dxw* attributes to 53 indicate slit resolution:: 54 55 data.dxl = 0.117 56 57 See :func:`sasmodels.resolution.slit_resolution` for details. 58 59 SESANS data is more complicated; if your SESANS format is not supported by 60 SasView you need to define a number of attributes beyond *x*, *y*. For 61 example:: 62 63 SElength = np.linspace(0, 2400, 61) # [A] 64 data = np.ones_like(SElength) 65 err_data = np.ones_like(SElength)*0.03 66 67 class Source: 68 wavelength = 6 # [A] 69 wavelength_unit = "A" 70 class Sample: 71 zacceptance = 0.1 # [A^-1] 72 thickness = 0.2 # [cm] 73 74 class SESANSData1D: 75 #q_zmax = 0.23 # [A^-1] 76 lam = 0.2 # [nm] 77 x = SElength 78 y = data 79 dy = err_data 80 sample = Sample() 81 data = SESANSData1D() 82 83 x, y = ... # create or load sesans 84 data = smd.Data 85 86 The *data* module defines various data plotters as well. 87 88 Using sasmodels directly 89 ======================== 90 91 Once you have a computational kernel and a data object, you can evaluate 92 the model for various parameters using 93 :class:`sasmodels.direct_model.DirectModel`. The resulting object *f* 94 will be callable as *f(par=value, ...)*, returning the $I(q)$ for the $q$ 95 values in the data. For example:: 96 97 import numpy as np 98 from sasmodels.data import empty_data1D 99 from sasmodels.core import load_model 100 from sasmodels.direct_model import DirectModel 101 102 # 120 points logarithmically spaced from 0.005 to 0.2, with dq/q = 5% 103 q = np.logspace(np.log10(5e-3), np.log10(2e-1), 120) 104 data = empty_data1D(q, resolution=0.05) 105 kernel = load_model("ellipsoid) 106 f = DirectModel(data, kernel) 107 Iq = f(radius_polar=100) 108 109 Polydispersity information is set with special parameter names: 110 111 * *par_pd* for polydispersity width, $\Delta p/p$, 112 * *par_pd_n* for the number of points in the distribution, 113 * *par_pd_type* for the distribution type (as a string), and 114 * *par_pd_nsigmas* for the limits of the distribution. 115 116 Using sasmodels through the bumps optimizer 117 =========================================== 118 119 Like DirectModel, you can wrap data and a kernel in a *bumps* model with 21 120 class:`sasmodels.bumps_model.Model` and create an 22 class:`sasmodels.bump _model.Experiment` that you can fit with the *bumps*121 class:`sasmodels.bumps_model.Experiment` that you can fit with the *bumps* 23 122 interface. Here is an example from the *example* directory such as 24 123 *example/model.py*:: … … 75 174 SasViewCom bumps.cli example/model.py --preview 76 175 77 Using sasmodels directly 78 ======================== 79 80 Bumps has a notion of parameter boxes in which you can set and retrieve 81 values. Instead of using bumps, you can create a directly callable function 82 with :class:`sasmodels.direct_model.DirectModel`. The resulting object *f* 83 will be callable as *f(par=value, ...)*, returning the $I(q)$ for the $q$ 84 values in the data. Polydisperse parameters use the same naming conventions 85 as in the bumps model, with e.g., radius_pd being the polydispersity associated 86 with radius. 176 Calling the computation kernel 177 ============================== 87 178 88 179 Getting a simple function that you can call on a set of q values and return -
sasmodels/compare.py
raf7a97c rbd7630d 1288 1288 1289 1289 if opts['datafile'] is not None: 1290 data = load_data(os.path.expanduser(opts['datafile'])) 1290 data0 = load_data(os.path.expanduser(opts['datafile'])) 1291 data = data0, data0 1291 1292 else: 1292 1293 # Hack around the fact that make_data doesn't take a pair of resolutions -
sasmodels/data.py
r7e923c2 rc88f983 183 183 *x* is spin echo length and *y* is polarization (P/P0). 184 184 """ 185 isSesans = True 185 186 def __init__(self, **kw): 186 187 Data1D.__init__(self, **kw) … … 301 302 self.wavelength_unit = "A" 302 303 304 class Sample(object): 305 """ 306 Sample attributes. 307 """ 308 def __init__(self): 309 # type: () -> None 310 pass 303 311 304 312 def empty_data1D(q, resolution=0.0, L=0., dL=0.): -
sasmodels/direct_model.py
r1a8c11c r7b9e4dd 31 31 from . import resolution2d 32 32 from .details import make_kernel_args, dispersion_mesh 33 from .modelinfo import DEFAULT_BACKGROUND 33 34 34 35 # pylint: disable=unused-import … … 349 350 350 351 # Need to pull background out of resolution for multiple scattering 351 background = pars.get('background', 0.)352 background = pars.get('background', DEFAULT_BACKGROUND) 352 353 pars = pars.copy() 353 354 pars['background'] = 0. -
sasmodels/generate.py
r5399809 r2773c66 991 991 docs = model_info.docs if model_info.docs is not None else "" 992 992 docs = convert_section_titles_to_boldface(docs) 993 pars = make_partable(model_info.parameters.COMMON 994 + model_info.parameters.kernel_parameters) 993 if model_info.structure_factor: 994 pars = model_info.parameters.kernel_parameters 995 else: 996 pars = model_info.parameters.COMMON + model_info.parameters.kernel_parameters 997 partable = make_partable(pars) 995 998 subst = dict(id=model_info.id.replace('_', '-'), 996 999 name=model_info.name, 997 1000 title=model_info.title, 998 parameters=par s,1001 parameters=partable, 999 1002 returns=Sq_units if model_info.structure_factor else Iq_units, 1000 1003 docs=docs) -
sasmodels/kernel_iq.c
rc57ee9e r2773c66 192 192 QACRotation *rotation, 193 193 double qx, double qy, 194 double *qa _out, double *qc_out)194 double *qab_out, double *qc_out) 195 195 { 196 // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2 196 197 const double dqc = rotation->R31*qx + rotation->R32*qy; 197 // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2 198 const double dqa = sqrt(-dqc*dqc + qx*qx + qy*qy); 199 200 *qa_out = dqa; 198 const double dqab_sq = -dqc*dqc + qx*qx + qy*qy; 199 //*qab_out = sqrt(fabs(dqab_sq)); 200 *qab_out = dqab_sq > 0.0 ? sqrt(dqab_sq) : 0.0; 201 201 *qc_out = dqc; 202 202 } -
sasmodels/modelinfo.py
r5399809 r2773c66 45 45 # Note that scale and background cannot be coordinated parameters whose value 46 46 # depends on the some polydisperse parameter with the current implementation 47 DEFAULT_BACKGROUND = 1e-3 47 48 COMMON_PARAMETERS = [ 48 49 ("scale", "", 1, (0.0, np.inf), "", "Source intensity"), 49 ("background", "1/cm", 1e-3, (-np.inf, np.inf), "", "Source background"),50 ("background", "1/cm", DEFAULT_BACKGROUND, (-np.inf, np.inf), "", "Source background"), 50 51 ] 51 52 assert (len(COMMON_PARAMETERS) == 2 -
sasmodels/models/ellipsoid.py
rd277229 r2773c66 125 125 import numpy as np 126 126 from numpy import inf, sin, cos, pi 127 128 try: 129 from numpy import cbrt 130 except ImportError: 131 def cbrt(x): return x ** (1.0/3.0) 127 132 128 133 name = "ellipsoid" -
sasmodels/models/spinodal.py
r71b751d rc88f983 3 3 ---------- 4 4 5 This model calculates the SAS signal of a phase separating solution 6 under spinodal decomposition. The scattering intensity $I(q)$ is calculated as 5 This model calculates the SAS signal of a phase separating system 6 undergoing spinodal decomposition. The scattering intensity $I(q)$ is calculated 7 as 7 8 8 9 .. math:: 9 10 I(q) = I_{max}\frac{(1+\gamma/2)x^2}{\gamma/2+x^{2+\gamma}}+B 10 11 11 where $x=q/q_0$ and $B$ is a flat background. The characteristic structure 12 length scales with the correlation peak at $q_0$. The exponent $\gamma$ is 13 equal to $d+1$ with d the dimensionality of the off-critical concentration 14 mixtures. A transition to $\gamma=2d$ is seen near the percolation threshold 15 into the critical concentration regime. 12 where $x=q/q_0$, $q_0$ is the peak position, $I_{max}$ is the intensity 13 at $q_0$ (parameterised as the $scale$ parameter), and $B$ is a flat 14 background. The spinodal wavelength is given by $2\pi/q_0$. 15 16 The exponent $\gamma$ is equal to $d+1$ for off-critical concentration 17 mixtures (smooth interfaces) and $2d$ for critical concentration mixtures 18 (entangled interfaces), where $d$ is the dimensionality (ie, 1, 2, 3) of the 19 system. Thus 2 <= $\gamma$ <= 6. A transition from $\gamma=d+1$ to $\gamma=2d$ 20 is expected near the percolation threshold. 21 22 As this function tends to zero as $q$ tends to zero, in practice it may be 23 necessary to combine it with another function describing the low-angle 24 scattering, or to simply omit the low-angle scattering from the fit. 16 25 17 26 References … … 22 31 Physica A 123,497 (1984). 23 32 24 Authorship and Verification 25 ---------------- ------------33 Revision History 34 ---------------- 26 35 27 * **Author:** Dirk Honecker **Date:** Oct 7, 2016 36 * **Author:** Dirk Honecker **Date:** Oct 7, 2016 37 * **Revised:** Steve King **Date:** Sep 7, 2018 28 38 """ 29 39 … … 34 44 title = "Spinodal decomposition model" 35 45 description = """\ 36 I(q) = scale ((1+gamma/2)x^2)/(gamma/2+x^(2+gamma))+background46 I(q) = Imax ((1+gamma/2)x^2)/(gamma/2+x^(2+gamma)) + background 37 47 38 48 List of default parameters: 39 scale = scaling 40 gamma = exponent 41 x = q/q_0 49 50 Imax = correlation peak intensity at q_0 51 background = incoherent background 52 gamma = exponent (see model documentation) 42 53 q_0 = correlation peak position [1/A] 43 background = Incoherent background""" 54 x = q/q_0""" 55 44 56 category = "shape-independent" 45 57 -
sasmodels/product.py
r33d7be3 r2773c66 42 42 pars = [] 43 43 if p_info.have_Fq: 44 par = parse_parameter("structure_factor_mode", "", 0, [["P*S","P*(1+beta*(S-1))"]], "", 45 "Structure factor calculation") 44 par = parse_parameter( 45 "structure_factor_mode", 46 "", 47 0, 48 [["P*S","P*(1+beta*(S-1))"]], 49 "", 50 "Structure factor calculation") 46 51 pars.append(par) 47 52 if p_info.effective_radius_type is not None: 48 par = parse_parameter("radius_effective_mode", "", 0, [["unconstrained"] + p_info.effective_radius_type], "", 49 "Effective radius calculation") 53 par = parse_parameter( 54 "radius_effective_mode", 55 "", 56 0, 57 [["unconstrained"] + p_info.effective_radius_type], 58 "", 59 "Effective radius calculation") 50 60 pars.append(par) 51 61 return pars … … 169 179 )) 170 180 171 def _intermediates_beta(F1, F2, S, scale, bg, effective_radius): 172 # type: (np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, float) -> OrderedDict[str, Union[np.ndarray, float]] 181 def _intermediates_beta(F1, # type: np.ndarray 182 F2, # type: np.ndarray 183 S, # type: np.ndarray 184 scale, # type: np.ndarray 185 bg, # type: np.ndarray 186 effective_radius # type: float 187 ): 188 # type: (...) -> OrderedDict[str, Union[np.ndarray, float]] 173 189 """ 174 190 Returns intermediate results for beta approximation-enabled product. The result may be an array or a float. -
sasmodels/sasview_model.py
raa44a6a r2773c66 543 543 # pylint: disable=unpacking-non-sequence 544 544 q, phi = x 545 return self.calculate_Iq([q*math.cos(phi)], [q*math.sin(phi)])[0] 546 else: 547 return self.calculate_Iq([x])[0] 545 result, _ = self.calculate_Iq([q*math.cos(phi)], [q*math.sin(phi)]) 546 return result[0] 547 else: 548 result, _ = self.calculate_Iq([x]) 549 return result[0] 548 550 549 551 … … 560 562 """ 561 563 if isinstance(x, (list, tuple)): 562 return self.calculate_Iq([x[0]], [x[1]])[0] 563 else: 564 return self.calculate_Iq([x])[0] 564 result, _ = self.calculate_Iq([x[0]], [x[1]]) 565 return result[0] 566 else: 567 result, _ = self.calculate_Iq([x]) 568 return result[0] 565 569 566 570 def evalDistribution(self, qdist): … … 596 600 # Check whether we have a list of ndarrays [qx,qy] 597 601 qx, qy = qdist 598 return self.calculate_Iq(qx, qy) 602 result, _ = self.calculate_Iq(qx, qy) 603 return result 599 604 600 605 elif isinstance(qdist, np.ndarray): 601 606 # We have a simple 1D distribution of q-values 602 return self.calculate_Iq(qdist) 607 result, _ = self.calculate_Iq(qdist) 608 return result 603 609 604 610 else: … … 640 646 if composition and composition[0] == 'product': # only P*S for now 641 647 with calculation_lock: 642 self._calculate_Iq(qx) 643 # for compatibility with sasview 4.3 644 results = self._intermediate_results() 645 return results["P(Q)"], results["S(Q)"] 648 _, lazy_results = self._calculate_Iq(qx) 649 # for compatibility with sasview 4.x 650 results = lazy_results() 651 pq_data = results.get("P(Q)") 652 sq_data = results.get("S(Q)") 653 return pq_data, sq_data 646 654 else: 647 655 return None 648 656 649 def calculate_Iq(self, qx, qy=None): 650 # type: (Sequence[float], Optional[Sequence[float]]) -> np.ndarray 657 def calculate_Iq(self, 658 qx, # type: Sequence[float] 659 qy=None # type: Optional[Sequence[float]] 660 ): 661 # type: (...) -> Tuple[np.ndarray, Callable[[], collections.OrderedDict[str, np.ndarray]]] 651 662 """ 652 663 Calculate Iq for one set of q with the current parameters. … … 656 667 This should NOT be used for fitting since it copies the *q* vectors 657 668 to the card for each evaluation. 669 670 The returned tuple contains the scattering intensity followed by a 671 callable which returns a dictionary of intermediate data from 672 ProductKernel. 658 673 """ 659 674 ## uncomment the following when trying to debug the uncoordinated calls … … 688 703 result = calculator(call_details, values, cutoff=self.cutoff, 689 704 magnetic=is_magnetic) 705 lazy_results = getattr(calculator, 'results', 706 lambda: collections.OrderedDict()) 690 707 #print("result", result) 691 self._intermediate_results = getattr(calculator, 'results', None) 708 692 709 calculator.release() 693 710 #self._model.release() 694 return result 711 712 return result, lazy_results 695 713 696 714 def calculate_ER(self): … … 881 899 q = np.linspace(-0.35, 0.35, 500) 882 900 qx, qy = np.meshgrid(q, q) 883 result = model.calculate_Iq(qx.flatten(), qy.flatten())901 result, _ = model.calculate_Iq(qx.flatten(), qy.flatten()) 884 902 result = result.reshape(qx.shape) 885 903
Note: See TracChangeset
for help on using the changeset viewer.