Changeset e526a9d in sasmodels


Ignore:
Timestamp:
Feb 4, 2018 9:40:01 PM (7 years ago)
Author:
GitHub <noreply@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
aadec17
Parents:
032646d (diff), 72e41a0 (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 Butler <butlerpd@…> (02/04/18 21:40:01)
git-committer:
GitHub <noreply@…> (02/04/18 21:40:01)
Message:

Merge pull request #59 from SasView?/boltzmann

add BoltzmannDistribution? and UniformDistribution?

Since Paul K approved the code and the only issue found by R. Heenan and P. Butler has been resolved by sasview PR 141 we should be able to merge this now.

Files:
3 added
10 edited
2 moved

Legend:

Unmodified
Added
Removed

  • doc/guide/pd/polydispersity.rst

    reda8b30 r92d330fd  
    4242calculations are generally more robust with more data points or more angles. 
    4343 
    44 The following five distribution functions are provided: 
     44The following distribution functions are provided: 
    4545 
    4646*  *Rectangular Distribution* 
     47*  *Uniform Distribution* 
    4748*  *Gaussian Distribution* 
    4849*  *Lognormal Distribution* 
    4950*  *Schulz Distribution* 
    5051*  *Array Distribution* 
     52*  *Boltzmann Distribution* 
    5153 
    5254These are all implemented as *number-average* distributions. 
    … …  
    8587    Rectangular distribution. 
    8688 
     89 
     90 
     91Uniform Distribution 
     92^^^^^^^^^^^^^^^^^^^^^^^^ 
     93 
     94The Uniform Distribution is defined as 
     95 
     96    .. math:: 
     97 
     98        f(x) = \frac{1}{\text{Norm}} 
     99        \begin{cases} 
     100          1 & \text{for } |x - \bar x| \leq \sigma \\ 
     101          0 & \text{for } |x - \bar x| > \sigma 
     102        \end{cases} 
     103 
     104    where $\bar x$ is the mean of the distribution, $\sigma$ is the half-width, and 
     105    *Norm* is a normalization factor which is determined during the numerical 
     106    calculation. 
     107 
     108    Note that the polydispersity is given by 
     109 
     110    .. math:: \text{PD} = \sigma / \bar x 
     111 
     112    .. figure:: pd_uniform.jpg 
     113 
     114        Uniform distribution. 
     115 
     116The value $N_\sigma$ is ignored for this distribution. 
     117 
    87118.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
    88119 
    … …  
    183214^^^^^^^^^^^^^^^^^^ 
    184215 
    185 This user-definable distribution should be given as as a simple ASCII text 
     216This user-definable distribution should be given as a simple ASCII text 
    186217file where the array is defined by two columns of numbers: $x$ and $f(x)$. 
    187218The $f(x)$ will be normalized to 1 during the computation. 
    … …  
    202233given for the model will have no affect, and will be ignored when computing 
    203234the average.  This means that any parameter with an array distribution will 
    204 not be fittable. 
     235not be fitable. 
     236 
     237.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 
     238 
     239Boltzmann Distribution 
     240^^^^^^^^^^^^^^^^^^^^^^ 
     241 
     242The Boltzmann Distribution is defined as 
     243 
     244.. math:: 
     245 
     246    f(x) = \frac{1}{\text{Norm}} 
     247           \exp\left(-\frac{ | x - \bar x | }{\sigma}\right) 
     248 
     249where $\bar x$ is the mean of the distribution and *Norm* is a normalization 
     250factor which is determined during the numerical calculation. 
     251The width is defined as 
     252 
     253.. math:: \sigma=\frac{k T}{E} 
     254 
     255which is the inverse Boltzmann factor, 
     256where $k$ is the Boltzmann constant, $T$ the temperature in Kelvin and $E$ a 
     257characteristic energy per particle. 
     258 
     259.. figure:: pd_boltzmann.jpg 
     260 
     261    Boltzmann distribution. 
    205262 
    206263.. ZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZZ 

  • sasmodels/weights.py

    r2d81cfe r3d58247  
    9797        return x, px 
    9898 
     99class UniformDispersion(Dispersion): 
     100    r""" 
     101    Uniform dispersion, with width $\sigma$. 
     102 
     103    .. math:: 
     104 
     105        w = 1 
     106    """ 
     107    type = "uniform" 
     108    default = dict(npts=35, width=0, nsigmas=None) 
     109    def _weights(self, center, sigma, lb, ub): 
     110        x = np.linspace(center-sigma, center+sigma, self.npts) 
     111        x = x[(x >= lb) & (x <= ub)] 
     112        return x, np.ones_like(x) 
    99113 
    100114class RectangleDispersion(Dispersion): 
    … …  
    107121    """ 
    108122    type = "rectangle" 
    109     default = dict(npts=35, width=0, nsigmas=1.70325) 
    110     def _weights(self, center, sigma, lb, ub): 
    111         x = self._linspace(center, sigma, lb, ub) 
    112         x = x[np.fabs(x-center) <= np.fabs(sigma)*sqrt(3.0)] 
    113         return x, np.ones_like(x) 
    114  
     123    default = dict(npts=35, width=0, nsigmas=1.73205) 
     124    def _weights(self, center, sigma, lb, ub): 
     125         x = self._linspace(center, sigma, lb, ub) 
     126         x = x[np.fabs(x-center) <= np.fabs(sigma)*sqrt(3.0)] 
     127         return x, np.ones_like(x) 
    115128 
    116129class LogNormalDispersion(Dispersion): 
    … …  
    190203        return x, px 
    191204 
     205class BoltzmannDispersion(Dispersion): 
     206    r""" 
     207    Boltzmann dispersion, with $\sigma=k T/E$. 
     208 
     209    .. math:: 
     210 
     211        w = \exp\left( -|x - c|/\sigma\right) 
     212    """ 
     213    type = "boltzmann" 
     214    default = dict(npts=35, width=0, nsigmas=3) 
     215    def _weights(self, center, sigma, lb, ub): 
     216        x = self._linspace(center, sigma, lb, ub) 
     217        px = np.exp(-np.abs(x-center) / np.abs(sigma)) 
     218        return x, px 
    192219 
    193220# dispersion name -> disperser lookup table. 
    … …  
    196223MODELS = OrderedDict((d.type, d) for d in ( 
    197224    RectangleDispersion, 
     225    UniformDispersion, 
    198226    ArrayDispersion, 
    199227    LogNormalDispersion, 
    200228    GaussianDispersion, 
    201229    SchulzDispersion, 
     230    BoltzmannDispersion 
    202231)) 
    203232 

  • appveyor.yml

    r1258e32 rfefc185  
    44  # /E:ON and /V:ON options are not enabled in the batch script interpreter 
    55  # See: http://stackoverflow.com/a/13751649/163740 
    6   CMD_IN_ENV: "cmd /E:ON /V:ON /C obvci_appveyor_python_build_env.cmd" 
     6  # 2018-01-19 PAK: probably irrelevant now that we are using tinycc rather than msvc 
     7  #CMD_IN_ENV: "cmd /E:ON /V:ON /C obvci_appveyor_python_build_env.cmd" 
    78 
    89  # Workaround for https://github.com/conda/conda-build/issues/636 
    … …  
    2223    - TARGET_ARCH: "x86" 
    2324      CONDA_PY: "35" 
    24       PY_CONDITION: "python >=3.5" 
     25      PY_CONDITION: "python >=3.5,<3.6" 
    2526      CONDA_INSTALL_LOCN: "C:\\Miniconda35" 
    2627    - TARGET_ARCH: "x64" 
    … …  
    3435    - TARGET_ARCH: "x64" 
    3536      CONDA_PY: "35" 
    36       PY_CONDITION: "python >=3.5" 
     37      PY_CONDITION: "python >=3.5,<3.6" 
    3738      CONDA_INSTALL_LOCN: "C:\\Miniconda35-x64" 
    3839 
    … …  
    4344 
    4445install: 
    45     # Set the CONDA_NPY, although it has no impact on the actual build. We need this because of a test within conda-build. 
     46    # Set the CONDA_NPY, although it has no impact on the actual build.  
     47    # We need this because of a test within conda-build. 
    4648    - cmd: set CONDA_NPY=19 
    4749 
    4850    # Remove cygwin (and therefore the git that comes with it). 
    49     - cmd: rmdir C:\cygwin /s /q 
     51    # 2018-01-19 PAK: irrelevant since we already pulled the repo 
     52    #- cmd: rmdir C:\cygwin /s /q 
     53 
     54    # Set the conda path; would be nice to do this 
     55    - cmd: path %CONDA_INSTALL_LOCN%\Scripts;%PATH% 
    5056 
    5157    # Use the pre-installed Miniconda for the desired arch 
    … …  
    5662    # so that we can update it. Then we remove it so that 
    5763    # we can do a proper activation. 
    58     - cmd: set "OLDPATH=%PATH%" 
    59     - cmd: set "PATH=%CONDA_INSTALL_LOCN%\\Scripts;%CONDA_INSTALL_LOCN%\\Library\\bin;%PATH%" 
    6064    - cmd: conda update --yes --quiet conda python 
    61     - cmd: set "PATH=%OLDPATH%" 
    6265    - cmd: call %CONDA_INSTALL_LOCN%\Scripts\activate.bat 
    63  
    64     - cmd: conda config --add channels conda-forge 
    65     - cmd: conda config --set show_channel_urls true 
    66     - cmd: conda update --yes --quiet conda 
    67     - cmd: conda install --yes --quiet obvious-ci 
    68     - cmd: conda install --yes --quiet numpy toolchain scipy cython cffi 
     66    #- cmd: conda config --add channels conda-forge 
     67    #- cmd: conda config --set show_channel_urls true 
     68    #- cmd: conda install --yes --quiet obvious-ci 
     69    # 2018-01-19 PAK: skipping toolchain cython and cffi (these were for pyopencl?) 
     70    - cmd: conda install --yes --quiet numpy scipy 
     71    # 2018-01-19 PAK: skipping pyopencl; this would be needed for deploy but maybe not for test 
    6972    #- cmd: conda install --yes --channel conda-forge pyopencl 
     73    # 2018-01-19 PAK: 3rd party packages might need msvc, so %CMD_IN_ENV% may be needed for pip 
    7074    - cmd: pip install bumps unittest-xml-reporting tinycc 
    7175 
    7276build_script: 
    7377    # Build the project 
    74     - "%CMD_IN_ENV% python setup.py build" 
     78    # 2018-01-19 PAK: maybe need one of this if using msvc? 
     79    #- "%CMD_IN_ENV% python setup.py build" 
     80    #- cmd /E:ON /V:ON /C obvci_appveyor_python_build_env.cmd python setup.py build 
     81    - python setup.py build 
    7582 
    7683test_script: 
    7784    # Run the project tests 
    78     - "%CMD_IN_ENV% python -m sasmodels.model_test dll all" 
     85    # 2018-01-19 PAK: maybe need one of this if using msvc? 
     86    #- "%CMD_IN_ENV% python -m sasmodels.model_test dll all" 
     87    #- cmd /E:ON /V:ON /C obvci_appveyor_python_build_env.cmd python -m sasmoels.model_test dll all 
     88    - python -m sasmodels.model_test dll all 

  • doc/developer/overview.rst

    r2a7e20e r0d5a655  
    190190*sasmodels.generate.PROJECTION*, with the default *PROJECTION=1* for 
    191191equirectangular and *PROJECTION=2* for sinusoidal.  The more complicated 
    192 Guyou and Postel projections are not implemented. See explore.jitter.draw_mesh 
     192Guyou and Postel projections are not implemented. See jitter.draw_mesh 
    193193for details. 
    194194 
    … …  
    274274code. 
    275275 
    276 *explore/jitter.py* is for exploring different options for handling 
    277 orientation and orientation dispersity.  It uses *explore/guyou.py* to 
     276*sasmodels/jitter.py* is for exploring different options for handling 
     277orientation and orientation dispersity.  It uses *sasmodels/guyou.py* to 
    278278generate the Guyou projection. 
    279279 

  • doc/genapi.py

    r706f466 r0d5a655  
    6969    ('exception', 'Annotate exceptions'), 
    7070    ('generate', 'Model parser'), 
     71    ('jitter', 'Orientation explorer'), 
     72    ('guyou', 'Guyou map projection'), 
    7173    ('kernel', 'Evaluator type definitions'), 
    7274    ('kernelcl', 'OpenCL model evaluator'), 

  • doc/guide/orientation/orientation.rst

    r5fb0634 r0d5a655  
    5252yaw angle $\theta$ about the $b$-axis and pitch angle $\phi$ about the 
    5353$a$-axis. 
     54 
     55You can explore the view and jitter angles interactively using 
     56:func:`sasmodels.jitter.run`.  Enter the following into the python 
     57interpreter:: 
     58 
     59    from sasmodels import jitter 
     60    jitter.run() 
    5461 
    5562More formally, starting with axes $a$-$b$-$c$ of the particle aligned 

  • explore/realspace.py

    r8fb2a94 r032646d  
    33import time 
    44from copy import copy 
     5import os 
     6import argparse 
     7from collections import OrderedDict 
    58 
    69import numpy as np 
    710from numpy import pi, radians, sin, cos, sqrt 
    8 from numpy.random import poisson, uniform 
     11from numpy.random import poisson, uniform, randn, rand 
    912from numpy.polynomial.legendre import leggauss 
    1013from scipy.integrate import simps 
    1114from scipy.special import j1 as J1 
     15 
     16try: 
     17    import numba 
     18    USE_NUMBA = True 
     19except ImportError: 
     20    USE_NUMBA = False 
    1221 
    1322# Definition of rotation matrices comes from wikipedia: 
    … …  
    8291        raise NotImplementedError() 
    8392 
     93    def dims(self): 
     94        # type: () -> float, float, float 
     95        raise NotImplementedError() 
     96 
    8497    def rotate(self, theta, phi, psi): 
    8598        self.rotation = rotation(theta, phi, psi) * self.rotation 
    … …  
    116129                     for s2 in shapes] 
    117130        self.r_max = max(distances + [s.r_max for s in shapes]) 
    118  
    119     def volume(self): 
    120         return sum(shape.volume() for shape in self.shapes) 
     131        self.volume = sum(shape.volume for shape in self.shapes) 
    121132 
    122133    def sample(self, density): 
    … …  
    133144        self._scale = np.array([a/2, b/2, c/2])[None, :] 
    134145        self.r_max = sqrt(a**2 + b**2 + c**2) 
    135  
    136     def volume(self): 
    137         return self.a*self.b*self.c 
     146        self.dims = a, b, c 
     147        self.volume = a*b*c 
    138148 
    139149    def sample(self, density): 
    140         num_points = poisson(density*self.a*self.b*self.c) 
     150        num_points = poisson(density*self.volume) 
    141151        points = self._scale*uniform(-1, 1, size=(num_points, 3)) 
    142152        values = self.value.repeat(points.shape[0]) 
    … …  
    152162        self._scale = np.array([ra, rb, length/2])[None, :] 
    153163        self.r_max = sqrt(4*max(ra, rb)**2 + length**2) 
    154  
    155     def volume(self): 
    156         return pi*self.ra*self.rb*self.length 
     164        self.dims = 2*ra, 2*rb, length 
     165        self.volume = pi*ra*rb*length 
    157166 
    158167    def sample(self, density): 
    159         # density of the bounding box 
     168        # randomly sample from a box of side length 2*r, excluding anything 
     169        # not in the cylinder 
    160170        num_points = poisson(density*4*self.ra*self.rb*self.length) 
    161171        points = uniform(-1, 1, size=(num_points, 3)) 
    … …  
    174184        self._scale = np.array([ra, rb, rc])[None, :] 
    175185        self.r_max = 2*max(ra, rb, rc) 
    176  
    177     def volume(self): 
    178         return 4*pi/3 * self.ra * self.rb * self.rc 
     186        self.dims = 2*ra, 2*rb, 2*rc 
     187        self.volume = 4*pi/3 * ra * rb * rc 
    179188 
    180189    def sample(self, density): 
    … …  
    194203        self.rotate(*orientation) 
    195204        self.shift(*center) 
     205        helix_length = helix_pitch * tube_length/sqrt(helix_radius**2 + helix_pitch**2) 
     206        total_radius = self.helix_radius + self.tube_radius 
    196207        self.helix_radius, self.helix_pitch = helix_radius, helix_pitch 
    197208        self.tube_radius, self.tube_length = tube_radius, tube_length 
    198         helix_length = helix_pitch * tube_length/sqrt(helix_radius**2 + helix_pitch**2) 
    199         self.r_max = sqrt((2*helix_radius + 2*tube_radius)*2 
    200                           + (helix_length + 2*tube_radius)**2) 
    201  
    202     def volume(self): 
     209        self.r_max = sqrt(4*total_radius + (helix_length + 2*tube_radius)**2) 
     210        self.dims = 2*total_radius, 2*total_radius, helix_length 
    203211        # small tube radius approximation; for larger tubes need to account 
    204212        # for the fact that the inner length is much shorter than the outer 
    205213        # length 
    206         return pi*self.tube_radius**2*self.tube_length 
     214        self.volume = pi*self.tube_radius**2*self.tube_length 
    207215 
    208216    def points(self, density): 
    … …  
    244252        return values, self._adjust(points) 
    245253 
    246 NUMBA = False 
    247 if NUMBA: 
     254def csbox(a=10, b=20, c=30, da=1, db=2, dc=3, slda=1, sldb=2, sldc=3, sld_core=4): 
     255    core = Box(a, b, c, sld_core) 
     256    side_a = Box(da, b, c, slda, center=((a+da)/2, 0, 0)) 
     257    side_b = Box(a, db, c, sldb, center=(0, (b+db)/2, 0)) 
     258    side_c = Box(a, b, dc, sldc, center=(0, 0, (c+dc)/2)) 
     259    side_a2 = copy(side_a).shift(-a-da, 0, 0) 
     260    side_b2 = copy(side_b).shift(0, -b-db, 0) 
     261    side_c2 = copy(side_c).shift(0, 0, -c-dc) 
     262    shape = Composite((core, side_a, side_b, side_c, side_a2, side_b2, side_c2)) 
     263    shape.dims = 2*da+a, 2*db+b, 2*dc+c 
     264    return shape 
     265 
     266def _Iqxy(values, x, y, z, qa, qb, qc): 
     267    """I(q) = |sum V(r) rho(r) e^(1j q.r)|^2 / sum V(r)""" 
     268    Iq = [abs(np.sum(values*np.exp(1j*(qa_k*x + qb_k*y + qc_k*z))))**2 
     269            for qa_k, qb_k, qc_k in zip(qa.flat, qb.flat, qc.flat)] 
     270    return Iq 
     271 
     272if USE_NUMBA: 
     273    # Override simple numpy solution with numba if available 
    248274    from numba import njit 
    249275    @njit("f8[:](f8[:],f8[:],f8[:],f8[:],f8[:],f8[:],f8[:])") 
    … …  
    252278        for j in range(len(Iq)): 
    253279            total = 0. + 0j 
    254             for k in range(len(Iq)): 
     280            for k in range(len(values)): 
    255281                total += values[k]*np.exp(1j*(qa[j]*x[k] + qb[j]*y[k] + qc[j]*z[k])) 
    256282            Iq[j] = abs(total)**2 
    … …  
    263289    values = rho*volume 
    264290    x, y, z = points.T 
     291    values, x, y, z, qa, qb, qc = [np.asarray(v, 'd') 
     292                                   for v in (values, x, y, z, qa, qb, qc)] 
    265293 
    266294    # I(q) = |sum V(r) rho(r) e^(1j q.r)|^2 / sum V(r) 
    267     if NUMBA: 
    268         Iq = _Iqxy(values, x, y, z, qa.flatten(), qb.flatten(), qc.flatten()) 
    269     else: 
    270         Iq = [abs(np.sum(values*np.exp(1j*(qa_k*x + qb_k*y + qc_k*z))))**2 
    271               for qa_k, qb_k, qc_k in zip(qa.flat, qb.flat, qc.flat)] 
     295    Iq = _Iqxy(values, x, y, z, qa.flatten(), qb.flatten(), qc.flatten()) 
    272296    return np.asarray(Iq).reshape(qx.shape) / np.sum(volume) 
    273297 
    … …  
    338362""" 
    339363 
    340 if NUMBA: 
     364if USE_NUMBA: 
     365    # Override simple numpy solution with numba if available 
    341366    @njit("f8[:](f8[:], f8[:], f8[:,:])") 
    342     def _calc_Pr_uniform_jit(r, rho, points): 
     367    def _calc_Pr_uniform(r, rho, points): 
    343368        dr = r[0] 
    344369        n_max = len(r) 
    … …  
    362387    # P(r) with uniform steps in r is 3x faster; check if we are uniform 
    363388    # before continuing 
     389    r, rho, points = [np.asarray(v, 'd') for v in (r, rho, points)] 
    364390    if np.max(np.abs(np.diff(r) - r[0])) > r[0]*0.01: 
    365391        Pr = _calc_Pr_nonuniform(r, rho, points) 
    366392    else: 
    367         if NUMBA: 
    368             Pr = _calc_Pr_uniform_jit(r, rho, points) 
    369         else: 
    370             Pr = _calc_Pr_uniform(r, rho, points) 
     393        Pr = _calc_Pr_uniform(r, rho, points) 
    371394    return Pr / Pr.max() 
    372395 
    … …  
    399422    return edges 
    400423 
     424# -------------- plotters ---------------- 
    401425def plot_calc(r, Pr, q, Iq, theory=None): 
    402426    import matplotlib.pyplot as plt 
    … …  
    410434    plt.ylabel('Iq') 
    411435    if theory is not None: 
    412         plt.loglog(theory[0], theory[1], '-', label='analytic') 
     436        plt.loglog(theory[0], theory[1]/theory[1][0], '-', label='analytic') 
    413437        plt.legend() 
    414438 
    … …  
    444468    ax.autoscale(True) 
    445469 
    446 def check_shape(shape, fn=None): 
    447     rho_solvent = 0 
    448     q = np.logspace(-3, 0, 200) 
    449     r = shape.r_bins(q, r_step=0.01) 
    450     sampling_density = 6*5000 / shape.volume() 
    451     rho, points = shape.sample(sampling_density) 
    452     t0 = time.time() 
    453     Pr = calc_Pr(r, rho-rho_solvent, points) 
    454     print("calc Pr time", time.time() - t0) 
    455     Iq = calc_Iq(q, r, Pr) 
    456     theory = (q, fn(q)) if fn is not None else None 
    457  
    458     import pylab 
    459     #plot_points(rho, points); pylab.figure() 
    460     plot_calc(r, Pr, q, Iq, theory=theory) 
    461     pylab.show() 
    462  
    463 def check_shape_2d(shape, fn=None, view=(0, 0, 0)): 
    464     rho_solvent = 0 
    465     nq, qmax = 100, 1.0 
    466     qx = np.linspace(0.0, qmax, nq) 
    467     qy = np.linspace(0.0, qmax, nq) 
    468     Qx, Qy = np.meshgrid(qx, qy) 
    469     sampling_density = 50000 / shape.volume() 
    470     #t0 = time.time() 
    471     rho, points = shape.sample(sampling_density) 
    472     #print("sample time", time.time() - t0) 
    473     t0 = time.time() 
    474     Iqxy = calc_Iqxy(Qx, Qy, rho, points, view=view) 
    475     print("calc time", time.time() - t0) 
    476     theory = fn(Qx, Qy) if fn is not None else None 
    477     Iqxy += 0.001 * Iqxy.max() 
    478     if theory is not None: 
    479         theory += 0.001 * theory.max() 
    480  
    481     import pylab 
    482     #plot_points(rho, points); pylab.figure() 
    483     plot_calc_2d(qx, qy, Iqxy, theory=theory) 
    484     pylab.show() 
    485  
     470# ----------- Analytic models -------------- 
    486471def sas_sinx_x(x): 
    487472    with np.errstate(all='ignore'): 
    … …  
    510495    for k, qk in enumerate(q): 
    511496        qab, qc = qk*sin_alpha, qk*cos_alpha 
    512         Fq = sas_2J1x_x(qab*radius) * j0(qc*length/2) 
     497        Fq = sas_2J1x_x(qab*radius) * sas_sinx_x(qc*length/2) 
    513498        Iq[k] = np.sum(w*Fq**2) 
    514     Iq = Iq/Iq[0] 
     499    Iq = Iq 
    515500    return Iq 
    516501 
    517502def cylinder_Iqxy(qx, qy, radius, length, view=(0, 0, 0)): 
    518503    qa, qb, qc = invert_view(qx, qy, view) 
    519     qab = np.sqrt(qa**2 + qb**2) 
    520     Fq = sas_2J1x_x(qab*radius) * j0(qc*length/2) 
     504    qab = sqrt(qa**2 + qb**2) 
     505    Fq = sas_2J1x_x(qab*radius) * sas_sinx_x(qc*length/2) 
    521506    Iq = Fq**2 
    522507    return Iq.reshape(qx.shape) 
    … …  
    524509def sphere_Iq(q, radius): 
    525510    Iq = sas_3j1x_x(q*radius)**2 
    526     return Iq/Iq[0] 
     511    return Iq 
     512 
     513def box_Iq(q, a, b, c): 
     514    z, w = leggauss(76) 
     515    outer_sum = np.zeros_like(q) 
     516    for cos_alpha, outer_w in zip((z+1)/2, w): 
     517        sin_alpha = sqrt(1.0-cos_alpha*cos_alpha) 
     518        qc = q*cos_alpha 
     519        siC = c*sas_sinx_x(c*qc/2) 
     520        inner_sum = np.zeros_like(q) 
     521        for beta, inner_w in zip((z + 1)*pi/4, w): 
     522            qa, qb = q*sin_alpha*sin(beta), q*sin_alpha*cos(beta) 
     523            siA = a*sas_sinx_x(a*qa/2) 
     524            siB = b*sas_sinx_x(b*qb/2) 
     525            Fq = siA*siB*siC 
     526            inner_sum += inner_w * Fq**2 
     527        outer_sum += outer_w * inner_sum 
     528    Iq = outer_sum / 4  # = outer*um*zm*8.0/(4.0*M_PI) 
     529    return Iq 
     530 
     531def box_Iqxy(qx, qy, a, b, c, view=(0, 0, 0)): 
     532    qa, qb, qc = invert_view(qx, qy, view) 
     533    sia = sas_sinx_x(qa*a/2) 
     534    sib = sas_sinx_x(qb*b/2) 
     535    sic = sas_sinx_x(qc*c/2) 
     536    Fq = sia*sib*sic 
     537    Iq = Fq**2 
     538    return Iq.reshape(qx.shape) 
    527539 
    528540def csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core): 
    … …  
    539551        sin_alpha = sqrt(1.0-cos_alpha*cos_alpha) 
    540552        qc = q*cos_alpha 
    541         siC = c*j0(c*qc/2) 
    542         siCt = tC*j0(tC*qc/2) 
     553        siC = c*sas_sinx_x(c*qc/2) 
     554        siCt = tC*sas_sinx_x(tC*qc/2) 
    543555        inner_sum = np.zeros_like(q) 
    544556        for beta, inner_w in zip((z + 1)*pi/4, w): 
    545557            qa, qb = q*sin_alpha*sin(beta), q*sin_alpha*cos(beta) 
    546             siA = a*j0(a*qa/2) 
    547             siB = b*j0(b*qb/2) 
    548             siAt = tA*j0(tA*qa/2) 
    549             siBt = tB*j0(tB*qb/2) 
     558            siA = a*sas_sinx_x(a*qa/2) 
     559            siB = b*sas_sinx_x(b*qb/2) 
     560            siAt = tA*sas_sinx_x(tA*qa/2) 
     561            siBt = tB*sas_sinx_x(tB*qb/2) 
    550562            if overlapping: 
    551563                Fq = (dr0*siA*siB*siC 
    … …  
    584596    return Iq.reshape(qx.shape) 
    585597 
    586 def check_cylinder(radius=25, length=125, rho=2.): 
     598# --------- Test cases ----------- 
     599 
     600def build_cylinder(radius=25, length=125, rho=2.): 
    587601    shape = EllipticalCylinder(radius, radius, length, rho) 
    588     fn = lambda q: cylinder_Iq(q, radius, length) 
    589     check_shape(shape, fn) 
    590  
    591 def check_cylinder_2d(radius=25, length=125, rho=2., view=(0, 0, 0)): 
    592     shape = EllipticalCylinder(radius, radius, length, rho) 
    593     fn = lambda qx, qy, view=view: cylinder_Iqxy(qx, qy, radius, length, view=view) 
    594     check_shape_2d(shape, fn, view=view) 
    595  
    596 def check_cylinder_2d_lattice(radius=25, length=125, rho=2., 
    597                               view=(0, 0, 0)): 
    598     nx, dx = 1, 2*radius 
    599     ny, dy = 30, 2*radius 
    600     nz, dz = 30, length 
    601     dx, dy, dz = 2*dx, 2*dy, 2*dz 
    602     def center(*args): 
    603         sigma = 0.333 
    604         space = 2 
    605         return [(space*n+np.random.randn()*sigma)*x for n, x in args] 
    606     shapes = [EllipticalCylinder(radius, radius, length, rho, 
    607                                  #center=(ix*dx, iy*dy, iz*dz) 
    608                                  orientation=np.random.randn(3)*0, 
    609                                  center=center((ix, dx), (iy, dy), (iz, dz)) 
    610                                 ) 
     602    fn = lambda q: cylinder_Iq(q, radius, length)*rho**2 
     603    fn_xy = lambda qx, qy, view: cylinder_Iqxy(qx, qy, radius, length, view=view)*rho**2 
     604    return shape, fn, fn_xy 
     605 
     606def build_sphere(radius=125, rho=2): 
     607    shape = TriaxialEllipsoid(radius, radius, radius, rho) 
     608    fn = lambda q: sphere_Iq(q, radius)*rho**2 
     609    fn_xy = lambda qx, qy, view: sphere_Iq(np.sqrt(qx**2+qy**2), radius)*rho**2 
     610    return shape, fn, fn_xy 
     611 
     612def build_box(a=10, b=20, c=30, rho=2.): 
     613    shape = Box(a, b, c, rho) 
     614    fn = lambda q: box_Iq(q, a, b, c)*rho**2 
     615    fn_xy = lambda qx, qy, view: box_Iqxy(qx, qy, a, b, c, view=view)*rho**2 
     616    return shape, fn, fn_xy 
     617 
     618def build_csbox(a=10, b=20, c=30, da=1, db=2, dc=3, slda=1, sldb=2, sldc=3, sld_core=4): 
     619    shape = csbox(a, b, c, da, db, dc, slda, sldb, sldc, sld_core) 
     620    fn = lambda q: csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core) 
     621    fn_xy = lambda qx, qy, view: csbox_Iqxy(qx, qy, a, b, c, da, db, dc, 
     622                                            slda, sldb, sldc, sld_core, view=view) 
     623    return shape, fn, fn_xy 
     624 
     625def build_cubic_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2, 
     626                  shuffle=0, rotate=0): 
     627    a, b, c = shape.dims 
     628    shapes = [copy(shape) 
     629              .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 
     630                     (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 
     631                     (iz+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 
     632              .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 
    611633              for ix in range(nx) 
    612634              for iy in range(ny) 
    613635              for iz in range(nz)] 
    614     shape = Composite(shapes) 
    615     fn = lambda qx, qy, view=view: cylinder_Iqxy(qx, qy, radius, length, view=view) 
    616     check_shape_2d(shape, fn, view=view) 
    617  
    618 def check_sphere(radius=125, rho=2): 
    619     shape = TriaxialEllipsoid(radius, radius, radius, rho) 
    620     fn = lambda q: sphere_Iq(q, radius) 
    621     check_shape(shape, fn) 
    622  
    623 def check_csbox(a=10, b=20, c=30, da=1, db=2, dc=3, slda=1, sldb=2, sldc=3, sld_core=4): 
    624     core = Box(a, b, c, sld_core) 
    625     side_a = Box(da, b, c, slda, center=((a+da)/2, 0, 0)) 
    626     side_b = Box(a, db, c, sldb, center=(0, (b+db)/2, 0)) 
    627     side_c = Box(a, b, dc, sldc, center=(0, 0, (c+dc)/2)) 
    628     side_a2 = copy(side_a).shift(-a-da, 0, 0) 
    629     side_b2 = copy(side_b).shift(0, -b-db, 0) 
    630     side_c2 = copy(side_c).shift(0, 0, -c-dc) 
    631     shape = Composite((core, side_a, side_b, side_c, side_a2, side_b2, side_c2)) 
    632     def fn(q): 
    633         return csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core) 
    634     #check_shape(shape, fn) 
    635  
    636     view = (20, 30, 40) 
    637     def fn_xy(qx, qy): 
    638         return csbox_Iqxy(qx, qy, a, b, c, da, db, dc, 
    639                           slda, sldb, sldc, sld_core, view=view) 
    640     check_shape_2d(shape, fn_xy, view=view) 
     636    lattice = Composite(shapes) 
     637    return lattice 
     638 
     639 
     640SHAPE_FUNCTIONS = OrderedDict([ 
     641    ("cylinder", build_cylinder), 
     642    ("sphere", build_sphere), 
     643    ("box", build_box), 
     644    ("csbox", build_csbox), 
     645]) 
     646SHAPES = list(SHAPE_FUNCTIONS.keys()) 
     647 
     648def check_shape(title, shape, fn=None, show_points=False, 
     649                mesh=100, qmax=1.0, r_step=0.01, samples=5000): 
     650    rho_solvent = 0 
     651    qmin = qmax/100. 
     652    q = np.logspace(np.log10(qmin), np.log10(qmax), mesh) 
     653    r = shape.r_bins(q, r_step=r_step) 
     654    sampling_density = samples / shape.volume 
     655    rho, points = shape.sample(sampling_density) 
     656    t0 = time.time() 
     657    Pr = calc_Pr(r, rho-rho_solvent, points) 
     658    print("calc Pr time", time.time() - t0) 
     659    Iq = calc_Iq(q, r, Pr) 
     660    theory = (q, fn(q)) if fn is not None else None 
     661 
     662    import pylab 
     663    if show_points: 
     664         plot_points(rho, points); pylab.figure() 
     665    plot_calc(r, Pr, q, Iq, theory=theory) 
     666    pylab.gcf().canvas.set_window_title(title) 
     667    pylab.show() 
     668 
     669def check_shape_2d(title, shape, fn=None, view=(0, 0, 0), show_points=False, 
     670                   mesh=100, qmax=1.0, samples=5000): 
     671    rho_solvent = 0 
     672    qx = np.linspace(0.0, qmax, mesh) 
     673    qy = np.linspace(0.0, qmax, mesh) 
     674    Qx, Qy = np.meshgrid(qx, qy) 
     675    sampling_density = samples / shape.volume 
     676    t0 = time.time() 
     677    rho, points = shape.sample(sampling_density) 
     678    print("point generation time", time.time() - t0) 
     679    t0 = time.time() 
     680    Iqxy = calc_Iqxy(Qx, Qy, rho, points, view=view) 
     681    print("calc Iqxy time", time.time() - t0) 
     682    theory = fn(Qx, Qy, view) if fn is not None else None 
     683    Iqxy += 0.001 * Iqxy.max() 
     684    if theory is not None: 
     685        theory += 0.001 * theory.max() 
     686 
     687    import pylab 
     688    if show_points: 
     689        plot_points(rho, points); pylab.figure() 
     690    plot_calc_2d(qx, qy, Iqxy, theory=theory) 
     691    pylab.gcf().canvas.set_window_title(title) 
     692    pylab.show() 
     693 
     694def main(): 
     695    parser = argparse.ArgumentParser( 
     696        description="Compute scattering from realspace sampling", 
     697        formatter_class=argparse.ArgumentDefaultsHelpFormatter, 
     698        ) 
     699    parser.add_argument('-d', '--dim', type=int, default=1, 
     700                        help='dimension 1 or 2') 
     701    parser.add_argument('-m', '--mesh', type=int, default=100, 
     702                        help='number of mesh points') 
     703    parser.add_argument('-s', '--samples', type=int, default=5000, 
     704                        help="number of sample points") 
     705    parser.add_argument('-q', '--qmax', type=float, default=0.5, 
     706                        help='max q') 
     707    parser.add_argument('-v', '--view', type=str, default='0,0,0', 
     708                        help='theta,phi,psi angles') 
     709    parser.add_argument('-n', '--lattice', type=str, default='1,1,1', 
     710                        help='lattice size') 
     711    parser.add_argument('-z', '--spacing', type=str, default='2,2,2', 
     712                        help='lattice spacing') 
     713    parser.add_argument('-r', '--rotate', type=float, default=0., 
     714                        help="rotation relative to lattice, gaussian < 30 degrees, uniform otherwise") 
     715    parser.add_argument('-w', '--shuffle', type=float, default=0., 
     716                        help="position relative to lattice, gaussian < 0.3, uniform otherwise") 
     717    parser.add_argument('-p', '--plot', action='store_true', 
     718                        help='plot points') 
     719    parser.add_argument('shape', choices=SHAPES, nargs='?', default=SHAPES[0], 
     720                        help='oriented shape') 
     721    parser.add_argument('pars', type=str, nargs='*', help='shape parameters') 
     722    opts = parser.parse_args() 
     723    pars = {key: float(value) for p in opts.pars for key, value in [p.split('=')]} 
     724    nx, ny, nz = [int(v) for v in opts.lattice.split(',')] 
     725    dx, dy, dz = [float(v) for v in opts.spacing.split(',')] 
     726    shuffle, rotate = opts.shuffle, opts.rotate 
     727    shape, fn, fn_xy = SHAPE_FUNCTIONS[opts.shape](**pars) 
     728    if nx > 1 or ny > 1 or nz > 1: 
     729        shape = build_cubic_lattice(shape, nx, ny, nz, dx, dy, dz, shuffle, rotate) 
     730    title = "%s(%s)" % (opts.shape, " ".join(opts.pars)) 
     731    if opts.dim == 1: 
     732        check_shape(title, shape, fn, show_points=opts.plot, 
     733                    mesh=opts.mesh, qmax=opts.qmax, samples=opts.samples) 
     734    else: 
     735        view = tuple(float(v) for v in opts.view.split(',')) 
     736        check_shape_2d(title, shape, fn_xy, view=view, show_points=opts.plot, 
     737                       mesh=opts.mesh, qmax=opts.qmax, samples=opts.samples) 
     738 
    641739 
    642740if __name__ == "__main__": 
    643     check_cylinder(radius=10, length=40) 
    644     #check_cylinder_2d(radius=10, length=40, view=(90,30,0)) 
    645     #check_cylinder_2d_lattice(radius=10, length=50, view=(90,30,0)) 
    646     #check_sphere() 
    647     #check_csbox() 
    648     #check_csbox(da=100, db=200, dc=300) 
     741    main() 

  • sasmodels/core.py

    r2d81cfe r7a516d0  
    148148    used with functions in generate to build the docs or extract model info. 
    149149    """ 
    150     if '@' in model_string: 
    151         parts = model_string.split('@') 
    152         if len(parts) != 2: 
    153             raise ValueError("Use P@S to apply a structure factor S to model P") 
    154         P_info, Q_info = [load_model_info(part) for part in parts] 
    155         return product.make_product_info(P_info, Q_info) 
    156  
    157     product_parts = [] 
    158     addition_parts = [] 
    159  
    160     addition_parts_names = model_string.split('+') 
    161     if len(addition_parts_names) >= 2: 
    162         addition_parts = [load_model_info(part) for part in addition_parts_names] 
    163     elif len(addition_parts_names) == 1: 
    164         product_parts_names = model_string.split('*') 
    165         if len(product_parts_names) >= 2: 
    166             product_parts = [load_model_info(part) for part in product_parts_names] 
    167         elif len(product_parts_names) == 1: 
    168             if "custom." in product_parts_names[0]: 
    169                 # Extract ModelName from "custom.ModelName" 
    170                 pattern = "custom.([A-Za-z0-9_-]+)" 
    171                 result = re.match(pattern, product_parts_names[0]) 
    172                 if result is None: 
    173                     raise ValueError("Model name in invalid format: " + product_parts_names[0]) 
    174                 model_name = result.group(1) 
    175                 # Use ModelName to find the path to the custom model file 
    176                 model_path = joinpath(CUSTOM_MODEL_PATH, model_name + ".py") 
    177                 if not os.path.isfile(model_path): 
    178                     raise ValueError("The model file {} doesn't exist".format(model_path)) 
    179                 kernel_module = custom.load_custom_kernel_module(model_path) 
    180                 return modelinfo.make_model_info(kernel_module) 
    181             # Model is a core model 
    182             kernel_module = generate.load_kernel_module(product_parts_names[0]) 
    183             return modelinfo.make_model_info(kernel_module) 
    184  
    185     model = None 
    186     if len(product_parts) > 1: 
    187         model = mixture.make_mixture_info(product_parts, operation='*') 
    188     if len(addition_parts) > 1: 
    189         if model is not None: 
    190             addition_parts.append(model) 
    191         model = mixture.make_mixture_info(addition_parts, operation='+') 
    192     return model 
     150    if "+" in model_string: 
     151        parts = [load_model_info(part) 
     152                 for part in model_string.split("+")] 
     153        return mixture.make_mixture_info(parts, operation='+') 
     154    elif "*" in model_string: 
     155        parts = [load_model_info(part) 
     156                 for part in model_string.split("*")] 
     157        return mixture.make_mixture_info(parts, operation='*') 
     158    elif "@" in model_string: 
     159        p_info, q_info = [load_model_info(part) 
     160                          for part in model_string.split("@")] 
     161        return product.make_product_info(p_info, q_info) 
     162    # We are now dealing with a pure model 
     163    elif "custom." in model_string: 
     164        pattern = "custom.([A-Za-z0-9_-]+)" 
     165        result = re.match(pattern, model_string) 
     166        if result is None: 
     167            raise ValueError("Model name in invalid format: " + model_string) 
     168        model_name = result.group(1) 
     169        # Use ModelName to find the path to the custom model file 
     170        model_path = joinpath(CUSTOM_MODEL_PATH, model_name + ".py") 
     171        if not os.path.isfile(model_path): 
     172            raise ValueError("The model file {} doesn't exist".format(model_path)) 
     173        kernel_module = custom.load_custom_kernel_module(model_path) 
     174        return modelinfo.make_model_info(kernel_module) 
     175    kernel_module = generate.load_kernel_module(model_string) 
     176    return modelinfo.make_model_info(kernel_module) 
    193177 
    194178 
    … …  
    336320    print("\n".join(list_models(kind))) 
    337321 
     322def test_load(): 
     323    # type: () -> None 
     324    """Check that model load works""" 
     325    def test_models(fst, snd): 
     326        """Confirm that two models produce the same parameters""" 
     327        fst = load_model(fst) 
     328        snd = load_model(snd) 
     329        # Remove the upper case characters frin the parameters, since 
     330        # they lasndel the order of events and we specfically are 
     331        # changin that aspect 
     332        fst = [[x for x in p.name if x == x.lower()] for p in fst.info.parameters.kernel_parameters] 
     333        snd = [[x for x in p.name if x == x.lower()] for p in snd.info.parameters.kernel_parameters] 
     334        assert sorted(fst) == sorted(snd), "{} != {}".format(fst, snd) 
     335 
     336 
     337    test_models( 
     338        "cylinder+sphere", 
     339        "sphere+cylinder") 
     340    test_models( 
     341        "cylinder*sphere", 
     342        "sphere*cylinder") 
     343    test_models( 
     344        "cylinder@hardsphere*sphere", 
     345        "sphere*cylinder@hardsphere") 
     346    test_models( 
     347        "barbell+sphere*cylinder@hardsphere", 
     348        "sphere*cylinder@hardsphere+barbell") 
     349    test_models( 
     350        "barbell+cylinder@hardsphere*sphere", 
     351        "cylinder@hardsphere*sphere+barbell") 
     352    test_models( 
     353        "barbell+sphere*cylinder@hardsphere", 
     354        "barbell+cylinder@hardsphere*sphere") 
     355    test_models( 
     356        "sphere*cylinder@hardsphere+barbell", 
     357        "cylinder@hardsphere*sphere+barbell") 
     358    test_models( 
     359        "barbell+sphere*cylinder@hardsphere", 
     360        "cylinder@hardsphere*sphere+barbell") 
     361    test_models( 
     362        "barbell+cylinder@hardsphere*sphere", 
     363        "sphere*cylinder@hardsphere+barbell") 
     364 
     365    #Test the the model produces the parameters that we would expect 
     366    model = load_model("cylinder@hardsphere*sphere") 
     367    actual = [p.name for p in model.info.parameters.kernel_parameters] 
     368    target = ("sld sld_solvent radius length theta phi volfraction" 
     369              " A_sld A_sld_solvent A_radius").split() 
     370    assert target == actual, "%s != %s"%(target, actual) 
     371 
    338372if __name__ == "__main__": 
    339373    list_models_main() 

  • sasmodels/jitter.py

    r8cfb486 r0d5a655  
    11#!/usr/bin/env python 
    22""" 
    3 Application to explore the difference between sasview 3.x orientation 
    4 dispersity and possible replacement algorithms. 
     3Jitter Explorer 
     4=============== 
     5 
     6Application to explore orientation angle and angular dispersity. 
    57""" 
    68from __future__ import division, print_function 
    … …  
    124126def draw_parallelepiped(ax, size, view, jitter, steps=None, alpha=1): 
    125127    """Draw a parallelepiped.""" 
    126     a,b,c = size 
     128    a, b, c = size 
    127129    x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1]) 
    128130    y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1]) 
    … …  
    141143    x, y, z = transform_xyz(view, jitter, x, y, z) 
    142144    ax.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha) 
     145 
     146    # Draw pink face on box. 
     147    # Since I can't control face color, instead draw a thin box situated just 
     148    # in front of the "a+" face. 
     149    if 1: 
     150        x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1]) 
     151        y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1]) 
     152        z = c*np.array([ 1, 1, 1, 1,-1,-1,-1,-1]) 
     153        x, y, z = transform_xyz(view, jitter, abs(x)*1.05, y, z) 
     154        ax.plot_trisurf(x, y, triangles=tri, Z=z, color=[1,0.6,0.6], alpha=alpha) 
    143155 
    144156    draw_labels(ax, view, jitter, [ 
    … …  
    608620    Show an interactive orientation and jitter demo. 
    609621 
    610     *model_name* is one of the models available in :func:`select_model`. 
     622    *model_name* is one of: sphere, cylinder, ellipsoid, parallelepiped, 
     623    triaxial_ellipsoid, or bcc_paracrystal. 
     624 
     625    *size* gives the dimensions (a, b, c) of the shape. 
     626 
     627    *dist* is the type of dispersition: gaussian, rectangle, or uniform. 
     628 
     629    *mesh* is the number of points in the dispersion mesh. 
     630 
     631    *projection* is the map projection to use for the mesh: equirectangular, 
     632    sinusoidal, guyou, azimuthal_equidistance, or azimuthal_equal_area. 
    611633    """ 
    612634    # projection number according to 1-order position in list, but 

  • sasmodels/mixture.py

    r2d81cfe recf895e  
    9494            # If model is a sum model, each constituent model gets its own scale parameter 
    9595            scale_prefix = prefix 
    96             if prefix == '' and part.operation == '*': 
     96            if prefix == '' and getattr(part, "operation", '') == '*': 
    9797                # `part` is a composition product model. Find the prefixes of 
    98                 # it's parameters to form a new prefix for the scale, eg: 
    99                 # a model with A*B*C will have ABC_scale 
     98                # it's parameters to form a new prefix for the scale. 
     99                # For example, a model with A*B*C will have ABC_scale. 
    100100                sub_prefixes = [] 
    101101                for param in part.parameters.kernel_parameters: 
    … …  
    233233        return self 
    234234 
    235     def next(self): 
     235    def __next__(self): 
    236236        # type: () -> Tuple[List[Callable], CallDetails, np.ndarray] 
    237237        if self.part_num >= len(self.parts): 
    … …  
    251251 
    252252        return kernel, call_details, values 
     253 
     254    # CRUFT: py2 support 
     255    next = __next__ 
    253256 
    254257    def _part_details(self, info, par_index): 

  • sasmodels/modelinfo.py

    r108e70e r95498a3  
    6969        processed.append(parse_parameter(*p)) 
    7070    partable = ParameterTable(processed) 
     71    partable.check_angles() 
    7172    return partable 
    7273 
    … …  
    421422        # type: (List[Parameter]) -> None 
    422423        self.kernel_parameters = parameters 
    423         self._check_angles() 
    424424        self._set_vector_lengths() 
    425425 
    … …  
    471471        self.pd_2d = set(p.name for p in self.call_parameters if p.polydisperse) 
    472472 
    473     def _check_angles(self): 
     473    def check_angles(self): 
     474        """ 
     475        Check that orientation angles are theta, phi and possibly psi. 
     476        """ 
    474477        theta = phi = psi = -1 
    475478        for k, p in enumerate(self.kernel_parameters): 
    … …  
    494497            if psi >= 0 and psi != phi+1: 
    495498                raise TypeError("psi must follow phi") 
    496             #if (psi >= 0 and psi != last_par) or (psi < 0 and phi != last_par): 
    497             #    raise TypeError("orientation parameters must appear at the " 
    498             #                    "end of the parameter table") 
     499            if (psi >= 0 and psi != last_par) or (psi < 0 and phi != last_par): 
     500                raise TypeError("orientation parameters must appear at the " 
     501                                "end of the parameter table") 
    499502        elif theta >= 0 or phi >= 0 or psi >= 0: 
    500503            raise TypeError("oriented shapes must have both theta and phi and maybe psi") 
Note: See TracChangeset for help on using the changeset viewer.