Changeset a3412a6 in sasmodels
- Timestamp:
- Mar 6, 2019 6:08:08 PM (6 years ago)
- Branches:
- ticket_1156
- Children:
- bd91e8f
- Parents:
- f752b9b (diff), 9150036 (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:
-
- 15 edited
Legend:
- Unmodified
- Added
- Removed
-
doc/guide/plugin.rst
raa8c6e0 r9150036 272 272 structure factor to account for interactions between particles. See 273 273 `Form_Factors`_ for more details. 274 275 **model_info = ...** lets you define a model directly, for example, by 276 loading and modifying existing models. This is done implicitly by 277 :func:`sasmodels.core.load_model_info`, which can create a mixture model 278 from a pair of existing models. For example:: 279 280 from sasmodels.core import load_model_info 281 model_info = load_model_info('sphere+cylinder') 282 283 See :class:`sasmodels.modelinfo.ModelInfo` for details about the model 284 attributes that are defined. 274 285 275 286 Model Parameters … … 894 905 - \frac{\sin(x)}{x}\left(\frac{1}{x} - \frac{3!}{x^3} + \frac{5!}{x^5} - \frac{7!}{x^7}\right) 895 906 896 For small arguments 907 For small arguments, 897 908 898 909 .. math:: -
example/multiscatfit.py
r49d1f8b8 r2c4a190 15 15 16 16 # Show the model without fitting 17 PYTHONPATH=..:../ explore:../../bumps:../../sasview/src python multiscatfit.py17 PYTHONPATH=..:../../bumps:../../sasview/src python multiscatfit.py 18 18 19 19 # Run the fit 20 PYTHONPATH=..:../ explore:../../bumps:../../sasview/src ../../bumps/run.py \20 PYTHONPATH=..:../../bumps:../../sasview/src ../../bumps/run.py \ 21 21 multiscatfit.py --store=/tmp/t1 22 22 … … 55 55 ) 56 56 57 # Tie the model to the data 58 M = Experiment(data=data, model=model) 59 60 # Stack mulitple scattering on top of the existing resolution function. 61 M.resolution = MultipleScattering(resolution=M.resolution, probability=0.) 62 57 63 # SET THE FITTING PARAMETERS 58 64 model.radius_polar.range(15, 3000) … … 65 71 model.scale.range(0, 0.1) 66 72 67 # Mulitple scattering probability parameter 68 # HACK: the probability is stuffed in as an extra parameter to the experiment. 69 probability = Parameter(name="probability", value=0.0) 70 probability.range(0.0, 0.9) 73 # The multiple scattering probability parameter is in the resolution function 74 # instead of the scattering function, so access it through M.resolution 75 M.scattering_probability.range(0.0, 0.9) 71 76 72 M = Experiment(data=data, model=model, extra_pars={'probability': probability}) 73 74 # Stack mulitple scattering on top of the existing resolution function. 75 # Because resolution functions in sasview don't have fitting parameters, 76 # we instead allow the multiple scattering calculator to take a function 77 # instead of a probability. This function returns the current value of 78 # the parameter. ** THIS IS TEMPORARY ** when multiple scattering is 79 # properly integrated into sasmodels and sasview, its fittable parameter 80 # will be treated like the model parameters. 81 M.resolution = MultipleScattering(resolution=M.resolution, 82 probability=lambda: probability.value, 83 ) 84 M._kernel_inputs = M.resolution.q_calc 77 # Let bumps know that we are fitting this experiment 85 78 problem = FitProblem(M) 86 79 -
sasmodels/__init__.py
ra1ec908 r37f38ff 14 14 defining new models. 15 15 """ 16 __version__ = "0.9 8"16 __version__ = "0.99" 17 17 18 18 def data_files(): -
sasmodels/bumps_model.py
r49d1f8b8 r2c4a190 35 35 # when bumps is not on the path. 36 36 from bumps.names import Parameter # type: ignore 37 from bumps.parameter import Reference # type: ignore 37 38 except ImportError: 38 39 pass … … 139 140 def __init__(self, data, model, cutoff=1e-5, name=None, extra_pars=None): 140 141 # type: (Data, Model, float) -> None 142 # Allow resolution function to define fittable parameters. We do this 143 # by creating reference parameters within the resolution object rather 144 # than modifying the object itself to use bumps parameters. We need 145 # to reset the parameters each time the object has changed. These 146 # additional parameters need to be returned from the fitting engine. 147 # To make them available to the user, they are added as top-level 148 # attributes to the experiment object. The only change to the 149 # resolution function is that it needs an optional 'fittable' attribute 150 # which maps the internal name to the user visible name for the 151 # for the parameter. 152 self._resolution = None 153 self._resolution_pars = {} 141 154 # remember inputs so we can inspect from outside 142 155 self.name = data.filename if name is None else name … … 145 158 self._interpret_data(data, model.sasmodel) 146 159 self._cache = {} 160 # CRUFT: no longer need extra parameters 161 # Multiple scattering probability is now retrieved directly from the 162 # multiple scattering resolution function. 147 163 self.extra_pars = extra_pars 148 164 … … 162 178 return len(self.Iq) 163 179 180 @property 181 def resolution(self): 182 return self._resolution 183 184 @resolution.setter 185 def resolution(self, value): 186 self._resolution = value 187 188 # Remove old resolution fitting parameters from experiment 189 for name in self._resolution_pars: 190 delattr(self, name) 191 192 # Create new resolution fitting parameters 193 res_pars = getattr(self._resolution, 'fittable', {}) 194 self._resolution_pars = { 195 name: Reference(self._resolution, refname, name=name) 196 for refname, name in res_pars.items() 197 } 198 199 # Add new resolution fitting parameters as experiment attributes 200 for name, ref in self._resolution_pars.items(): 201 setattr(self, name, ref) 202 164 203 def parameters(self): 165 204 # type: () -> Dict[str, Parameter] … … 168 207 """ 169 208 pars = self.model.parameters() 170 if self.extra_pars :209 if self.extra_pars is not None: 171 210 pars.update(self.extra_pars) 211 pars.update(self._resolution_pars) 172 212 return pars 173 213 -
sasmodels/direct_model.py
rc1799d3 r9150036 224 224 else: 225 225 Iq, dIq = None, None 226 #self._theory = np.zeros_like(q)227 q_vectors = [res.q_calc]228 226 elif self.data_type == 'Iqxy': 229 227 #if not model.info.parameters.has_2d: … … 242 240 res = resolution2d.Pinhole2D(data=data, index=index, 243 241 nsigma=3.0, accuracy=accuracy) 244 #self._theory = np.zeros_like(self.Iq)245 q_vectors = res.q_calc246 242 elif self.data_type == 'Iq': 247 243 index = (data.x >= data.qmin) & (data.x <= data.qmax) … … 268 264 else: 269 265 res = resolution.Perfect1D(data.x[index]) 270 271 #self._theory = np.zeros_like(self.Iq)272 q_vectors = [res.q_calc]273 266 elif self.data_type == 'Iq-oriented': 274 267 index = (data.x >= data.qmin) & (data.x <= data.qmax) … … 286 279 qx_width=data.dxw[index], 287 280 qy_width=data.dxl[index]) 288 q_vectors = res.q_calc289 281 else: 290 282 raise ValueError("Unknown data type") # never gets here … … 292 284 # Remember function inputs so we can delay loading the function and 293 285 # so we can save/restore state 294 self._kernel_inputs = q_vectors295 286 self._kernel = None 296 287 self.Iq, self.dIq, self.index = Iq, dIq, index … … 329 320 # type: (ParameterSet, float) -> np.ndarray 330 321 if self._kernel is None: 331 self._kernel = self._model.make_kernel(self._kernel_inputs) 322 # TODO: change interfaces so that resolution returns kernel inputs 323 # Maybe have resolution always return a tuple, or maybe have 324 # make_kernel accept either an ndarray or a pair of ndarrays. 325 kernel_inputs = self.resolution.q_calc 326 if isinstance(kernel_inputs, np.ndarray): 327 kernel_inputs = (kernel_inputs,) 328 self._kernel = self._model.make_kernel(kernel_inputs) 332 329 333 330 # Need to pull background out of resolution for multiple scattering -
sasmodels/multiscat.py
rb3703f5 r2c4a190 342 342 343 343 *probability* is related to the expected number of scattering 344 events in the sample $\lambda$ as $p = 1 = e^{-\lambda}$. As a 345 hack to allow probability to be a fitted parameter, the "value" 346 can be a function that takes no parameters and returns the current 347 value of the probability. *coverage* determines how many scattering 348 steps to consider. The default is 0.99, which sets $n$ such that 349 $1 \ldots n$ covers 99% of the Poisson probability mass function. 344 events in the sample $\lambda$ as $p = 1 - e^{-\lambda}$. 345 *coverage* determines how many scattering steps to consider. The 346 default is 0.99, which sets $n$ such that $1 \ldots n$ covers 99% 347 of the Poisson probability mass function. 350 348 351 349 *is2d* is True then 2D scattering is used, otherwise it accepts … … 399 397 self.qmin = qmin 400 398 self.nq = nq 401 self.probability = probability399 self.probability = 0. if probability is None else probability 402 400 self.coverage = coverage 403 401 self.is2d = is2d … … 456 454 self.Iqxy = None # type: np.ndarray 457 455 456 # Label probability as a fittable parameter, and give its external name 457 # Note that the external name must be a valid python identifier, since 458 # is will be set as an experiment attribute. 459 self.fittable = {'probability': 'scattering_probability'} 460 458 461 def apply(self, theory): 459 462 if self.is2d: … … 463 466 Iq_calc = Iq_calc.reshape(self.nq, self.nq) 464 467 468 # CRUFT: don't need probability as a function anymore 465 469 probability = self.probability() if callable(self.probability) else self.probability 466 470 coverage = self.coverage -
sasmodels/sasview_model.py
ra8a1d48 r9150036 25 25 from . import core 26 26 from . import custom 27 from . import kernelcl 27 28 from . import product 28 29 from . import generate … … 30 31 from . import modelinfo 31 32 from .details import make_kernel_args, dispersion_mesh 33 from .kernelcl import reset_environment 32 34 33 35 # pylint: disable=unused-import … … 68 70 #: has changed since we last reloaded. 69 71 _CACHED_MODULE = {} # type: Dict[str, "module"] 72 73 def reset_environment(): 74 # type: () -> None 75 """ 76 Clear the compute engine context so that the GUI can change devices. 77 78 This removes all compiled kernels, even those that are active on fit 79 pages, but they will be restored the next time they are needed. 80 """ 81 kernelcl.reset_environment() 82 for model in MODELS.values(): 83 model._model = None 70 84 71 85 def find_model(modelname): … … 696 710 def _calculate_Iq(self, qx, qy=None): 697 711 if self._model is None: 698 self._model = core.build_model(self._model_info) 712 # Only need one copy of the compiled kernel regardless of how many 713 # times it is used, so store it in the class. Also, to reset the 714 # compute engine, need to clear out all existing compiled kernels, 715 # which is much easier to do if we store them in the class. 716 self.__class__._model = core.build_model(self._model_info) 699 717 if qy is not None: 700 718 q_vectors = [np.asarray(qx), np.asarray(qy)] -
explore/realspace.py
r362a66f r5778279 9 9 import numpy as np 10 10 from numpy import pi, radians, sin, cos, sqrt 11 from numpy.random import poisson, uniform, randn, rand 11 from numpy.random import poisson, uniform, randn, rand, randint 12 12 from numpy.polynomial.legendre import leggauss 13 13 from scipy.integrate import simps … … 78 78 79 79 80 I3 = np.matrix([[1., 0, 0], [0, 1, 0], [0, 0, 1]]) 81 80 82 class Shape: 81 rotation = np.matrix([[1., 0, 0], [0, 1, 0], [0, 0, 1]])83 rotation = I3 82 84 center = np.array([0., 0., 0.])[:, None] 83 85 r_max = None 86 lattice_size = np.array((1, 1, 1)) 87 lattice_spacing = np.array((1., 1., 1.)) 88 lattice_distortion = 0.0 89 lattice_rotation = 0.0 90 lattice_type = "" 84 91 85 92 def volume(self): … … 96 103 97 104 def rotate(self, theta, phi, psi): 98 self.rotation = rotation(theta, phi, psi) * self.rotation 105 if theta != 0. or phi != 0. or psi != 0.: 106 self.rotation = rotation(theta, phi, psi) * self.rotation 99 107 return self 100 108 … … 103 111 return self 104 112 113 def lattice(self, size=(1, 1, 1), spacing=(2, 2, 2), type="sc", 114 distortion=0.0, rotation=0.0): 115 self.lattice_size = np.asarray(size, 'i') 116 self.lattice_spacing = np.asarray(spacing, 'd') 117 self.lattice_type = type 118 self.lattice_distortion = distortion 119 self.lattice_rotation = rotation 120 105 121 def _adjust(self, points): 106 points = np.asarray(self.rotation * np.matrix(points.T)) + self.center 122 if self.rotation is I3: 123 points = points.T + self.center 124 else: 125 points = np.asarray(self.rotation * np.matrix(points.T)) + self.center 126 if self.lattice_type: 127 points = self._apply_lattice(points) 107 128 return points.T 108 129 109 def r_bins(self, q, over_sampling=1, r_step=0.): 110 r_max = min(2 * pi / q[0], self.r_max) 130 def r_bins(self, q, over_sampling=10, r_step=0.): 131 if self.lattice_type: 132 r_max = np.sqrt(np.sum(self.lattice_size*self.lattice_spacing*self.dims)**2)/2 133 else: 134 r_max = self.r_max 135 #r_max = min(2 * pi / q[0], r_max) 111 136 if r_step == 0.: 112 137 r_step = 2 * pi / q[-1] / over_sampling 113 138 #r_step = 0.01 114 139 return np.arange(r_step, r_max, r_step) 140 141 def _apply_lattice(self, points): 142 """Spread points to different lattice positions""" 143 size = self.lattice_size 144 spacing = self.lattice_spacing 145 shuffle = self.lattice_distortion 146 rotate = self.lattice_rotation 147 lattice = self.lattice_type 148 149 if rotate != 0: 150 # To vectorize the rotations we will need to unwrap the matrix multiply 151 raise NotImplementedError("don't handle rotations yet") 152 153 # Determine the number of lattice points in the lattice 154 shapes_per_cell = 2 if lattice == "bcc" else 4 if lattice == "fcc" else 1 155 number_of_lattice_points = np.prod(size) * shapes_per_cell 156 157 # For each point in the original shape, figure out which lattice point 158 # to translate it to. This is both cell index (i*ny*nz + j*nz + k) as 159 # well as the point in the cell (corner, body center or face center). 160 nsamples = points.shape[1] 161 lattice_point = randint(number_of_lattice_points, size=nsamples) 162 163 # Translate the cell index into the i,j,k coordinates of the senter 164 cell_index = lattice_point // shapes_per_cell 165 center = np.vstack((cell_index//(size[1]*size[2]), 166 (cell_index%(size[1]*size[2]))//size[2], 167 cell_index%size[2])) 168 center = np.asarray(center, dtype='d') 169 if lattice == "bcc": 170 center[:, lattice_point % shapes_per_cell == 1] += [[0.5], [0.5], [0.5]] 171 elif lattice == "fcc": 172 center[:, lattice_point % shapes_per_cell == 1] += [[0.0], [0.5], [0.5]] 173 center[:, lattice_point % shapes_per_cell == 2] += [[0.5], [0.0], [0.5]] 174 center[:, lattice_point % shapes_per_cell == 3] += [[0.5], [0.5], [0.0]] 175 176 # Each lattice point has its own displacement from the ideal position. 177 # Not checking that shapes do not overlap if displacement is too large. 178 offset = shuffle*(randn(3, number_of_lattice_points) if shuffle < 0.3 179 else rand(3, number_of_lattice_points)) 180 center += offset[:, cell_index] 181 182 # Each lattice point has its own rotation. Rotate the point prior to 183 # applying any displacement. 184 # rotation = rotate*(randn(size=(shapes, 3)) if shuffle < 30 else rand(size=(nsamples, 3))) 185 # for k in shapes: points[k] = rotation[k]*points[k] 186 points += center*(np.array([spacing])*np.array(self.dims)).T 187 return points 115 188 116 189 class Composite(Shape): … … 669 742 Iq = 100 * np.ones_like(qx) 670 743 data = Data2D(x=qx, y=qy, z=Iq, dx=None, dy=None, dz=np.sqrt(Iq)) 671 data.x_bins = qx[0, :]672 data.y_bins = qy[:, 0]744 data.x_bins = qx[0, :] 745 data.y_bins = qy[:, 0] 673 746 data.filename = "fake data" 674 747 … … 695 768 return shape, fn, fn_xy 696 769 697 def build_sphere(radius=125, rho=2): 770 DEFAULT_SPHERE_RADIUS = 125 771 DEFAULT_SPHERE_CONTRAST = 2 772 def build_sphere(radius=DEFAULT_SPHERE_RADIUS, rho=DEFAULT_SPHERE_CONTRAST): 698 773 shape = TriaxialEllipsoid(radius, radius, radius, rho) 699 774 fn = lambda q: sphere_Iq(q, radius)*rho**2 … … 751 826 return shape, fn, fn_xy 752 827 753 def build_ cubic_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2,754 shuffle=0, rotate=0):828 def build_sc_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2, 829 shuffle=0, rotate=0): 755 830 a, b, c = shape.dims 756 shapes= [copy(shape)831 corners= [copy(shape) 757 832 .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 758 833 (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, … … 762 837 for iy in range(ny) 763 838 for iz in range(nz)] 764 lattice = Composite( shapes)839 lattice = Composite(corners) 765 840 return lattice 766 841 842 def build_bcc_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2, 843 shuffle=0, rotate=0): 844 a, b, c = shape.dims 845 corners = [copy(shape) 846 .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 847 (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 848 (iz+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 849 .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 850 for ix in range(nx) 851 for iy in range(ny) 852 for iz in range(nz)] 853 centers = [copy(shape) 854 .shift((ix+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 855 (iy+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 856 (iz+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 857 .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 858 for ix in range(nx) 859 for iy in range(ny) 860 for iz in range(nz)] 861 lattice = Composite(corners + centers) 862 return lattice 863 864 def build_fcc_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2, 865 shuffle=0, rotate=0): 866 a, b, c = shape.dims 867 corners = [copy(shape) 868 .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 869 (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 870 (iz+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 871 .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 872 for ix in range(nx) 873 for iy in range(ny) 874 for iz in range(nz)] 875 faces_a = [copy(shape) 876 .shift((ix+0.0+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 877 (iy+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 878 (iz+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 879 .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 880 for ix in range(nx) 881 for iy in range(ny) 882 for iz in range(nz)] 883 faces_b = [copy(shape) 884 .shift((ix+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 885 (iy+0.0+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 886 (iz+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 887 .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 888 for ix in range(nx) 889 for iy in range(ny) 890 for iz in range(nz)] 891 faces_c = [copy(shape) 892 .shift((ix+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 893 (iy+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 894 (iz+0.0+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 895 .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 896 for ix in range(nx) 897 for iy in range(ny) 898 for iz in range(nz)] 899 lattice = Composite(corners + faces_a + faces_b + faces_c) 900 return lattice 767 901 768 902 SHAPE_FUNCTIONS = OrderedDict([ … … 775 909 ]) 776 910 SHAPES = list(SHAPE_FUNCTIONS.keys()) 911 LATTICE_FUNCTIONS = OrderedDict([ 912 ("sc", build_sc_lattice), 913 ("bcc", build_bcc_lattice), 914 ("fcc", build_fcc_lattice), 915 ]) 916 LATTICE_TYPES = list(LATTICE_FUNCTIONS.keys()) 777 917 778 918 def check_shape(title, shape, fn=None, show_points=False, … … 783 923 r = shape.r_bins(q, r_step=r_step) 784 924 sampling_density = samples / shape.volume 925 print("sampling points") 785 926 rho, points = shape.sample(sampling_density) 927 print("calculating Pr") 786 928 t0 = time.time() 787 929 Pr = calc_Pr(r, rho-rho_solvent, points) … … 792 934 import pylab 793 935 if show_points: 794 936 plot_points(rho, points); pylab.figure() 795 937 plot_calc(r, Pr, q, Iq, theory=theory, title=title) 796 938 pylab.gcf().canvas.set_window_title(title) … … 806 948 Qx, Qy = np.meshgrid(qx, qy) 807 949 sampling_density = samples / shape.volume 950 print("sampling points") 808 951 t0 = time.time() 809 952 rho, points = shape.sample(sampling_density) … … 844 987 help='lattice size') 845 988 parser.add_argument('-z', '--spacing', type=str, default='2,2,2', 846 help='lattice spacing') 989 help='lattice spacing (relative to shape)') 990 parser.add_argument('-t', '--type', choices=LATTICE_TYPES, 991 default=LATTICE_TYPES[0], 992 help='lattice type') 847 993 parser.add_argument('-r', '--rotate', type=float, default=0., 848 994 help="rotation relative to lattice, gaussian < 30 degrees, uniform otherwise") … … 858 1004 nx, ny, nz = [int(v) for v in opts.lattice.split(',')] 859 1005 dx, dy, dz = [float(v) for v in opts.spacing.split(',')] 860 shuffle, rotate= opts.shuffle, opts.rotate1006 distortion, rotation = opts.shuffle, opts.rotate 861 1007 shape, fn, fn_xy = SHAPE_FUNCTIONS[opts.shape](**pars) 862 if nx > 1 or ny > 1 or nz > 1: 863 shape = build_cubic_lattice(shape, nx, ny, nz, dx, dy, dz, shuffle, rotate) 1008 view = tuple(float(v) for v in opts.view.split(',')) 1009 # If comparing a sphere in a cubic lattice, compare against the 1010 # corresponding paracrystalline model. 1011 if opts.shape == "sphere" and dx == dy == dz and nx*ny*nz > 1: 1012 radius = pars.get('radius', DEFAULT_SPHERE_RADIUS) 1013 model_name = opts.type + "_paracrystal" 1014 model_pars = { 1015 "scale": 1., 1016 "background": 0., 1017 "lattice_spacing": 2*radius*dx, 1018 "lattice_distortion": distortion, 1019 "radius": radius, 1020 "sld": pars.get('rho', DEFAULT_SPHERE_CONTRAST), 1021 "sld_solvent": 0., 1022 "theta": view[0], 1023 "phi": view[1], 1024 "psi": view[2], 1025 } 1026 fn, fn_xy = wrap_sasmodel(model_name, **model_pars) 1027 if nx*ny*nz > 1: 1028 if rotation != 0: 1029 print("building %s lattice"%opts.type) 1030 build_lattice = LATTICE_FUNCTIONS[opts.type] 1031 shape = build_lattice(shape, nx, ny, nz, dx, dy, dz, 1032 distortion, rotation) 1033 else: 1034 shape.lattice(size=(nx, ny, nz), spacing=(dx, dy, dz), 1035 type=opts.type, 1036 rotation=rotation, distortion=distortion) 1037 864 1038 title = "%s(%s)" % (opts.shape, " ".join(opts.pars)) 865 1039 if opts.dim == 1: … … 867 1041 mesh=opts.mesh, qmax=opts.qmax, samples=opts.samples) 868 1042 else: 869 view = tuple(float(v) for v in opts.view.split(','))870 1043 check_shape_2d(title, shape, fn_xy, view=view, show_points=opts.plot, 871 1044 mesh=opts.mesh, qmax=opts.qmax, samples=opts.samples) -
sasmodels/convert.py
r610ef23 rb3f4831 166 166 oldpars = _hand_convert_3_1_2_to_4_1(name, oldpars) 167 167 if version < (4, 2, 0): 168 oldpars = _hand_convert_4_1_to_4_2(name, oldpars) 168 169 oldpars = _rename_magnetic_pars(oldpars) 170 return oldpars 171 172 def _hand_convert_4_1_to_4_2(name, oldpars): 173 if name in ('bcc_paracrystal', 'fcc_paracrystal', 'sc_paracrystal'): 174 oldpars['lattice_spacing'] = oldpars.pop('dnn') 175 oldpars['lattice_distortion'] = oldpars.pop('d_factor') 169 176 return oldpars 170 177 -
sasmodels/models/bcc_paracrystal.c
r108e70e r642046e 1 1 static double 2 bcc_Zq(double qa, double qb, double qc, double dnn, double d_factor)2 bcc_Zq(double qa, double qb, double qc, double lattice_spacing, double lattice_distortion) 3 3 { 4 4 // Equations from Matsuoka 26-27-28, multiplied by |q| … … 17 17 // => exp(a)^2 - 2 cos(d ak) exp(a) + 1) 18 18 // => (exp(a) - 2 cos(d ak)) * exp(a) + 1 19 const double arg = -0.5*square( dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);19 const double arg = -0.5*square(lattice_spacing*lattice_distortion)*(a1*a1 + a2*a2 + a3*a3); 20 20 const double exp_arg = exp(arg); 21 21 const double Zq = -cube(expm1(2.0*arg)) 22 / ( ((exp_arg - 2.0*cos( dnn*a1))*exp_arg + 1.0)23 * ((exp_arg - 2.0*cos( dnn*a2))*exp_arg + 1.0)24 * ((exp_arg - 2.0*cos( dnn*a3))*exp_arg + 1.0));22 / ( ((exp_arg - 2.0*cos(lattice_spacing*a1))*exp_arg + 1.0) 23 * ((exp_arg - 2.0*cos(lattice_spacing*a2))*exp_arg + 1.0) 24 * ((exp_arg - 2.0*cos(lattice_spacing*a3))*exp_arg + 1.0)); 25 25 26 26 #elif 0 … … 36 36 // = tanh(-a) / [1 - cos(d a_k)/cosh(-a)] 37 37 // 38 const double arg = 0.5*square( dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);38 const double arg = 0.5*square(lattice_spacing*lattice_distortion)*(a1*a1 + a2*a2 + a3*a3); 39 39 const double sinh_qd = sinh(arg); 40 40 const double cosh_qd = cosh(arg); 41 const double Zq = sinh_qd/(cosh_qd - cos( dnn*a1))42 * sinh_qd/(cosh_qd - cos( dnn*a2))43 * sinh_qd/(cosh_qd - cos( dnn*a3));41 const double Zq = sinh_qd/(cosh_qd - cos(lattice_spacing*a1)) 42 * sinh_qd/(cosh_qd - cos(lattice_spacing*a2)) 43 * sinh_qd/(cosh_qd - cos(lattice_spacing*a3)); 44 44 #else 45 const double arg = 0.5*square( dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);45 const double arg = 0.5*square(lattice_spacing*lattice_distortion)*(a1*a1 + a2*a2 + a3*a3); 46 46 const double tanh_qd = tanh(arg); 47 47 const double cosh_qd = cosh(arg); 48 const double Zq = tanh_qd/(1.0 - cos( dnn*a1)/cosh_qd)49 * tanh_qd/(1.0 - cos( dnn*a2)/cosh_qd)50 * tanh_qd/(1.0 - cos( dnn*a3)/cosh_qd);48 const double Zq = tanh_qd/(1.0 - cos(lattice_spacing*a1)/cosh_qd) 49 * tanh_qd/(1.0 - cos(lattice_spacing*a2)/cosh_qd) 50 * tanh_qd/(1.0 - cos(lattice_spacing*a3)/cosh_qd); 51 51 #endif 52 52 … … 57 57 // occupied volume fraction calculated from lattice symmetry and sphere radius 58 58 static double 59 bcc_volume_fraction(double radius, double dnn)59 bcc_volume_fraction(double radius, double lattice_spacing) 60 60 { 61 return 2.0*sphere_volume( sqrt(0.75)*radius/dnn);61 return 2.0*sphere_volume(radius/lattice_spacing); 62 62 } 63 63 … … 69 69 70 70 71 static double Iq(double q, double dnn,72 double d_factor, double radius,71 static double Iq(double q, double lattice_spacing, 72 double lattice_distortion, double radius, 73 73 double sld, double solvent_sld) 74 74 { … … 94 94 const double qa = qab*cos_phi; 95 95 const double qb = qab*sin_phi; 96 const double form = bcc_Zq(qa, qb, qc, dnn, d_factor);96 const double form = bcc_Zq(qa, qb, qc, lattice_spacing, lattice_distortion); 97 97 inner_sum += GAUSS_W[j] * form; 98 98 } … … 103 103 const double Zq = outer_sum/(4.0*M_PI); 104 104 const double Pq = sphere_form(q, radius, sld, solvent_sld); 105 return bcc_volume_fraction(radius, dnn) * Pq * Zq;105 return bcc_volume_fraction(radius, lattice_spacing) * Pq * Zq; 106 106 } 107 107 108 108 109 109 static double Iqabc(double qa, double qb, double qc, 110 double dnn, double d_factor, double radius,110 double lattice_spacing, double lattice_distortion, double radius, 111 111 double sld, double solvent_sld) 112 112 { 113 113 const double q = sqrt(qa*qa + qb*qb + qc*qc); 114 const double Zq = bcc_Zq(qa, qb, qc, dnn, d_factor);114 const double Zq = bcc_Zq(qa, qb, qc, lattice_spacing, lattice_distortion); 115 115 const double Pq = sphere_form(q, radius, sld, solvent_sld); 116 return bcc_volume_fraction(radius, dnn) * Pq * Zq;116 return bcc_volume_fraction(radius, lattice_spacing) * Pq * Zq; 117 117 } -
sasmodels/models/bcc_paracrystal.py
rda7b26b rb3f4831 34 34 .. math:: 35 35 36 V_\text{lattice} = \frac{ 16\pi}{3} \frac{R^3}{\left(D\sqrt{2}\right)^3}36 V_\text{lattice} = \frac{8\pi}{3} \frac{R^3}{\left(2D/\sqrt{3}\right)^3} 37 37 38 38 … … 104 104 105 105 * **Author:** NIST IGOR/DANSE **Date:** pre 2010 106 * **Last Modified by:** Paul Butler **Date:** September 29, 2016107 * **Last Reviewed by:** Richard Heenan **Date:** March 21, 2016106 * **Last Modified by:** Paul Butler **Date:** September 16, 2018 107 * **Last Reviewed by:** Paul Butler **Date:** September 16, 2018 108 108 """ 109 109 … … 127 127 # pylint: disable=bad-whitespace, line-too-long 128 128 # ["name", "units", default, [lower, upper], "type","description" ], 129 parameters = [[" dnn", "Ang", 220, [-inf, inf], "", "Nearest neighbour distance"],130 [" d_factor", "", 0.06, [-inf, inf], "", "Paracrystal distortion factor"],129 parameters = [["lattice_spacing", "Ang", 220, [-inf, inf], "", "Lattice spacing"], 130 ["lattice_distortion", "", 0.06, [-inf, inf], "", "Paracrystal distortion factor"], 131 131 ["radius", "Ang", 40, [0, inf], "volume", "Particle radius"], 132 132 ["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld", "Particle scattering length density"], … … 149 149 # in this range 'cuz its easy. 150 150 radius = 10**np.random.uniform(1.3, 4) 151 d_factor= 10**np.random.uniform(-2, -0.7) # sigma_d in 0.01-0.7152 dnn_fraction = np.random.beta(a=10, b=1)153 dnn = radius*4/np.sqrt(3)/dnn_fraction151 lattice_distortion = 10**np.random.uniform(-2, -0.7) # sigma_d in 0.01-0.7 152 lattice_spacing_fraction = np.random.beta(a=10, b=1) 153 lattice_spacing = radius*4/np.sqrt(3)/lattice_spacing_fraction 154 154 pars = dict( 155 155 #sld=1, sld_solvent=0, scale=1, background=1e-32, 156 dnn=dnn,157 d_factor=d_factor,156 lattice_spacing=lattice_spacing, 157 lattice_distortion=lattice_distortion, 158 158 radius=radius, 159 159 ) -
sasmodels/models/fcc_paracrystal.c
r71b751d rcc8b183 1 1 static double 2 fcc_Zq(double qa, double qb, double qc, double dnn, double d_factor)2 fcc_Zq(double qa, double qb, double qc, double lattice_spacing, double lattice_distortion) 3 3 { 4 4 // Equations from Matsuoka 17-18-19, multiplied by |q| … … 16 16 // => exp(a)^2 - 2 cos(d ak) exp(a) + 1) 17 17 // => (exp(a) - 2 cos(d ak)) * exp(a) + 1 18 const double arg = -0.5*square( dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);18 const double arg = -0.5*square(lattice_spacing*lattice_distortion)*(a1*a1 + a2*a2 + a3*a3); 19 19 const double exp_arg = exp(arg); 20 20 const double Zq = -cube(expm1(2.0*arg)) 21 / ( ((exp_arg - 2.0*cos( dnn*a1))*exp_arg + 1.0)22 * ((exp_arg - 2.0*cos( dnn*a2))*exp_arg + 1.0)23 * ((exp_arg - 2.0*cos( dnn*a3))*exp_arg + 1.0));21 / ( ((exp_arg - 2.0*cos(lattice_spacing*a1))*exp_arg + 1.0) 22 * ((exp_arg - 2.0*cos(lattice_spacing*a2))*exp_arg + 1.0) 23 * ((exp_arg - 2.0*cos(lattice_spacing*a3))*exp_arg + 1.0)); 24 24 25 25 return Zq; … … 29 29 // occupied volume fraction calculated from lattice symmetry and sphere radius 30 30 static double 31 fcc_volume_fraction(double radius, double dnn)31 fcc_volume_fraction(double radius, double lattice_spacing) 32 32 { 33 return 4.0*sphere_volume( M_SQRT1_2*radius/dnn);33 return 4.0*sphere_volume(radius/lattice_spacing); 34 34 } 35 35 … … 41 41 42 42 43 static double Iq(double q, double dnn,44 double d_factor, double radius,43 static double Iq(double q, double lattice_spacing, 44 double lattice_distortion, double radius, 45 45 double sld, double solvent_sld) 46 46 { … … 66 66 const double qa = qab*cos_phi; 67 67 const double qb = qab*sin_phi; 68 const double form = fcc_Zq(qa, qb, qc, dnn, d_factor);68 const double form = fcc_Zq(qa, qb, qc, lattice_spacing, lattice_distortion); 69 69 inner_sum += GAUSS_W[j] * form; 70 70 } … … 76 76 const double Pq = sphere_form(q, radius, sld, solvent_sld); 77 77 78 return fcc_volume_fraction(radius, dnn) * Pq * Zq;78 return fcc_volume_fraction(radius, lattice_spacing) * Pq * Zq; 79 79 } 80 80 81 81 static double Iqabc(double qa, double qb, double qc, 82 double dnn, double d_factor, double radius,82 double lattice_spacing, double lattice_distortion, double radius, 83 83 double sld, double solvent_sld) 84 84 { 85 85 const double q = sqrt(qa*qa + qb*qb + qc*qc); 86 86 const double Pq = sphere_form(q, radius, sld, solvent_sld); 87 const double Zq = fcc_Zq(qa, qb, qc, dnn, d_factor);88 return fcc_volume_fraction(radius, dnn) * Pq * Zq;87 const double Zq = fcc_Zq(qa, qb, qc, lattice_spacing, lattice_distortion); 88 return fcc_volume_fraction(radius, lattice_spacing) * Pq * Zq; 89 89 } -
sasmodels/models/fcc_paracrystal.py
rda7b26b rb3f4831 3 3 #note - calculation requires double precision 4 4 r""" 5 .. warning:: This model and this model description are under review following 6 concerns raised by SasView users. If you need to use this model, 7 please email help@sasview.org for the latest situation. *The 5 .. warning:: This model and this model description are under review following 6 concerns raised by SasView users. If you need to use this model, 7 please email help@sasview.org for the latest situation. *The 8 8 SasView Developers. September 2018.* 9 9 … … 100 100 101 101 Authorship and Verification 102 --------------------------- 102 ---------------------------- 103 103 104 104 * **Author:** NIST IGOR/DANSE **Date:** pre 2010 105 * **Last Modified by:** Paul Butler **Date:** September 29, 2016106 * **Last Reviewed by:** Richard Heenan **Date:** March 21, 2016105 * **Last Modified by:** Paul Butler **Date:** September 16, 2018 106 * **Last Reviewed by:** Paul Butler **Date:** September 16, 2018 107 107 """ 108 108 … … 123 123 # pylint: disable=bad-whitespace, line-too-long 124 124 # ["name", "units", default, [lower, upper], "type","description"], 125 parameters = [[" dnn", "Ang", 220, [-inf, inf], "", "Nearest neighbour distance"],126 [" d_factor", "", 0.06, [-inf, inf], "", "Paracrystal distortion factor"],125 parameters = [["lattice_spacing", "Ang", 220, [-inf, inf], "", "Lattice spacing"], 126 ["lattice_distortion", "", 0.06, [-inf, inf], "", "Paracrystal distortion factor"], 127 127 ["radius", "Ang", 40, [0, inf], "volume", "Particle radius"], 128 128 ["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld", "Particle scattering length density"], … … 139 139 # copied from bcc_paracrystal 140 140 radius = 10**np.random.uniform(1.3, 4) 141 d_factor= 10**np.random.uniform(-2, -0.7) # sigma_d in 0.01-0.7142 dnn_fraction = np.random.beta(a=10, b=1)143 dnn = radius*4/np.sqrt(2)/dnn_fraction141 lattice_distortion = 10**np.random.uniform(-2, -0.7) # sigma_d in 0.01-0.7 142 lattice_spacing_fraction = np.random.beta(a=10, b=1) 143 lattice_spacing = radius*4/np.sqrt(2)/lattice_spacing_fraction 144 144 pars = dict( 145 145 #sld=1, sld_solvent=0, scale=1, background=1e-32, 146 dnn=dnn,147 d_factor=d_factor,146 latice_spacing=lattice_spacing, 147 lattice_distortion=d_factor, 148 148 radius=radius, 149 149 ) -
sasmodels/models/sc_paracrystal.c
r71b751d rcc8b183 1 1 static double 2 sc_Zq(double qa, double qb, double qc, double dnn, double d_factor)2 sc_Zq(double qa, double qb, double qc, double lattice_spacing, double lattice_distortion) 3 3 { 4 4 // Equations from Matsuoka 9-10-11, multiplied by |q| … … 16 16 // => exp(a)^2 - 2 cos(d ak) exp(a) + 1) 17 17 // => (exp(a) - 2 cos(d ak)) * exp(a) + 1 18 const double arg = -0.5*square( dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);18 const double arg = -0.5*square(lattice_spacing*lattice_distortion)*(a1*a1 + a2*a2 + a3*a3); 19 19 const double exp_arg = exp(arg); 20 20 const double Zq = -cube(expm1(2.0*arg)) 21 / ( ((exp_arg - 2.0*cos( dnn*a1))*exp_arg + 1.0)22 * ((exp_arg - 2.0*cos( dnn*a2))*exp_arg + 1.0)23 * ((exp_arg - 2.0*cos( dnn*a3))*exp_arg + 1.0));21 / ( ((exp_arg - 2.0*cos(lattice_spacing*a1))*exp_arg + 1.0) 22 * ((exp_arg - 2.0*cos(lattice_spacing*a2))*exp_arg + 1.0) 23 * ((exp_arg - 2.0*cos(lattice_spacing*a3))*exp_arg + 1.0)); 24 24 25 25 return Zq; … … 28 28 // occupied volume fraction calculated from lattice symmetry and sphere radius 29 29 static double 30 sc_volume_fraction(double radius, double dnn)30 sc_volume_fraction(double radius, double lattice_spacing) 31 31 { 32 return sphere_volume(radius/ dnn);32 return sphere_volume(radius/lattice_spacing); 33 33 } 34 34 … … 41 41 42 42 static double 43 Iq(double q, double dnn,44 double d_factor, double radius,43 Iq(double q, double lattice_spacing, 44 double lattice_distortion, double radius, 45 45 double sld, double solvent_sld) 46 46 { … … 67 67 const double qa = qab*cos_phi; 68 68 const double qb = qab*sin_phi; 69 const double form = sc_Zq(qa, qb, qc, dnn, d_factor);69 const double form = sc_Zq(qa, qb, qc, lattice_spacing, lattice_distortion); 70 70 inner_sum += GAUSS_W[j] * form; 71 71 } … … 77 77 const double Pq = sphere_form(q, radius, sld, solvent_sld); 78 78 79 return sc_volume_fraction(radius, dnn) * Pq * Zq;79 return sc_volume_fraction(radius, lattice_spacing) * Pq * Zq; 80 80 } 81 81 82 82 static double 83 83 Iqabc(double qa, double qb, double qc, 84 double dnn, double d_factor, double radius,84 double lattice_spacing, double lattice_distortion, double radius, 85 85 double sld, double solvent_sld) 86 86 { 87 87 const double q = sqrt(qa*qa + qb*qb + qc*qc); 88 88 const double Pq = sphere_form(q, radius, sld, solvent_sld); 89 const double Zq = sc_Zq(qa, qb, qc, dnn, d_factor);90 return sc_volume_fraction(radius, dnn) * Pq * Zq;89 const double Zq = sc_Zq(qa, qb, qc, lattice_spacing, lattice_distortion); 90 return sc_volume_fraction(radius, lattice_spacing) * Pq * Zq; 91 91 } -
sasmodels/models/sc_paracrystal.py
rda7b26b rb3f4831 1 1 r""" 2 .. warning:: This model and this model description are under review following 3 concerns raised by SasView users. If you need to use this model, 4 please email help@sasview.org for the latest situation. *The 2 .. warning:: This model and this model description are under review following 3 concerns raised by SasView users. If you need to use this model, 4 please email help@sasview.org for the latest situation. *The 5 5 SasView Developers. September 2018.* 6 6 7 7 Definition 8 8 ---------- … … 104 104 105 105 Authorship and Verification 106 --------------------------- 106 ---------------------------- 107 107 108 108 * **Author:** NIST IGOR/DANSE **Date:** pre 2010 109 * **Last Modified by:** Paul Butler **Date:** September 29, 2016110 * **Last Reviewed by:** Richard Heenan **Date:** March 21, 2016109 * **Last Modified by:** Paul Butler **Date:** September 16, 2018 110 * **Last Reviewed by:** Paul Butler **Date:** September 16, 2018 111 111 """ 112 112 … … 129 129 scale: volume fraction of spheres 130 130 bkg:background, R: radius of sphere 131 dnn: Nearest neighbor distance132 d_factor: Paracrystal distortion factor131 lattice_spacing: Nearest neighbor distance 132 lattice_distortion: Paracrystal distortion factor 133 133 radius: radius of the spheres 134 134 sldSph: SLD of the sphere … … 139 139 # pylint: disable=bad-whitespace, line-too-long 140 140 # ["name", "units", default, [lower, upper], "type","description"], 141 parameters = [[" dnn", "Ang", 220.0, [0.0, inf], "", "Nearest neighbor distance"],142 [" d_factor", "", 0.06, [-inf, inf], "","Paracrystal distortion factor"],141 parameters = [["lattice_spacing", "Ang", 220.0, [0.0, inf], "", "Lattice spacing"], 142 ["lattice_distortion", "", 0.06, [-inf, inf], "", "Paracrystal distortion factor"], 143 143 ["radius", "Ang", 40.0, [0.0, inf], "volume", "Radius of sphere"], 144 144 ["sld", "1e-6/Ang^2", 3.0, [0.0, inf], "sld", "Sphere scattering length density"], … … 155 155 # copied from bcc_paracrystal 156 156 radius = 10**np.random.uniform(1.3, 4) 157 d_factor= 10**np.random.uniform(-2, -0.7) # sigma_d in 0.01-0.7158 dnn_fraction = np.random.beta(a=10, b=1)159 dnn = radius*4/np.sqrt(4)/dnn_fraction157 lattice_distortion = 10**np.random.uniform(-2, -0.7) # sigma_d in 0.01-0.7 158 lattice_spacing_fraction = np.random.beta(a=10, b=1) 159 lattice_spacing = radius*4/np.sqrt(4)/lattice_spacing_fraction 160 160 pars = dict( 161 161 #sld=1, sld_solvent=0, scale=1, background=1e-32, 162 dnn=dnn,163 d_factor=d_factor,162 lattice_spacing=lattice_spacing, 163 lattice_distortion=lattice_distortion, 164 164 radius=radius, 165 165 )
Note: See TracChangeset
for help on using the changeset viewer.