Changeset 0c86b79 in sasmodels

Ignore:
Timestamp:
Mar 31, 2019 7:59:48 AM (7 months ago)
Parents:
1dd78a2b (diff), 225bf94 (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.
git-author:
Paul Kienzle <pkienzle@…> (03/31/19 07:59:48)
git-committer:
Message:
Files:
144 edited

Unmodified
Removed
• .gitignore

 re9ed2de *.so *.obj *.o /doc/_build/ /doc/api/ /.pydevproject /.idea .vscode /sasmodels.egg-info/ /example/Fit_*/
• sasmodels/modelinfo.py

 ra34b811 from os.path import abspath, basename, splitext import inspect import logging import numpy as np  # type: ignore from . import autoc # Optional typing TestCondition = Tuple[ParameterSetUser, TestInput, TestValue] # pylint: enable=unused-import logger = logging.getLogger(__name__) # If MAX_PD changes, need to change the loop macros in kernel_iq.c info.profile = getattr(kernel_module, 'profile', None) # type: ignore info.sesans = getattr(kernel_module, 'sesans', None) # type: ignore info.random = getattr(kernel_module, 'random', None) info.hidden = getattr(kernel_module, 'hidden', None) # type: ignore info.lineno = {} _find_source_lines(info, kernel_module) if getattr(kernel_module, 'py2c', False): try: warnings = autoc.convert(info, kernel_module) except Exception as exc: warnings = [str(exc)] if warnings: warnings.append("while converting %s from C to python"%name) if len(warnings) > 2: warnings = "\n".join(warnings) else: warnings = " ".join(warnings) logger.warn(warnings) # Default single and opencl to True for C models.  Python models have callable Iq. # Needs to come after autoc.convert since the Iq symbol may have been # converted from python to C info.opencl = getattr(kernel_module, 'opencl', not callable(info.Iq)) info.single = getattr(kernel_module, 'single', not callable(info.Iq)) info.random = getattr(kernel_module, 'random', None) info.hidden = getattr(kernel_module, 'hidden', None) # type: ignore # Set control flag for explicitly set parameters, e.g., in the RPA model. control = getattr(kernel_module, 'control', None) raise ValueError("oriented python models not supported") info.lineno = {} _find_source_lines(info, kernel_module) return info
• doc/guide/index.rst

 rda5536f pd/polydispersity.rst resolution.rst plugin.rst fitting_sq.rst magnetism/magnetism.rst orientation/orientation.rst sesans/sans_to_sesans.rst sesans/sesans_fitting.rst plugin.rst scripting.rst refs.rst
• doc/guide/plugin.rst

 r9150036 **Note: The order of the parameters in the definition will be the order of the parameters in the user interface and the order of the parameters in Fq(), Iq(), Iqac(), Iqabc(), form_volume() and shell_volume(). Iqac(), Iqabc(), radius_effective(), form_volume() and shell_volume(). And** *scale* **and** *background* **parameters are implicit to all models, so they do not need to be included in the parameter table.** can take arbitrary values, even for integer parameters, so your model should round the incoming parameter value to the nearest integer inside your model you should round to the nearest integer.  In C code, you can do this using:: you should round to the nearest integer.  In C code, you can do this using: .. code-block:: c static double ............. .. note:: Pure python models do not yet support direct computation of  or $^2$. Neither do they support orientational distributions or magnetism (use C models if these are required). For pure python models, define the *Iq* function:: Models should define *form_volume(par1, par2, ...)* where the parameter list includes the *volume* parameters in order.  This is used for a weighted volume normalization so that scattering is on an absolute scale.  If *form_volume* is not defined, then the default *form_volume = 1.0* will be used. volume normalization so that scattering is on an absolute scale.  For solid shapes, the *I(q)* function should use *form_volume* squared as its scale factor.  If *form_volume* is not defined, then the default *form_volume = 1.0* will be used. Hollow shapes, where the volume fraction of particle corresponds to the material in the shell rather than the volume enclosed by the shape, must also define a *shell_volume(par1, par2, ...)* function.  The parameters are the same as for *form_volume*.  The *I(q)* calculation should use *shell_volume* squared as its scale factor for the volume normalization. The structure factor calculation needs *form_volume* in order to properly scale the volume fraction parameter, so both functions are required for hollow shapes. Note: Pure python models do not yet support direct computation of the average of $F(q)$ and $F^2(q)$. are the same as for *form_volume*.  Here the *I(q)* function should use *shell_volume* squared instead of *form_volume* squared so that the scale parameter corresponds to the volume fraction of material in the sample. The structure factor calculation needs the volume fraction of the filled shapes for its calculation, so the volume fraction parameter in the model is automatically scaled by *form_volume/shell_volume* prior to calling the structure factor. Embedded C Models """ This expands into the equivalent C code:: This expands into the equivalent C code: .. code-block:: c double Iq(double q, double par1, double par2, ...); includes only the volume parameters. *form_volume* defines the volume of the shell for hollow shapes. As in *shell_volume* defines the volume of the shell for hollow shapes. As in python models, it includes only the volume parameters. Rather than returning NAN from Iq, you must define the *INVALID(v)*.  The *v* parameter lets you access all the parameters in the model using *v.par1*, *v.par2*, etc. For example:: *v.par1*, *v.par2*, etc. For example: .. code-block:: c #define INVALID(v) (v.bell_radius < v.radius) Instead of defining the *Iq* function, models can define *Fq* as something like:: something like: .. code-block:: c double Fq(double q, double *F1, double *F2, double par1, double par2, ...); laboratory frame and beam travelling along $-z$. The oriented C model is called using *Iqabc(qa, qb, qc, par1, par2, ...)* where The oriented C model (oriented pure Python models are not supported) is called using *Iqabc(qa, qb, qc, par1, par2, ...)* where *par1*, etc. are the parameters to the model.  If the shape is rotationally symmetric about *c* then *psi* is not needed, and the model is called Using the $z, w$ values for Gauss-Legendre integration in "lib/gauss76.c", the numerical integration is then:: numerical integration is then: .. code-block:: c double outer_sum = 0.0; to compute the proper magnetism and orientation, which you can implement using *Iqxy(qx, qy, par1, par2, ...)*. **Note: Magnetism is not supported in pure Python models.** Special Functions memory, and wrong answers computed. The conclusion from a very long and strange debugging session was that any arrays that you declare in your model should be a multiple of four. For example:: model should be a multiple of four. For example: .. code-block:: c double Iq(q, p1, p2, ...) structure factor is the *hardsphere* interaction, which uses the effective radius of the form factor as an input to the structure factor model.  The effective radius is the average radius of the form averaged over all the polydispersity values. :: def ER(radius, thickness): """Effective radius of a core-shell sphere.""" return radius + thickness Now consider the *core_shell_sphere*, which has a simple effective radius equal to the radius of the core plus the thickness of the shell, as shown above. Given polydispersity over *(r1, r2, ..., rm)* in radius and *(t1, t2, ..., tn)* in thickness, *ER* is called with a mesh grid covering all possible combinations of radius and thickness. That is, *radius* is *(r1, r2, ..., rm, r1, r2, ..., rm, ...)* and *thickness* is *(t1, t1, ... t1, t2, t2, ..., t2, ...)*. The *ER* function returns one effective radius for each combination. The effective radius calculator weights each of these according to the polydispersity distributions and calls the structure factor with the average *ER*. :: def VR(radius, thickness): """Sphere and shell volumes for a core-shell sphere.""" whole = 4.0/3.0 * pi * (radius + thickness)**3 core = 4.0/3.0 * pi * radius**3 return whole, whole - core Core-shell type models have an additional volume ratio which scales the structure factor.  The *VR* function returns the volume of the whole sphere and the volume of the shell. Like *ER*, there is one return value for each point in the mesh grid. *NOTE: we may be removing or modifying this feature soon. As of the time of writing, core-shell sphere returns (1., 1.) for VR, giving a volume ratio of 1.0.* factor model.  The effective radius is the weighted average over all values of the shape in polydisperse systems. There can be many notions of effective radius, depending on the shape.  For a sphere it is clearly just the radius, but for an ellipsoid of revolution we provide average curvature, equivalent sphere radius, minimum radius and maximum radius.  These options are listed as *radius_effective_modes* in the python model defintion, and must be computed by the *radius_effective* function which takes the *radius_effective_mode* parameter as an integer, along with the various model parameters.  Unlike normal C/Python arrays, the first mode is 1, the second is 2, etc.  Mode 0 indicates that the effective radius from the shape is to be ignored in favour of the the effective radius parameter in the structure factor model. Consider the core-shell sphere, which defines the following effective radius modes in the python model:: radius_effective_modes = [ "outer radius", "core radius", ] and the following function in the C-file for the model: .. code-block:: c static double radius_effective(int mode, double radius, double thickness) { switch (mode) { case 0: return radius + thickness; case 1: return radius; default: return 0.; } } static double form_volume(double radius, double thickness) { return M_4PI_3 * cube(radius + thickness); } Given polydispersity over *(r1, r2, ..., rm)* in radius and *(t1, t2, ..., tn)* in thickness, *radius_effective* is called over a mesh grid covering all possible combinations of radius and thickness, with a single *(ri, tj)* pair in each call. The weights each of these results according to the polydispersity distributions and calls the structure factor with the average effective radius.  Similarly, for *form_volume*. Hollow models have an additional volume ratio which is needed to scale the structure factor.  The structure factor uses the volume fraction of the filled particles as part of its density estimate, but the scale factor for the scattering intensity (as non-solvent volume fraction / volume) is determined by the shell volume only.  Therefore the *shell_volume* function is needed to compute the form:shell volume ratio, which then scales the *volfraction* parameter prior to calling the structure factor calculator. In the case of a hollow sphere, this would be: .. code-block:: c static double shell_volume(double radius, double thickness) { double whole = M_4PI_3 * cube(radius + thickness); double core = M_4PI_3 * cube(radius); return whole - core; } If *shell_volume* is not present, then *form_volume* and *shell_volume* are assumed to be equal, and the shape is considered solid. Unit Tests and a check that the model runs. If you are not using sasmodels from SasView, skip this step. Recommended Testing ................... **NB: For now, this more detailed testing is only possible if you have a SasView build environment available!** If the model compiles and runs, you can next run the unit tests that | 2016-10-25 Steve King | 2017-05-07 Paul Kienzle - Moved from sasview to sasmodels docs | 2019-03-28 Paul Kienzle - Update docs for radius_effective and shell_volume
• doc/guide/resolution.rst

 r0db85af .. sm_help.rst .. This is a port of the original SasView html help file to ReSTructured text .. by S King, ISIS, during SasView CodeCamp-III in Feb 2015. .. resolution.rst .. This is a port of the original SasView html help file sm_help to ReSTructured .. text by S King, ISIS, during SasView CodeCamp-III in Feb 2015. resolution contribution into a model calculation/simulation (which by definition will be exact) to make it more representative of what has been measured experimentally - a process called *smearing*. Sasmodels does the latter. experimentally - a process called *smearing*. The Sasmodels component of SasView does the latter. Both smearing and desmearing rely on functions to describe the resolution for the instrument and stored with the data file.  If not, they will need to be set manually before fitting. .. note:: Problems may be encountered if the data set loaded by SasView is a concatenation of SANS data from several detector distances where, of course, the worst Q resolution is next to the beam stop at each detector distance. (This will also be noticeable in the residuals plot where there will be poor overlap). SasView sensibly orders all the input data points by increasing Q for nicer-looking plots, however, the dQ data can then vary considerably from point to point. If 'Use dQ data' smearing is selected then spikes may appear in the model fits, whereas if 'None' or 'Custom Pinhole Smear' are selected the fits look normal. In such instances, possible solutions are to simply remove the data with poor Q resolution from the shorter detector distances, or to fit the data from different detector distances simultaneously.

• explore/beta/sasfit_compare.py

 r119073a #radius_effective=12.59921049894873, ) target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars) target = sasmodels_theory(q, model, radius_effective_mode=0, structure_factor_mode=1, **pars) actual = ellipsoid_pe(q, norm='sasview', **pars) title = " ".join(("sasmodels", model, pd_type))

 rdb8756e # Run tests STATUS=0 python -m sasmodels.model_test opencl_and_dll all || STATUS=$? python -m sasmodels.model_test -v all || STATUS=$? python -m sasmodels.resolution || STATUS=$? • extra/pylint.rc  r6e68289 locally-enabled, locally-disabled, #no-else-return, #invalid-name, fixme, too-many-lines, too-many-locals, too-many-branches, too-many-statements, #too-many-instance-attributes, #too-many-return-statements, #no-self-use, [REPORTS] good-names=_, input, id, h, i, j, k, n, p, v, w, x, y, z, q, qx, qy, qz, Q, Qx, Qy, Qz, dt, dx, dy, dz, dq, h, i, j, k, n, p, u, v, w, x, y, z, r, dr, q, qx, qy, qz, Q, Qx, Qy, Qz, nx, ny, nz, nq, nr, dt, dx, dy, dz, dq, S, P, a, b, c, Iq, dIq, Iqxy, Iq_calc, ER, call_ER, VR, call_VR, # Regular expression matching correct variable names variable-rgx=[a-z_][a-z0-9_{2D}{1D}]{2,30}$ variable-rgx=[a-z_][a-z0-9_{2D}{1D}]{1,30}$# Naming hint for variable names # Regular expression matching correct argument names argument-rgx=[a-z_][a-z0-9_]{2,30}$ argument-rgx=[a-z_][a-z0-9_]{1,30}$# Naming hint for argument names # (useful for modules/projects where namespaces are manipulated during runtime # and thus existing member attributes cannot be deduced by static analysis ignored-modules=bumps,sas,numpy,numpy.random,scipy,scipy.special ignored-modules=bumps,sas,numpy,numpy.random,scipy,scipy.special,pyopencl,pycuda,matplotlib.cm,ipyvolume,ipywidget # List of classes names for which member attributes should not be checked • sasmodels/bumps_model.py  r2c4a190 from .kernel import KernelModel from .modelinfo import ModelInfo from .resolution import Resolution Data = Union[Data1D, Data2D] except ImportError: @property def resolution(self): # type: () -> Union[None, Resolution] """ :class:sasmodels.Resolution applied to the data, if any. """ return self._resolution @resolution.setter def resolution(self, value): # type: (Resolution) -> None """ :class:sasmodels.Resolution applied to the data, if any. """ self._resolution = value Save the model parameters and data into a file. Not Implemented. Not Implemented except for sesans fits. """ if self.data_type == "sesans": • sasmodels/compare.py  rc1799d3 set_integration_size(model_info, ngauss) if (dtype != "default" and not dtype.endswith('!') if (dtype != "default" and not dtype.endswith('!') and not (kernelcl.use_opencl() or kernelcuda.use_cuda())): raise RuntimeError("OpenCL not available " + kernelcl.OPENCL_ERROR) opts['pars'] = parse_pars(opts, maxdim=maxdim) if opts['pars'] is None: return return limits result = run_models(opts, verbose=True) if opts['plot']: args = dict((k, v) for k, v in args.items() if "_pd" not in k and ":" not in k and k not in ("background", "scale", "theta", "phi", "psi")) and ":" not in k and k not in ("background", "scale", "theta", "phi", "psi")) args = args.copy() # work with trimmed data, not the full set sorted_err = np.sort(abs(err.compressed())) if len(sorted_err) == 0: if sorted_err.size == 0: print(label + " no valid values") return • sasmodels/convert.py  r21c93c3 def revert_name(model_info): """Translate model name back to the name used in SasView 3.x""" oldname, _ = CONVERSION_TABLE.get(model_info.id, [None, {}]) return oldname def test_backward_forward(): """ Test conversion of model parameters from 4.x to 3.x and back. """ from .core import list_models L = lambda name: _check_one(name, seed=1) • sasmodels/core.py  rd92182f def test_composite_order(): """ Check that mixture models produce the same result independent of ordder. """ def test_models(fst, snd): """Confirm that two models produce the same parameters""" def build_test(first, second): """Construct pair model test""" test = lambda description: test_models(first, second) description = first + " vs. " + second # type: () -> None """Check that model load works""" from .product import RADIUS_ID, VOLFRAC_ID, STRUCTURE_MODE_ID, RADIUS_MODE_ID #Test the the model produces the parameters that we would expect model = load_model("cylinder@hardsphere*sphere") actual = [p.name for p in model.info.parameters.kernel_parameters] target = ("sld sld_solvent radius length theta phi" " radius_effective volfraction " " structure_factor_mode radius_effective_mode" " A_sld A_sld_solvent A_radius").split() target = ["sld", "sld_solvent", "radius", "length", "theta", "phi", RADIUS_ID, VOLFRAC_ID, STRUCTURE_MODE_ID, RADIUS_MODE_ID, "A_sld", "A_sld_solvent", "A_radius"] assert target == actual, "%s != %s"%(target, actual) def list_models_main(): # type: () -> None # type: () -> int """ Run list_models as a main program. See :func:list_models for the try: models = list_models(kind) except Exception as exc: print("\n".join(models)) except Exception: print(list_models.__doc__) return 1 print("\n".join(list_models(kind))) return 0 if __name__ == "__main__": • sasmodels/data.py  rc88f983 mdata = masked_array(data.y, data.mask.copy()) mdata[~np.isfinite(mdata)] = masked if view is 'log': if view == 'log': mdata[mdata <= 0] = masked plt.errorbar(data.x, scale*mdata, yerr=data.dy, fmt='.') mtheory = masked_array(theory) mtheory[~np.isfinite(mtheory)] = masked if view is 'log': if view == 'log': mtheory[mtheory <= 0] = masked plt.plot(theory_x, scale*mtheory, '-') • sasmodels/details.py  r885753a npars = kernel.info.parameters.npars nvalues = kernel.info.parameters.nvalues scalars = [value for value, _dispersity, _weight in mesh] scalars = [value for value, dispersity, weight in mesh] # skipping scale and background when building values and weights _values, dispersity, weights = zip(*mesh[2:npars+2]) if npars else ((), (), ()) #weights = correct_theta_weights(kernel.info.parameters, dispersity, weights) length = np.array([len(w) for w in weights]) value, dispersity, weight = zip(*mesh[2:npars+2]) if npars else ((), (), ()) #weight = correct_theta_weights(kernel.info.parameters, dispersity, weight) length = np.array([len(w) for w in weight]) offset = np.cumsum(np.hstack((0, length))) call_details = make_details(kernel.info, length, offset[:-1], offset[-1]) data_len = nvalues + 2*sum(len(v) for v in dispersity) extra = (32 - data_len%32)%32 data = np.hstack((scalars,) + dispersity + weights + ZEROS[:extra]) data = np.hstack((scalars,) + dispersity + weight + ZEROS[:extra]) data = data.astype(kernel.dtype) is_magnetic = convert_magnetism(kernel.info.parameters, data) • sasmodels/direct_model.py  r9150036 from . import resolution2d from .details import make_kernel_args, dispersion_mesh from .modelinfo import DEFAULT_BACKGROUND from .product import RADIUS_MODE_ID """ mesh = get_mesh(calculator.info, pars, dim=calculator.dim, mono=mono) #print("pars", list(zip(*mesh))[0]) #print("in call_kernel: pars:", list(zip(*mesh))[0]) call_details, values, is_magnetic = make_kernel_args(calculator, mesh) #print("values:", values) #print("in call_kernel: values:", values) return calculator(call_details, values, cutoff, is_magnetic) Use parameter *radius_effective_mode* to select the effective radius calculation. calculation to use amongst the *radius_effective_modes* list given in the model. """ R_eff_type = int(pars.pop(RADIUS_MODE_ID, 1.0)) mesh = get_mesh(calculator.info, pars, dim=calculator.dim, mono=mono) #print("pars", list(zip(*mesh))[0]) #print("in call_Fq: pars", list(zip(*mesh))[0]) call_details, values, is_magnetic = make_kernel_args(calculator, mesh) #print("values:", values) #print("in call_Fq: values:", values) return calculator.Fq(call_details, values, cutoff, is_magnetic, R_eff_type) active = lambda name: True #print("pars",[p.id for p in parameters.call_parameters]) #print("in get_mesh: pars:",[p.id for p in parameters.call_parameters]) mesh = [_get_par_weights(p, values, active(p.name)) for p in parameters.call_parameters] • sasmodels/generate.py  ra8a1d48 which are hollow. *effective_radius(mode, p1, p2, ...)* returns the effective radius of *radius_effective(mode, p1, p2, ...)* returns the effective radius of the form with particular dimensions. Mode determines the type of effective radius returned, with mode=1 for equivalent volume. def get_data_path(external_dir, target_file): """ Search for the target file relative in the installed application. Search first in the location of the generate module in case we are running directly from the distribution. Search next to the python executable for windows installs. Search in the ../Resources directory next to the executable for Mac OS/X installs. """ path = abspath(dirname(__file__)) if exists(joinpath(path, target_file)): # # NOTE: there is an RST_PROLOG at the end of this file which is NOT # used for the bundled documentation. Still as long as we are defining the macros # in two places any new addition should define the macro in both places. # used for the bundled documentation. Still as long as we are defining the # macros in two places any new addition should define the macro in both places. RST_UNITS = { "Ang": "|Ang|", def test_tag_float(): """check that floating point constants are properly identified and tagged with 'f'""" """Check that floating point constants are identified and tagged with 'f'""" cases = """ return False _SHELL_VOLUME_PATTERN = re.compile(r"(^|\s)double\s+shell_volume[(]", flags=re.MULTILINE) _SHELL_VOLUME_PATTERN = re.compile(r"(^|\s)double\s+shell_volume[(]", flags=re.MULTILINE) def contains_shell_volume(source): # type: (List[str]) -> bool logger.warn("oriented shapes should define Iqac or Iqabc") else: raise ValueError("Expected function Iqac or Iqabc for oriented shape") raise ValueError("Expected Iqac or Iqabc for oriented shape") # Define the parameter table for p in partable.kernel_parameters)) # Define the function calls call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) 0.0" call_radius_effective = "#define CALL_RADIUS_EFFECTIVE(_mode, _v) 0.0" if partable.form_volume_parameters: refs = _call_pars("_v.", partable.form_volume_parameters) if is_hollow: call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = form_volume(%s); _shell = shell_volume(%s); } while (0)"%((",".join(refs),)*2) call_volume = ( "#define CALL_VOLUME(_form, _shell, _v) " "do { _form = form_volume(%s); _shell = shell_volume(%s); } " "while (0)") % ((",".join(refs),)*2) else: call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = _shell = form_volume(%s); } while (0)"%(",".join(refs)) if model_info.effective_radius_type: call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) effective_radius(_mode, %s)"%(",".join(refs)) call_volume = ( "#define CALL_VOLUME(_form, _shell, _v) " "do { _form = _shell = form_volume(%s); } " "while (0)") % (",".join(refs)) if model_info.radius_effective_modes: call_radius_effective = ( "#define CALL_RADIUS_EFFECTIVE(_mode, _v) " "radius_effective(_mode, %s)") % (",".join(refs)) else: # Model doesn't have volume. We could make the kernel run a little # faster by not using/transferring the volume normalizations, but # the ifdef's reduce readability more than is worthwhile. call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = _shell = 1.0; } while (0)" call_volume = ( "#define CALL_VOLUME(_form, _shell, _v) " "do { _form = _shell = 1.0; } while (0)") source.append(call_volume) source.append(call_effective_radius) source.append(call_radius_effective) model_refs = _call_pars("_v.", partable.iq_parameters) source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars)) source.append("#define PROJECTION %d"%PROJECTION) # TODO: allow mixed python/opencl kernels? ocl = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name) dll = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name) result = { 'dll': '\n'.join(source+dll[0]+dll[1]+dll[2]), 'opencl': '\n'.join(source+ocl[0]+ocl[1]+ocl[2]), } wrappers = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name) code = '\n'.join(source + wrappers[0] + wrappers[1] + wrappers[2]) # Note: Identical code for dll and opencl. This may not always be the case # so leave them as separate entries in the returned value. result = {'dll': code, 'opencl': code} return result • sasmodels/gengauss.py  rff31782 #!/usr/bin/env python """ Generate the Gauss-Legendre integration points and save them as a C file. """ from __future__ import division, print_function import sys import numpy as np def gengauss(n, path): """ Save the Gauss-Legendre integration points for length *n* into file *path*. """ z, w = leggauss(n) • sasmodels/guyou.py  rb3703f5 def ellipticJi(u, v, m): """Returns [sn, cn, dn](u + iv|m).""" scalar = np.isscalar(u) and np.isscalar(v) and np.isscalar(m) u, v, m = np.broadcast_arrays(u, v, m) ] # calculate F(phi+ipsi|m). # see Abramowitz and Stegun, 17.4.11. def ellipticFi(phi, psi, m): """Returns F(phi+ipsi|m). See Abramowitz and Stegun, 17.4.11.""" if np.any(phi == 0): scalar = np.isscalar(phi) and np.isscalar(psi) and np.isscalar(m) plt.plot(x, y, 'g') for long in range(-limit, limit+1, step): x, y = guyou_invert(scale*long, scale*long_line) for longitude in range(-limit, limit+1, step): x, y = guyou_invert(scale*longitude, scale*long_line) plt.plot(x, y, 'b') plt.xlabel('longitude') • sasmodels/jitter.py  r7d97437 #!/usr/bin/env python # -*- coding: utf-8 -*- """ Jitter Explorer Application to explore orientation angle and angular dispersity. From the command line:: # Show docs python -m sasmodels.jitter --help # Guyou projection jitter, uniform over 20 degree theta and 10 in phi python -m sasmodels.jitter --projection=guyou --dist=uniform --jitter=20,10,0 From a jupyter cell:: import ipyvolume as ipv from sasmodels import jitter import importlib; importlib.reload(jitter) jitter.set_plotter("ipv") size = (10, 40, 100) view = (20, 0, 0) #size = (15, 15, 100) #view = (60, 60, 0) dview = (0, 0, 0) #dview = (5, 5, 0) #dview = (15, 180, 0) #dview = (180, 15, 0) projection = 'equirectangular' #projection = 'azimuthal_equidistance' #projection = 'guyou' #projection = 'sinusoidal' #projection = 'azimuthal_equal_area' dist = 'uniform' #dist = 'gaussian' jitter.run(size=size, view=view, jitter=dview, dist=dist, projection=projection) #filename = projection+('_theta' if dview[0] == 180 else '_phi' if dview[1] == 180 else '') #ipv.savefig(filename+'.png') """ from __future__ import division, print_function import argparse try: # CRUFT: travis-ci does not support mpl_toolkits.mplot3d import mpl_toolkits.mplot3d # Adds projection='3d' option to subplot except ImportError: pass import matplotlib as mpl import matplotlib.pyplot as plt from matplotlib.widgets import Slider import numpy as np from numpy import pi, cos, sin, sqrt, exp, degrees, radians def draw_beam(axes, view=(0, 0)): from numpy import pi, cos, sin, sqrt, exp, log, degrees, radians, arccos, arctan2 # Too many complaints about variable names from pylint: # a, b, c, u, v, x, y, z, dx, dy, dz, px, py, pz, R, Rx, Ry, Rz, ... # pylint: disable=invalid-name def draw_beam(axes, view=(0, 0), alpha=0.5, steps=25): """ Draw the beam going from source at (0, 0, 1) to detector at (0, 0, -1) """ #axes.plot([0,0],[0,0],[1,-1]) #axes.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8) steps = 25 u = np.linspace(0, 2 * np.pi, steps) v = np.linspace(-1, 1, steps) #axes.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=alpha) u = np.linspace(0, 2 * pi, steps) v = np.linspace(-1, 1, 2) r = 0.02 x = r*np.outer(np.cos(u), np.ones_like(v)) y = r*np.outer(np.sin(u), np.ones_like(v)) x = r*np.outer(cos(u), np.ones_like(v)) y = r*np.outer(sin(u), np.ones_like(v)) z = 1.3*np.outer(np.ones_like(u), v) points = Rz(phi)*Ry(theta)*points x, y, z = [v.reshape(shape) for v in points] axes.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5) axes.plot_surface(x, y, z, color='yellow', alpha=alpha) # TODO: draw endcaps on beam ## Drawing tiny balls on the end will work #draw_sphere(axes, radius=0.02, center=(0, 0, 1.3), color='yellow', alpha=alpha) #draw_sphere(axes, radius=0.02, center=(0, 0, -1.3), color='yellow', alpha=alpha) ## The following does not work #triangles = [(0, i+1, i+2) for i in range(steps-2)] #x_cap, y_cap = x[:, 0], y[:, 0] #for z_cap in z[:, 0], z[:, -1]: # axes.plot_trisurf(x_cap, y_cap, z_cap, triangles, # color='yellow', alpha=alpha) def draw_ellipsoid(axes, size, view, jitter, steps=25, alpha=1): """Draw an ellipsoid.""" a, b, c = size u = np.linspace(0, 2 * np.pi, steps) v = np.linspace(0, np.pi, steps) x = a*np.outer(np.cos(u), np.sin(v)) y = b*np.outer(np.sin(u), np.sin(v)) z = c*np.outer(np.ones_like(u), np.cos(v)) u = np.linspace(0, 2 * pi, steps) v = np.linspace(0, pi, steps) x = a*np.outer(cos(u), sin(v)) y = b*np.outer(sin(u), sin(v)) z = c*np.outer(np.ones_like(u), cos(v)) x, y, z = transform_xyz(view, jitter, x, y, z) axes.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha) axes.plot_surface(x, y, z, color='w', alpha=alpha) draw_labels(axes, view, jitter, [ def draw_sc(axes, size, view, jitter, steps=None, alpha=1): """Draw points for simple cubic paracrystal""" # pylint: disable=unused-argument atoms = _build_sc() _draw_crystal(axes, size, view, jitter, atoms=atoms) def draw_fcc(axes, size, view, jitter, steps=None, alpha=1): """Draw points for face-centered cubic paracrystal""" # pylint: disable=unused-argument # Build the simple cubic crystal atoms = _build_sc() def draw_bcc(axes, size, view, jitter, steps=None, alpha=1): """Draw points for body-centered cubic paracrystal""" # pylint: disable=unused-argument # Build the simple cubic crystal atoms = _build_sc() return atoms def draw_parallelepiped(axes, size, view, jitter, steps=None, alpha=1): """Draw a parallelepiped.""" def draw_box(axes, size, view): """Draw a wireframe box at a particular view.""" a, b, c = size x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1]) y = b*np.array([+1, +1, -1, -1, +1, +1, -1, -1]) z = c*np.array([+1, +1, +1, +1, -1, -1, -1, -1]) x, y, z = transform_xyz(view, None, x, y, z) def _draw(i, j): axes.plot([x[i], x[j]], [y[i], y[j]], [z[i], z[j]], color='black') _draw(0, 1) _draw(0, 2) _draw(0, 3) _draw(7, 4) _draw(7, 5) _draw(7, 6) def draw_parallelepiped(axes, size, view, jitter, steps=None, color=(0.6, 1.0, 0.6), alpha=1): """Draw a parallelepiped surface, with view and jitter.""" # pylint: disable=unused-argument a, b, c = size x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1]) x, y, z = transform_xyz(view, jitter, x, y, z) axes.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha) # Draw pink face on box. axes.plot_trisurf(x, y, triangles=tri, Z=z, color=color, alpha=alpha, linewidth=0) # Colour the c+ face of the box. # Since I can't control face color, instead draw a thin box situated just # in front of the "c+" face. Use the c face so that rotations about psi # rotate that face. if 1: if 0: # pylint: disable=using-constant-test color = (1, 0.6, 0.6) # pink x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1]) y = b*np.array([+1, +1, -1, -1, +1, +1, -1, -1]) z = c*np.array([+1, +1, +1, +1, -1, -1, -1, -1]) x, y, z = transform_xyz(view, jitter, x, y, abs(z)+0.001) axes.plot_trisurf(x, y, triangles=tri, Z=z, color=[1, 0.6, 0.6], alpha=alpha) axes.plot_trisurf(x, y, triangles=tri, Z=z, color=color, alpha=alpha) draw_labels(axes, view, jitter, [ ]) def draw_sphere(axes, radius=10., steps=100): def draw_sphere(axes, radius=1.0, steps=25, center=(0, 0, 0), color='w', alpha=1.): """Draw a sphere""" u = np.linspace(0, 2 * np.pi, steps) v = np.linspace(0, np.pi, steps) x = radius * np.outer(np.cos(u), np.sin(v)) y = radius * np.outer(np.sin(u), np.sin(v)) z = radius * np.outer(np.ones(np.size(u)), np.cos(v)) axes.plot_surface(x, y, z, rstride=4, cstride=4, color='w') def draw_jitter(axes, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0), draw_shape=draw_parallelepiped): u = np.linspace(0, 2 * pi, steps) v = np.linspace(0, pi, steps) x = radius * np.outer(cos(u), sin(v)) + center[0] y = radius * np.outer(sin(u), sin(v)) + center[1] z = radius * np.outer(np.ones(np.size(u)), cos(v)) + center[2] axes.plot_surface(x, y, z, color=color, alpha=alpha) #axes.plot_wireframe(x, y, z) def draw_axes(axes, origin=(-1, -1, -1), length=(2, 2, 2)): """Draw wireframe axes lines, with given origin and length""" x, y, z = origin dx, dy, dz = length axes.plot([x, x+dx], [y, y], [z, z], color='black') axes.plot([x, x], [y, y+dy], [z, z], color='black') axes.plot([x, x], [y, y], [z, z+dz], color='black') def draw_person_on_sphere(axes, view, height=0.5, radius=1.0): """ Draw a person on the surface of a sphere. *view* indicates (latitude, longitude, orientation) """ limb_offset = height * 0.05 head_radius = height * 0.10 head_height = height - head_radius neck_length = head_radius * 0.50 shoulder_height = height - 2*head_radius - neck_length torso_length = shoulder_height * 0.55 torso_radius = torso_length * 0.30 leg_length = shoulder_height - torso_length arm_length = torso_length * 0.90 def _draw_part(y, z): x = np.zeros_like(y) xp, yp, zp = transform_xyz(view, None, x, y, z + radius) axes.plot(xp, yp, zp, color='k') # circle for head u = np.linspace(0, 2 * pi, 40) y = head_radius * cos(u) z = head_radius * sin(u) + head_height _draw_part(y, z) # rectangle for body y = np.array([-torso_radius, torso_radius, torso_radius, -torso_radius, -torso_radius]) z = np.array([0., 0, torso_length, torso_length, 0]) + leg_length _draw_part(y, z) # arms y = np.array([-torso_radius - limb_offset, -torso_radius - limb_offset, -torso_radius]) z = np.array([shoulder_height - arm_length, shoulder_height, shoulder_height]) _draw_part(y, z) _draw_part(-y, z) # pylint: disable=invalid-unary-operand-type # legs y = np.array([-torso_radius + limb_offset, -torso_radius + limb_offset]) z = np.array([0, leg_length]) _draw_part(y, z) _draw_part(-y, z) # pylint: disable=invalid-unary-operand-type limits = [-radius-height, radius+height] axes.set_xlim(limits) axes.set_ylim(limits) axes.set_zlim(limits) axes.set_axis_off() def draw_jitter(axes, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0), draw_shape=draw_parallelepiped, projection='equirectangular', alpha=0.8, views=None): """ Represent jitter as a set of shapes at different orientations. """ project, project_weight = get_projection(projection) # set max diagonal to 0.95 scale = 0.95/sqrt(sum(v**2 for v in size)) size = tuple(scale*v for v in size) #np.random.seed(10) #cloud = np.random.randn(10,3) cloud = [ [-1, -1, -1], [-1, -1, +0], [-1, -1, +1], [-1, +0, -1], [-1, +0, +0], [-1, +0, +1], [-1, +1, -1], [-1, +1, +0], [-1, +1, +1], [+0, -1, -1], [+0, -1, +0], [+0, -1, +1], [+0, +0, -1], [+0, +0, +0], [+0, +0, +1], [+0, +1, -1], [+0, +1, +0], [+0, +1, +1], [+1, -1, -1], [+1, -1, +0], [+1, -1, +1], [+1, +0, -1], [+1, +0, +0], [+1, +0, +1], [+1, +1, -1], [+1, +1, +0], [+1, +1, +1], ] dtheta, dphi, dpsi = jitter if dtheta == 0: cloud = [v for v in cloud if v[0] == 0] if dphi == 0: cloud = [v for v in cloud if v[1] == 0] if dpsi == 0: cloud = [v for v in cloud if v[2] == 0] draw_shape(axes, size, view, [0, 0, 0], steps=100, alpha=0.8) scale = {'gaussian':1, 'rectangle':1/sqrt(3), 'uniform':1/3}[dist] for point in cloud: delta = [scale*dtheta*point[0], scale*dphi*point[1], scale*dpsi*point[2]] draw_shape(axes, size, view, delta, alpha=0.8) base = {'gaussian':3, 'rectangle':sqrt(3), 'uniform':1}[dist] def _steps(delta): if views is None: n = max(3, min(25, 2*int(base*delta/5))) else: n = views return base*delta*np.linspace(-1, 1, n) if delta > 0 else [0.] for theta in _steps(dtheta): for phi in _steps(dphi): for psi in _steps(dpsi): w = project_weight(theta, phi, 1.0, 1.0) if w > 0: dview = project(theta, phi, psi) draw_shape(axes, size, view, dview, alpha=alpha) for v in 'xyz': a, b, c = size lim = np.sqrt(a**2 + b**2 + c**2) lim = sqrt(a**2 + b**2 + c**2) getattr(axes, 'set_'+v+'lim')([-lim, lim]) getattr(axes, v+'axis').label.set_text(v) #getattr(axes, v+'axis').label.set_text(v) PROJECTIONS = [ 'azimuthal_equal_area', ] def draw_mesh(axes, view, jitter, radius=1.2, n=11, dist='gaussian', projection='equirectangular'): """ Draw the dispersion mesh showing the theta-phi orientations at which the model will be evaluated. def get_projection(projection): """ jitter projections Should allow free movement in theta, but phi is distorted. """ # pylint: disable=unused-argument # TODO: try Kent distribution instead of a gaussian warped by projection dist_x = np.linspace(-1, 1, n) weights = np.ones_like(dist_x) if dist == 'gaussian': dist_x *= 3 weights = exp(-0.5*dist_x**2) elif dist == 'rectangle': # Note: uses sasmodels ridiculous definition of rectangle width dist_x *= sqrt(3) elif dist == 'uniform': pass else: raise ValueError("expected dist to be gaussian, rectangle or uniform") if projection == 'equirectangular': #define PROJECTION 1 def _rotate(theta_i, phi_j): return Rx(phi_j)*Ry(theta_i) def _project(theta_i, phi_j, psi): latitude, longitude = theta_i, phi_j return latitude, longitude, psi, 'xyz' #return Rx(phi_j)*Ry(theta_i) def _weight(theta_i, phi_j, w_i, w_j): return w_i*w_j*abs(cos(radians(theta_i))) elif projection == 'sinusoidal': #define PROJECTION 2 def _rotate(theta_i, phi_j): def _project(theta_i, phi_j, psi): latitude = theta_i scale = cos(radians(latitude)) longitude = phi_j/scale if abs(phi_j) < abs(scale)*180 else 0 #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) return Rx(longitude)*Ry(latitude) def _weight(theta_i, phi_j, w_i, w_j): return latitude, longitude, psi, 'xyz' #return Rx(longitude)*Ry(latitude) def _project(theta_i, phi_j, w_i, w_j): latitude = theta_i scale = cos(radians(latitude)) return active*w_i*w_j elif projection == 'guyou': #define PROJECTION 3 (eventually?) def _rotate(theta_i, phi_j): def _project(theta_i, phi_j, psi): from .guyou import guyou_invert #latitude, longitude = guyou_invert([theta_i], [phi_j]) longitude, latitude = guyou_invert([phi_j], [theta_i]) return Rx(longitude[0])*Ry(latitude[0]) return latitude, longitude, psi, 'xyz' #return Rx(longitude[0])*Ry(latitude[0]) def _weight(theta_i, phi_j, w_i, w_j): return w_i*w_j elif projection == 'azimuthal_equidistance': # Note: Rz Ry, not Rx Ry def _rotate(theta_i, phi_j): elif projection == 'azimuthal_equidistance': # Note that calculates angles for Rz Ry rather than Rx Ry def _project(theta_i, phi_j, psi): latitude = sqrt(theta_i**2 + phi_j**2) longitude = degrees(np.arctan2(phi_j, theta_i)) longitude = degrees(arctan2(phi_j, theta_i)) #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) return Rz(longitude)*Ry(latitude) return latitude, longitude, psi-longitude, 'zyz' #R = Rz(longitude)*Ry(latitude)*Rz(psi) #return R_to_xyz(R) #return Rz(longitude)*Ry(latitude) def _weight(theta_i, phi_j, w_i, w_j): # Weighting for each point comes from the integral: return weight*w_i*w_j if latitude < 180 else 0 elif projection == 'azimuthal_equal_area': def _rotate(theta_i, phi_j): # Note that calculates angles for Rz Ry rather than Rx Ry def _project(theta_i, phi_j, psi): radius = min(1, sqrt(theta_i**2 + phi_j**2)/180) latitude = 180-degrees(2*np.arccos(radius)) longitude = degrees(np.arctan2(phi_j, theta_i)) latitude = 180-degrees(2*arccos(radius)) longitude = degrees(arctan2(phi_j, theta_i)) #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) return Rz(longitude)*Ry(latitude) return latitude, longitude, psi, 'zyz' #R = Rz(longitude)*Ry(latitude)*Rz(psi) #return R_to_xyz(R) #return Rz(longitude)*Ry(latitude) def _weight(theta_i, phi_j, w_i, w_j): latitude = sqrt(theta_i**2 + phi_j**2) raise ValueError("unknown projection %r"%projection) return _project, _weight def R_to_xyz(R): """ Return phi, theta, psi Tait-Bryan angles corresponding to the given rotation matrix. Extracting Euler Angles from a Rotation Matrix Mike Day, Insomniac Games https://d3cw3dd2w32x2b.cloudfront.net/wp-content/uploads/2012/07/euler-angles1.pdf Based on: Shoemakeâs "Euler Angle Conversion", Graphics Gems IV, pp. 222-229 """ phi = arctan2(R[1, 2], R[2, 2]) theta = arctan2(-R[0, 2], sqrt(R[0, 0]**2 + R[0, 1]**2)) psi = arctan2(R[0, 1], R[0, 0]) return degrees(phi), degrees(theta), degrees(psi) def draw_mesh(axes, view, jitter, radius=1.2, n=11, dist='gaussian', projection='equirectangular'): """ Draw the dispersion mesh showing the theta-phi orientations at which the model will be evaluated. """ _project, _weight = get_projection(projection) def _rotate(theta, phi, z): dview = _project(theta, phi, 0.) if dview[3] == 'zyz': return Rz(dview[1])*Ry(dview[0])*z else: # dview[3] == 'xyz': return Rx(dview[1])*Ry(dview[0])*z dist_x = np.linspace(-1, 1, n) weights = np.ones_like(dist_x) if dist == 'gaussian': dist_x *= 3 weights = exp(-0.5*dist_x**2) elif dist == 'rectangle': # Note: uses sasmodels ridiculous definition of rectangle width dist_x *= sqrt(3) elif dist == 'uniform': pass else: raise ValueError("expected dist to be gaussian, rectangle or uniform") # mesh in theta, phi formed by rotating z dtheta, dphi, dpsi = jitter dtheta, dphi, dpsi = jitter # pylint: disable=unused-variable z = np.matrix([[0], [0], [radius]]) points = np.hstack([_rotate(theta_i, phi_j)*z points = np.hstack([_rotate(theta_i, phi_j, z) for theta_i in dtheta*dist_x for phi_j in dphi*dist_x]) Points are stored in a 3 x n numpy matrix, not a numpy array or tuple. """ dtheta, dphi, dpsi = jitter points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points if jitter is None: return points # Hack to deal with the fact that azimuthal_equidistance uses euler angles if len(jitter) == 4: dtheta, dphi, dpsi, _ = jitter points = Rz(dphi)*Ry(dtheta)*Rz(dpsi)*points else: dtheta, dphi, dpsi = jitter points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points return points """ theta, phi, psi = view points = Rz(phi)*Ry(theta)*Rz(psi)*points points = Rz(phi)*Ry(theta)*Rz(psi)*points # viewing angle #points = Rz(psi)*Ry(pi/2-theta)*Rz(phi)*points # 1-D integration angles #points = Rx(phi)*Ry(theta)*Rz(psi)*points # angular dispersion angle return points def orient_relative_to_beam_quaternion(view, points): """ Apply the view transform to a set of points. Points are stored in a 3 x n numpy matrix, not a numpy array or tuple. This variant uses quaternions rather than rotation matrices for the computation. It works but it is not used because it doesn't solve any problems. The challenge of mapping theta/phi/psi to SO(3) does not disappear by calculating the transform differently. """ theta, phi, psi = view x, y, z = [1, 0, 0], [0, 1, 0], [0, 0, 1] q = Quaternion(1, [0, 0, 0]) ## Compose a rotation about the three axes by rotating ## the unit vectors before applying the rotation. #q = Quaternion.from_angle_axis(theta, q.rot(x)) * q #q = Quaternion.from_angle_axis(phi, q.rot(y)) * q #q = Quaternion.from_angle_axis(psi, q.rot(z)) * q ## The above turns out to be equivalent to reversing ## the order of application, so ignore it and use below. q = q * Quaternion.from_angle_axis(theta, x) q = q * Quaternion.from_angle_axis(phi, y) q = q * Quaternion.from_angle_axis(psi, z) ## Reverse the order by post-multiply rather than pre-multiply #q = Quaternion.from_angle_axis(theta, x) * q #q = Quaternion.from_angle_axis(phi, y) * q #q = Quaternion.from_angle_axis(psi, z) * q #print("axes psi", q.rot(np.matrix([x, y, z]).T)) return q.rot(points) #orient_relative_to_beam = orient_relative_to_beam_quaternion # === Quaterion class definition === BEGIN # Simple stand-alone quaternion class # Note: this code works but isn't unused since quaternions didn't solve the # representation problem. Leave it here in case we want to revisit this later. #import numpy as np class Quaternion(object): r""" Quaternion(w, r) = w + ir[0] + jr[1] + kr[2] Quaternion.from_angle_axis(theta, r) for a rotation of angle theta about an axis oriented toward the direction r. This defines a unit quaternion, normalizing$r$to the unit vector$\hat r$, and setting quaternion$Q = \cos \theta + \sin \theta \hat r$Quaternion objects can be multiplied, which applies a rotation about the given axis, allowing composition of rotations without risk of gimbal lock. The resulting quaternion is applied to a set of points using *Q.rot(v)*. """ def __init__(self, w, r): self.w = w self.r = np.asarray(r, dtype='d') @staticmethod def from_angle_axis(theta, r): """Build quaternion as rotation theta about axis r""" theta = np.radians(theta)/2 r = np.asarray(r) w = np.cos(theta) r = np.sin(theta)*r/np.dot(r, r) return Quaternion(w, r) def __mul__(self, other): """Multiply quaterions""" if isinstance(other, Quaternion): w = self.w*other.w - np.dot(self.r, other.r) r = self.w*other.r + other.w*self.r + np.cross(self.r, other.r) return Quaternion(w, r) raise NotImplementedError("Quaternion * non-quaternion not implemented") def rot(self, v): """Transform point *v* by quaternion""" v = np.asarray(v).T use_transpose = (v.shape[-1] != 3) if use_transpose: v = v.T v = v + np.cross(2*self.r, np.cross(self.r, v) + self.w*v) #v = v + 2*self.w*np.cross(self.r, v) + np.cross(2*self.r, np.cross(self.r, v)) if use_transpose: v = v.T return v.T def conj(self): """Conjugate quaternion""" return Quaternion(self.w, -self.r) def inv(self): """Inverse quaternion""" return self.conj()/self.norm()**2 def norm(self): """Quaternion length""" return np.sqrt(self.w**2 + np.sum(self.r**2)) def __str__(self): return "%g%+gi%+gj%+gk"%(self.w, self.r[0], self.r[1], self.r[2]) def test_qrot(): """Quaternion checks""" # Define rotation of 60 degrees around an axis in y-z that is 60 degrees # from y. The rotation axis is determined by rotating the point [0, 1, 0] # about x. ax = Quaternion.from_angle_axis(60, [1, 0, 0]).rot([0, 1, 0]) q = Quaternion.from_angle_axis(60, ax) # Set the point to be rotated, and its expected rotated position. p = [1, -1, 2] target = [(10+4*np.sqrt(3))/8, (1+2*np.sqrt(3))/8, (14-3*np.sqrt(3))/8] #print(q, q.rot(p) - target) assert max(abs(q.rot(p) - target)) < 1e-14 #test_qrot() #import sys; sys.exit() # === Quaterion class definition === END # translate between number of dimension of dispersity and the number of or the top of the range, depending on whether *mode* is 'central' or 'top'. """ if portion == 1.0: return data.min(), data.max() elif mode == 'central': data = np.sort(data.flatten()) offset = int(portion*len(data)/2 + 0.5) return data[offset], data[-offset] elif mode == 'top': data = np.sort(data.flatten()) offset = int(portion*len(data) + 0.5) return data[offset], data[-1] if portion < 1.0: if mode == 'central': data = np.sort(data.flatten()) offset = int(portion*len(data)/2 + 0.5) return data[offset], data[-offset] if mode == 'top': data = np.sort(data.flatten()) offset = int(portion*len(data) + 0.5) return data[offset], data[-1] # Default: full range return data.min(), data.max() def draw_scattering(calculator, axes, view, jitter, dist='gaussian'): # compute the pattern qx, qy = calculator._data.x_bins, calculator._data.y_bins qx, qy = calculator.qxqy Iqxy = calculator(**pars).reshape(len(qx), len(qy)) # scale it and draw it Iqxy = np.log(Iqxy) Iqxy = log(Iqxy) if calculator.limits: # use limits from orientation (0,0,0) vmin = vmax*10**-7 #vmin, vmax = clipped_range(Iqxy, portion=portion, mode='top') #vmin, vmax = Iqxy.min(), Iqxy.max() #print("range",(vmin,vmax)) #qx, qy = np.meshgrid(qx, qy) if 0: if 0: # pylint: disable=using-constant-test level = np.asarray(255*(Iqxy - vmin)/(vmax - vmin), 'i') level[level < 0] = 0 from matplotlib import pylab as plt colors = plt.get_cmap()(level) axes.plot_surface(qx, qy, -1.1, rstride=1, cstride=1, facecolors=colors) elif 1: #from matplotlib import cm #colors = cm.coolwarm(level) #colors = cm.gist_yarg(level) #colors = cm.Wistia(level) colors[level <= 0, 3] = 0. # set floor to transparent x, y = np.meshgrid(qx/qx.max(), qy/qy.max()) axes.plot_surface(x, y, -1.1*np.ones_like(x), facecolors=colors) elif 1: # pylint: disable=using-constant-test axes.contourf(qx/qx.max(), qy/qy.max(), Iqxy, zdir='z', offset=-1.1, levels=np.linspace(vmin, vmax, 24)) calculator = DirectModel(data, model) # Remember the data axes so we can plot the results calculator.qxqy = (q, q) # stuff the values for non-orientation parameters into the calculator calculator.pars = pars.copy() # under rotation or angular dispersion Iqxy = calculator(theta=0, phi=0, psi=0, **calculator.pars) Iqxy = np.log(Iqxy) Iqxy = log(Iqxy) vmin, vmax = clipped_range(Iqxy, 0.95, mode='top') calculator.limits = vmin, vmax+1 } def run(model_name='parallelepiped', size=(10, 40, 100), view=(0, 0, 0), jitter=(0, 0, 0), dist='gaussian', mesh=30, projection='equirectangular'): *size* gives the dimensions (a, b, c) of the shape. *view* gives the initial view (theta, phi, psi) of the shape. *view* gives the initial jitter (dtheta, dphi, dpsi) of the shape. *dist* is the type of dispersition: gaussian, rectangle, or uniform. calculator, size = select_calculator(model_name, n=150, size=size) draw_shape = DRAW_SHAPES.get(model_name, draw_parallelepiped) #draw_shape = draw_fcc ## uncomment to set an independent the colour range for every view calculator.limits = None ## initial view #theta, dtheta = 70., 10. #phi, dphi = -45., 3. #psi, dpsi = -45., 3. theta, phi, psi = 0, 0, 0 dtheta, dphi, dpsi = 0, 0, 0 PLOT_ENGINE(calculator, draw_shape, size, view, jitter, dist, mesh, projection) def _mpl_plot(calculator, draw_shape, size, view, jitter, dist, mesh, projection): # Note: travis-ci does not support mpl_toolkits.mplot3d, but this shouldn't be # an issue since we are lazy-loading the package on a path that isn't tested. # Importing mplot3d adds projection='3d' option to subplot import mpl_toolkits.mplot3d # pylint: disable=unused-variable import matplotlib as mpl import matplotlib.pyplot as plt from matplotlib.widgets import Slider ## create the plot window ## add control widgets to plot axes_theta = plt.axes([0.1, 0.15, 0.45, 0.04], **props) axes_phi = plt.axes([0.1, 0.1, 0.45, 0.04], **props) axes_psi = plt.axes([0.1, 0.05, 0.45, 0.04], **props) stheta = Slider(axes_theta, 'Theta', -90, 90, valinit=theta) sphi = Slider(axes_phi, 'Phi', -180, 180, valinit=phi) spsi = Slider(axes_psi, 'Psi', -180, 180, valinit=psi) axes_dtheta = plt.axes([0.75, 0.15, 0.15, 0.04], **props) axes_dphi = plt.axes([0.75, 0.1, 0.15, 0.04], **props) axes_dpsi = plt.axes([0.75, 0.05, 0.15, 0.04], **props) axes_theta = plt.axes([0.05, 0.15, 0.50, 0.04], **props) axes_phi = plt.axes([0.05, 0.10, 0.50, 0.04], **props) axes_psi = plt.axes([0.05, 0.05, 0.50, 0.04], **props) stheta = Slider(axes_theta, u'Îž', -90, 90, valinit=0) sphi = Slider(axes_phi, u'Ï', -180, 180, valinit=0) spsi = Slider(axes_psi, u'Ï', -180, 180, valinit=0) axes_dtheta = plt.axes([0.70, 0.15, 0.20, 0.04], **props) axes_dphi = plt.axes([0.70, 0.1, 0.20, 0.04], **props) axes_dpsi = plt.axes([0.70, 0.05, 0.20, 0.04], **props) # Note: using ridiculous definition of rectangle distribution, whose width # in sasmodels is sqrt(3) times the given width. Divide by sqrt(3) to keep # the maximum width to 90. dlimit = DIST_LIMITS[dist] sdtheta = Slider(axes_dtheta, 'dTheta', 0, 2*dlimit, valinit=dtheta) sdphi = Slider(axes_dphi, 'dPhi', 0, 2*dlimit, valinit=dphi) sdpsi = Slider(axes_dpsi, 'dPsi', 0, 2*dlimit, valinit=dpsi) sdtheta = Slider(axes_dtheta, u'ÎÎž', 0, 2*dlimit, valinit=0) sdphi = Slider(axes_dphi, u'ÎÏ', 0, 2*dlimit, valinit=0) sdpsi = Slider(axes_dpsi, u'ÎÏ', 0, 2*dlimit, valinit=0) ## initial view and jitter theta, phi, psi = view stheta.set_val(theta) sphi.set_val(phi) spsi.set_val(psi) dtheta, dphi, dpsi = jitter sdtheta.set_val(dtheta) sdphi.set_val(dphi) sdpsi.set_val(dpsi) ## callback to draw the new view def update(val, axis=None): def _update(val, axis=None): # pylint: disable=unused-argument view = stheta.val, sphi.val, spsi.val jitter = sdtheta.val, sdphi.val, sdpsi.val # set small jitter as 0 if multiple pd dims dims = sum(v > 0 for v in jitter) limit = [0, 0.5, 5][dims] limit = [0, 0.5, 5, 5][dims] jitter = [0 if v < limit else v for v in jitter] axes.cla() draw_beam(axes, (0, 0)) draw_jitter(axes, view, jitter, dist=dist, size=size, draw_shape=draw_shape) #draw_jitter(axes, view, (0,0,0)) ## Visualize as person on globe #draw_sphere(axes, radius=0.5) #draw_person_on_sphere(axes, view, radius=0.5) ## Move beam instead of shape #draw_beam(axes, -view[:2]) #draw_jitter(axes, (0,0,0), (0,0,0), views=3) ## Move shape and draw scattering draw_beam(axes, (0, 0), alpha=1.) #draw_person_on_sphere(axes, view, radius=1.2, height=0.5) draw_jitter(axes, view, jitter, dist=dist, size=size, alpha=1., draw_shape=draw_shape, projection=projection, views=3) draw_mesh(axes, view, jitter, dist=dist, n=mesh, projection=projection) draw_scattering(calculator, axes, view, jitter, dist=dist) plt.gcf().canvas.draw() ## bind control widgets to view updater stheta.on_changed(lambda v: update(v, 'theta')) sphi.on_changed(lambda v: update(v, 'phi')) spsi.on_changed(lambda v: update(v, 'psi')) sdtheta.on_changed(lambda v: update(v, 'dtheta')) sdphi.on_changed(lambda v: update(v, 'dphi')) sdpsi.on_changed(lambda v: update(v, 'dpsi')) stheta.on_changed(lambda v: _update(v, 'theta')) sphi.on_changed(lambda v: _update(v, 'phi')) spsi.on_changed(lambda v: _update(v, 'psi')) sdtheta.on_changed(lambda v: _update(v, 'dtheta')) sdphi.on_changed(lambda v: _update(v, 'dphi')) sdpsi.on_changed(lambda v: _update(v, 'dpsi')) ## initialize view update(None, 'phi') _update(None, 'phi') ## go interactive plt.show() def map_colors(z, kw): """ Process matplotlib-style colour arguments. Pulls 'cmap', 'alpha', 'vmin', and 'vmax' from th *kw* dictionary, setting the *kw['color']* to an RGB array. These are ignored if 'c' or 'color' are set inside *kw*. """ from matplotlib import cm cmap = kw.pop('cmap', cm.coolwarm) alpha = kw.pop('alpha', None) vmin = kw.pop('vmin', z.min()) vmax = kw.pop('vmax', z.max()) c = kw.pop('c', None) color = kw.pop('color', c) if color is None: znorm = ((z - vmin) / (vmax - vmin)).clip(0, 1) color = cmap(znorm) elif isinstance(color, np.ndarray) and color.shape == z.shape: color = cmap(color) if alpha is None: if isinstance(color, np.ndarray): color = color[..., :3] else: color[..., 3] = alpha kw['color'] = color def make_vec(*args): """Turn all elements of *args* into numpy arrays""" #return [np.asarray(v, 'd').flatten() for v in args] return [np.asarray(v, 'd') for v in args] def make_image(z, kw): """Convert numpy array *z* into a *PIL* RGB image.""" import PIL.Image from matplotlib import cm cmap = kw.pop('cmap', cm.coolwarm) znorm = (z-z.min())/z.ptp() c = cmap(znorm) c = c[..., :3] rgb = np.asarray(c*255, 'u1') image = PIL.Image.fromarray(rgb, mode='RGB') return image _IPV_MARKERS = { 'o': 'sphere', } _IPV_COLORS = { 'w': 'white', 'k': 'black', 'c': 'cyan', 'm': 'magenta', 'y': 'yellow', 'r': 'red', 'g': 'green', 'b': 'blue', } def _ipv_fix_color(kw): alpha = kw.pop('alpha', None) color = kw.get('color', None) if isinstance(color, str): color = _IPV_COLORS.get(color, color) kw['color'] = color if alpha is not None: color = kw['color'] #TODO: convert color to [r, g, b, a] if not already if isinstance(color, (tuple, list)): if len(color) == 3: color = (color[0], color[1], color[2], alpha) else: color = (color[0], color[1], color[2], alpha*color[3]) color = np.array(color) if isinstance(color, np.ndarray) and color.shape[-1] == 4: color[..., 3] = alpha kw['color'] = color def _ipv_set_transparency(kw, obj): color = kw.get('color', None) if (isinstance(color, np.ndarray) and color.shape[-1] == 4 and (color[..., 3] != 1.0).any()): obj.material.transparent = True obj.material.side = "FrontSide" def ipv_axes(): """ Build a matplotlib style Axes interface for ipyvolume """ import ipyvolume as ipv class Axes(object): """ Matplotlib Axes3D style interface to ipyvolume renderer. """ # pylint: disable=no-self-use,no-init # transparency can be achieved by setting the following: # mesh.color = [r, g, b, alpha] # mesh.material.transparent = True # mesh.material.side = "FrontSide" # smooth(ish) rotation can be achieved by setting: # slide.continuous_update = True # figure.animation = 0. # mesh.material.x = x # mesh.material.y = y # mesh.material.z = z # maybe need to synchronize update of x/y/z to avoid shimmy when moving def plot(self, x, y, z, **kw): """mpl style plot interface for ipyvolume""" _ipv_fix_color(kw) x, y, z = make_vec(x, y, z) ipv.plot(x, y, z, **kw) def plot_surface(self, x, y, z, **kw): """mpl style plot_surface interface for ipyvolume""" facecolors = kw.pop('facecolors', None) if facecolors is not None: kw['color'] = facecolors _ipv_fix_color(kw) x, y, z = make_vec(x, y, z) h = ipv.plot_surface(x, y, z, **kw) _ipv_set_transparency(kw, h) #h.material.side = "DoubleSide" return h def plot_trisurf(self, x, y, triangles=None, Z=None, **kw): """mpl style plot_trisurf interface for ipyvolume""" kw.pop('linewidth', None) _ipv_fix_color(kw) x, y, z = make_vec(x, y, Z) if triangles is not None: triangles = np.asarray(triangles) h = ipv.plot_trisurf(x, y, z, triangles=triangles, **kw) _ipv_set_transparency(kw, h) return h def scatter(self, x, y, z, **kw): """mpl style scatter interface for ipyvolume""" x, y, z = make_vec(x, y, z) map_colors(z, kw) marker = kw.get('marker', None) kw['marker'] = _IPV_MARKERS.get(marker, marker) h = ipv.scatter(x, y, z, **kw) _ipv_set_transparency(kw, h) return h def contourf(self, x, y, v, zdir='z', offset=0, levels=None, **kw): """mpl style contourf interface for ipyvolume""" # pylint: disable=unused-argument # Don't use contour for now (although we might want to later) self.pcolor(x, y, v, zdir='z', offset=offset, **kw) def pcolor(self, x, y, v, zdir='z', offset=0, **kw): """mpl style pcolor interface for ipyvolume""" # pylint: disable=unused-argument x, y, v = make_vec(x, y, v) image = make_image(v, kw) xmin, xmax = x.min(), x.max() ymin, ymax = y.min(), y.max() x = np.array([[xmin, xmax], [xmin, xmax]]) y = np.array([[ymin, ymin], [ymax, ymax]]) z = x*0 + offset u = np.array([[0., 1], [0, 1]]) v = np.array([[0., 0], [1, 1]]) h = ipv.plot_mesh(x, y, z, u=u, v=v, texture=image, wireframe=False) _ipv_set_transparency(kw, h) h.material.side = "DoubleSide" return h def text(self, *args, **kw): """mpl style text interface for ipyvolume""" pass def set_xlim(self, limits): """mpl style set_xlim interface for ipyvolume""" ipv.xlim(*limits) def set_ylim(self, limits): """mpl style set_ylim interface for ipyvolume""" ipv.ylim(*limits) def set_zlim(self, limits): """mpl style set_zlim interface for ipyvolume""" ipv.zlim(*limits) def set_axes_on(self): """mpl style set_axes_on interface for ipyvolume""" ipv.style.axis_on() def set_axis_off(self): """mpl style set_axes_off interface for ipyvolume""" ipv.style.axes_off() return Axes() def _ipv_plot(calculator, draw_shape, size, view, jitter, dist, mesh, projection): from IPython.display import display import ipywidgets as widgets import ipyvolume as ipv axes = ipv_axes() def _draw(view, jitter): camera = ipv.gcf().camera #print(ipv.gcf().__dict__.keys()) #print(dir(ipv.gcf())) ipv.figure(animation=0.) # no animation when updating object mesh # set small jitter as 0 if multiple pd dims dims = sum(v > 0 for v in jitter) limit = [0, 0.5, 5, 5][dims] jitter = [0 if v < limit else v for v in jitter] ## Visualize as person on globe #draw_beam(axes, (0, 0)) #draw_sphere(axes, radius=0.5) #draw_person_on_sphere(axes, view, radius=0.5) ## Move beam instead of shape #draw_beam(axes, view=(-view[0], -view[1])) #draw_jitter(axes, view=(0,0,0), jitter=(0,0,0)) ## Move shape and draw scattering draw_beam(axes, (0, 0), steps=25) draw_jitter(axes, view, jitter, dist=dist, size=size, alpha=1.0, draw_shape=draw_shape, projection=projection) draw_mesh(axes, view, jitter, dist=dist, n=mesh, radius=0.95, projection=projection) draw_scattering(calculator, axes, view, jitter, dist=dist) draw_axes(axes, origin=(-1, -1, -1.1)) ipv.style.box_off() ipv.style.axes_off() ipv.xyzlabel(" ", " ", " ") ipv.gcf().camera = camera ipv.show() trange, prange = (-180., 180., 1.), (-180., 180., 1.) dtrange, dprange = (0., 180., 1.), (0., 180., 1.) ## Super simple interfaca, but uses non-ascii variable namese # Îž Ï Ï ÎÎž ÎÏ ÎÏ #def update(**kw): # view = kw['Îž'], kw['Ï'], kw['Ï'] # jitter = kw['ÎÎž'], kw['ÎÏ'], kw['ÎÏ'] # draw(view, jitter) #widgets.interact(update, Îž=trange, Ï=prange, Ï=prange, ÎÎž=dtrange, ÎÏ=dprange, ÎÏ=dprange) def _update(theta, phi, psi, dtheta, dphi, dpsi): _draw(view=(theta, phi, psi), jitter=(dtheta, dphi, dpsi)) def _slider(name, slice, init=0.): return widgets.FloatSlider( value=init, min=slice[0], max=slice[1], step=slice[2], description=name, disabled=False, #continuous_update=True, continuous_update=False, orientation='horizontal', readout=True, readout_format='.1f', ) theta = _slider(u'Îž', trange, view[0]) phi = _slider(u'Ï', prange, view[1]) psi = _slider(u'Ï', prange, view[2]) dtheta = _slider(u'ÎÎž', dtrange, jitter[0]) dphi = _slider(u'ÎÏ', dprange, jitter[1]) dpsi = _slider(u'ÎÏ', dprange, jitter[2]) fields = { 'theta': theta, 'phi': phi, 'psi': psi, 'dtheta': dtheta, 'dphi': dphi, 'dpsi': dpsi, } ui = widgets.HBox([ widgets.VBox([theta, phi, psi]), widgets.VBox([dtheta, dphi, dpsi]) ]) out = widgets.interactive_output(_update, fields) display(ui, out) _ENGINES = { "matplotlib": _mpl_plot, "mpl": _mpl_plot, #"plotly": _plotly_plot, "ipvolume": _ipv_plot, "ipv": _ipv_plot, } PLOT_ENGINE = _ENGINES["matplotlib"] def set_plotter(name): """ Setting the plotting engine to matplotlib/ipyvolume or equivalently mpl/ipv. """ global PLOT_ENGINE PLOT_ENGINE = _ENGINES[name] def main(): """ Command line interface to the jitter viewer. """ parser = argparse.ArgumentParser( description="Display jitter", parser.add_argument('-s', '--size', type=str, default='10,40,100', help='a,b,c lengths') parser.add_argument('-v', '--view', type=str, default='0,0,0', help='initial view angles') parser.add_argument('-j', '--jitter', type=str, default='0,0,0', help='initial angular dispersion') parser.add_argument('-d', '--distribution', choices=DISTRIBUTIONS, default=DISTRIBUTIONS[0], help='oriented shape') opts = parser.parse_args() size = tuple(int(v) for v in opts.size.split(',')) run(opts.shape, size=size, size = tuple(float(v) for v in opts.size.split(',')) view = tuple(float(v) for v in opts.view.split(',')) jitter = tuple(float(v) for v in opts.jitter.split(',')) run(opts.shape, size=size, view=view, jitter=jitter, mesh=opts.mesh, dist=opts.distribution, projection=opts.projection) • sasmodels/kernel.py  r36a2418 class KernelModel(object): """ Model definition for the compute engine. """ info = None # type: ModelInfo dtype = None # type: np.dtype def make_kernel(self, q_vectors): # type: (List[np.ndarray]) -> "Kernel" """ Instantiate a kernel for evaluating the model at *q_vectors*. """ raise NotImplementedError("need to implement make_kernel") def release(self): # type: () -> None """ Free resources associated with the kernel. """ pass class Kernel(object): #: kernel dimension, either "1d" or "2d" """ Instantiated model for the compute engine, applied to a particular *q*. Subclasses should define :meth:_call_kernel to evaluate the kernel over its inputs. """ #: Kernel dimension, either "1d" or "2d". dim = None # type: str #: Model info. info = None # type: ModelInfo results = None # type: List[np.ndarray] #: Numerical precision for the computation. dtype = None # type: np.dtype #: q values at which the kernel is to be evaluated q_input = None # type: Any #: Place to hold result of :meth:_call_kernel for subclass. result = None # type: np.ndarray def Iq(self, call_details, values, cutoff, magnetic): built into the model, and do not need an additional scale. """ _, F2, _, shell_volume, _ = self.Fq(call_details, values, cutoff, magnetic, effective_radius_type=0) _, F2, _, shell_volume, _ = self.Fq(call_details, values, cutoff, magnetic, radius_effective_mode=0) combined_scale = values[0]/shell_volume background = values[1] __call__ = Iq def Fq(self, call_details, values, cutoff, magnetic, effective_radius_type=0): def Fq(self, call_details, values, cutoff, magnetic, radius_effective_mode=0): # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray r""" .. math:: I(q) = \text{scale} * P (1 + ^2/ (S - 1)) + \text{background} I(q) = \text{scale} P (1 + ^2/ (S - 1)) + \text{background} = \text{scale}/ ( + ^2 (S - 1)) + \text{background} .. math:: P(q) = \frac{\sum w_k F^2(q, x_k) / \sum w_k}{\sum w_k V_k / \sum w_k} P(q)=\frac{\sum w_k F^2(q, x_k) / \sum w_k}{\sum w_k V_k / \sum w_k} The form factor itself is scaled by volume and contrast to compute the volume fraction of the particles. The model can have several different ways to compute effective radius, with the *effective_radius_type* parameter used to select amongst them. The *radius_effective_mode* parameter used to select amongst them. The volume fraction of particles should be determined from the total volume fraction of the form, not just the shell volume fraction. hollow and solid shapes. """ self._call_kernel(call_details, values, cutoff, magnetic, effective_radius_type) self._call_kernel(call_details, values, cutoff, magnetic, radius_effective_mode) #print("returned",self.q_input.q, self.result) nout = 2 if self.info.have_Fq and self.dim == '1d' else 1 form_volume = self.result[nout*self.q_input.nq + 1]/total_weight shell_volume = self.result[nout*self.q_input.nq + 2]/total_weight effective_radius = self.result[nout*self.q_input.nq + 3]/total_weight radius_effective = self.result[nout*self.q_input.nq + 3]/total_weight if shell_volume == 0.: shell_volume = 1. F1 = self.result[1:nout*self.q_input.nq:nout]/total_weight if nout == 2 else None F1 = (self.result[1:nout*self.q_input.nq:nout]/total_weight if nout == 2 else None) F2 = self.result[0:nout*self.q_input.nq:nout]/total_weight return F1, F2, effective_radius, shell_volume, form_volume/shell_volume return F1, F2, radius_effective, shell_volume, form_volume/shell_volume def release(self): # type: () -> None """ Free resources associated with the kernel instance. """ pass def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): def _call_kernel(self, call_details, values, cutoff, magnetic, radius_effective_mode): # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray """ • sasmodels/kernel_iq.c  r12f4c19 // parameters in the parameter table. // CALL_VOLUME(form, shell, table) : assign form and shell values // CALL_EFFECTIVE_RADIUS(type, table) : call the R_eff function // CALL_RADIUS_EFFECTIVE(mode, table) : call the R_eff function // CALL_IQ(q, table) : call the Iq function for 1D calcs. // CALL_IQ_A(q, table) : call the Iq function with |q| for 2D data. static void set_spin_weights(double in_spin, double out_spin, double weight[6]) { double norm; in_spin = clip(in_spin, 0.0, 1.0); out_spin = clip(out_spin, 0.0, 1.0); // However, since the weights are applied to the final intensity and // are not interned inside the I(q) function, we want the full // weight and not the square root. Any function using // set_spin_weights as part of calculating an amplitude will need to // manually take that square root, but there is currently no such // function. weight[0] = (1.0-in_spin) * (1.0-out_spin); // dd weight[1] = (1.0-in_spin) * out_spin; // du weight[2] = in_spin * (1.0-out_spin); // ud weight[3] = in_spin * out_spin; // uu // weight and not the square root. Anyway no function will ever use // set_spin_weights as part of calculating an amplitude, as the weights are // related to polarisation efficiency of the instrument. The weights serve to // construct various magnet scattering cross sections, which are linear combinations // of the spin-resolved cross sections. The polarisation efficiency e_in and e_out // are parameters ranging from 0.5 (unpolarised) beam to 1 (perfect optics). // For in_spin or out_spin <0.5 one assumes a CS, where the spin is reversed/flipped // with respect to the initial supermirror polariser. The actual polarisation efficiency // in this case is however e_in/out = 1-in/out_spin. if (out_spin < 0.5){norm=1-out_spin;} else{norm=out_spin;} // The norm is needed to make sure that the scattering cross sections are //correctly weighted, such that the sum of spin-resolved measurements adds up to // the unpolarised or half-polarised scattering cross section. No intensity weighting // needed on the incoming polariser side (assuming that a user), has normalised // to the incoming flux with polariser in for SANSPOl and unpolarised beam, respectively. weight[0] = (1.0-in_spin) * (1.0-out_spin) / norm; // dd weight[1] = (1.0-in_spin) * out_spin / norm; // du weight[2] = in_spin * (1.0-out_spin) / norm; // ud weight[3] = in_spin * out_spin / norm; // uu weight[4] = weight[1]; // du.imag weight[5] = weight[2]; // ud.imag switch (xs) { default: // keep compiler happy; condition ensures xs in [0,1,2,3] case 0: // uu => sld - D M_perpx case 0: // dd => sld - D M_perpx return sld - px*perp; case 1: // ud.real => -D M_perpy case 1: // du.real => -D M_perpy return py*perp; case 2: // du.real => -D M_perpy case 2: // ud.real => -D M_perpy return py*perp; case 3: // dd => sld + D M_perpx case 3: // uu => sld + D M_perpx return sld + px*perp; } pglobal double *result, // nq+1 return values, again with padding const double cutoff, // cutoff in the dispersity weight product int32_t effective_radius_type // which effective radius to compute int32_t radius_effective_mode // which effective radius to compute ) { weighted_form += weight * form; weighted_shell += weight * shell; if (effective_radius_type != 0) { weighted_radius += weight * CALL_EFFECTIVE_RADIUS(effective_radius_type, local_values.table); if (radius_effective_mode != 0) { weighted_radius += weight * CALL_RADIUS_EFFECTIVE(radius_effective_mode, local_values.table); } BUILD_ROTATION(); • sasmodels/kernelcl.py  r0d26e91 from __future__ import print_function import sys import os import warnings from time import perf_counter as clock except ImportError: # CRUFT: python < 3.3 import sys if sys.platform.count("darwin") > 0: from time import time as clock # CRUFT: pyopencl < 2017.1 (as of June 2016 needs quotes around include path). def quote_path(v): # type: (str) -> str """ Quote the path if it is not already quoted. def fix_pyopencl_include(): # type: (None) -> None """ Monkey patch pyopencl to allow spaces in include file path. """ import pyopencl as cl if hasattr(cl, '_DEFAULT_INCLUDE_OPTIONS'): cl._DEFAULT_INCLUDE_OPTIONS = [ quote_path(v) for v in cl._DEFAULT_INCLUDE_OPTIONS # pylint: disable=protected-access import pyopencl if hasattr(pyopencl, '_DEFAULT_INCLUDE_OPTIONS'): pyopencl._DEFAULT_INCLUDE_OPTIONS = [ quote_path(v) for v in pyopencl._DEFAULT_INCLUDE_OPTIONS ] def use_opencl(): # type: () -> bool """Return True if OpenCL is the default computational engine""" sas_opencl = os.environ.get("SAS_OPENCL", "OpenCL").lower() return HAVE_OPENCL and sas_opencl != "none" and not sas_opencl.startswith("cuda") ENV = None def reset_environment(): """ Call to create a new OpenCL context, such as after a change to SAS_OPENCL. # type: () -> "GpuEnvironment" """ Return a new OpenCL context, such as after a change to SAS_OPENCL. """ global ENV ENV = GpuEnvironment() if use_opencl() else None return ENV def environment(): return [cl.create_some_context(interactive=False)] except Exception as exc: # TODO: Should warnings instead be put into logging.warn? warnings.warn(str(exc)) warnings.warn("pyopencl.create_some_context() failed. The " "environment variable 'SAS_OPENCL' or 'PYOPENCL_CTX' might " "not be set correctly") warnings.warn( "pyopencl.create_some_context() failed. The environment " "variable 'SAS_OPENCL' or 'PYOPENCL_CTX' might not be set " "correctly") return _get_default_context() def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): radius_effective_mode): # type: (CallDetails, np.ndarray, float, bool, int) -> np.ndarray env = environment() self._result_b, # Result storage. self._as_dtype(cutoff), # Probability cutoff. np.uint32(effective_radius_type), # R_eff mode. np.uint32(radius_effective_mode), # R_eff mode. ] • sasmodels/kernelcuda.py  rfa26e78 import os import warnings import logging import time def use_cuda(): # type: None -> bool """Returns True if CUDA is the default compute engine.""" sas_opencl = os.environ.get("SAS_OPENCL", "CUDA").lower() return HAVE_CUDA and sas_opencl.startswith("cuda") if not HAVE_CUDA: raise RuntimeError("CUDA startup failed with ***" + CUDA_ERROR + "***; using C compiler instead") + CUDA_ERROR + "***; using C compiler instead") reset_environment() if ENV is None: """ for match in FUNCTION_PATTERN.finditer(source): print(match.group('qualifiers').replace('\n',r'\n'), match.group('function'), '(') print(match.group('qualifiers').replace('\n', r'\n'), match.group('function'), '(') return source def release(self): """Free the CUDA device associated with this context.""" if self.context is not None: self.context.pop() def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): radius_effective_mode): # type: (CallDetails, np.ndarray, float, bool, int) -> np.ndarray self._result_b, # Result storage. self._as_dtype(cutoff), # Probability cutoff. np.uint32(effective_radius_type), # R_eff mode. np.uint32(radius_effective_mode), # R_eff mode. ] grid = partition(self.q_input.nq) • sasmodels/kerneldll.py  r3199b17 def compile(source, output): def compile_model(source, output): # type: (str, str) -> None """ with os.fdopen(system_fd, "w") as file_handle: file_handle.write(source) compile(source=filename, output=dll) compile_model(source=filename, output=dll) # Comment the following to keep the generated C file. # Note: If there is a syntax error then compile raises an error def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): radius_effective_mode): # type: (CallDetails, np.ndarray, float, bool, int) -> np.ndarray self.result.ctypes.data, # Result storage. self._as_dtype(cutoff), # Probability cutoff. effective_radius_type, # R_eff mode. radius_effective_mode, # R_eff mode. ] • sasmodels/kernelpy.py  r3199b17 from numpy import cbrt except ImportError: def cbrt(x): return x ** (1.0/3.0) def cbrt(x): """Return cubed root of x.""" return x ** (1.0/3.0) from .generate import F64 self.info = model_info self.dtype = np.dtype('d') logger.info("make python model " + self.info.name) logger.info("make python model %s", self.info.name) def make_kernel(self, q_vectors): """Instantiate the python kernel with input *q_vectors*""" q_input = PyInput(q_vectors, dtype=F64) return PyKernel(self.info, q_input) volume = model_info.form_volume shell = model_info.shell_volume radius = model_info.effective_radius radius = model_info.radius_effective self._volume = ((lambda: (shell(*volume_args), volume(*volume_args))) if shell and volume else (lambda: [volume(*volume_args)]*2) if volume else (lambda mode: 1.0)) def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): def _call_kernel(self, call_details, values, cutoff, magnetic, radius_effective_mode): # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray if magnetic: #print("Calling python kernel") #call_details.show(values) radius = ((lambda: 0.0) if effective_radius_type == 0 else (lambda: self._radius(effective_radius_type))) radius = ((lambda: 0.0) if radius_effective_mode == 0 else (lambda: self._radius(radius_effective_mode))) self.result = _loops( self._parameter_vector, self._form, self._volume, radius, • sasmodels/list_pars.py  r2d81cfe from .compare import columnize def find_pars(type=None): def find_pars(kind=None): """ Find all parameters in all models. model_info = load_model_info(name) for p in model_info.parameters.kernel_parameters: if type is None or p.type == type: if kind is None or p.type == kind: partable.setdefault(p.name, []) partable[p.name].append(name) return partable def list_pars(names_only=True, type=None): def list_pars(names_only=True, kind=None): """ Print all parameters in all models. occurs in. """ partable = find_pars(type) partable = find_pars(kind) if names_only: print(columnize(list(sorted(partable.keys())))) help="list models which use this argument") parser.add_argument( 'type', default="any", nargs='?', 'kind', default="any", nargs='?', metavar="volume|orientation|sld|none|any", choices=['volume', 'orientation', 'sld', None, 'any'], type=lambda v: None if v == 'any' else '' if v == 'none' else v, help="only list arguments of the given type") type=lambda v: None if v == 'any' else None if v == 'none' else v, help="only list arguments of the given kind") args = parser.parse_args() list_pars(names_only=not args.verbose, type=args.type) list_pars(names_only=not args.verbose, kind=args.kind) if __name__ == "__main__": • sasmodels/mixture.py  r39a06c9 def random(): """Random set of model parameters for mixture model""" combined_pars = {} for k, part in enumerate(parts): class MixtureModel(KernelModel): """ Model definition for mixture of models. """ def __init__(self, model_info, parts): # type: (ModelInfo, List[KernelModel]) -> None kernels = [part.make_kernel(q_vectors) for part in self.parts] return MixtureKernel(self.info, kernels) make_kernel.__doc__ = KernelModel.make_kernel.__doc__ def release(self): # type: () -> None """ Free resources associated with the model. """ """Free resources associated with the model.""" for part in self.parts: part.release() release.__doc__ = KernelModel.release.__doc__ class MixtureKernel(Kernel): """ Instantiated kernel for mixture of models. """ def __init__(self, model_info, kernels): # type: (ModelInfo, List[Kernel]) -> None self.results = [] # type: List[np.ndarray] def __call__(self, call_details, values, cutoff, magnetic): def Iq(self, call_details, values, cutoff, magnetic): # type: (CallDetails, np.ndarray, np.ndarry, float, bool) -> np.ndarray scale, background = values[0:2] # remember the parts for plotting later self.results = [] # type: List[np.ndarray] parts = MixtureParts(self.info, self.kernels, call_details, values) parts = _MixtureParts(self.info, self.kernels, call_details, values) for kernel, kernel_details, kernel_values in parts: #print("calling kernel", kernel.info.name) return scale*total + background Iq.__doc__ = Kernel.Iq.__doc__ __call__ = Iq def release(self): # type: () -> None """Free resources associated with the kernel.""" for k in self.kernels: k.release() class MixtureParts(object): # Note: _MixtureParts doesn't implement iteration correctly, and only allows # a single iterator to be active at once. It doesn't matter in this case # since _MixtureParts is only used in one place, but it is not clean style. class _MixtureParts(object): """ Mixture component iterator. """ def __init__(self, model_info, kernels, call_details, values): # type: (ModelInfo, List[Kernel], CallDetails, np.ndarray) -> None self.values = values self.spin_index = model_info.parameters.npars + 2 # The following are redefined by __iter__, but set them here so that # lint complains a little less. self.part_num = -1 self.par_index = -1 self.mag_index = -1 #call_details.show(values) • sasmodels/model_test.py  rd92182f import numpy as np # type: ignore from . import core from .core import list_models, load_model_info, build_model from .direct_model import call_kernel, call_Fq test_method_name = "test_%s_python" % model_info.id test = ModelTestCase(test_name, model_info, test_method_name, platform="dll", # so that dtype="double", stash=stash) test_method_name, platform="dll", # so that dtype="double", stash=stash) suite.addTest(test) else: # kernel implemented in C test_method_name = "test_%s_dll" % model_info.id test = ModelTestCase(test_name, model_info, test_method_name, platform="dll", dtype="double", stash=stash) test_method_name, platform="dll", dtype="double", stash=stash) suite.addTest(test) # presence of *single=False* in the model file. test = ModelTestCase(test_name, model_info, test_method_name, platform="ocl", dtype=None, stash=stash) test_method_name, platform="ocl", dtype=None, stash=stash) #print("defining", test_name) suite.addTest(test) # presence of *single=False* in the model file. test = ModelTestCase(test_name, model_info, test_method_name, platform="cuda", dtype=None, stash=stash) test_method_name, platform="cuda", dtype=None, stash=stash) #print("defining", test_name) suite.addTest(test) s_name = pars.pop('@S') ps_test = [pars] + list(test[1:]) #print("PS TEST PARAMS!!!",ps_test) # build the P@S model s_info = load_model_info(s_name) platform=self.platform) # run the tests #self.info = ps_model.info #print("SELF.INFO PARAMS!!!",[p.id for p in self.info.parameters.call_parameters]) #print("PS MODEL PARAMETERS:",[p.id for p in ps_model.info.parameters.call_parameters]) results.append(self.run_one(ps_model, ps_test)) """Run a single test case.""" user_pars, x, y = test[:3] pars = expand_pars(self.info.parameters, user_pars) invalid = invalid_pars(self.info.parameters, pars) #print("PS MODEL PARAMETERS:",[p.id for p in model.info.parameters.call_parameters]) pars = expand_pars(model.info.parameters, user_pars) invalid = invalid_pars(model.info.parameters, pars) if invalid: raise ValueError("Unknown parameters in test: " + ", ".join(invalid)) else: y1 = y y2 = test[3] if not isinstance(test[3], list) else [test[3]] F1, F2, R_eff, volume, volume_ratio = call_Fq(kernel, pars) if F1 is not None: # F1 is none for models with Iq instead of Fq self._check_vectors(x, y1, F1, 'F') self._check_vectors(x, y2, F2, 'F^2') y2 = test[3] if isinstance(test[3], list) else [test[3]] F, Fsq, R_eff, volume, volume_ratio = call_Fq(kernel, pars) if F is not None: # F is none for models with Iq instead of Fq self._check_vectors(x, y1, F, 'F') self._check_vectors(x, y2, Fsq, 'F^2') self._check_scalar(test[4], R_eff, 'R_eff') self._check_scalar(test[5], volume, 'volume') self._check_scalar(test[6], volume_ratio, 'form:shell ratio') return F2 return Fsq def _check_scalar(self, target, actual, name): if target is None: # smoke test --- make sure it runs and produces a value self.assertTrue(not np.isnan(actual), 'invalid %s: %s' % (name, actual)) elif np.isnan(target): # make sure nans match self.assertTrue(np.isnan(actual), '%s: expected:%s; actual:%s' % (name, target, actual)) else: # is_near does not work for infinite values, so also test # for exact values. self.assertTrue(target == actual or is_near(target, actual, 5), '%s: expected:%s; actual:%s' % (name, target, actual)) self.assertTrue(is_near(target, actual, 5), '%s: expected:%s; actual:%s' % (name, target, actual)) def _check_vectors(self, x, target, actual, name='I'): '%s(...) returned wrong length'%name) for xi, yi, actual_yi in zip(x, target, actual): if yi is None: # smoke test --- make sure it runs and produces a value self.assertTrue(not np.isnan(actual_yi), 'invalid %s(%s): %s' % (name, xi, actual_yi)) elif np.isnan(yi): # make sure nans match self.assertTrue(np.isnan(actual_yi), '%s(%s): expected:%s; actual:%s' % (name, xi, yi, actual_yi)) else: # is_near does not work for infinite values, so also test # for exact values. self.assertTrue(yi == actual_yi or is_near(yi, actual_yi, 5), '%s(%s); expected:%s; actual:%s' % (name, xi, yi, actual_yi)) self.assertTrue(is_near(yi, actual_yi, 5), '%s(%s): expected:%s; actual:%s' % (name, xi, target, actual)) return ModelTestCase invalid = [] for par in sorted(pars.keys()): # special handling of R_eff mode, which is not a usual parameter # Ignore the R_eff mode parameter when checking for valid parameters. # It is an allowed parameter for a model even though it does not exist # in the parameter table. The call_Fq() function pops it from the # parameter list and sends it directly to kernel.Fq(). if par == product.RADIUS_MODE_ID: continue """ Returns true if *actual* is within *digits* significant digits of *target*. """ import math shift = 10**math.ceil(math.log10(abs(target))) return abs(target-actual)/shift < 1.5*10**-digits *taget* zero and inf should match *actual* zero and inf. If you want to accept eps for zero, choose a value such as 1e-10, which must match up to +/- 1e-15 when *digits* is the default value of 5. If *target* is None, then just make sure that *actual* is not NaN. If *target* is NaN, make sure *actual* is NaN. """ if target is None: # target is None => actual cannot be NaN return not np.isnan(actual) elif target == 0.: # target is 0. => actual must be 0. # Note: if small values are allowed, then use maybe test zero against eps instead? return actual == 0. elif np.isfinite(target): shift = np.ceil(np.log10(abs(target))) return abs(target-actual) < 1.5*10**(shift-digits) elif target == actual: # target is inf => actual must be inf of same sign return True else: # target is NaN => actual must be NaN return np.isnan(target) == np.isnan(actual) # CRUFT: old interface; should be deprecated and removed def run_one(model_name): # type: (str) -> str """ [Deprecated] Run the tests associated with *model_name*. Use the following instead:: succss, output = check_model(load_model_info(model_name)) """ # msg = "use check_model(model_info) rather than run_one(model_name)" # warnings.warn(msg, category=DeprecationWarning, stacklevel=2) return output success, output = check_model(model_info) _, output = check_model(model_info) return output def check_model(model_info): # type: (ModelInfo) -> str # type: (ModelInfo) -> Tuple[bool, str] """ Run the tests for a single model, capturing the output. loaders.append('cuda') tests = make_suite(loaders, ['all']) def build_test(test): def _build_test(test): # In order for nosetest to show the test name, wrap the test.run_all # instance in function that takes the test name as a parameter which for test in tests: yield build_test(test) yield _build_test(test) "models will be tested. See core.list_models() for " "names of other groups, such as 'py' or 'single'.") args, models = parser.parse_known_args() if args.engine == "opencl": opts = parser.parse_args() if opts.engine == "opencl": if not use_opencl(): print("opencl is not available") return 1 loaders = ['opencl'] elif args.engine == "dll": elif opts.engine == "dll": loaders = ["dll"] elif args.engine == "cuda": elif opts.engine == "cuda": if not use_cuda(): print("cuda is not available") return 1 loaders = ['cuda'] elif args.engine == "all": elif opts.engine == "all": loaders = ['dll'] if use_opencl(): loaders.append('cuda') else: print("unknown engine " + args.engine) print("unknown engine " + opts.engine) return 1 runner = TestRunner(verbosity=args.verbose, **test_args) result = runner.run(make_suite(loaders, args.models)) runner = TestRunner(verbosity=opts.verbose, **test_args) result = runner.run(make_suite(loaders, opts.models)) return 1 if result.failures or result.errors else 0 • sasmodels/models/_spherepy.py  r304c775 to the output of the software provided by the NIST (Kline, 2006). References ---------- John Wiley and Sons, New York, (1955) Source ------ _spherepy.py _ sphere.c _ Authorship and Verification ---------------------------- * **Author: P Kienzle** * **Last Modified by:** * **Last Reviewed by:** S King and P Parker **Date:** 2013/09/09 and 2014/01/06 * **Source added by :** Steve King **Date:** March 25, 2019 """ def form_volume(radius): """Calculate volume for sphere""" return 1.333333333333333 * pi * radius ** 3 def effective_radius(mode, radius): return radius def radius_effective(mode, radius): """Calculate R_eff for sphere""" return radius if mode else 0. def Iq(q, sld, sld_solvent, radius): """Calculate I(q) for sphere""" #print "q",q #print "sld,r",sld,sld_solvent,radius • sasmodels/models/adsorbed_layer.py  r2d81cfe Layers*, *Macromol. Symp.*, 190 (2002) 33-42. Source ------ adsorbed_layer.py _ Authorship and Verification ---------------------------- * **Last Modified by:** Paul Kienzle **Date:** April 14, 2016 * **Last Reviewed by:** Steve King **Date:** March 18, 2016 * **Source added by :** Steve King **Date:** March 25, 2019 """ def Iq(q, second_moment, adsorbed_amount, density_shell, radius, volfraction, sld_shell, sld_solvent): """Return I(q) for adsorbed layer model.""" with errstate(divide='ignore'): aa = ((sld_shell - sld_solvent)/density_shell * adsorbed_amount) / q def random(): """Return a random parameter set for the model.""" # only care about the value of second_moment: # curve = scale * e**(-second_moment^2 q^2)/q^2 • sasmodels/models/barbell.c  r99658f6 static double effective_radius(int mode, double radius_bell, double radius, double length) radius_effective(int mode, double radius_bell, double radius, double length) { switch (mode) { • sasmodels/models/barbell.py  r99658f6 ---------- .. [#] H Kaya, *J. Appl. Cryst.*, 37 (2004) 37 223-230 .. [#] H Kaya, *J. Appl. Cryst.*, 37 (2004) 223-230 .. [#] H Kaya and N R deSouza, *J. Appl. Cryst.*, 37 (2004) 508-509 (addenda and errata) L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949). .. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659 Source ------ barbell.py _ barbell.c _ Authorship and Verification * **Last Modified by:** Paul Butler **Date:** March 20, 2016 * **Last Reviewed by:** Richard Heenan **Date:** January 4, 2017 * **Source added by :** Steve King **Date:** March 25, 2019 """ source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "barbell.c"] have_Fq = True effective_radius_type = [ "equivalent cylinder excluded volume","equivalent volume sphere", "radius", "half length", "half total length", have_Fq = True radius_effective_modes = [ "equivalent cylinder excluded volume", "equivalent volume sphere", "radius", "half length", "half total length", ] def random(): """Return a random parameter set for the model.""" # TODO: increase volume range once problem with bell radius is fixed # The issue is that bell radii of more than about 200 fail at high q • sasmodels/models/bcc_paracrystal.py  rda7b26b r""" .. warning:: This model and this model description are under review following concerns raised by SasView users. If you need to use this model, please email help@sasview.org for the latest situation. *The .. warning:: This model and this model description are under review following concerns raised by SasView users. If you need to use this model, please email help@sasview.org for the latest situation. *The SasView Developers. September 2018.* (Corrections to FCC and BCC lattice structure calculation) Source ------ bcc_paracrystal.py _ bcc_paracrystal.c _ Authorship and Verification --------------------------- * **Last Modified by:** Paul Butler **Date:** September 29, 2016 * **Last Reviewed by:** Richard Heenan **Date:** March 21, 2016 * **Source added by :** Steve King **Date:** March 25, 2019 """ def random(): """Return a random parameter set for the model.""" # Define lattice spacing as a multiple of the particle radius # using the formulat a = 4 r/sqrt(3). Systems which are ordered • sasmodels/models/be_polyelectrolyte.py  rca77fc1 constituting the polymer monomer and the solvent molecules, respectively. -$v_p$and$v_s$are the partial molar volume of the polymer and the -$v_p$and$v_s$are the partial molar volume of the polymer and the solvent, respectively. -$C_s$is the concentration of monovalent salt(1/|Ang^3| - internally converted from mol/L). -$\alpha$is the degree of ionization (the ratio of charged monomers to the total -$\alpha$is the degree of ionization (the ratio of charged monomers to the total number of monomers) dimensionally useful units of 1/|Ang^3|, only the converted version of the polymer concentration was actually being used in the calculation while the unconverted salt concentration (still in apparent units of mol/L) was being used. This was carried through to Sasmodels as used for SasView 4.1 (though the line of code converting the salt concentration to the new units was removed somewhere along the line). Simple dimensional analysis of the calculation shows that the converted salt concentration should be used, which the original code suggests was the intention, so this has now been corrected (for SasView 4.2). unconverted salt concentration (still in apparent units of mol/L) was being used. This was carried through to Sasmodels as used for SasView 4.1 (though the line of code converting the salt concentration to the new units was removed somewhere along the line). Simple dimensional analysis of the calculation shows that the converted salt concentration should be used, which the original code suggests was the intention, so this has now been corrected (for SasView 4.2). Once better validation has been performed this note will be removed. .. [#] E Raphael, J F Joanny, *Europhysics Letters*, 11 (1990) 179 Source ------ be_polyelectrolyte.py _ Authorship and Verification ---------------------------- * **Last Modified by:** Paul Butler **Date:** September 25, 2018 * **Last Reviewed by:** Paul Butler **Date:** September 25, 2018 * **Source added by :** Steve King **Date:** March 25, 2019 """ :params: see parameter table :return: 1-D form factor for polyelectrolytes in low salt parameter names, units, default values, and behavior (volume, sld etc) are defined in the parameter table. The concentrations are converted from def random(): """Return a random parameter set for the model.""" # TODO: review random be_polyelectrolyte model generation pars = dict( • sasmodels/models/binary_hard_sphere.py  r2d81cfe .. [#] S R Kline, *J Appl. Cryst.*, 39 (2006) 895 Source ------ binary_hard_sphere.py _ binary_hard_sphere.c _ Authorship and Verification ---------------------------- * **Last Modified by:** Paul Butler **Date:** March 20, 2016 * **Last Reviewed by:** Paul Butler **Date:** March 20, 2016 * **Source added by :** Steve King **Date:** March 25, 2019 """ def random(): """Return a random parameter set for the model.""" # TODO: binary_hard_sphere fails at low qr radius_lg = 10**np.random.uniform(2, 4.7) • sasmodels/models/broad_peak.py  r2d81cfe None. Source ------ broad_peak.py _ Authorship and Verification ---------------------------- * **Last Modified by:** Paul kienle **Date:** July 24, 2016 * **Last Reviewed by:** Richard Heenan **Date:** March 21, 2016 * **Source added by :** Steve King **Date:** March 25, 2019 """ def random(): """Return a random parameter set for the model.""" pars = dict( scale=1, • sasmodels/models/capped_cylinder.c  r99658f6 static double effective_radius(int mode, double radius, double radius_cap, double length) radius_effective(int mode, double radius, double radius_cap, double length) { switch (mode) { • sasmodels/models/capped_cylinder.py  r99658f6 .. [#] H Kaya and N-R deSouza, *J. Appl. Cryst.*, 37 (2004) 508-509 (addenda and errata) L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949). .. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659 Source ------ capped_cylinder.py _ capped_cylinder.c _ Authorship and Verification * **Last Modified by:** Paul Butler **Date:** September 30, 2016 * **Last Reviewed by:** Richard Heenan **Date:** January 4, 2017 * **Source added by :** Steve King **Date:** March 25, 2019 """ source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "capped_cylinder.c"] have_Fq = True effective_radius_type = [ "equivalent cylinder excluded volume", "equivalent volume sphere", "radius", "half length", "half total length", radius_effective_modes = [ "equivalent cylinder excluded volume", "equivalent volume sphere", "radius", "half length", "half total length", ] def random(): """Return a random parameter set for the model.""" # TODO: increase volume range once problem with bell radius is fixed # The issue is that bell radii of more than about 200 fail at high q • sasmodels/models/core_multi_shell.c  rd42dd4a static double effective_radius(int mode, double core_radius, double fp_n, double thickness[]) radius_effective(int mode, double core_radius, double fp_n, double thickness[]) { switch (mode) { • sasmodels/models/core_multi_shell.py  ree60aa7 Neutron Scattering*, Plenum Press, New York, 1987. Source ------ core_multi_shell.py _ core_multi_shell.c _ Authorship and Verification ---------------------------- * **Last Modified by:** Paul Kienzle **Date:** September 12, 2016 * **Last Reviewed by:** Paul Kienzle **Date:** September 12, 2016 * **Source added by :** Steve King **Date:** March 25, 2019 """ from __future__ import division source = ["lib/sas_3j1x_x.c", "core_multi_shell.c"] have_Fq = True effective_radius_type = ["outer radius", "core radius"] radius_effective_modes = ["outer radius", "core radius"] def random(): """Return a random parameter set for the model.""" num_shells = np.minimum(np.random.poisson(3)+1, 10) total_radius = 10**np.random.uniform(1.7, 4) • sasmodels/models/core_shell_bicelle.c  r99658f6 static double effective_radius(int mode, double radius, double thick_rim, double thick_face, double length) radius_effective(int mode, double radius, double thick_rim, double thick_face, double length) { switch (mode) { • sasmodels/models/core_shell_bicelle.py  r99658f6 from Proquest _ L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949). .. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659 Source ------ core_shell_bicelle.py _ core_shell_bicelle.c _ Authorship and Verification * **Last Modified by:** Paul Butler **Date:** September 30, 2016 * **Last Reviewed by:** Richard Heenan **Date:** January 4, 2017 * **Source added by :** Steve King **Date:** March 25, 2019 """ "core_shell_bicelle.c"] have_Fq = True effective_radius_type = [ "excluded volume","equivalent volume sphere", "outer rim radius", radius_effective_modes = [ "excluded volume", "equivalent volume sphere", "outer rim radius", "half outer thickness", "half diagonal", ] def random(): """Return a random parameter set for the model.""" pars = dict( radius=10**np.random.uniform(1.3, 3), • sasmodels/models/core_shell_bicelle_elliptical.c  r99658f6 static double effective_radius(int mode, double r_minor, double x_core, double thick_rim, double thick_face, double length) radius_effective(int mode, double r_minor, double x_core, double thick_rim, double thick_face, double length) { switch (mode) { • sasmodels/models/core_shell_bicelle_elliptical.py  r99658f6 ---------- .. [#] L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949). .. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659 Source ------ core_shell_bicelle_elliptical.py _ core_shell_bicelle_elliptical.c _ Authorship and Verification * **Last Modified by:** Richard Heenan **Date:** December 14, 2016 * **Last Reviewed by:** Paul Kienzle **Date:** Feb 28, 2018 * **Source added by :** Steve King **Date:** March 25, 2019 """ "core_shell_bicelle_elliptical.c"] have_Fq = True effective_radius_type = [ "equivalent cylinder excluded volume", "equivalent volume sphere", "outer rim average radius", "outer rim min radius", radius_effective_modes = [ "equivalent cylinder excluded volume", "equivalent volume sphere", "outer rim average radius", "outer rim min radius", "outer max radius", "half outer thickness", "half diagonal", ] def random(): """Return a random parameter set for the model.""" outer_major = 10**np.random.uniform(1, 4.7) outer_minor = 10**np.random.uniform(1, 4.7) • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c  r99658f6 static double effective_radius(int mode, double r_minor, double x_core, double thick_rim, double thick_face, double length) radius_effective(int mode, double r_minor, double x_core, double thick_rim, double thick_face, double length) { switch (mode) { • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py  r99658f6 ---------- .. [#] L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949). .. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659 Source ------ core_shell_bicelle_elliptical_belt_rough.py _ core_shell_bicelle_elliptical_belt_rough.c _ Authorship and Verification * **Last Modified by:** Richard Heenan new 2d orientation **Date:** October 5, 2017 * **Last Reviewed by:** Richard Heenan 2d calc seems agree with 1d **Date:** Nov 2, 2017 * **Source added by :** Steve King **Date:** March 25, 2019 """ "core_shell_bicelle_elliptical_belt_rough.c"] have_Fq = True effective_radius_type = [ "equivalent cylinder excluded volume", "equivalent volume sphere", "outer rim average radius", "outer rim min radius", radius_effective_modes = [ "equivalent cylinder excluded volume", "equivalent volume sphere", "outer rim average radius", "outer rim min radius", "outer max radius", "half outer thickness", "half diagonal", ] # TODO: No random() for core-shell bicelle elliptical belt rough demo = dict(scale=1, background=0, #[{'radius': 30.0, 'x_core': 3.0, 'thick_rim':8.0, 'thick_face':14.0, 'length':50.0}, 'VR', 1], [{'radius': 30.0, 'x_core': 3.0, 'thick_rim':8.0, 'thick_face':14.0, 'length':50.0, 'sld_core':4.0, 'sld_face':7.0, 'sld_rim':1.0, 'sld_solvent':6.0, 'background':0.0}, [{'radius': 30.0, 'x_core': 3.0, 'thick_rim': 8.0, 'thick_face': 14.0, 'length': 50.0, 'sld_core': 4.0, 'sld_face': 7.0, 'sld_rim': 1.0, 'sld_solvent': 6.0, 'background': 0.0}, 0.015, 189.328], #[{'theta':80., 'phi':10.}, (qx, qy), 7.88866563001 ], • sasmodels/models/core_shell_cylinder.c  r99658f6 static double effective_radius(int mode, double radius, double thickness, double length) radius_effective(int mode, double radius, double thickness, double length) { switch (mode) { • sasmodels/models/core_shell_cylinder.py  r99658f6 density of the solvent, and *background* is the background level. The outer radius of the shell is given by$R+T$and the total length of the outer shell is given by$L+2T$.$J1$is the first order Bessel function. shell is given by$L+2T$.$J_1$is the first order Bessel function. .. _core-shell-cylinder-geometry: 1445-1452 .. [#kline] S R Kline, *J Appl. Cryst.*, 39 (2006) 895 L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949). .. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659 Source ------ core_shell_cylinder.py _ core_shell_cylinder.c _ Authorship and Verification * **Last Modified by:** Paul Kienzle **Date:** Aug 8, 2016 * **Last Reviewed by:** Richard Heenan **Date:** March 18, 2016 * **Source added by :** Steve King **Date:** March 25, 2019 """ source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "core_shell_cylinder.c"] have_Fq = True effective_radius_type = [ radius_effective_modes = [ "excluded volume", "equivalent volume sphere", "outer radius", "half outer length", "half min outer dimension", "half max outer dimension", "half outer diagonal", def random(): """Return a random parameter set for the model.""" outer_radius = 10**np.random.uniform(1, 4.7) # Use a distribution with a preference for thin shell or thin core • sasmodels/models/core_shell_ellipsoid.c  r99658f6 static double effective_radius(int mode, double radius_equat_core, double x_core, double thick_shell, double x_polar_shell) radius_effective(int mode, double radius_equat_core, double x_core, double thick_shell, double x_polar_shell) { const double radius_equat_tot = radius_equat_core + thick_shell; • sasmodels/models/core_shell_ellipsoid.py  r99658f6 ---------- Parameters for this model are the core axial ratio X and a shell thickness, which are more often what we would like to determine and makes the model better behaved, particularly when polydispersity is applied than the four independent radii used in the original parameterization of this model. Parameters for this model are the core axial ratio$X_{core}$and a shell thickness$t_{shell}$, which are more often what we would like to determine and make the model better behaved, particularly when polydispersity is applied, than the four independent radii used in the original parameterization of this model. the poles, of a prolate ellipsoid. When *X_core < 1* the core is oblate; when *X_core > 1* it is prolate. *X_core = 1* is a spherical core. For a fixed shell thickness *XpolarShell = 1*, to scale the shell thickness pro-rata with the radius set or constrain *XpolarShell = X_core*. When including an$S(q)$, the radius in$S(q)$is calculated to be that of a sphere with the same 2nd virial coefficient of the outer surface of the ellipsoid. This may have some undesirable effects if the aspect ratio of the ellipsoid is large (ie, if$X << 1$or$X >> 1$), when the$S(q)$- which assumes spheres - will not in any case be valid. Generating a custom product model will enable separate effective volume fraction and effective radius in the$S(q)$. When$X_{core}$< 1 the core is oblate; when$X_{core}$> 1 it is prolate.$X_{core}$= 1 is a spherical core. For a fixed shell thickness$X_{polar shell}$= 1, to scale$t_{shell}$pro-rata with the radius set or constrain$X_{polar shell}$=$X_{core}$. .. note:: When including an$S(q)$, the radius in$S(q)$is calculated to be that of a sphere with the same 2nd virial coefficient of the outer surface of the ellipsoid. This may have some undesirable effects if the aspect ratio of the ellipsoid is large (ie, if$X << 1$or$X >> 1$), when the$S(q)$- which assumes spheres - will not in any case be valid. Generating a custom product model will enable separate effective volume fraction and effective radius in the$S(q). If SAS data are in absolute units, and the SLDs are correct, then scale should where .. In following equation SK changed radius\_equat\_core to R_e .. math:: :nowrap: \begin{align*} F(q,\alpha) = &f(q,radius\_equat\_core,radius\_equat\_core.x\_core,\alpha) \\ &+ f(q,radius\_equat\_core + thick\_shell, radius\_equat\_core.x\_core + thick\_shell.x\_polar\_shell,\alpha) F(q,\alpha) = &f(q,R_e,R_e.x_{core},\alpha) \\ &+ f(q,R_e + t_{shell}, R_e.x_{core} + t_{shell}.x_{polar shell},\alpha) \end{align*}V = (4/3)\pi R_pR_e^2$is the volume of the ellipsoid ,$R_p$is the polar radius along the rotational axis of the ellipsoid,$R_e$is the equatorial radius perpendicular to the rotational axis of the ellipsoid and$\Delta \rho$(contrast) is the scattering length density difference, either$(sld\_core - sld\_shell)$or$(sld\_shell - sld\_solvent)$. equatorial radius perpendicular to the rotational axis of the ellipsoid,$t_{shell}$is the thickness of the shell at the equator, and$\Delta \rho$(the contrast) is the scattering length density difference, either$(\rho_{core} - \rho_{shell})$or$(\rho_{shell} - \rho_{solvent})$. For randomly oriented particles: ---------- see for example: Kotlarchyk, M.; Chen, S.-H. J. Chem. Phys., 1983, 79, 2461. Berr, S. J. Phys. Chem., 1987, 91, 4760. .. [#] Kotlarchyk, M.; Chen, S.-H. *J. Chem. Phys.*, 1983, 79, 2461 .. [#] Berr, S. *J. Phys. Chem.*, 1987, 91, 4760 Source ------ core_shell_ellipsoid.py _ core_shell_ellipsoid.c _ Authorship and Verification * **Author:** NIST IGOR/DANSE **Date:** pre 2010 * **Last Modified by:** Richard Heenan (reparametrised model) **Date:** 2015 * **Last Reviewed by:** Richard Heenan **Date:** October 6, 2016 * **Last Reviewed by:** Steve King **Date:** March 27, 2019 * **Source added by :** Steve King **Date:** March 25, 2019 """ source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "core_shell_ellipsoid.c"] have_Fq = True effective_radius_type = [ "average outer curvature", "equivalent volume sphere", radius_effective_modes = [ "average outer curvature", "equivalent volume sphere", "min outer radius", "max outer radius", ] def random(): """Return a random parameter set for the model.""" volume = 10**np.random.uniform(5, 12) outer_polar = 10**np.random.uniform(1.3, 4) • sasmodels/models/core_shell_parallelepiped.c  r99658f6 static double effective_radius(int mode, double length_a, double length_b, double length_c, radius_effective(int mode, double length_a, double length_b, double length_c, double thick_rim_a, double thick_rim_b, double thick_rim_c) { • sasmodels/models/core_shell_parallelepiped.py  r99658f6 from Proquest _ L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949). .. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659 Source ------ core_shell_parallelepiped.py _ core_shell_parallelepiped.c _ Authorship and Verification * **Last Reviewed by:** Paul Butler **Date:** May 24, 2018 - documentation updated * **Source added by :** Steve King **Date:** March 25, 2019 """ import numpy as np from numpy import pi, inf, sqrt, cos, sin from numpy import inf name = "core_shell_parallelepiped" source = ["lib/gauss76.c", "core_shell_parallelepiped.c"] have_Fq = True effective_radius_type = [ "equivalent cylinder excluded volume", radius_effective_modes = [ "equivalent cylinder excluded volume", "equivalent volume sphere", "half outer length_a", "half outer length_b", "half outer length_c", def random(): """Return a random parameter set for the model.""" outer = 10**np.random.uniform(1, 4.7, size=3) thick = np.random.beta(0.5, 0.5, size=3)*(outer-2) + 1 # 2d values not tested against other codes or models if 0: # pak: model rewrite; need to update tests qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) qx, qy = 0.2 * np.cos(np.pi/6.), 0.2 * np.sin(np.pi/6.) tests = [[{}, 0.2, 0.533149288477], [{}, [0.2], [0.533149288477]], • sasmodels/models/core_shell_sphere.c  rd42dd4a static double effective_radius(int mode, double radius, double thickness) radius_effective(int mode, double radius, double thickness) { switch (mode) { • sasmodels/models/core_shell_sphere.py  r304c775 effective radius for$S(Q)$when$P(Q) \cdot S(Q)$is applied. References ---------- A Guinier and G Fournet, *Small-Angle Scattering of X-Rays*, John Wiley and Sons, New York, (1955) Validation ---------- the output of the software provided by NIST (Kline, 2006). Figure 1 shows a comparison of the output of our model and the output of the NIST software. References ---------- .. [#] A Guinier and G Fournet, *Small-Angle Scattering of X-Rays*, John Wiley and Sons, New York, (1955) Source ------ core_shell_sphere.py _ core_shell_sphere.c _ Authorship and Verification ---------------------------- * **Author:** * **Last Modified by:** * **Last Reviewed by:** * **Source added by :** Steve King **Date:** March 25, 2019 """ source = ["lib/sas_3j1x_x.c", "lib/core_shell.c", "core_shell_sphere.c"] have_Fq = True effective_radius_type = ["outer radius", "core radius"] radius_effective_modes = ["outer radius", "core radius"] demo = dict(scale=1, background=0, radius=60, thickness=10, def random(): """Return a random parameter set for the model.""" outer_radius = 10**np.random.uniform(1.3, 4.3) # Use a distribution with a preference for thin shell or thin core • sasmodels/models/correlation_length.py  ra807206 ---------- B Hammouda, D L Ho and S R Kline, Insight into Clustering in Poly(ethylene oxide) Solutions, Macromolecules, 37 (2004) 6932-6937 .. [#] B Hammouda, D L Ho and S R Kline, Insight into Clustering in Poly(ethylene oxide) Solutions, Macromolecules, 37 (2004) 6932-6937 Source ------ correlation_length.py _ Authorship and Verification ---------------------------- * **Author:** * **Last Modified by:** * **Last Reviewed by:** * **Source added by :** Steve King **Date:** March 25, 2019 """ • sasmodels/models/cylinder.c  r99658f6 static double effective_radius(int mode, double radius, double length) radius_effective(int mode, double radius, double length) { switch (mode) { • sasmodels/models/cylinder.py  r99658f6 ---------- J. S. Pedersen, Adv. Colloid Interface Sci. 70, 171-210 (1997). G. Fournet, Bull. Soc. Fr. Mineral. Cristallogr. 74, 39-113 (1951). L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949). .. [#] J. Pedersen, *Adv. Colloid Interface Sci.*, 70 (1997) 171-210 .. [#] G. Fournet, *Bull. Soc. Fr. Mineral. Cristallogr.*, 74 (1951) 39-113 .. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659 Source ------ cylinder.py _ cylinder.c _ Authorship and Verification ---------------------------- * **Author:** * **Last Modified by:** * **Last Reviewed by:** * **Source added by :** Steve King **Date:** March 25, 2019 """ source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "cylinder.c"] have_Fq = True effective_radius_type = [ radius_effective_modes = [ "excluded volume", "equivalent volume sphere", "radius", "half length", "half min dimension", "half max dimension", "half diagonal", def random(): """Return a random parameter set for the model.""" volume = 10**np.random.uniform(5, 12) length = 10**np.random.uniform(-2, 2)*volume**0.333 phi_pd=10, phi_pd_n=5) # pylint: disable=bad-whitespace, line-too-long # Test 1-D and 2-D models qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) # After redefinition of angles, find new tests values. Was 10 10 in old coords theta, phi = 80.1534480601659, 10.1510817110481 # (10, 10) in sasview 3.x tests = [ [{}, 0.2, 0.042761386790780453], [{}, [0.2], [0.042761386790780453]], # new coords [{'theta':80.1534480601659, 'phi':10.1510817110481}, (qx, qy), 0.03514647218513852], [{'theta':80.1534480601659, 'phi':10.1510817110481}, [(qx, qy)], [0.03514647218513852]], # old coords #[{'theta':10.0, 'phi':10.0}, (qx, qy), 0.03514647218513852], #[{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.03514647218513852]], [{'theta': theta, 'phi': phi}, (qx, qy), 0.03514647218513852], [{'theta': theta, 'phi': phi}, [(qx, qy)], [0.03514647218513852]], ] del qx, qy # not necessary to delete, but cleaner # Default radius and length radius, length = parameters[2][2], parameters[3][2] tests.extend([ ({'radius_effective_mode': 0}, 0.1, None, None, 0., pi*radius**2*length, 1.0), ({'radius_effective_mode': 1}, 0.1, None, None, 0.5*(0.75*radius*(2.0*radius*length + (radius + length)*(pi*radius + length)))**(1./3.), None, None), ({'radius_effective_mode': 2}, 0.1, None, None, (0.75*radius**2*length)**(1./3.), None, None), ({'radius_effective_mode': 3}, 0.1, None, None, radius, None, None), ({'radius_effective_mode': 4}, 0.1, None, None, length/2., None, None), ({'radius_effective_mode': 5}, 0.1, None, None, min(radius, length/2.), None, None), ({'radius_effective_mode': 6}, 0.1, None, None, max(radius, length/2.), None, None), ({'radius_effective_mode': 7}, 0.1, None, None, np.sqrt(4*radius**2 + length**2)/2., None, None), ]) del radius, length # pylint: enable=bad-whitespace, line-too-long del qx, qy, theta, phi # not necessary to delete, but cleaner def _extend_with_reff_tests(radius, length): """Test R_eff and form volume calculations""" # V and Vr are the same for each R_eff mode V = pi*radius**2*length # shell volume = form volume for solid objects Vr = 1.0 # form:shell volume ratio # Use test value for I(0.2) from above to check Fsq value. Need to # remove scale and background before testing. q = 0.2 scale, background = V, 0.001 Fsq = (0.042761386790780453 - background)*scale F = None # Need target value for # Various values for R_eff, depending on mode r_effs = [ 0., 0.5*(0.75*radius*(2.0*radius*length + (radius + length)*(pi*radius + length)))**(1./3.), (0.75*radius**2*length)**(1./3.), radius, length/2., min(radius, length/2.), max(radius, length/2.), np.sqrt(4*radius**2 + length**2)/2., ] tests.extend([ ({'radius_effective_mode': 0}, q, F, Fsq, r_effs[0], V, Vr), ({'radius_effective_mode': 1}, q, F, Fsq, r_effs[1], V, Vr), ({'radius_effective_mode': 2}, q, F, Fsq, r_effs[2], V, Vr), ({'radius_effective_mode': 3}, q, F, Fsq, r_effs[3], V, Vr), ({'radius_effective_mode': 4}, q, F, Fsq, r_effs[4], V, Vr), ({'radius_effective_mode': 5}, q, F, Fsq, r_effs[5], V, Vr), ({'radius_effective_mode': 6}, q, F, Fsq, r_effs[6], V, Vr), ({'radius_effective_mode': 7}, q, F, Fsq, r_effs[7], V, Vr), ]) # Test Reff and volume with default model parameters _extend_with_reff_tests(parameters[2][2], parameters[3][2]) del _extend_with_reff_tests # ADDED by: RKH ON: 18Mar2016 renamed sld's etc • sasmodels/models/dab.py  r304c775 ---------- P Debye, H R Anderson, H Brumberger, *Scattering by an Inhomogeneous Solid. II. The Correlation Function and its Application*, *J. Appl. Phys.*, 28(6) (1957) 679 .. [#] P Debye, H R Anderson, H Brumberger, *Scattering by an Inhomogeneous Solid. II. The Correlation Function and its Application*, *J. Appl. Phys.*, 28(6) (1957) 679 .. [#] P Debye, A M Bueche, *Scattering by an Inhomogeneous Solid*, *J. Appl. Phys.*, 20 (1949) 518 P Debye, A M Bueche, *Scattering by an Inhomogeneous Solid*, *J. Appl. Phys.*, 20 (1949) 518 Source ------ *2013/09/09 - Description reviewed by King, S and Parker, P.* dab.py _ Authorship and Verification ---------------------------- * **Author:** * **Last Modified by:** * **Last Reviewed by:** Steve King & Peter Parker **Date:** September 09, 2013 * **Source added by :** Steve King **Date:** March 25, 2019 """ def random(): """Return a random parameter set for the model.""" pars = dict( scale=10**np.random.uniform(1, 4), cor_length=10**np.random.uniform(0.3, 3), #background = 0, # background = 0, ) pars['scale'] /= pars['cor_length']**3 • sasmodels/models/ellipsoid.c  r99658f6 static double effective_radius(int mode, double radius_polar, double radius_equatorial) radius_effective(int mode, double radius_polar, double radius_equatorial) { switch (mode) { • sasmodels/models/ellipsoid.py  r99658f6 The$\theta$and$\phi$parameters are not used for the 1D output. Validation ---------- L A Feigin and D I Svergun. *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, Plenum Press, New York, 1987. A. Isihara. J. Chem. Phys. 18(1950) 1446-1449 .. [#] L A Feigin and D I Svergun. *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, Plenum Press, New York, 1987 .. [#] A. Isihara. *J. Chem. Phys.*, 18 (1950) 1446-1449 Source ------ ellipsoid.py _ ellipsoid.c _ Authorship and Verification * **Converted to sasmodels by:** Helen Park **Date:** July 9, 2014 * **Last Modified by:** Paul Kienzle **Date:** March 22, 2017 * **Source added by :** Steve King **Date:** March 25, 2019 """ from __future__ import division import numpy as np from numpy import inf, sin, cos, pi try: from numpy import cbrt except ImportError: def cbrt(x): return x ** (1.0/3.0) name = "ellipsoid" source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "ellipsoid.c"] have_Fq = True effective_radius_type = [ radius_effective_modes = [ "average curvature", "equivalent volume sphere", "min radius", "max radius", ] def random(): """Return a random parameter set for the model.""" volume = 10**np.random.uniform(5, 12) radius_polar = 10**np.random.uniform(1.3, 4) • sasmodels/models/elliptical_cylinder.c  r99658f6 static double effective_radius(int mode, double radius_minor, double r_ratio, double length) radius_effective(int mode, double radius_minor, double r_ratio, double length) { switch (mode) { • sasmodels/models/elliptical_cylinder.py  r99658f6 ---------- L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, Plenum, New York, (1987) [see table 3.4] L. Onsager, Ann. New York Acad. Sci. 51, 627-659 (1949). .. [#] L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, Plenum, New York, (1987) [see table 3.4] .. [#] L. Onsager, *Ann. New York Acad. Sci.*, 51 (1949) 627-659 Source ------ elliptical_cylinder.py _ elliptical_cylinder.c _ Authorship and Verification * **Last Modified by:** * **Last Reviewed by:** Richard Heenan - corrected equation in docs **Date:** December 21, 2016 * **Source added by :** Steve King **Date:** March 25, 2019 """ source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "elliptical_cylinder.c"] have_Fq = True effective_radius_type = [ "equivalent cylinder excluded volume", "equivalent volume sphere", "average radius", "min radius", "max radius", radius_effective_modes = [ "equivalent cylinder excluded volume", "equivalent volume sphere", "average radius", "min radius", "max radius", "equivalent circular cross-section", "half length", "half min dimension", "half max dimension", "half diagonal", def random(): """Return a random parameter set for the model.""" # V = pi * radius_major * radius_minor * length; volume = 10**np.random.uniform(3, 9) length = 10**np.random.uniform(1, 3) axis_ratio = 10**np.random.uniform(0, 2) radius_minor = np.sqrt(volume/length/axis_ratio) radius_minor = sqrt(volume/length/axis_ratio) volfrac = 10**np.random.uniform(-4, -2) pars = dict( tests = [ # [{'radius_minor': 20.0, 'axis_ratio': 1.5, 'length':400.0}, 'ER', 79.89245454155024], # [{'radius_minor': 20.0, 'axis_ratio': 1.2, 'length':300.0}, 'VR', 1], #[{'radius_minor': 20.0, 'axis_ratio': 1.5, 'length':400.0}, 'ER', 79.89245454155024], #[{'radius_minor': 20.0, 'axis_ratio': 1.2, 'length':300.0}, 'VR', 1], # The SasView test result was 0.00169, with a background of 0.001 • sasmodels/models/fcc_paracrystal.py  rda7b26b #note - calculation requires double precision r""" .. warning:: This model and this model description are under review following concerns raised by SasView users. If you need to use this model, please email help@sasview.org for the latest situation. *The .. warning:: This model and this model description are under review following concerns raised by SasView users. If you need to use this model, please email help@sasview.org for the latest situation. *The SasView Developers. September 2018.* (Corrections to FCC and BCC lattice structure calculation) Source ------ fcc_paracrystal.py _ fcc_paracrystal.c _ Authorship and Verification --------------------------- * **Last Modified by:** Paul Butler **Date:** September 29, 2016 * **Last Reviewed by:** Richard Heenan **Date:** March 21, 2016 * **Source added by :** Steve King **Date:** March 25, 2019 """ def random(): """Return a random parameter set for the model.""" # copied from bcc_paracrystal radius = 10**np.random.uniform(1.3, 4) • sasmodels/models/flexible_cylinder.py  r2d81cfe The Kuhn length$(b = 2*l_p)$is also used to describe the stiffness of a chain. The returned value is in units of$cm^{-1}$, on absolute scale. In the parameters, the sld and sld\_solvent represent the SLD of the cylinder and solvent respectively. Our model uses the form factor calculations implemented in a c-library provided by the NIST Center for Neutron Research (Kline, 2006). From the reference: Our model uses the form factor calculations in reference [1] as implemented in a c-library provided by the NIST Center for Neutron Research (Kline, 2006). This states: 'Method 3 With Excluded Volume' is used. See equations (13,26-27) in the original reference for the details. .. note:: There are several typos in the original reference that have been corrected by WRC [2]. Details of the corrections are in the reference below. Most notably - Equation (13): the term$(1 - w(QR))$should swap position with$w(QR)$- Equations (23) and (24) are incorrect; WRC has entered these into Mathematica and solved analytically. The results were then converted to code. - Equation (27) should be$q0 = max(a3/(Rg^2)^{1/2},3)$instead of$max(a3*b(Rg^2)^{1/2},3)$- The scattering function is negative for a range of parameter values and q-values that are experimentally accessible. A correction function has been added to give the proper behavior. **This is a model with complex behaviour depending on the ratio of**$L/b$**and the reader is strongly encouraged to read reference [1] before use.** .. note:: There are several typos in the original reference that have been corrected by WRC [2]. Details of the corrections are in the reference below. Most notably - Equation (13): the term$(1 - w(QR))$should swap position with$w(QR)$- Equations (23) and (24) are incorrect; WRC has entered these into Mathematica and solved analytically. The results were then converted to code. - Equation (27) should be$q0 = max(a3/(Rg^2)^{1/2},3)$instead of$max(a3*b(Rg^2)^{1/2},3)$- The scattering function is negative for a range of parameter values and q-values that are experimentally accessible. A correction function has been added to give the proper behavior. **This is a model with complex behaviour depending on the ratio of**$L/b$**and the reader is strongly encouraged to read reference [1] before use.** References ---------- J S Pedersen and P Schurtenberger. *Scattering functions of semiflexible polymers with and without excluded volume effects.* Macromolecules, 29 (1996) 7602-7612 .. [#] J S Pedersen and P Schurtenberger. *Scattering functions of semiflexible polymers with and without excluded volume effects.* Macromolecules, 29 (1996) 7602-7612 Correction of the formula can be found in W R Chen, P D Butler and L J Magid, *Incorporating Intermicellar Interactions in the Fitting of SANS Data from Cationic Wormlike Micelles.* Langmuir, 22(15) 2006 6539-6548 .. [#] W R Chen, P D Butler and L J Magid, *Incorporating Intermicellar Interactions in the Fitting of SANS Data from Cationic Wormlike Micelles.* Langmuir, 22(15) 2006 6539-6548 Source ------ flexible_cylinder.py _ flexible_cylinder.c _ wrc_cyl.c _ Authorship and Verification ---------------------------- * **Author:** * **Last Modified by:** * **Last Reviewed by:** Steve King **Date:** March 26, 2019 * **Source added by :** Steve King **Date:** March 25, 2019 """ name = "flexible_cylinder" title = "Flexible cylinder where the form factor is normalized by the volume" \ title = "Flexible cylinder where the form factor is normalized by the volume " \ "of the cylinder." description = """Note : scale and contrast = (sld - sld_solvent) are both def random(): """Return a random parameter set for the model.""" length = 10**np.random.uniform(2, 6) radius = 10**np.random.uniform(1, 3) • sasmodels/models/flexible_cylinder_elliptical.py  r2d81cfe The non-negligible diameter of the cylinder is included by accounting for excluded volume interactions within the walk of a single cylinder. **Inter-cylinder interactions are NOT provided for.** The form factor is normalized by the particle volume such that ----------- The function calculated in a similar way to that for the flexible_cylinder model from the reference given below using the author's "Method 3 With Excluded Volume". The function is calculated in a similar way to that for the :ref:flexible-cylinder model in reference [1] below using the author's "Method 3 With Excluded Volume". The model is a parameterization of simulations of a discrete representation of the worm-like chain model of Kratky and Porod applied in the pseudo-continuous There are several typos in the original reference that have been corrected by WRC. Details of the corrections are in the reference below. Most notably by WRC [2]. Details of the corrections are in the reference below. Most notably - Equation (13): the term$(1 - w(QR))$should swap position with$w(QR)$code. - Equation (27) should be$q0 = max(a3/sqrt(RgSquare),3)$instead of$max(a3*b/sqrt(RgSquare),3)$- Equation (27) should be$q0 = max(a3/(Rg^2)^{1/2},3)$instead of$max(a3*b(Rg^2)^{1/2},3)$- The scattering function is negative for a range of parameter values and The cross section of the cylinder is elliptical, with minor radius$a$. The major radius is larger, so of course, **the axis ratio (parameter 5) must be The major radius is larger, so of course, **the axis_ratio must be greater than one.** Simple constraints should be applied during curve fitting to maintain this inequality. The returned value is in units of$cm^{-1}$, on absolute scale. In the parameters, the$sld$and$sld\_solvent$represent the SLD of the these parameters must be held fixed during model fitting. **No inter-cylinder interference effects are included in this calculation.** **This is a model with complex behaviour depending on the ratio of**$L/b\$ **and the reader is strongly encouraged to read reference [1] before use.** References ---------- J S Pedersen and P Schurtenberger. *Scattering functions of semiflexible polymers with and without excluded volume effects.* Macromolecules, 29 (1996) 7602-7612 .. [#] J S Pedersen and P Schurtenberger. *Scattering functions of semiflexible polymers with and without excluded volume effects.* Macromolecules, 29 (1996) 7602-7612 Correction of the formula can be found in W R Chen, P D Butler and L J Magid, *Incorporating Intermicellar Interactions in the Fitting of SANS Data from Cationic Wormlike Micelles.* Langmuir, 22(15) 2006 6539-6548 .. [#] W R Chen, P D Butler and L J Magid, *Incorporating Intermicellar Interactions in the Fitting of SANS Data from Cationic Wormlike Micelles.* Langmuir, 22(15) 2006 6539-6548 Source ------ flexible_cylinder_elliptical.py _ flexible_cylinder_elliptical.c _ wrc_cyl.c _ Authorship and Verification ---------------------------- * **Author:** * **Last Modified by:** Richard Heenan **Date:** December, 2016 * **Last Reviewed by:** Steve King **Date:** March 26, 2019 * **Source added by :** Steve King **Date:** March 25, 2019 """ def random(): """Return a random parameter set for the model.""" length = 10**np.random.uniform(2, 6) radius = 10**np.random.uniform(1, 3)
 r2d81cfe .. [#] J Teixeira, *J. Appl. Cryst.*, 21 (1988) 781-785 Source ------ fractal.py _ fractal.c _ Authorship and Verification ---------------------------- * **Last Modified by:** Paul Butler **Date:** March 12, 2017 * **Last Reviewed by:** Paul Butler **Date:** March 12, 2017 * **Source added by :** Steve King **Date:** March 25, 2019 """ from __future__ import division def random(): """Return a random parameter set for the model.""" radius = 10**np.random.uniform(0.7, 4) #radius = 5
 r2cc8aa2 .. [#Kline]  S R Kline, *J Appl. Cryst.*, 39 (2006) 895 Source ------ fractal_core_shell.py _ fractal_core_shell.c _ Authorship and Verification ---------------------------- * **Last Modified by:** Paul Butler and Paul Kienzle **Date:** November 27, 2016 * **Last Reviewed by:** Paul Butler and Paul Kienzle **Date:** November 27, 2016 * **Source added by :** Steve King **Date:** March 25, 2019 """ import numpy as np from numpy import pi, inf from numpy import inf name = "fractal_core_shell" def random(): """Return a random parameter set for the model.""" outer_radius = 10**np.random.uniform(0.7, 4) # Use a distribution with a preference for thin shell or thin core cor_length=100.0) # TODO: why is there an ER function here? def ER(radius, thickness): """ #tests = [[{'radius': 20.0, 'thickness': 10.0}, 'ER', 30.0], tests = [ #         # At some point the SasView 3.x test result was deemed incorrect. The #following tests were verified against NIST IGOR macros ver 7.850. #NOTE: NIST macros do only provide for a polydispers core (no option #for a poly shell or for a monodisperse core.  The results seemed #extremely sensitive to the core PD, varying non monotonically all #the way to a PD of 1e-6. From 1e-6 to 1e-9 no changes in the #results were observed and the values below were taken using PD=1e-9. #Non-monotonically = I(0.001)=188 to 140 to 177 back to 160 etc. &n