Changeset bd91e8f in sasmodels


Ignore:
Timestamp:
Mar 30, 2019 4:59:53 PM (8 months ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
ticket_1156
Children:
345cc8f
Parents:
a3412a6 (diff), e15a822 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' into ticket_1156

Files:
2 added
150 edited

Legend:

Unmodified
Added
Removed
  • .travis.yml

    r5c36bf1 rf3767c2  
    99  - os: linux 
    1010    env: 
    11     - PY=3.6 
     11    - PY=3.7 
    1212  - os: osx 
    1313    language: generic 
     
    1717    language: generic 
    1818    env: 
    19     - PY=3.5 
     19    - PY=3.7 
    2020branches: 
    2121  only: 
  • doc/guide/index.rst

    rda5536f rbc69321  
    1212   pd/polydispersity.rst 
    1313   resolution.rst 
     14   plugin.rst 
     15   fitting_sq.rst 
    1416   magnetism/magnetism.rst 
    1517   orientation/orientation.rst 
    1618   sesans/sans_to_sesans.rst 
    1719   sesans/sesans_fitting.rst 
    18    plugin.rst 
    1920   scripting.rst 
    2021   refs.rst 
  • doc/guide/plugin.rst

    r9150036 re15a822  
    303303**Note: The order of the parameters in the definition will be the order of the 
    304304parameters in the user interface and the order of the parameters in Fq(), Iq(), 
    305 Iqac(), Iqabc(), form_volume() and shell_volume(). 
     305Iqac(), Iqabc(), radius_effective(), form_volume() and shell_volume(). 
    306306And** *scale* **and** *background* **parameters are implicit to all models, 
    307307so they do not need to be included in the parameter table.** 
     
    387387can take arbitrary values, even for integer parameters, so your model should 
    388388round the incoming parameter value to the nearest integer inside your model 
    389 you should round to the nearest integer.  In C code, you can do this using:: 
     389you should round to the nearest integer.  In C code, you can do this using: 
     390 
     391.. code-block:: c 
    390392 
    391393    static double 
     
    454456............. 
    455457 
     458.. note:: 
     459 
     460   Pure python models do not yet support direct computation of $<F(Q)^2>$ or 
     461   $<F(Q)>^2$. Neither do they support orientational distributions or magnetism 
     462   (use C models if these are required). 
     463 
    456464For pure python models, define the *Iq* function:: 
    457465 
     
    499507Models should define *form_volume(par1, par2, ...)* where the parameter 
    500508list includes the *volume* parameters in order.  This is used for a weighted 
    501 volume normalization so that scattering is on an absolute scale.  If 
    502 *form_volume* is not defined, then the default *form_volume = 1.0* will be 
    503 used. 
     509volume normalization so that scattering is on an absolute scale.  For 
     510solid shapes, the *I(q)* function should use *form_volume* squared 
     511as its scale factor.  If *form_volume* is not defined, then the default 
     512*form_volume = 1.0* will be used. 
    504513 
    505514Hollow shapes, where the volume fraction of particle corresponds to the 
    506515material in the shell rather than the volume enclosed by the shape, must 
    507516also define a *shell_volume(par1, par2, ...)* function.  The parameters 
    508 are the same as for *form_volume*.  The *I(q)* calculation should use 
    509 *shell_volume* squared as its scale factor for the volume normalization. 
    510 The structure factor calculation needs *form_volume* in order to properly 
    511 scale the volume fraction parameter, so both functions are required for 
    512 hollow shapes. 
    513  
    514 Note: Pure python models do not yet support direct computation of the 
    515 average of $F(q)$ and $F^2(q)$. 
     517are the same as for *form_volume*.  Here the *I(q)* function should use 
     518*shell_volume* squared instead of *form_volume* squared so that the scale 
     519parameter corresponds to the volume fraction of material in the sample. 
     520The structure factor calculation needs the volume fraction of the filled 
     521shapes for its calculation, so the volume fraction parameter in the model 
     522is automatically scaled by *form_volume/shell_volume* prior to calling the 
     523structure factor. 
    516524 
    517525Embedded C Models 
     
    524532    """ 
    525533 
    526 This expands into the equivalent C code:: 
     534This expands into the equivalent C code: 
     535 
     536.. code-block:: c 
    527537 
    528538    double Iq(double q, double par1, double par2, ...); 
     
    535545includes only the volume parameters. 
    536546 
    537 *form_volume* defines the volume of the shell for hollow shapes. As in 
     547*shell_volume* defines the volume of the shell for hollow shapes. As in 
    538548python models, it includes only the volume parameters. 
    539549 
     
    567577Rather than returning NAN from Iq, you must define the *INVALID(v)*.  The 
    568578*v* parameter lets you access all the parameters in the model using 
    569 *v.par1*, *v.par2*, etc. For example:: 
     579*v.par1*, *v.par2*, etc. For example: 
     580 
     581.. code-block:: c 
    570582 
    571583    #define INVALID(v) (v.bell_radius < v.radius) 
     
    582594 
    583595Instead of defining the *Iq* function, models can define *Fq* as 
    584 something like:: 
     596something like: 
     597 
     598.. code-block:: c 
    585599 
    586600    double Fq(double q, double *F1, double *F2, double par1, double par2, ...); 
     
    614628laboratory frame and beam travelling along $-z$. 
    615629 
    616 The oriented C model is called using *Iqabc(qa, qb, qc, par1, par2, ...)* where 
     630The oriented C model (oriented pure Python models are not supported) 
     631is called using *Iqabc(qa, qb, qc, par1, par2, ...)* where 
    617632*par1*, etc. are the parameters to the model.  If the shape is rotationally 
    618633symmetric about *c* then *psi* is not needed, and the model is called 
     
    642657 
    643658Using the $z, w$ values for Gauss-Legendre integration in "lib/gauss76.c", the 
    644 numerical integration is then:: 
     659numerical integration is then: 
     660 
     661.. code-block:: c 
    645662 
    646663    double outer_sum = 0.0; 
     
    718735to compute the proper magnetism and orientation, which you can implement 
    719736using *Iqxy(qx, qy, par1, par2, ...)*. 
     737 
     738**Note: Magnetism is not supported in pure Python models.** 
    720739 
    721740Special Functions 
     
    9831002  memory, and wrong answers computed. The conclusion from a very long and 
    9841003  strange debugging session was that any arrays that you declare in your 
    985   model should be a multiple of four. For example:: 
     1004  model should be a multiple of four. For example: 
     1005 
     1006  .. code-block:: c 
    9861007 
    9871008      double Iq(q, p1, p2, ...) 
     
    10151036structure factor is the *hardsphere* interaction, which 
    10161037uses the effective radius of the form factor as an input to the structure 
    1017 factor model.  The effective radius is the average radius of the 
    1018 form averaged over all the polydispersity values. 
    1019  
    1020 :: 
    1021  
    1022     def ER(radius, thickness): 
    1023         """Effective radius of a core-shell sphere.""" 
    1024         return radius + thickness 
    1025  
    1026 Now consider the *core_shell_sphere*, which has a simple effective radius 
    1027 equal to the radius of the core plus the thickness of the shell, as 
    1028 shown above. Given polydispersity over *(r1, r2, ..., rm)* in radius and 
    1029 *(t1, t2, ..., tn)* in thickness, *ER* is called with a mesh 
    1030 grid covering all possible combinations of radius and thickness. 
    1031 That is, *radius* is *(r1, r2, ..., rm, r1, r2, ..., rm, ...)* 
    1032 and *thickness* is *(t1, t1, ... t1, t2, t2, ..., t2, ...)*. 
    1033 The *ER* function returns one effective radius for each combination. 
    1034 The effective radius calculator weights each of these according to 
    1035 the polydispersity distributions and calls the structure factor 
    1036 with the average *ER*. 
    1037  
    1038 :: 
    1039  
    1040     def VR(radius, thickness): 
    1041         """Sphere and shell volumes for a core-shell sphere.""" 
    1042         whole = 4.0/3.0 * pi * (radius + thickness)**3 
    1043         core = 4.0/3.0 * pi * radius**3 
    1044         return whole, whole - core 
    1045  
    1046 Core-shell type models have an additional volume ratio which scales 
    1047 the structure factor.  The *VR* function returns the volume of 
    1048 the whole sphere and the volume of the shell. Like *ER*, there is 
    1049 one return value for each point in the mesh grid. 
    1050  
    1051 *NOTE: we may be removing or modifying this feature soon. As of the 
    1052 time of writing, core-shell sphere returns (1., 1.) for VR, giving a volume 
    1053 ratio of 1.0.* 
     1038factor model.  The effective radius is the weighted average over all 
     1039values of the shape in polydisperse systems. 
     1040 
     1041There can be many notions of effective radius, depending on the shape.  For 
     1042a sphere it is clearly just the radius, but for an ellipsoid of revolution 
     1043we provide average curvature, equivalent sphere radius, minimum radius and 
     1044maximum radius.  These options are listed as *radius_effective_modes* in 
     1045the python model defintion, and must be computed by the *radius_effective* 
     1046function which takes the *radius_effective_mode* parameter as an integer, 
     1047along with the various model parameters.  Unlike normal C/Python arrays, 
     1048the first mode is 1, the second is 2, etc.  Mode 0 indicates that the 
     1049effective radius from the shape is to be ignored in favour of the the 
     1050effective radius parameter in the structure factor model. 
     1051 
     1052 
     1053Consider the core-shell sphere, which defines the following effective radius 
     1054modes in the python model:: 
     1055 
     1056    radius_effective_modes = [ 
     1057        "outer radius", 
     1058        "core radius", 
     1059    ] 
     1060 
     1061and the following function in the C-file for the model: 
     1062 
     1063.. code-block:: c 
     1064 
     1065    static double 
     1066    radius_effective(int mode, double radius, double thickness) 
     1067    { 
     1068        switch (mode) { 
     1069            case 0: return radius + thickness; 
     1070            case 1: return radius; 
     1071            default: return 0.; 
     1072        } 
     1073    } 
     1074 
     1075    static double 
     1076    form_volume(double radius, double thickness) 
     1077    { 
     1078        return M_4PI_3 * cube(radius + thickness); 
     1079    } 
     1080 
     1081Given polydispersity over *(r1, r2, ..., rm)* in radius and *(t1, t2, ..., tn)* 
     1082in thickness, *radius_effective* is called over a mesh grid covering all 
     1083possible combinations of radius and thickness, with a single *(ri, tj)* pair 
     1084in each call. The weights each of these results according to the 
     1085polydispersity distributions and calls the structure factor with the average 
     1086effective radius.  Similarly, for *form_volume*. 
     1087 
     1088Hollow models have an additional volume ratio which is needed to scale the 
     1089structure factor.  The structure factor uses the volume fraction of the filled 
     1090particles as part of its density estimate, but the scale factor for the 
     1091scattering intensity (as non-solvent volume fraction / volume) is determined 
     1092by the shell volume only.  Therefore the *shell_volume* function is 
     1093needed to compute the form:shell volume ratio, which then scales the 
     1094*volfraction* parameter prior to calling the structure factor calculator. 
     1095In the case of a hollow sphere, this would be: 
     1096 
     1097.. code-block:: c 
     1098 
     1099    static double 
     1100    shell_volume(double radius, double thickness) 
     1101    { 
     1102        double whole = M_4PI_3 * cube(radius + thickness); 
     1103        double core = M_4PI_3 * cube(radius); 
     1104        return whole - core; 
     1105    } 
     1106 
     1107If *shell_volume* is not present, then *form_volume* and *shell_volume* are 
     1108assumed to be equal, and the shape is considered solid. 
    10541109 
    10551110Unit Tests 
     
    11081163and a check that the model runs. 
    11091164 
    1110 If you are not using sasmodels from SasView, skip this step. 
    1111  
    11121165Recommended Testing 
    11131166................... 
     1167 
     1168**NB: For now, this more detailed testing is only possible if you have a 
     1169SasView build environment available!** 
    11141170 
    11151171If the model compiles and runs, you can next run the unit tests that 
     
    12481304| 2016-10-25 Steve King 
    12491305| 2017-05-07 Paul Kienzle - Moved from sasview to sasmodels docs 
     1306| 2019-03-28 Paul Kienzle - Update docs for radius_effective and shell_volume 
  • doc/guide/resolution.rst

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

    rc462169 r66dbbfb  
    2727import os, sys 
    2828sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.realpath(__file__)))) 
     29import warnings 
    2930 
    3031import numpy as np 
     
    3637 
    3738import sasmodels.special as sp 
    38  
    39 # Need to parse shape early since it determines the kernel function 
    40 # that will be used for the various integrators 
    41 shape = 'parallelepiped' if len(sys.argv) < 2 else sys.argv[1] 
    42 Qstr = '0.005' if len(sys.argv) < 3 else sys.argv[2] 
    4339 
    4440DTYPE = 'd' 
     
    201197    return norm, Fq 
    202198 
    203 if shape == 'sphere': 
    204     RADIUS = 50  # integer for the sake of mpf 
    205     NORM, KERNEL = make_sphere(radius=RADIUS) 
    206     NORM_MP, KERNEL_MP = make_sphere(radius=RADIUS, env=MPenv) 
    207 elif shape == 'cylinder': 
    208     #RADIUS, LENGTH = 10, 100000 
    209     RADIUS, LENGTH = 10, 300  # integer for the sake of mpf 
    210     NORM, KERNEL = make_cylinder(radius=RADIUS, length=LENGTH) 
    211     NORM_MP, KERNEL_MP = make_cylinder(radius=RADIUS, length=LENGTH, env=MPenv) 
    212 elif shape == 'triaxial_ellipsoid': 
    213     #A, B, C = 4450, 14000, 47 
    214     A, B, C = 445, 140, 47  # integer for the sake of mpf 
    215     NORM, KERNEL = make_triaxial_ellipsoid(A, B, C) 
    216     NORM_MP, KERNEL_MP = make_triaxial_ellipsoid(A, B, C, env=MPenv) 
    217 elif shape == 'parallelepiped': 
    218     #A, B, C = 4450, 14000, 47 
    219     A, B, C = 445, 140, 47  # integer for the sake of mpf 
    220     NORM, KERNEL = make_parallelepiped(A, B, C) 
    221     NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv) 
    222 elif shape == 'core_shell_parallelepiped': 
    223     #A, B, C = 4450, 14000, 47 
    224     #A, B, C = 445, 140, 47  # integer for the sake of mpf 
    225     A, B, C = 114, 1380, 6800 
    226     DA, DB, DC = 21, 58, 2300 
    227     SLDA, SLDB, SLDC = "5", "-0.3", "11.5" 
    228     ## default parameters from sasmodels 
    229     #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 400,75,35,10,10,10,2,4,2 
    230     ## swap A-B-C to C-B-A 
    231     #A, B, C, DA, DB, DC, SLDA, SLDB, SLDC = C, B, A, DC, DB, DA, SLDC, SLDB, SLDA 
    232     #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 10,20,30,100,200,300,1,2,3 
    233     #SLD_SOLVENT,CONTRAST = 0, 4 
    234     if 1: # C shortest 
    235         B, C = C, B 
    236         DB, DC = DC, DB 
    237         SLDB, SLDC = SLDC, SLDB 
    238     elif 0: # C longest 
    239         A, C = C, A 
    240         DA, DC = DC, DA 
    241         SLDA, SLDC = SLDC, SLDA 
    242     #NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 
    243     NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 
    244     NORM_MP, KERNEL_MP = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC, env=MPenv) 
    245 elif shape == 'paracrystal': 
    246     LATTICE = 'bcc' 
    247     #LATTICE = 'fcc' 
    248     #LATTICE = 'sc' 
    249     DNN, D_FACTOR = 220, '0.06'  # mpmath needs to initialize floats from string 
    250     RADIUS = 40  # integer for the sake of mpf 
    251     NORM, KERNEL = make_paracrystal( 
    252         radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE) 
    253     NORM_MP, KERNEL_MP = make_paracrystal( 
    254         radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE, env=MPenv) 
    255 else: 
    256     raise ValueError("Unknown shape %r"%shape) 
     199NORM = 1.0  # type: float 
     200KERNEL = None  # type: CALLABLE[[ndarray, ndarray, ndarray], ndarray] 
     201NORM_MP = 1  # type: mpf 
     202KERNEL = None  # type: CALLABLE[[mpf, mpf, mpf], mpf] 
     203 
     204SHAPES = [ 
     205    'sphere', 
     206    'cylinder', 
     207    'triaxial_ellipsoid', 
     208    'parallelepiped', 
     209    'core_shell_parallelepiped', 
     210    'fcc_paracrystal', 
     211    'bcc_paracrystal', 
     212    'sc_paracrystal', 
     213] 
     214def build_shape(shape, **pars): 
     215    global NORM, KERNEL 
     216    global NORM_MP, KERNEL_MP 
     217 
     218    # Note: using integer or string defaults for the sake of mpf 
     219    if shape == 'sphere': 
     220        RADIUS = pars.get('radius', 50) 
     221        NORM, KERNEL = make_sphere(radius=RADIUS) 
     222        NORM_MP, KERNEL_MP = make_sphere(radius=RADIUS, env=MPenv) 
     223    elif shape == 'cylinder': 
     224        #RADIUS, LENGTH = 10, 100000 
     225        RADIUS = pars.get('radius', 10) 
     226        LENGTH = pars.get('radius', 300) 
     227        NORM, KERNEL = make_cylinder(radius=RADIUS, length=LENGTH) 
     228        NORM_MP, KERNEL_MP = make_cylinder(radius=RADIUS, length=LENGTH, env=MPenv) 
     229    elif shape == 'triaxial_ellipsoid': 
     230        #A, B, C = 4450, 14000, 47 
     231        A = pars.get('a', 445) 
     232        B = pars.get('b', 140) 
     233        C = pars.get('c', 47) 
     234        NORM, KERNEL = make_triaxial_ellipsoid(A, B, C) 
     235        NORM_MP, KERNEL_MP = make_triaxial_ellipsoid(A, B, C, env=MPenv) 
     236    elif shape == 'parallelepiped': 
     237        #A, B, C = 4450, 14000, 47 
     238        A = pars.get('a', 445) 
     239        B = pars.get('b', 140) 
     240        C = pars.get('c', 47) 
     241        NORM, KERNEL = make_parallelepiped(A, B, C) 
     242        NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv) 
     243    elif shape == 'core_shell_parallelepiped': 
     244        #A, B, C = 4450, 14000, 47 
     245        #A, B, C = 445, 140, 47  # integer for the sake of mpf 
     246        A = pars.get('a', 114) 
     247        B = pars.get('b', 1380) 
     248        C = pars.get('c', 6800) 
     249        DA = pars.get('da', 21) 
     250        DB = pars.get('db', 58) 
     251        DC = pars.get('dc', 2300) 
     252        SLDA = pars.get('slda', "5") 
     253        SLDB = pars.get('sldb', "-0.3") 
     254        SLDC = pars.get('sldc', "11.5") 
     255        ## default parameters from sasmodels 
     256        #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 400,75,35,10,10,10,2,4,2 
     257        ## swap A-B-C to C-B-A 
     258        #A, B, C, DA, DB, DC, SLDA, SLDB, SLDC = C, B, A, DC, DB, DA, SLDC, SLDB, SLDA 
     259        #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 10,20,30,100,200,300,1,2,3 
     260        #SLD_SOLVENT,CONTRAST = 0, 4 
     261        if 1: # C shortest 
     262            B, C = C, B 
     263            DB, DC = DC, DB 
     264            SLDB, SLDC = SLDC, SLDB 
     265        elif 0: # C longest 
     266            A, C = C, A 
     267            DA, DC = DC, DA 
     268            SLDA, SLDC = SLDC, SLDA 
     269        #NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 
     270        NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 
     271        NORM_MP, KERNEL_MP = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC, env=MPenv) 
     272    elif shape.endswith('paracrystal'): 
     273        LATTICE, _ = shape.split('_') 
     274        DNN = pars.get('dnn', 220) 
     275        D_FACTOR = pars.get('d_factor', '0.06') 
     276        RADIUS = pars.get('radius', 40) 
     277        NORM, KERNEL = make_paracrystal( 
     278            radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE) 
     279        NORM_MP, KERNEL_MP = make_paracrystal( 
     280            radius=RADIUS, dnn=DNN, d_factor=D_FACTOR, lattice=LATTICE, env=MPenv) 
     281    else: 
     282        raise ValueError("Unknown shape %r"%shape) 
    257283 
    258284# Note: hardcoded in mp_quad 
     
    273299 
    274300# 2D integration functions 
    275 def mp_quad_2d(q, shape): 
     301def mp_quad_2d(q): 
    276302    evals = [0] 
    277303    def integrand(theta, phi): 
     
    393419    print("simpson-%d"%n, n**2, simps(simps(Zq, dx=dx), dx=dy)*SCALE/(4*pi)) 
    394420    print("romb-%d"%n, n**2, romb(romb(Zq, dx=dx), dx=dy)*SCALE/(4*pi)) 
     421 
     422def quadpy_method(q, rule): 
     423    """ 
     424    Use *rule*="name:index" where name and index are chosen from below. 
     425 
     426    Available rule names and the corresponding indices:: 
     427 
     428        AlbrechtCollatz: [1-5] 
     429        BazantOh: 9, 11, 13 
     430        HeoXu: 13, 15, 17, 19-[1-2], 21-[1-6], 23-[1-3], 25-[1-2], 27-[1-3], 
     431            29, 31, 33, 35, 37, 39-[1-2] 
     432        FliegeMaier: 4, 9, 16, 25 
     433        Lebedev: 3[a-c], 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 35, 
     434            41, 47, 53, 59, 65, 71, 77 83, 89, 95, 101, 107, 113, 119, 125, 131 
     435        McLaren: [1-10] 
     436        Stroud: U3 3-1, U3 5-[1-5], U3 7-[1-2], U3 8-1, U3 9-[1-3], 
     437            U3 11-[1-3], U3 14-1 
     438    """ 
     439    try: 
     440        import quadpy 
     441    except ImportError: 
     442        warnings.warn("use 'pip install quadpy' to enable quadpy.sphere tests") 
     443        return 
     444 
     445    from quadpy.sphere import (AlbrechtCollatz, BazantOh, HeoXu, 
     446        FliegeMaier, Lebedev, McLaren, Stroud, integrate_spherical) 
     447    RULES = { 
     448        'AlbrechtCollatz': AlbrechtCollatz, 
     449        'BazantOh': BazantOh, 
     450        'HeoXu': HeoXu, 
     451        'FliegeMaier': FliegeMaier, 
     452        'Lebedev': Lebedev, 
     453        'McLaren': McLaren, 
     454        'Stroud': Stroud, 
     455    } 
     456    int_index = 'AlbrechtCollatz', 'McLaren' 
     457 
     458    rule_name, rule_index = rule.split(':') 
     459    index = int(rule_index) if rule_name in int_index else rule_index 
     460    rule_obj = RULES[rule_name](index) 
     461    fn = lambda azimuthal, polar: kernel_2d(q=q, theta=polar, phi=azimuthal) 
     462    Iq = integrate_spherical(fn, rule=rule_obj)/(4*pi) 
     463    print("%s degree=%d points=%s => %.15g" 
     464          % (rule, rule_obj.degree, len(rule_obj.points), Iq)) 
    395465 
    396466def plot_2d(q, n=300): 
     
    414484    pylab.show() 
    415485 
    416 def main(Qstr): 
    417     Q = float(Qstr) 
    418     if shape == 'sphere': 
     486def main(): 
     487    import argparse 
     488 
     489    parser = argparse.ArgumentParser( 
     490        description="asymmetric integration explorer", 
     491        formatter_class=argparse.ArgumentDefaultsHelpFormatter, 
     492        ) 
     493    parser.add_argument('-s', '--shape', choices=SHAPES, 
     494                        default='parallelepiped', 
     495                        help='oriented shape') 
     496    parser.add_argument('-q', '--q_value', type=str, default='0.005', 
     497                        help='Q value to evaluate') 
     498    parser.add_argument('pars', type=str, nargs='*', default=[], 
     499                        help='p=val for p in shape parameters') 
     500    opts = parser.parse_args() 
     501    pars = {k: v for par in opts.pars for k, v in [par.split('=')]} 
     502    build_shape(opts.shape, **pars) 
     503 
     504    Q = float(opts.q_value) 
     505    if opts.shape == 'sphere': 
    419506        print("exact", NORM*sp.sas_3j1x_x(Q*RADIUS)**2) 
    420     print("gauss-20", *gauss_quad_2d(Q, n=20)) 
    421     print("gauss-76", *gauss_quad_2d(Q, n=76)) 
    422     print("gauss-150", *gauss_quad_2d(Q, n=150)) 
    423     print("gauss-500", *gauss_quad_2d(Q, n=500)) 
    424     print("gauss-1025", *gauss_quad_2d(Q, n=1025)) 
    425     print("gauss-2049", *gauss_quad_2d(Q, n=2049)) 
    426     print("gauss-20 usub", *gauss_quad_usub(Q, n=20)) 
    427     print("gauss-76 usub", *gauss_quad_usub(Q, n=76)) 
    428     print("gauss-150 usub", *gauss_quad_usub(Q, n=150)) 
     507 
     508    # Methods from quadpy, if quadpy is available 
     509    #  AlbrechtCollatz: [1-5] 
     510    #  BazantOh: 9, 11, 13 
     511    #  HeoXu: 13, 15, 17, 19-[1-2], 21-[1-6], 23-[1-3], 25-[1-2], 27-[1-3], 
     512    #     29, 31, 33, 35, 37, 39-[1-2] 
     513    #  FliegeMaier: 4, 9, 16, 25 
     514    #  Lebedev: 3[a-c], 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 35, 
     515    #     41, 47, 53, 59, 65, 71, 77 83, 89, 95, 101, 107, 113, 119, 125, 131 
     516    #  McLaren: [1-10] 
     517    #  Stroud: U3 3-1, U3 5-[1-5], U3 7-[1-2], U3 8-1, U3 9-[1-3], 
     518    #     U3 11-[1-3], U3 14-1 
     519    quadpy_method(Q, "AlbrechtCollatz:5") 
     520    quadpy_method(Q, "HeoXu:39-2") 
     521    quadpy_method(Q, "FliegeMaier:25") 
     522    quadpy_method(Q, "Lebedev:19") 
     523    quadpy_method(Q, "Lebedev:131") 
     524    quadpy_method(Q, "McLaren:10") 
     525    quadpy_method(Q, "Stroud:U3 14-1") 
     526 
     527    print("gauss-20 points=%d => %.15g" % gauss_quad_2d(Q, n=20)) 
     528    print("gauss-76 points=%d => %.15g" % gauss_quad_2d(Q, n=76)) 
     529    print("gauss-150 points=%d => %.15g" % gauss_quad_2d(Q, n=150)) 
     530    print("gauss-500 points=%d => %.15g" % gauss_quad_2d(Q, n=500)) 
     531    print("gauss-1025 points=%d => %.15g" % gauss_quad_2d(Q, n=1025)) 
     532    print("gauss-2049 points=%d => %.15g" % gauss_quad_2d(Q, n=2049)) 
     533    print("gauss-20 usub points=%d => %.15g" % gauss_quad_usub(Q, n=20)) 
     534    print("gauss-76 usub points=%d => %.15g" % gauss_quad_usub(Q, n=76)) 
     535    print("gauss-150 usub points=%d => %.15g" % gauss_quad_usub(Q, n=150)) 
     536 
    429537    #gridded_2d(Q, n=2**8+1) 
    430538    gridded_2d(Q, n=2**10+1) 
    431539    #gridded_2d(Q, n=2**12+1) 
    432540    #gridded_2d(Q, n=2**15+1) 
    433     if shape not in ('paracrystal', 'core_shell_parallelepiped'): 
    434         # adaptive forms on models for which the calculations are fast enough 
     541    # adaptive forms on models for which the calculations are fast enough 
     542    SLOW_SHAPES = { 
     543        'fcc_paracrystal', 'bcc_paracrystal', 'sc_paracrystal', 
     544        'core_shell_parallelepiped', 
     545    } 
     546    if opts.shape not in SLOW_SHAPES: 
    435547        print("dblquad", *scipy_dblquad_2d(Q)) 
    436548        print("semi-romberg-100", *semi_romberg_2d(Q, n=100)) 
    437549        print("romberg", *scipy_romberg_2d(Q)) 
    438550        with mp.workprec(100): 
    439             print("mpmath", *mp_quad_2d(mp.mpf(Qstr), shape)) 
     551            print("mpmath", *mp_quad_2d(mp.mpf(opts.q_value))) 
    440552    plot_2d(Q, n=200) 
    441553 
    442554if __name__ == "__main__": 
    443     main(Qstr) 
     555    main() 
  • explore/beta/sasfit_compare.py

    r119073a ra34b811  
    408408        #radius_effective=12.59921049894873, 
    409409        ) 
    410     target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars) 
     410    target = sasmodels_theory(q, model, radius_effective_mode=0, structure_factor_mode=1, **pars) 
    411411    actual = ellipsoid_pe(q, norm='sasview', **pars) 
    412412    title = " ".join(("sasmodels", model, pd_type)) 
  • explore/beta/sasfit_compare_new.py

    r2586ab72 ra34b811  
    2121#Used to calculate F(q) for the cylinder, sphere, ellipsoid models 
    2222# RKH adding vesicle and hollow_cylinder to test sasview special cases of ER and VR 
    23 # BEWARE there are issues here if you run this from python3 (i.e. qt5 version of sasview), so do "source activate sasview" 
     23# There were issues here from python3 (i.e. qt5 version of sasview),fixed by Paul K (else do "source activate sasview") 
    2424def sas_sinx_x(x): 
    2525    with np.errstate(all='ignore'): 
     
    275275    if radius_effective is None: 
    276276        radius_effective = radius 
    277         print("radius_effective  ",radius_effective) 
    278277        average_volume = total_volume/total_weight 
    279278    if norm == 'sasfit': 
     
    285284        F2 *= 1e-12 
    286285        IQD = F2/average_volume*1e8*volfraction 
     286    print("in sphere_r : radius_effective  ",radius_effective,"  volfraction",volfraction) 
    287287    SQ = hardsphere_simple(q, radius_effective, volfraction) 
    288288    beta = F1**2/F2 
     
    332332    if radius_effective is None: 
    333333        radius_effective = radius_eff/total_weight 
    334         print("radius_effective  ",radius_effective) 
     334        print("in vesicle_pe : radius_effective  ",radius_effective) 
    335335    if norm == 'sasfit': 
    336336        IQD = F2 
     
    348348        F2 *= 1e-12 
    349349        IQD = F2/average_volume*1e8*volfraction 
    350     # RKH SORT THIS OUT WHEN HAVE NEW MODEL 
    351     vfff=volfraction*(30./20.)**3 
    352     print("radius_effective, fudged volfaction for S(Q) ",radius_effective,vfff) 
     350    # RKH SORT THIS OUT WHEN HAVE NEW MODEL  Vtot/Vshell = (R+T)^3/ ( (R+T)^3 -R^3 ) 
     351    vfff=volfraction*( 1.0/(1.0 - ( (radius/(radius+thickness))**3 ))) 
     352    print("in vesicle_pe : radius_effective, fudged volfraction for S(Q) ",radius_effective,vfff) 
    353353    SQ = hardsphere_simple(q, radius_effective, vfff) 
    354354    beta = F1**2/F2 
     
    359359# 
    360360#polydispersity for hollow_cylinder 
    361 def hollow_cylinder_pe(q,radius=20, thickness=10, sld=4, sld_solvent=1, volfraction=0.15, 
    362         radius_pd=0.1, thickness_pd=0.2, length_pd=0.05, radius_pd_type="gaussian", thickness_pd_type="gaussian", 
    363         radius_effective=None, background=0, scale=1, 
     361def hollow_cylinder_pe(q,radius=20, thickness=10, length=80, sld=4, sld_solvent=1, volfraction=0.15, 
     362        radius_pd=0.1, thickness_pd=0.2, length_pd=0.05, radius_pd_type="gaussian", length_pd_type="gaussian", 
     363        thickness_pd_type="gaussian", radius_effective=None, background=0, scale=1, 
    364364                 norm='sasview'): 
    365365    if norm not in ['sasview', 'sasfit', 'yun']: 
     
    374374        T_val, T_prob = schulz_distribution(thickness, thickness_pd*thickness, 0, inf) 
    375375    if length_pd_type == 'gaussian': 
    376         T_val, T_prob = gaussian_distribution(length, length_pd*length, 0, inf) 
     376        L_val, L_prob = gaussian_distribution(length, length_pd*length, 0, inf) 
    377377    elif length_pd_type == 'schulz': 
    378378        L_val, L_prob = schulz_distribution(length, length_pd*length, 0, inf) 
     
    402402    if radius_effective is None: 
    403403        radius_effective = radius_eff/total_weight 
    404         print("radius_effective  ",radius_effective) 
     404        print("in hollow_cylinder : radius_effective  ",radius_effective) 
    405405    if norm == 'sasfit': 
    406406        IQD = F2 
     
    420420# RKH SORT THIS OUT WHEN HAVE NEW MODEL 
    421421    vfff=volfraction 
    422     print("radius_effective, volfaction for S(Q) ",radius_effective,vfff) 
     422    print("in hollow_cylinder : radius_effective, volfaction for S(Q) ",radius_effective,vfff) 
    423423    SQ = hardsphere_simple(q, radius_effective, vfff) 
    424424    beta = F1**2/F2 
     
    548548    #Iq = None#= Sq = None 
    549549    r = dict(I._kernel.results()) 
    550     #print(r) 
     550    print(r) 
    551551    return Theory(Q=q, F1=None, F2=None, P=Pq, S=None, I=None, Seff=r["S_eff(Q)"], Ibeta=Iq) 
    552552 
     
    591591        radius_pd=.1, radius_pd_type=pd_type, 
    592592        volfraction=0.15, 
    593         radius_effective=12.59921049894873,  # equivalent average sphere radius 
     593        radius_effective=20  # equivalent average sphere radius 
     594#       NOTE sasview computes its own radius_effective in "target" (the print(r) at end of sasmodels_theory will reveal its value), 
     595#       this one is only used locally for "actual" 
    594596        ) 
    595     target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars) 
     597    target = sasmodels_theory(q, model, radius_effective_mode=0, structure_factor_mode=1, **pars) 
    596598    actual = sphere_r(q, norm='sasview', **pars) 
    597599    title = " ".join(("sasmodels", model, pd_type)) 
     
    604606    q = np.logspace(-5, 0, 250) 
    605607    model = 'vesicle' 
    606     print("F2 seems OK, but big issues with S(Q), wo work in progress") 
    607 # NOTE parameers for vesicle model are soon to be changed to remove "volfraction" 
     608    print("F2 seems OK, but big issues with S(Q), so work in progress") 
     609# NOTE parameters for vesicle model are soon to be changed to remove "volfraction" (though still there in 5.0) 
    608610    pars = dict( 
    609611        radius=20, thickness=10, sld=4, sld_solvent=1, volfraction=0.03, 
    610612        background=0, 
    611         radius_pd=0.1, thickness_pd=0.0, radius_pd_type=pd_type, thickness_pd_type=pd_type 
    612         #,radius_effective=12.59921049894873,  # equivalent average sphere radius 
    613         ) 
    614     target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars) 
     613        radius_pd=0.01, thickness_pd=0.01, radius_pd_type=pd_type, thickness_pd_type=pd_type, 
     614        radius_effective=30. ) 
     615        # equivalent average sphere radius for local "actual" to match what sasview uses, use None to compute average outer radius here, 
     616 
     617    target = sasmodels_theory(q, model, radius_effective_mode=0, structure_factor_mode=1, **pars) 
    615618    actual = vesicle_pe(q, norm='sasview', **pars) 
    616619    title = " ".join(("sasmodels", model, pd_type)) 
     
    629632        radius=20, thickness=10, length=80, sld=4, sld_solvent=1, 
    630633        background=0, 
    631         radius_pd=0.1, thickness_pd=0.0, length_pd=0.0, radius_pd_type=pd_type, thickness_pd_type=pd_type, length_pd_type=pd_type 
    632         #,radius_effective=12.59921049894873,  # equivalent average sphere radius 
    633         ) 
    634     target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars) 
     634        radius_pd=0.1, thickness_pd=0.0, length_pd=0.0, radius_pd_type=pd_type, thickness_pd_type=pd_type, length_pd_type=pd_type, 
     635        radius_effective=40.687) 
     636        # equivalent average sphere radius for local "actual" to match what sasview uses 
     637    target = sasmodels_theory(q, model, radius_effective_mode=0, structure_factor_mode=1, **pars) 
    635638    actual = hollow_cylinder_pe(q, norm='sasview', **pars) 
    636639# RKH monodisp was OK,    actual = hollow_cylinder_theta(q,radius=20, thickness=10, length=80, sld=4, sld_solvent=1 ) 
     
    651654        volfraction=0.15, 
    652655        radius_effective=270.7543927018, 
    653         #radius_effective=12.59921049894873, 
     656# if change radius_effective to some other value, the S(Q) from sasview does not agree 
    654657        ) 
    655     target = sasmodels_theory(q, model, effective_radius_mode=0, structure_factor_mode=1, **pars) 
     658    target = sasmodels_theory(q, model, radius_effective_mode=0, structure_factor_mode=1, **pars) 
    656659    actual = ellipsoid_pe(q, norm='sasview', **pars) 
    657660# RKH test       actual = ellipsoid_theta(q, radius_polar=20, radius_equatorial=400, sld=4, sld_solvent=1, volfraction=0.15, radius_effective=270.) 
     
    751754    } 
    752755 
    753     Q, IQ = load_sasfit(data_file('richard_test.txt')) 
    754     Q, IQSD = load_sasfit(data_file('richard_test2.txt')) 
    755     Q, IQBD = load_sasfit(data_file('richard_test3.txt')) 
     756    Q, IQ = load_sasfit(data_file('sasfit_sphere_schulz_IQD.txt')) 
     757    Q, IQSD = load_sasfit(data_file('sasfit_sphere_schulz_IQSD.txt')) 
     758    Q, IQBD = load_sasfit(data_file('sasfit_sphere_schulz_IQBD.txt')) 
    756759    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD) 
    757760    actual = sphere_r(Q, norm="sasfit", **pars) 
     
    772775    } 
    773776 
    774     Q, IQ = load_sasfit(data_file('richard_test4.txt')) 
    775     Q, IQSD = load_sasfit(data_file('richard_test5.txt')) 
    776     Q, IQBD = load_sasfit(data_file('richard_test6.txt')) 
     777    Q, IQ = load_sasfit(data_file('sasfit_ellipsoid_schulz_IQD.txt')) 
     778    Q, IQSD = load_sasfit(data_file('sasfit_ellipsoid_schulz_IQSD.txt')) 
     779    Q, IQBD = load_sasfit(data_file('sasfit_ellipsoid_schulz_IQBD.txt')) 
    777780    target = Theory(Q=Q, F1=None, F2=None, P=IQ, S=None, I=IQSD, Seff=None, Ibeta=IQBD) 
    778781    actual = ellipsoid_pe(Q, norm="sasfit", **pars) 
  • extra/build_linux.sh

    rdb8756e r0d362b7  
    1414# Run tests 
    1515STATUS=0 
    16 python -m sasmodels.model_test opencl_and_dll all || STATUS=$? 
     16python -m sasmodels.model_test -v all || STATUS=$? 
    1717python -m sasmodels.resolution || STATUS=$? 
    1818 
  • extra/pylint.rc

    r6e68289 rdf20715  
    6767    locally-enabled, 
    6868    locally-disabled, 
     69    #no-else-return, 
     70    #invalid-name, 
     71    fixme, 
     72    too-many-lines, 
     73    too-many-locals, 
     74    too-many-branches, 
     75    too-many-statements, 
     76    #too-many-instance-attributes, 
     77    #too-many-return-statements, 
     78    #no-self-use, 
    6979 
    7080[REPORTS] 
     
    109119good-names=_, 
    110120    input, id, 
    111     h, i, j, k, n, p, v, w, x, y, z, 
    112     q, qx, qy, qz, Q, Qx, Qy, Qz, 
    113     dt, dx, dy, dz, dq, 
     121    h, i, j, k, n, p, u, v, w, x, y, z, r, dr, 
     122    q, qx, qy, qz, Q, Qx, Qy, Qz, nx, ny, nz, nq, nr, 
     123    dt, dx, dy, dz, dq, S, P, a, b, c, 
    114124    Iq, dIq, Iqxy, Iq_calc, 
    115125    ER, call_ER, VR, call_VR, 
     
    133143 
    134144# Regular expression matching correct variable names 
    135 variable-rgx=[a-z_][a-z0-9_{2D}{1D}]{2,30}$ 
     145variable-rgx=[a-z_][a-z0-9_{2D}{1D}]{1,30}$ 
    136146 
    137147# Naming hint for variable names 
     
    151161 
    152162# Regular expression matching correct argument names 
    153 argument-rgx=[a-z_][a-z0-9_]{2,30}$ 
     163argument-rgx=[a-z_][a-z0-9_]{1,30}$ 
    154164 
    155165# Naming hint for argument names 
     
    280290# (useful for modules/projects where namespaces are manipulated during runtime 
    281291# and thus existing member attributes cannot be deduced by static analysis 
    282 ignored-modules=bumps,sas,numpy,numpy.random,scipy,scipy.special 
     292ignored-modules=bumps,sas,numpy,numpy.random,scipy,scipy.special,pyopencl,pycuda,matplotlib.cm,ipyvolume,ipywidget 
    283293 
    284294# List of classes names for which member attributes should not be checked 
  • sasmodels/bumps_model.py

    r2c4a190 rb297ba9  
    2626    from .kernel import KernelModel 
    2727    from .modelinfo import ModelInfo 
     28    from .resolution import Resolution 
    2829    Data = Union[Data1D, Data2D] 
    2930except ImportError: 
     
    180181    @property 
    181182    def resolution(self): 
     183        # type: () -> Union[None, Resolution] 
     184        """ 
     185        :class:`sasmodels.Resolution` applied to the data, if any. 
     186        """ 
    182187        return self._resolution 
    183188 
    184189    @resolution.setter 
    185190    def resolution(self, value): 
     191        # type: (Resolution) -> None 
     192        """ 
     193        :class:`sasmodels.Resolution` applied to the data, if any. 
     194        """ 
    186195        self._resolution = value 
    187196 
     
    270279        Save the model parameters and data into a file. 
    271280 
    272         Not Implemented. 
     281        Not Implemented except for sesans fits. 
    273282        """ 
    274283        if self.data_type == "sesans": 
  • sasmodels/compare.py

    rc1799d3 rb297ba9  
    725725        set_integration_size(model_info, ngauss) 
    726726 
    727     if (dtype != "default" and not dtype.endswith('!')  
     727    if (dtype != "default" and not dtype.endswith('!') 
    728728            and not (kernelcl.use_opencl() or kernelcuda.use_cuda())): 
    729729        raise RuntimeError("OpenCL not available " + kernelcl.OPENCL_ERROR) 
     
    772772        opts['pars'] = parse_pars(opts, maxdim=maxdim) 
    773773        if opts['pars'] is None: 
    774             return 
     774            return limits 
    775775        result = run_models(opts, verbose=True) 
    776776        if opts['plot']: 
     
    822822    args = dict((k, v) for k, v in args.items() 
    823823                if "_pd" not in k 
    824                    and ":" not in k 
    825                    and k not in ("background", "scale", "theta", "phi", "psi")) 
     824                and ":" not in k 
     825                and k not in ("background", "scale", "theta", "phi", "psi")) 
    826826    args = args.copy() 
    827827 
     
    896896    # work with trimmed data, not the full set 
    897897    sorted_err = np.sort(abs(err.compressed())) 
    898     if len(sorted_err) == 0: 
     898    if sorted_err.size == 0: 
    899899        print(label + "  no valid values") 
    900900        return 
  • sasmodels/convert.py

    rb3f4831 rbd91e8f  
    105105                    translation[newid+str(k)] = oldid+str(k) 
    106106    # Remove control parameter from the result 
    107     if model_info.control: 
    108         translation[model_info.control] = "CONTROL" 
     107    control_pars = [p.id for p in model_info.parameters.kernel_parameters 
     108                    if p.is_control] 
     109    if control_pars: 
     110        control_id = control_pars[0] 
     111        translation[control_id] = "CONTROL" 
    109112    return translation 
    110113 
     
    388391 
    389392def revert_name(model_info): 
     393    """Translate model name back to the name used in SasView 3.x""" 
    390394    oldname, _ = CONVERSION_TABLE.get(model_info.id, [None, {}]) 
    391395    return oldname 
     
    657661 
    658662def test_backward_forward(): 
     663    """ 
     664    Test conversion of model parameters from 4.x to 3.x and back. 
     665    """ 
    659666    from .core import list_models 
    660667    L = lambda name: _check_one(name, seed=1) 
  • sasmodels/core.py

    r07646b6 r9562dd2  
    6464        * all: all models 
    6565        * py: python models only 
    66         * c: compiled models only 
    67         * single: models which support single precision 
    68         * double: models which require double precision 
    69         * opencl: controls if OpenCL is supperessed 
    70         * 1d: models which are 1D only, or 2D using abs(q) 
    71         * 2d: models which can be 2D 
    72         * magnetic: models with an sld 
    73         * nommagnetic: models without an sld 
     66        * c: c models only 
     67        * single: c models which support single precision 
     68        * double: c models which require double precision 
     69        * opencl: c models which run in opencl 
     70        * dll: c models which do not run in opencl 
     71        * 1d: models without orientation 
     72        * 2d: models with orientation 
     73        * magnetic: models supporting magnetic sld 
     74        * nommagnetic: models without magnetic parameter 
    7475 
    7576    For multiple conditions, combine with plus.  For example, *c+single+2d* 
     
    9596    info = load_model_info(name) 
    9697    pars = info.parameters.kernel_parameters 
    97     if kind == "py" and callable(info.Iq): 
    98         return True 
    99     elif kind == "c" and not callable(info.Iq): 
    100         return True 
    101     elif kind == "double" and not info.single: 
    102         return True 
    103     elif kind == "single" and info.single: 
    104         return True 
    105     elif kind == "opencl" and info.opencl: 
    106         return True 
    107     elif kind == "2d" and any(p.type == 'orientation' for p in pars): 
    108         return True 
    109     elif kind == "1d" and all(p.type != 'orientation' for p in pars): 
    110         return True 
    111     elif kind == "magnetic" and any(p.type == 'sld' for p in pars): 
    112         return True 
    113     elif kind == "nonmagnetic" and any(p.type != 'sld' for p in pars): 
    114         return True 
     98    # TODO: may be adding Fq to the list at some point 
     99    is_pure_py = callable(info.Iq) 
     100    if kind == "py": 
     101        return is_pure_py 
     102    elif kind == "c": 
     103        return not is_pure_py 
     104    elif kind == "double": 
     105        return not info.single and not is_pure_py 
     106    elif kind == "single": 
     107        return info.single and not is_pure_py 
     108    elif kind == "opencl": 
     109        return info.opencl 
     110    elif kind == "dll": 
     111        return not info.opencl and not is_pure_py 
     112    elif kind == "2d": 
     113        return any(p.type == 'orientation' for p in pars) 
     114    elif kind == "1d": 
     115        return all(p.type != 'orientation' for p in pars) 
     116    elif kind == "magnetic": 
     117        return any(p.type == 'sld' for p in pars) 
     118    elif kind == "nonmagnetic": 
     119        return not any(p.type == 'sld' for p in pars) 
    115120    return False 
    116121 
     
    317322    return numpy_dtype, fast, platform 
    318323 
    319 def list_models_main(): 
    320     # type: () -> None 
    321     """ 
    322     Run list_models as a main program.  See :func:`list_models` for the 
    323     kinds of models that can be requested on the command line. 
    324     """ 
    325     import sys 
    326     kind = sys.argv[1] if len(sys.argv) > 1 else "all" 
    327     print("\n".join(list_models(kind))) 
    328  
    329324def test_composite_order(): 
     325    """ 
     326    Check that mixture models produce the same result independent of ordder. 
     327    """ 
    330328    def test_models(fst, snd): 
    331329        """Confirm that two models produce the same parameters""" 
     
    342340 
    343341    def build_test(first, second): 
     342        """Construct pair model test""" 
    344343        test = lambda description: test_models(first, second) 
    345344        description = first + " vs. " + second 
     
    377376    # type: () -> None 
    378377    """Check that model load works""" 
     378    from .product import RADIUS_ID, VOLFRAC_ID, STRUCTURE_MODE_ID, RADIUS_MODE_ID 
    379379    #Test the the model produces the parameters that we would expect 
    380380    model = load_model("cylinder@hardsphere*sphere") 
    381381    actual = [p.name for p in model.info.parameters.kernel_parameters] 
    382     target = ("sld sld_solvent radius length theta phi" 
    383               " radius_effective volfraction " 
    384               " structure_factor_mode radius_effective_mode" 
    385               " A_sld A_sld_solvent A_radius").split() 
     382    target = ["sld", "sld_solvent", "radius", "length", "theta", "phi", 
     383              RADIUS_ID, VOLFRAC_ID, STRUCTURE_MODE_ID, RADIUS_MODE_ID, 
     384              "A_sld", "A_sld_solvent", "A_radius"] 
    386385    assert target == actual, "%s != %s"%(target, actual) 
     386 
     387def list_models_main(): 
     388    # type: () -> int 
     389    """ 
     390    Run list_models as a main program.  See :func:`list_models` for the 
     391    kinds of models that can be requested on the command line. 
     392    """ 
     393    import sys 
     394    kind = sys.argv[1] if len(sys.argv) > 1 else "all" 
     395    try: 
     396        models = list_models(kind) 
     397        print("\n".join(models)) 
     398    except Exception: 
     399        print(list_models.__doc__) 
     400        return 1 
     401    return 0 
    387402 
    388403if __name__ == "__main__": 
  • sasmodels/data.py

    rc88f983 rb297ba9  
    500500            mdata = masked_array(data.y, data.mask.copy()) 
    501501            mdata[~np.isfinite(mdata)] = masked 
    502             if view is 'log': 
     502            if view == 'log': 
    503503                mdata[mdata <= 0] = masked 
    504504            plt.errorbar(data.x, scale*mdata, yerr=data.dy, fmt='.') 
     
    514514            mtheory = masked_array(theory) 
    515515            mtheory[~np.isfinite(mtheory)] = masked 
    516             if view is 'log': 
     516            if view == 'log': 
    517517                mtheory[mtheory <= 0] = masked 
    518518            plt.plot(theory_x, scale*mtheory, '-') 
  • sasmodels/details.py

    r885753a rb297ba9  
    236236    npars = kernel.info.parameters.npars 
    237237    nvalues = kernel.info.parameters.nvalues 
    238     scalars = [value for value, _dispersity, _weight in mesh] 
     238    scalars = [value for value, dispersity, weight in mesh] 
    239239    # skipping scale and background when building values and weights 
    240     _values, dispersity, weights = zip(*mesh[2:npars+2]) if npars else ((), (), ()) 
    241     #weights = correct_theta_weights(kernel.info.parameters, dispersity, weights) 
    242     length = np.array([len(w) for w in weights]) 
     240    value, dispersity, weight = zip(*mesh[2:npars+2]) if npars else ((), (), ()) 
     241    #weight = correct_theta_weights(kernel.info.parameters, dispersity, weight) 
     242    length = np.array([len(w) for w in weight]) 
    243243    offset = np.cumsum(np.hstack((0, length))) 
    244244    call_details = make_details(kernel.info, length, offset[:-1], offset[-1]) 
     
    246246    data_len = nvalues + 2*sum(len(v) for v in dispersity) 
    247247    extra = (32 - data_len%32)%32 
    248     data = np.hstack((scalars,) + dispersity + weights + ZEROS[:extra]) 
     248    data = np.hstack((scalars,) + dispersity + weight + ZEROS[:extra]) 
    249249    data = data.astype(kernel.dtype) 
    250250    is_magnetic = convert_magnetism(kernel.info.parameters, data) 
  • sasmodels/direct_model.py

    r9150036 re220acc  
    3131from . import resolution2d 
    3232from .details import make_kernel_args, dispersion_mesh 
    33 from .modelinfo import DEFAULT_BACKGROUND 
    3433from .product import RADIUS_MODE_ID 
    3534 
     
    6059    """ 
    6160    mesh = get_mesh(calculator.info, pars, dim=calculator.dim, mono=mono) 
    62     #print("pars", list(zip(*mesh))[0]) 
     61    #print("in call_kernel: pars:", list(zip(*mesh))[0]) 
    6362    call_details, values, is_magnetic = make_kernel_args(calculator, mesh) 
    64     #print("values:", values) 
     63    #print("in call_kernel: values:", values) 
    6564    return calculator(call_details, values, cutoff, is_magnetic) 
    6665 
     
    7372 
    7473    Use parameter *radius_effective_mode* to select the effective radius 
    75     calculation. 
     74    calculation to use amongst the *radius_effective_modes* list given in the 
     75    model. 
    7676    """ 
    7777    R_eff_type = int(pars.pop(RADIUS_MODE_ID, 1.0)) 
    7878    mesh = get_mesh(calculator.info, pars, dim=calculator.dim, mono=mono) 
    79     #print("pars", list(zip(*mesh))[0]) 
     79    #print("in call_Fq: pars", list(zip(*mesh))[0]) 
    8080    call_details, values, is_magnetic = make_kernel_args(calculator, mesh) 
    81     #print("values:", values) 
     81    #print("in call_Fq: values:", values) 
    8282    return calculator.Fq(call_details, values, cutoff, is_magnetic, R_eff_type) 
    8383 
     
    119119        active = lambda name: True 
    120120 
    121     #print("pars",[p.id for p in parameters.call_parameters]) 
     121    #print("in get_mesh: pars:",[p.id for p in parameters.call_parameters]) 
    122122    mesh = [_get_par_weights(p, values, active(p.name)) 
    123123            for p in parameters.call_parameters] 
  • sasmodels/generate.py

    ra8a1d48 ra34b811  
    2727    which are hollow. 
    2828 
    29     *effective_radius(mode, p1, p2, ...)* returns the effective radius of 
     29    *radius_effective(mode, p1, p2, ...)* returns the effective radius of 
    3030    the form with particular dimensions.  Mode determines the type of 
    3131    effective radius returned, with mode=1 for equivalent volume. 
     
    190190 
    191191def get_data_path(external_dir, target_file): 
     192    """ 
     193    Search for the target file relative in the installed application. 
     194 
     195    Search first in the location of the generate module in case we are 
     196    running directly from the distribution.  Search next to the python 
     197    executable for windows installs.  Search in the ../Resources directory 
     198    next to the executable for Mac OS/X installs. 
     199    """ 
    192200    path = abspath(dirname(__file__)) 
    193201    if exists(joinpath(path, target_file)): 
     
    225233# 
    226234# NOTE: there is an RST_PROLOG at the end of this file which is NOT 
    227 # used for the bundled documentation. Still as long as we are defining the macros 
    228 # in two places any new addition should define the macro in both places. 
     235# used for the bundled documentation. Still as long as we are defining the 
     236# macros in two places any new addition should define the macro in both places. 
    229237RST_UNITS = { 
    230238    "Ang": "|Ang|", 
     
    530538 
    531539def test_tag_float(): 
    532     """check that floating point constants are properly identified and tagged with 'f'""" 
     540    """Check that floating point constants are identified and tagged with 'f'""" 
    533541 
    534542    cases = """ 
     
    707715    return False 
    708716 
    709 _SHELL_VOLUME_PATTERN = re.compile(r"(^|\s)double\s+shell_volume[(]", flags=re.MULTILINE) 
     717_SHELL_VOLUME_PATTERN = re.compile(r"(^|\s)double\s+shell_volume[(]", 
     718                                   flags=re.MULTILINE) 
    710719def contains_shell_volume(source): 
    711720    # type: (List[str]) -> bool 
     
    801810            logger.warn("oriented shapes should define Iqac or Iqabc") 
    802811        else: 
    803             raise ValueError("Expected function Iqac or Iqabc for oriented shape") 
     812            raise ValueError("Expected Iqac or Iqabc for oriented shape") 
    804813 
    805814    # Define the parameter table 
     
    811820                              for p in partable.kernel_parameters)) 
    812821    # Define the function calls 
    813     call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) 0.0" 
     822    call_radius_effective = "#define CALL_RADIUS_EFFECTIVE(_mode, _v) 0.0" 
    814823    if partable.form_volume_parameters: 
    815824        refs = _call_pars("_v.", partable.form_volume_parameters) 
    816825        if is_hollow: 
    817             call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = form_volume(%s); _shell = shell_volume(%s); } while (0)"%((",".join(refs),)*2) 
     826            call_volume = ( 
     827                "#define CALL_VOLUME(_form, _shell, _v) " 
     828                "do { _form = form_volume(%s); _shell = shell_volume(%s); } " 
     829                "while (0)") % ((",".join(refs),)*2) 
    818830        else: 
    819             call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = _shell = form_volume(%s); } while (0)"%(",".join(refs)) 
    820         if model_info.effective_radius_type: 
    821             call_effective_radius = "#define CALL_EFFECTIVE_RADIUS(_mode, _v) effective_radius(_mode, %s)"%(",".join(refs)) 
     831            call_volume = ( 
     832                "#define CALL_VOLUME(_form, _shell, _v) " 
     833                "do { _form = _shell = form_volume(%s); } " 
     834                "while (0)") % (",".join(refs)) 
     835        if model_info.radius_effective_modes: 
     836            call_radius_effective = ( 
     837                "#define CALL_RADIUS_EFFECTIVE(_mode, _v) " 
     838                "radius_effective(_mode, %s)") % (",".join(refs)) 
    822839    else: 
    823840        # Model doesn't have volume.  We could make the kernel run a little 
    824841        # faster by not using/transferring the volume normalizations, but 
    825842        # the ifdef's reduce readability more than is worthwhile. 
    826         call_volume = "#define CALL_VOLUME(_form, _shell, _v) do { _form = _shell = 1.0; } while (0)" 
     843        call_volume = ( 
     844            "#define CALL_VOLUME(_form, _shell, _v) " 
     845            "do { _form = _shell = 1.0; } while (0)") 
    827846    source.append(call_volume) 
    828     source.append(call_effective_radius) 
     847    source.append(call_radius_effective) 
    829848    model_refs = _call_pars("_v.", partable.iq_parameters) 
    830849 
     
    879898    source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars)) 
    880899    source.append("#define PROJECTION %d"%PROJECTION) 
    881     # TODO: allow mixed python/opencl kernels? 
    882     ocl = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name) 
    883     dll = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name) 
    884  
    885     result = { 
    886         'dll': '\n'.join(source+dll[0]+dll[1]+dll[2]), 
    887         'opencl': '\n'.join(source+ocl[0]+ocl[1]+ocl[2]), 
    888     } 
     900    wrappers = _kernels(kernel_code, call_iq, clear_iq, 
     901                        call_iqxy, clear_iqxy, model_info.name) 
     902    code = '\n'.join(source + wrappers[0] + wrappers[1] + wrappers[2]) 
     903 
     904    # Note: Identical code for dll and opencl.  This may not always be the case 
     905    # so leave them as separate entries in the returned value. 
     906    result = {'dll': code, 'opencl': code} 
    889907    return result 
    890908 
  • sasmodels/gengauss.py

    rff31782 rb297ba9  
    11#!/usr/bin/env python 
     2""" 
     3Generate the Gauss-Legendre integration points and save them as a C file. 
     4""" 
    25from __future__ import division, print_function 
    3  
    4 import sys 
    56 
    67import numpy as np 
     
    89 
    910def gengauss(n, path): 
     11    """ 
     12    Save the Gauss-Legendre integration points for length *n* into file *path*. 
     13    """ 
    1014    z, w = leggauss(n) 
    1115 
  • sasmodels/guyou.py

    rb3703f5 rb297ba9  
    6868 
    6969def ellipticJi(u, v, m): 
     70    """Returns [sn, cn, dn](u + iv|m).""" 
    7071    scalar = np.isscalar(u) and np.isscalar(v) and np.isscalar(m) 
    7172    u, v, m = np.broadcast_arrays(u, v, m) 
     
    100101    ] 
    101102 
    102 # calculate F(phi+ipsi|m). 
    103 # see Abramowitz and Stegun, 17.4.11. 
    104103def ellipticFi(phi, psi, m): 
     104    """Returns F(phi+ipsi|m). See Abramowitz and Stegun, 17.4.11.""" 
    105105    if np.any(phi == 0): 
    106106        scalar = np.isscalar(phi) and np.isscalar(psi) and np.isscalar(m) 
     
    207207        plt.plot(x, y, 'g') 
    208208 
    209     for long in range(-limit, limit+1, step): 
    210         x, y = guyou_invert(scale*long, scale*long_line) 
     209    for longitude in range(-limit, limit+1, step): 
     210        x, y = guyou_invert(scale*longitude, scale*long_line) 
    211211        plt.plot(x, y, 'b') 
    212212    plt.xlabel('longitude') 
  • sasmodels/jitter.py

    r7d97437 r4e28511  
    11#!/usr/bin/env python 
     2# -*- coding: utf-8 -*- 
    23""" 
    34Jitter Explorer 
     
    56 
    67Application to explore orientation angle and angular dispersity. 
     8 
     9From the command line:: 
     10 
     11    # Show docs 
     12    python -m sasmodels.jitter --help 
     13 
     14    # Guyou projection jitter, uniform over 20 degree theta and 10 in phi 
     15    python -m sasmodels.jitter --projection=guyou --dist=uniform --jitter=20,10,0 
     16 
     17From a jupyter cell:: 
     18 
     19    import ipyvolume as ipv 
     20    from sasmodels import jitter 
     21    import importlib; importlib.reload(jitter) 
     22    jitter.set_plotter("ipv") 
     23 
     24    size = (10, 40, 100) 
     25    view = (20, 0, 0) 
     26 
     27    #size = (15, 15, 100) 
     28    #view = (60, 60, 0) 
     29 
     30    dview = (0, 0, 0) 
     31    #dview = (5, 5, 0) 
     32    #dview = (15, 180, 0) 
     33    #dview = (180, 15, 0) 
     34 
     35    projection = 'equirectangular' 
     36    #projection = 'azimuthal_equidistance' 
     37    #projection = 'guyou' 
     38    #projection = 'sinusoidal' 
     39    #projection = 'azimuthal_equal_area' 
     40 
     41    dist = 'uniform' 
     42    #dist = 'gaussian' 
     43 
     44    jitter.run(size=size, view=view, jitter=dview, dist=dist, projection=projection) 
     45    #filename = projection+('_theta' if dview[0] == 180 else '_phi' if dview[1] == 180 else '') 
     46    #ipv.savefig(filename+'.png') 
    747""" 
    848from __future__ import division, print_function 
     
    1050import argparse 
    1151 
    12 try: # CRUFT: travis-ci does not support mpl_toolkits.mplot3d 
    13     import mpl_toolkits.mplot3d  # Adds projection='3d' option to subplot 
    14 except ImportError: 
    15     pass 
    16  
    17 import matplotlib as mpl 
    18 import matplotlib.pyplot as plt 
    19 from matplotlib.widgets import Slider 
    2052import numpy as np 
    21 from numpy import pi, cos, sin, sqrt, exp, degrees, radians 
    22  
    23 def draw_beam(axes, view=(0, 0)): 
     53from numpy import pi, cos, sin, sqrt, exp, log, degrees, radians, arccos, arctan2 
     54 
     55# Too many complaints about variable names from pylint: 
     56#    a, b, c, u, v, x, y, z, dx, dy, dz, px, py, pz, R, Rx, Ry, Rz, ... 
     57# pylint: disable=invalid-name 
     58 
     59def draw_beam(axes, view=(0, 0), alpha=0.5, steps=25): 
    2460    """ 
    2561    Draw the beam going from source at (0, 0, 1) to detector at (0, 0, -1) 
    2662    """ 
    2763    #axes.plot([0,0],[0,0],[1,-1]) 
    28     #axes.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=0.8) 
    29  
    30     steps = 25 
    31     u = np.linspace(0, 2 * np.pi, steps) 
    32     v = np.linspace(-1, 1, steps) 
     64    #axes.scatter([0]*100,[0]*100,np.linspace(1, -1, 100), alpha=alpha) 
     65 
     66    u = np.linspace(0, 2 * pi, steps) 
     67    v = np.linspace(-1, 1, 2) 
    3368 
    3469    r = 0.02 
    35     x = r*np.outer(np.cos(u), np.ones_like(v)) 
    36     y = r*np.outer(np.sin(u), np.ones_like(v)) 
     70    x = r*np.outer(cos(u), np.ones_like(v)) 
     71    y = r*np.outer(sin(u), np.ones_like(v)) 
    3772    z = 1.3*np.outer(np.ones_like(u), v) 
    3873 
     
    4277    points = Rz(phi)*Ry(theta)*points 
    4378    x, y, z = [v.reshape(shape) for v in points] 
    44  
    45     axes.plot_surface(x, y, z, rstride=4, cstride=4, color='y', alpha=0.5) 
     79    axes.plot_surface(x, y, z, color='yellow', alpha=alpha) 
     80 
     81    # TODO: draw endcaps on beam 
     82    ## Drawing tiny balls on the end will work 
     83    #draw_sphere(axes, radius=0.02, center=(0, 0, 1.3), color='yellow', alpha=alpha) 
     84    #draw_sphere(axes, radius=0.02, center=(0, 0, -1.3), color='yellow', alpha=alpha) 
     85    ## The following does not work 
     86    #triangles = [(0, i+1, i+2) for i in range(steps-2)] 
     87    #x_cap, y_cap = x[:, 0], y[:, 0] 
     88    #for z_cap in z[:, 0], z[:, -1]: 
     89    #    axes.plot_trisurf(x_cap, y_cap, z_cap, triangles, 
     90    #                      color='yellow', alpha=alpha) 
     91 
    4692 
    4793def draw_ellipsoid(axes, size, view, jitter, steps=25, alpha=1): 
    4894    """Draw an ellipsoid.""" 
    4995    a, b, c = size 
    50     u = np.linspace(0, 2 * np.pi, steps) 
    51     v = np.linspace(0, np.pi, steps) 
    52     x = a*np.outer(np.cos(u), np.sin(v)) 
    53     y = b*np.outer(np.sin(u), np.sin(v)) 
    54     z = c*np.outer(np.ones_like(u), np.cos(v)) 
     96    u = np.linspace(0, 2 * pi, steps) 
     97    v = np.linspace(0, pi, steps) 
     98    x = a*np.outer(cos(u), sin(v)) 
     99    y = b*np.outer(sin(u), sin(v)) 
     100    z = c*np.outer(np.ones_like(u), cos(v)) 
    55101    x, y, z = transform_xyz(view, jitter, x, y, z) 
    56102 
    57     axes.plot_surface(x, y, z, rstride=4, cstride=4, color='w', alpha=alpha) 
     103    axes.plot_surface(x, y, z, color='w', alpha=alpha) 
    58104 
    59105    draw_labels(axes, view, jitter, [ 
     
    68114def draw_sc(axes, size, view, jitter, steps=None, alpha=1): 
    69115    """Draw points for simple cubic paracrystal""" 
     116    # pylint: disable=unused-argument 
    70117    atoms = _build_sc() 
    71118    _draw_crystal(axes, size, view, jitter, atoms=atoms) 
     
    73120def draw_fcc(axes, size, view, jitter, steps=None, alpha=1): 
    74121    """Draw points for face-centered cubic paracrystal""" 
     122    # pylint: disable=unused-argument 
    75123    # Build the simple cubic crystal 
    76124    atoms = _build_sc() 
     
    88136def draw_bcc(axes, size, view, jitter, steps=None, alpha=1): 
    89137    """Draw points for body-centered cubic paracrystal""" 
     138    # pylint: disable=unused-argument 
    90139    # Build the simple cubic crystal 
    91140    atoms = _build_sc() 
     
    124173    return atoms 
    125174 
    126 def draw_parallelepiped(axes, size, view, jitter, steps=None, alpha=1): 
    127     """Draw a parallelepiped.""" 
     175def draw_box(axes, size, view): 
     176    """Draw a wireframe box at a particular view.""" 
     177    a, b, c = size 
     178    x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1]) 
     179    y = b*np.array([+1, +1, -1, -1, +1, +1, -1, -1]) 
     180    z = c*np.array([+1, +1, +1, +1, -1, -1, -1, -1]) 
     181    x, y, z = transform_xyz(view, None, x, y, z) 
     182    def _draw(i, j): 
     183        axes.plot([x[i], x[j]], [y[i], y[j]], [z[i], z[j]], color='black') 
     184    _draw(0, 1) 
     185    _draw(0, 2) 
     186    _draw(0, 3) 
     187    _draw(7, 4) 
     188    _draw(7, 5) 
     189    _draw(7, 6) 
     190 
     191def draw_parallelepiped(axes, size, view, jitter, steps=None, 
     192                        color=(0.6, 1.0, 0.6), alpha=1): 
     193    """Draw a parallelepiped surface, with view and jitter.""" 
     194    # pylint: disable=unused-argument 
    128195    a, b, c = size 
    129196    x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1]) 
     
    142209 
    143210    x, y, z = transform_xyz(view, jitter, x, y, z) 
    144     axes.plot_trisurf(x, y, triangles=tri, Z=z, color='w', alpha=alpha) 
    145  
    146     # Draw pink face on box. 
     211    axes.plot_trisurf(x, y, triangles=tri, Z=z, color=color, alpha=alpha, 
     212                      linewidth=0) 
     213 
     214    # Colour the c+ face of the box. 
    147215    # Since I can't control face color, instead draw a thin box situated just 
    148216    # in front of the "c+" face.  Use the c face so that rotations about psi 
    149217    # rotate that face. 
    150     if 1: 
     218    if 0: # pylint: disable=using-constant-test 
     219        color = (1, 0.6, 0.6)  # pink 
    151220        x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1]) 
    152221        y = b*np.array([+1, +1, -1, -1, +1, +1, -1, -1]) 
    153222        z = c*np.array([+1, +1, +1, +1, -1, -1, -1, -1]) 
    154223        x, y, z = transform_xyz(view, jitter, x, y, abs(z)+0.001) 
    155         axes.plot_trisurf(x, y, triangles=tri, Z=z, color=[1, 0.6, 0.6], alpha=alpha) 
     224        axes.plot_trisurf(x, y, triangles=tri, Z=z, color=color, alpha=alpha) 
    156225 
    157226    draw_labels(axes, view, jitter, [ 
     
    164233    ]) 
    165234 
    166 def draw_sphere(axes, radius=10., steps=100): 
     235def draw_sphere(axes, radius=1.0, steps=25, 
     236                center=(0, 0, 0), color='w', alpha=1.): 
    167237    """Draw a sphere""" 
    168     u = np.linspace(0, 2 * np.pi, steps) 
    169     v = np.linspace(0, np.pi, steps) 
    170  
    171     x = radius * np.outer(np.cos(u), np.sin(v)) 
    172     y = radius * np.outer(np.sin(u), np.sin(v)) 
    173     z = radius * np.outer(np.ones(np.size(u)), np.cos(v)) 
    174     axes.plot_surface(x, y, z, rstride=4, cstride=4, color='w') 
    175  
    176 def draw_jitter(axes, view, jitter, dist='gaussian', size=(0.1, 0.4, 1.0), 
    177                 draw_shape=draw_parallelepiped): 
     238    u = np.linspace(0, 2 * pi, steps) 
     239    v = np.linspace(0, pi, steps) 
     240 
     241    x = radius * np.outer(cos(u), sin(v)) + center[0] 
     242    y = radius * np.outer(sin(u), sin(v)) + center[1] 
     243    z = radius * np.outer(np.ones(np.size(u)), cos(v)) + center[2] 
     244    axes.plot_surface(x, y, z, color=color, alpha=alpha) 
     245    #axes.plot_wireframe(x, y, z) 
     246 
     247def draw_axes(axes, origin=(-1, -1, -1), length=(2, 2, 2)): 
     248    """Draw wireframe axes lines, with given origin and length""" 
     249    x, y, z = origin 
     250    dx, dy, dz = length 
     251    axes.plot([x, x+dx], [y, y], [z, z], color='black') 
     252    axes.plot([x, x], [y, y+dy], [z, z], color='black') 
     253    axes.plot([x, x], [y, y], [z, z+dz], color='black') 
     254 
     255def draw_person_on_sphere(axes, view, height=0.5, radius=1.0): 
     256    """ 
     257    Draw a person on the surface of a sphere. 
     258 
     259    *view* indicates (latitude, longitude, orientation) 
     260    """ 
     261    limb_offset = height * 0.05 
     262    head_radius = height * 0.10 
     263    head_height = height - head_radius 
     264    neck_length = head_radius * 0.50 
     265    shoulder_height = height - 2*head_radius - neck_length 
     266    torso_length = shoulder_height * 0.55 
     267    torso_radius = torso_length * 0.30 
     268    leg_length = shoulder_height - torso_length 
     269    arm_length = torso_length * 0.90 
     270 
     271    def _draw_part(y, z): 
     272        x = np.zeros_like(y) 
     273        xp, yp, zp = transform_xyz(view, None, x, y, z + radius) 
     274        axes.plot(xp, yp, zp, color='k') 
     275 
     276    # circle for head 
     277    u = np.linspace(0, 2 * pi, 40) 
     278    y = head_radius * cos(u) 
     279    z = head_radius * sin(u) + head_height 
     280    _draw_part(y, z) 
     281 
     282    # rectangle for body 
     283    y = np.array([-torso_radius, torso_radius, torso_radius, -torso_radius, -torso_radius]) 
     284    z = np.array([0., 0, torso_length, torso_length, 0]) + leg_length 
     285    _draw_part(y, z) 
     286 
     287    # arms 
     288    y = np.array([-torso_radius - limb_offset, -torso_radius - limb_offset, -torso_radius]) 
     289    z = np.array([shoulder_height - arm_length, shoulder_height, shoulder_height]) 
     290    _draw_part(y, z) 
     291    _draw_part(-y, z)  # pylint: disable=invalid-unary-operand-type 
     292 
     293    # legs 
     294    y = np.array([-torso_radius + limb_offset, -torso_radius + limb_offset]) 
     295    z = np.array([0, leg_length]) 
     296    _draw_part(y, z) 
     297    _draw_part(-y, z)  # pylint: disable=invalid-unary-operand-type 
     298 
     299    limits = [-radius-height, radius+height] 
     300    axes.set_xlim(limits) 
     301    axes.set_ylim(limits) 
     302    axes.set_zlim(limits) 
     303    axes.set_axis_off() 
     304 
     305def draw_jitter(axes, view, jitter, dist='gaussian', 
     306                size=(0.1, 0.4, 1.0), 
     307                draw_shape=draw_parallelepiped, 
     308                projection='equirectangular', 
     309                alpha=0.8, 
     310                views=None): 
    178311    """ 
    179312    Represent jitter as a set of shapes at different orientations. 
    180313    """ 
     314    project, project_weight = get_projection(projection) 
     315 
    181316    # set max diagonal to 0.95 
    182317    scale = 0.95/sqrt(sum(v**2 for v in size)) 
    183318    size = tuple(scale*v for v in size) 
    184319 
    185     #np.random.seed(10) 
    186     #cloud = np.random.randn(10,3) 
    187     cloud = [ 
    188         [-1, -1, -1], 
    189         [-1, -1, +0], 
    190         [-1, -1, +1], 
    191         [-1, +0, -1], 
    192         [-1, +0, +0], 
    193         [-1, +0, +1], 
    194         [-1, +1, -1], 
    195         [-1, +1, +0], 
    196         [-1, +1, +1], 
    197         [+0, -1, -1], 
    198         [+0, -1, +0], 
    199         [+0, -1, +1], 
    200         [+0, +0, -1], 
    201         [+0, +0, +0], 
    202         [+0, +0, +1], 
    203         [+0, +1, -1], 
    204         [+0, +1, +0], 
    205         [+0, +1, +1], 
    206         [+1, -1, -1], 
    207         [+1, -1, +0], 
    208         [+1, -1, +1], 
    209         [+1, +0, -1], 
    210         [+1, +0, +0], 
    211         [+1, +0, +1], 
    212         [+1, +1, -1], 
    213         [+1, +1, +0], 
    214         [+1, +1, +1], 
    215     ] 
    216320    dtheta, dphi, dpsi = jitter 
    217     if dtheta == 0: 
    218         cloud = [v for v in cloud if v[0] == 0] 
    219     if dphi == 0: 
    220         cloud = [v for v in cloud if v[1] == 0] 
    221     if dpsi == 0: 
    222         cloud = [v for v in cloud if v[2] == 0] 
    223     draw_shape(axes, size, view, [0, 0, 0], steps=100, alpha=0.8) 
    224     scale = {'gaussian':1, 'rectangle':1/sqrt(3), 'uniform':1/3}[dist] 
    225     for point in cloud: 
    226         delta = [scale*dtheta*point[0], scale*dphi*point[1], scale*dpsi*point[2]] 
    227         draw_shape(axes, size, view, delta, alpha=0.8) 
     321    base = {'gaussian':3, 'rectangle':sqrt(3), 'uniform':1}[dist] 
     322    def _steps(delta): 
     323        if views is None: 
     324            n = max(3, min(25, 2*int(base*delta/5))) 
     325        else: 
     326            n = views 
     327        return base*delta*np.linspace(-1, 1, n) if delta > 0 else [0.] 
     328    for theta in _steps(dtheta): 
     329        for phi in _steps(dphi): 
     330            for psi in _steps(dpsi): 
     331                w = project_weight(theta, phi, 1.0, 1.0) 
     332                if w > 0: 
     333                    dview = project(theta, phi, psi) 
     334                    draw_shape(axes, size, view, dview, alpha=alpha) 
    228335    for v in 'xyz': 
    229336        a, b, c = size 
    230         lim = np.sqrt(a**2 + b**2 + c**2) 
     337        lim = sqrt(a**2 + b**2 + c**2) 
    231338        getattr(axes, 'set_'+v+'lim')([-lim, lim]) 
    232         getattr(axes, v+'axis').label.set_text(v) 
     339        #getattr(axes, v+'axis').label.set_text(v) 
    233340 
    234341PROJECTIONS = [ 
     
    238345    'azimuthal_equal_area', 
    239346] 
    240 def draw_mesh(axes, view, jitter, radius=1.2, n=11, dist='gaussian', 
    241               projection='equirectangular'): 
    242     """ 
    243     Draw the dispersion mesh showing the theta-phi orientations at which 
    244     the model will be evaluated. 
    245  
     347def get_projection(projection): 
     348 
     349    """ 
    246350    jitter projections 
    247351    <https://en.wikipedia.org/wiki/List_of_map_projections> 
     
    297401        Should allow free movement in theta, but phi is distorted. 
    298402    """ 
     403    # pylint: disable=unused-argument 
    299404    # TODO: try Kent distribution instead of a gaussian warped by projection 
    300405 
    301     dist_x = np.linspace(-1, 1, n) 
    302     weights = np.ones_like(dist_x) 
    303     if dist == 'gaussian': 
    304         dist_x *= 3 
    305         weights = exp(-0.5*dist_x**2) 
    306     elif dist == 'rectangle': 
    307         # Note: uses sasmodels ridiculous definition of rectangle width 
    308         dist_x *= sqrt(3) 
    309     elif dist == 'uniform': 
    310         pass 
    311     else: 
    312         raise ValueError("expected dist to be gaussian, rectangle or uniform") 
    313  
    314406    if projection == 'equirectangular':  #define PROJECTION 1 
    315         def _rotate(theta_i, phi_j): 
    316             return Rx(phi_j)*Ry(theta_i) 
     407        def _project(theta_i, phi_j, psi): 
     408            latitude, longitude = theta_i, phi_j 
     409            return latitude, longitude, psi, 'xyz' 
     410            #return Rx(phi_j)*Ry(theta_i) 
    317411        def _weight(theta_i, phi_j, w_i, w_j): 
    318412            return w_i*w_j*abs(cos(radians(theta_i))) 
    319413    elif projection == 'sinusoidal':  #define PROJECTION 2 
    320         def _rotate(theta_i, phi_j): 
     414        def _project(theta_i, phi_j, psi): 
    321415            latitude = theta_i 
    322416            scale = cos(radians(latitude)) 
    323417            longitude = phi_j/scale if abs(phi_j) < abs(scale)*180 else 0 
    324418            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
    325             return Rx(longitude)*Ry(latitude) 
    326         def _weight(theta_i, phi_j, w_i, w_j): 
     419            return latitude, longitude, psi, 'xyz' 
     420            #return Rx(longitude)*Ry(latitude) 
     421        def _project(theta_i, phi_j, w_i, w_j): 
    327422            latitude = theta_i 
    328423            scale = cos(radians(latitude)) 
     
    330425            return active*w_i*w_j 
    331426    elif projection == 'guyou':  #define PROJECTION 3  (eventually?) 
    332         def _rotate(theta_i, phi_j): 
     427        def _project(theta_i, phi_j, psi): 
    333428            from .guyou import guyou_invert 
    334429            #latitude, longitude = guyou_invert([theta_i], [phi_j]) 
    335430            longitude, latitude = guyou_invert([phi_j], [theta_i]) 
    336             return Rx(longitude[0])*Ry(latitude[0]) 
     431            return latitude, longitude, psi, 'xyz' 
     432            #return Rx(longitude[0])*Ry(latitude[0]) 
    337433        def _weight(theta_i, phi_j, w_i, w_j): 
    338434            return w_i*w_j 
    339     elif projection == 'azimuthal_equidistance':  # Note: Rz Ry, not Rx Ry 
    340         def _rotate(theta_i, phi_j): 
     435    elif projection == 'azimuthal_equidistance': 
     436        # Note that calculates angles for Rz Ry rather than Rx Ry 
     437        def _project(theta_i, phi_j, psi): 
    341438            latitude = sqrt(theta_i**2 + phi_j**2) 
    342             longitude = degrees(np.arctan2(phi_j, theta_i)) 
     439            longitude = degrees(arctan2(phi_j, theta_i)) 
    343440            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
    344             return Rz(longitude)*Ry(latitude) 
     441            return latitude, longitude, psi-longitude, 'zyz' 
     442            #R = Rz(longitude)*Ry(latitude)*Rz(psi) 
     443            #return R_to_xyz(R) 
     444            #return Rz(longitude)*Ry(latitude) 
    345445        def _weight(theta_i, phi_j, w_i, w_j): 
    346446            # Weighting for each point comes from the integral: 
     
    376476            return weight*w_i*w_j if latitude < 180 else 0 
    377477    elif projection == 'azimuthal_equal_area': 
    378         def _rotate(theta_i, phi_j): 
     478        # Note that calculates angles for Rz Ry rather than Rx Ry 
     479        def _project(theta_i, phi_j, psi): 
    379480            radius = min(1, sqrt(theta_i**2 + phi_j**2)/180) 
    380             latitude = 180-degrees(2*np.arccos(radius)) 
    381             longitude = degrees(np.arctan2(phi_j, theta_i)) 
     481            latitude = 180-degrees(2*arccos(radius)) 
     482            longitude = degrees(arctan2(phi_j, theta_i)) 
    382483            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
    383             return Rz(longitude)*Ry(latitude) 
     484            return latitude, longitude, psi, 'zyz' 
     485            #R = Rz(longitude)*Ry(latitude)*Rz(psi) 
     486            #return R_to_xyz(R) 
     487            #return Rz(longitude)*Ry(latitude) 
    384488        def _weight(theta_i, phi_j, w_i, w_j): 
    385489            latitude = sqrt(theta_i**2 + phi_j**2) 
     
    389493        raise ValueError("unknown projection %r"%projection) 
    390494 
     495    return _project, _weight 
     496 
     497def R_to_xyz(R): 
     498    """ 
     499    Return phi, theta, psi Tait-Bryan angles corresponding to the given rotation matrix. 
     500 
     501    Extracting Euler Angles from a Rotation Matrix 
     502    Mike Day, Insomniac Games 
     503    https://d3cw3dd2w32x2b.cloudfront.net/wp-content/uploads/2012/07/euler-angles1.pdf 
     504    Based on: Shoemake’s "Euler Angle Conversion", Graphics Gems IV, pp.  222-229 
     505    """ 
     506    phi = arctan2(R[1, 2], R[2, 2]) 
     507    theta = arctan2(-R[0, 2], sqrt(R[0, 0]**2 + R[0, 1]**2)) 
     508    psi = arctan2(R[0, 1], R[0, 0]) 
     509    return degrees(phi), degrees(theta), degrees(psi) 
     510 
     511def draw_mesh(axes, view, jitter, radius=1.2, n=11, dist='gaussian', 
     512              projection='equirectangular'): 
     513    """ 
     514    Draw the dispersion mesh showing the theta-phi orientations at which 
     515    the model will be evaluated. 
     516    """ 
     517 
     518    _project, _weight = get_projection(projection) 
     519    def _rotate(theta, phi, z): 
     520        dview = _project(theta, phi, 0.) 
     521        if dview[3] == 'zyz': 
     522            return Rz(dview[1])*Ry(dview[0])*z 
     523        else:  # dview[3] == 'xyz': 
     524            return Rx(dview[1])*Ry(dview[0])*z 
     525 
     526 
     527    dist_x = np.linspace(-1, 1, n) 
     528    weights = np.ones_like(dist_x) 
     529    if dist == 'gaussian': 
     530        dist_x *= 3 
     531        weights = exp(-0.5*dist_x**2) 
     532    elif dist == 'rectangle': 
     533        # Note: uses sasmodels ridiculous definition of rectangle width 
     534        dist_x *= sqrt(3) 
     535    elif dist == 'uniform': 
     536        pass 
     537    else: 
     538        raise ValueError("expected dist to be gaussian, rectangle or uniform") 
     539 
    391540    # mesh in theta, phi formed by rotating z 
    392     dtheta, dphi, dpsi = jitter 
     541    dtheta, dphi, dpsi = jitter  # pylint: disable=unused-variable 
    393542    z = np.matrix([[0], [0], [radius]]) 
    394     points = np.hstack([_rotate(theta_i, phi_j)*z 
     543    points = np.hstack([_rotate(theta_i, phi_j, z) 
    395544                        for theta_i in dtheta*dist_x 
    396545                        for phi_j in dphi*dist_x]) 
     
    470619    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple. 
    471620    """ 
    472     dtheta, dphi, dpsi = jitter 
    473     points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points 
     621    if jitter is None: 
     622        return points 
     623    # Hack to deal with the fact that azimuthal_equidistance uses euler angles 
     624    if len(jitter) == 4: 
     625        dtheta, dphi, dpsi, _ = jitter 
     626        points = Rz(dphi)*Ry(dtheta)*Rz(dpsi)*points 
     627    else: 
     628        dtheta, dphi, dpsi = jitter 
     629        points = Rx(dphi)*Ry(dtheta)*Rz(dpsi)*points 
    474630    return points 
    475631 
     
    481637    """ 
    482638    theta, phi, psi = view 
    483     points = Rz(phi)*Ry(theta)*Rz(psi)*points 
     639    points = Rz(phi)*Ry(theta)*Rz(psi)*points # viewing angle 
     640    #points = Rz(psi)*Ry(pi/2-theta)*Rz(phi)*points # 1-D integration angles 
     641    #points = Rx(phi)*Ry(theta)*Rz(psi)*points  # angular dispersion angle 
    484642    return points 
     643 
     644def orient_relative_to_beam_quaternion(view, points): 
     645    """ 
     646    Apply the view transform to a set of points. 
     647 
     648    Points are stored in a 3 x n numpy matrix, not a numpy array or tuple. 
     649 
     650    This variant uses quaternions rather than rotation matrices for the 
     651    computation.  It works but it is not used because it doesn't solve 
     652    any problems.  The challenge of mapping theta/phi/psi to SO(3) does 
     653    not disappear by calculating the transform differently. 
     654    """ 
     655    theta, phi, psi = view 
     656    x, y, z = [1, 0, 0], [0, 1, 0], [0, 0, 1] 
     657    q = Quaternion(1, [0, 0, 0]) 
     658    ## Compose a rotation about the three axes by rotating 
     659    ## the unit vectors before applying the rotation. 
     660    #q = Quaternion.from_angle_axis(theta, q.rot(x)) * q 
     661    #q = Quaternion.from_angle_axis(phi, q.rot(y)) * q 
     662    #q = Quaternion.from_angle_axis(psi, q.rot(z)) * q 
     663    ## The above turns out to be equivalent to reversing 
     664    ## the order of application, so ignore it and use below. 
     665    q = q * Quaternion.from_angle_axis(theta, x) 
     666    q = q * Quaternion.from_angle_axis(phi, y) 
     667    q = q * Quaternion.from_angle_axis(psi, z) 
     668    ## Reverse the order by post-multiply rather than pre-multiply 
     669    #q = Quaternion.from_angle_axis(theta, x) * q 
     670    #q = Quaternion.from_angle_axis(phi, y) * q 
     671    #q = Quaternion.from_angle_axis(psi, z) * q 
     672    #print("axes psi", q.rot(np.matrix([x, y, z]).T)) 
     673    return q.rot(points) 
     674#orient_relative_to_beam = orient_relative_to_beam_quaternion 
     675 
     676# === Quaterion class definition === BEGIN 
     677# Simple stand-alone quaternion class 
     678 
     679# Note: this code works but isn't unused since quaternions didn't solve the 
     680# representation problem.  Leave it here in case we want to revisit this later. 
     681 
     682#import numpy as np 
     683class Quaternion(object): 
     684    r""" 
     685    Quaternion(w, r) = w + ir[0] + jr[1] + kr[2] 
     686 
     687    Quaternion.from_angle_axis(theta, r) for a rotation of angle theta about 
     688    an axis oriented toward the direction r.  This defines a unit quaternion, 
     689    normalizing $r$ to the unit vector $\hat r$, and setting quaternion 
     690    $Q = \cos \theta + \sin \theta \hat r$ 
     691 
     692    Quaternion objects can be multiplied, which applies a rotation about the 
     693    given axis, allowing composition of rotations without risk of gimbal lock. 
     694    The resulting quaternion is applied to a set of points using *Q.rot(v)*. 
     695    """ 
     696    def __init__(self, w, r): 
     697        self.w = w 
     698        self.r = np.asarray(r, dtype='d') 
     699 
     700    @staticmethod 
     701    def from_angle_axis(theta, r): 
     702        """Build quaternion as rotation theta about axis r""" 
     703        theta = np.radians(theta)/2 
     704        r = np.asarray(r) 
     705        w = np.cos(theta) 
     706        r = np.sin(theta)*r/np.dot(r, r) 
     707        return Quaternion(w, r) 
     708 
     709    def __mul__(self, other): 
     710        """Multiply quaterions""" 
     711        if isinstance(other, Quaternion): 
     712            w = self.w*other.w - np.dot(self.r, other.r) 
     713            r = self.w*other.r + other.w*self.r + np.cross(self.r, other.r) 
     714            return Quaternion(w, r) 
     715        raise NotImplementedError("Quaternion * non-quaternion not implemented") 
     716 
     717    def rot(self, v): 
     718        """Transform point *v* by quaternion""" 
     719        v = np.asarray(v).T 
     720        use_transpose = (v.shape[-1] != 3) 
     721        if use_transpose: 
     722            v = v.T 
     723        v = v + np.cross(2*self.r, np.cross(self.r, v) + self.w*v) 
     724        #v = v + 2*self.w*np.cross(self.r, v) + np.cross(2*self.r, np.cross(self.r, v)) 
     725        if use_transpose: 
     726            v = v.T 
     727        return v.T 
     728 
     729    def conj(self): 
     730        """Conjugate quaternion""" 
     731        return Quaternion(self.w, -self.r) 
     732 
     733    def inv(self): 
     734        """Inverse quaternion""" 
     735        return self.conj()/self.norm()**2 
     736 
     737    def norm(self): 
     738        """Quaternion length""" 
     739        return np.sqrt(self.w**2 + np.sum(self.r**2)) 
     740 
     741    def __str__(self): 
     742        return "%g%+gi%+gj%+gk"%(self.w, self.r[0], self.r[1], self.r[2]) 
     743 
     744def test_qrot(): 
     745    """Quaternion checks""" 
     746    # Define rotation of 60 degrees around an axis in y-z that is 60 degrees 
     747    # from y.  The rotation axis is determined by rotating the point [0, 1, 0] 
     748    # about x. 
     749    ax = Quaternion.from_angle_axis(60, [1, 0, 0]).rot([0, 1, 0]) 
     750    q = Quaternion.from_angle_axis(60, ax) 
     751    # Set the point to be rotated, and its expected rotated position. 
     752    p = [1, -1, 2] 
     753    target = [(10+4*np.sqrt(3))/8, (1+2*np.sqrt(3))/8, (14-3*np.sqrt(3))/8] 
     754    #print(q, q.rot(p) - target) 
     755    assert max(abs(q.rot(p) - target)) < 1e-14 
     756#test_qrot() 
     757#import sys; sys.exit() 
     758# === Quaterion class definition === END 
    485759 
    486760# translate between number of dimension of dispersity and the number of 
     
    504778    or the top of the range, depending on whether *mode* is 'central' or 'top'. 
    505779    """ 
    506     if portion == 1.0: 
    507         return data.min(), data.max() 
    508     elif mode == 'central': 
    509         data = np.sort(data.flatten()) 
    510         offset = int(portion*len(data)/2 + 0.5) 
    511         return data[offset], data[-offset] 
    512     elif mode == 'top': 
    513         data = np.sort(data.flatten()) 
    514         offset = int(portion*len(data) + 0.5) 
    515         return data[offset], data[-1] 
     780    if portion < 1.0: 
     781        if mode == 'central': 
     782            data = np.sort(data.flatten()) 
     783            offset = int(portion*len(data)/2 + 0.5) 
     784            return data[offset], data[-offset] 
     785        if mode == 'top': 
     786            data = np.sort(data.flatten()) 
     787            offset = int(portion*len(data) + 0.5) 
     788            return data[offset], data[-1] 
     789    # Default: full range 
     790    return data.min(), data.max() 
    516791 
    517792def draw_scattering(calculator, axes, view, jitter, dist='gaussian'): 
     
    544819 
    545820    # compute the pattern 
    546     qx, qy = calculator._data.x_bins, calculator._data.y_bins 
     821    qx, qy = calculator.qxqy 
    547822    Iqxy = calculator(**pars).reshape(len(qx), len(qy)) 
    548823 
    549824    # scale it and draw it 
    550     Iqxy = np.log(Iqxy) 
     825    Iqxy = log(Iqxy) 
    551826    if calculator.limits: 
    552827        # use limits from orientation (0,0,0) 
     
    556831        vmin = vmax*10**-7 
    557832        #vmin, vmax = clipped_range(Iqxy, portion=portion, mode='top') 
     833    #vmin, vmax = Iqxy.min(), Iqxy.max() 
    558834    #print("range",(vmin,vmax)) 
    559835    #qx, qy = np.meshgrid(qx, qy) 
    560     if 0: 
     836    if 0:  # pylint: disable=using-constant-test 
    561837        level = np.asarray(255*(Iqxy - vmin)/(vmax - vmin), 'i') 
    562838        level[level < 0] = 0 
     839        from matplotlib import pylab as plt 
    563840        colors = plt.get_cmap()(level) 
    564         axes.plot_surface(qx, qy, -1.1, rstride=1, cstride=1, facecolors=colors) 
    565     elif 1: 
     841        #from matplotlib import cm 
     842        #colors = cm.coolwarm(level) 
     843        #colors = cm.gist_yarg(level) 
     844        #colors = cm.Wistia(level) 
     845        colors[level <= 0, 3] = 0.  # set floor to transparent 
     846        x, y = np.meshgrid(qx/qx.max(), qy/qy.max()) 
     847        axes.plot_surface(x, y, -1.1*np.ones_like(x), facecolors=colors) 
     848    elif 1:  # pylint: disable=using-constant-test 
    566849        axes.contourf(qx/qx.max(), qy/qy.max(), Iqxy, zdir='z', offset=-1.1, 
    567850                      levels=np.linspace(vmin, vmax, 24)) 
     
    594877    calculator = DirectModel(data, model) 
    595878 
     879    # Remember the data axes so we can plot the results 
     880    calculator.qxqy = (q, q) 
     881 
    596882    # stuff the values for non-orientation parameters into the calculator 
    597883    calculator.pars = pars.copy() 
     
    601887    # under rotation or angular dispersion 
    602888    Iqxy = calculator(theta=0, phi=0, psi=0, **calculator.pars) 
    603     Iqxy = np.log(Iqxy) 
     889    Iqxy = log(Iqxy) 
    604890    vmin, vmax = clipped_range(Iqxy, 0.95, mode='top') 
    605891    calculator.limits = vmin, vmax+1 
     
    692978} 
    693979 
     980 
    694981def run(model_name='parallelepiped', size=(10, 40, 100), 
     982        view=(0, 0, 0), jitter=(0, 0, 0), 
    695983        dist='gaussian', mesh=30, 
    696984        projection='equirectangular'): 
     
    702990 
    703991    *size* gives the dimensions (a, b, c) of the shape. 
     992 
     993    *view* gives the initial view (theta, phi, psi) of the shape. 
     994 
     995    *view* gives the initial jitter (dtheta, dphi, dpsi) of the shape. 
    704996 
    705997    *dist* is the type of dispersition: gaussian, rectangle, or uniform. 
     
    7211013    calculator, size = select_calculator(model_name, n=150, size=size) 
    7221014    draw_shape = DRAW_SHAPES.get(model_name, draw_parallelepiped) 
     1015    #draw_shape = draw_fcc 
    7231016 
    7241017    ## uncomment to set an independent the colour range for every view 
     
    7261019    calculator.limits = None 
    7271020 
    728     ## initial view 
    729     #theta, dtheta = 70., 10. 
    730     #phi, dphi = -45., 3. 
    731     #psi, dpsi = -45., 3. 
    732     theta, phi, psi = 0, 0, 0 
    733     dtheta, dphi, dpsi = 0, 0, 0 
     1021    PLOT_ENGINE(calculator, draw_shape, size, view, jitter, dist, mesh, projection) 
     1022 
     1023def _mpl_plot(calculator, draw_shape, size, view, jitter, dist, mesh, projection): 
     1024    # Note: travis-ci does not support mpl_toolkits.mplot3d, but this shouldn't be 
     1025    # an issue since we are lazy-loading the package on a path that isn't tested. 
     1026    # Importing mplot3d adds projection='3d' option to subplot 
     1027    import mpl_toolkits.mplot3d  # pylint: disable=unused-variable 
     1028    import matplotlib as mpl 
     1029    import matplotlib.pyplot as plt 
     1030    from matplotlib.widgets import Slider 
    7341031 
    7351032    ## create the plot window 
     
    7521049 
    7531050    ## add control widgets to plot 
    754     axes_theta = plt.axes([0.1, 0.15, 0.45, 0.04], **props) 
    755     axes_phi = plt.axes([0.1, 0.1, 0.45, 0.04], **props) 
    756     axes_psi = plt.axes([0.1, 0.05, 0.45, 0.04], **props) 
    757     stheta = Slider(axes_theta, 'Theta', -90, 90, valinit=theta) 
    758     sphi = Slider(axes_phi, 'Phi', -180, 180, valinit=phi) 
    759     spsi = Slider(axes_psi, 'Psi', -180, 180, valinit=psi) 
    760  
    761     axes_dtheta = plt.axes([0.75, 0.15, 0.15, 0.04], **props) 
    762     axes_dphi = plt.axes([0.75, 0.1, 0.15, 0.04], **props) 
    763     axes_dpsi = plt.axes([0.75, 0.05, 0.15, 0.04], **props) 
     1051    axes_theta = plt.axes([0.05, 0.15, 0.50, 0.04], **props) 
     1052    axes_phi = plt.axes([0.05, 0.10, 0.50, 0.04], **props) 
     1053    axes_psi = plt.axes([0.05, 0.05, 0.50, 0.04], **props) 
     1054    stheta = Slider(axes_theta, u'Ξ', -90, 90, valinit=0) 
     1055    sphi = Slider(axes_phi, u'φ', -180, 180, valinit=0) 
     1056    spsi = Slider(axes_psi, u'ψ', -180, 180, valinit=0) 
     1057 
     1058    axes_dtheta = plt.axes([0.70, 0.15, 0.20, 0.04], **props) 
     1059    axes_dphi = plt.axes([0.70, 0.1, 0.20, 0.04], **props) 
     1060    axes_dpsi = plt.axes([0.70, 0.05, 0.20, 0.04], **props) 
     1061 
    7641062    # Note: using ridiculous definition of rectangle distribution, whose width 
    7651063    # in sasmodels is sqrt(3) times the given width.  Divide by sqrt(3) to keep 
    7661064    # the maximum width to 90. 
    7671065    dlimit = DIST_LIMITS[dist] 
    768     sdtheta = Slider(axes_dtheta, 'dTheta', 0, 2*dlimit, valinit=dtheta) 
    769     sdphi = Slider(axes_dphi, 'dPhi', 0, 2*dlimit, valinit=dphi) 
    770     sdpsi = Slider(axes_dpsi, 'dPsi', 0, 2*dlimit, valinit=dpsi) 
    771  
     1066    sdtheta = Slider(axes_dtheta, u'Δξ', 0, 2*dlimit, valinit=0) 
     1067    sdphi = Slider(axes_dphi, u'Δφ', 0, 2*dlimit, valinit=0) 
     1068    sdpsi = Slider(axes_dpsi, u'Δψ', 0, 2*dlimit, valinit=0) 
     1069 
     1070    ## initial view and jitter 
     1071    theta, phi, psi = view 
     1072    stheta.set_val(theta) 
     1073    sphi.set_val(phi) 
     1074    spsi.set_val(psi) 
     1075    dtheta, dphi, dpsi = jitter 
     1076    sdtheta.set_val(dtheta) 
     1077    sdphi.set_val(dphi) 
     1078    sdpsi.set_val(dpsi) 
    7721079 
    7731080    ## callback to draw the new view 
    774     def update(val, axis=None): 
     1081    def _update(val, axis=None): 
     1082        # pylint: disable=unused-argument 
    7751083        view = stheta.val, sphi.val, spsi.val 
    7761084        jitter = sdtheta.val, sdphi.val, sdpsi.val 
    7771085        # set small jitter as 0 if multiple pd dims 
    7781086        dims = sum(v > 0 for v in jitter) 
    779         limit = [0, 0.5, 5][dims] 
     1087        limit = [0, 0.5, 5, 5][dims] 
    7801088        jitter = [0 if v < limit else v for v in jitter] 
    7811089        axes.cla() 
    782         draw_beam(axes, (0, 0)) 
    783         draw_jitter(axes, view, jitter, dist=dist, size=size, draw_shape=draw_shape) 
    784         #draw_jitter(axes, view, (0,0,0)) 
     1090 
     1091        ## Visualize as person on globe 
     1092        #draw_sphere(axes, radius=0.5) 
     1093        #draw_person_on_sphere(axes, view, radius=0.5) 
     1094 
     1095        ## Move beam instead of shape 
     1096        #draw_beam(axes, -view[:2]) 
     1097        #draw_jitter(axes, (0,0,0), (0,0,0), views=3) 
     1098 
     1099        ## Move shape and draw scattering 
     1100        draw_beam(axes, (0, 0), alpha=1.) 
     1101        #draw_person_on_sphere(axes, view, radius=1.2, height=0.5) 
     1102        draw_jitter(axes, view, jitter, dist=dist, size=size, alpha=1., 
     1103                    draw_shape=draw_shape, projection=projection, views=3) 
    7851104        draw_mesh(axes, view, jitter, dist=dist, n=mesh, projection=projection) 
    7861105        draw_scattering(calculator, axes, view, jitter, dist=dist) 
     1106 
    7871107        plt.gcf().canvas.draw() 
    7881108 
    7891109    ## bind control widgets to view updater 
    790     stheta.on_changed(lambda v: update(v, 'theta')) 
    791     sphi.on_changed(lambda v: update(v, 'phi')) 
    792     spsi.on_changed(lambda v: update(v, 'psi')) 
    793     sdtheta.on_changed(lambda v: update(v, 'dtheta')) 
    794     sdphi.on_changed(lambda v: update(v, 'dphi')) 
    795     sdpsi.on_changed(lambda v: update(v, 'dpsi')) 
     1110    stheta.on_changed(lambda v: _update(v, 'theta')) 
     1111    sphi.on_changed(lambda v: _update(v, 'phi')) 
     1112    spsi.on_changed(lambda v: _update(v, 'psi')) 
     1113    sdtheta.on_changed(lambda v: _update(v, 'dtheta')) 
     1114    sdphi.on_changed(lambda v: _update(v, 'dphi')) 
     1115    sdpsi.on_changed(lambda v: _update(v, 'dpsi')) 
    7961116 
    7971117    ## initialize view 
    798     update(None, 'phi') 
     1118    _update(None, 'phi') 
    7991119 
    8001120    ## go interactive 
    8011121    plt.show() 
    8021122 
     1123 
     1124def map_colors(z, kw): 
     1125    """ 
     1126    Process matplotlib-style colour arguments. 
     1127 
     1128    Pulls 'cmap', 'alpha', 'vmin', and 'vmax' from th *kw* dictionary, setting 
     1129    the *kw['color']* to an RGB array.  These are ignored if 'c' or 'color' are 
     1130    set inside *kw*. 
     1131    """ 
     1132    from matplotlib import cm 
     1133 
     1134    cmap = kw.pop('cmap', cm.coolwarm) 
     1135    alpha = kw.pop('alpha', None) 
     1136    vmin = kw.pop('vmin', z.min()) 
     1137    vmax = kw.pop('vmax', z.max()) 
     1138    c = kw.pop('c', None) 
     1139    color = kw.pop('color', c) 
     1140    if color is None: 
     1141        znorm = ((z - vmin) / (vmax - vmin)).clip(0, 1) 
     1142        color = cmap(znorm) 
     1143    elif isinstance(color, np.ndarray) and color.shape == z.shape: 
     1144        color = cmap(color) 
     1145    if alpha is None: 
     1146        if isinstance(color, np.ndarray): 
     1147            color = color[..., :3] 
     1148    else: 
     1149        color[..., 3] = alpha 
     1150    kw['color'] = color 
     1151 
     1152def make_vec(*args): 
     1153    """Turn all elements of *args* into numpy arrays""" 
     1154    #return [np.asarray(v, 'd').flatten() for v in args] 
     1155    return [np.asarray(v, 'd') for v in args] 
     1156 
     1157def make_image(z, kw): 
     1158    """Convert numpy array *z* into a *PIL* RGB image.""" 
     1159    import PIL.Image 
     1160    from matplotlib import cm 
     1161 
     1162    cmap = kw.pop('cmap', cm.coolwarm) 
     1163 
     1164    znorm = (z-z.min())/z.ptp() 
     1165    c = cmap(znorm) 
     1166    c = c[..., :3] 
     1167    rgb = np.asarray(c*255, 'u1') 
     1168    image = PIL.Image.fromarray(rgb, mode='RGB') 
     1169    return image 
     1170 
     1171 
     1172_IPV_MARKERS = { 
     1173    'o': 'sphere', 
     1174} 
     1175_IPV_COLORS = { 
     1176    'w': 'white', 
     1177    'k': 'black', 
     1178    'c': 'cyan', 
     1179    'm': 'magenta', 
     1180    'y': 'yellow', 
     1181    'r': 'red', 
     1182    'g': 'green', 
     1183    'b': 'blue', 
     1184} 
     1185def _ipv_fix_color(kw): 
     1186    alpha = kw.pop('alpha', None) 
     1187    color = kw.get('color', None) 
     1188    if isinstance(color, str): 
     1189        color = _IPV_COLORS.get(color, color) 
     1190        kw['color'] = color 
     1191    if alpha is not None: 
     1192        color = kw['color'] 
     1193        #TODO: convert color to [r, g, b, a] if not already 
     1194        if isinstance(color, (tuple, list)): 
     1195            if len(color) == 3: 
     1196                color = (color[0], color[1], color[2], alpha) 
     1197            else: 
     1198                color = (color[0], color[1], color[2], alpha*color[3]) 
     1199            color = np.array(color) 
     1200        if isinstance(color, np.ndarray) and color.shape[-1] == 4: 
     1201            color[..., 3] = alpha 
     1202            kw['color'] = color 
     1203 
     1204def _ipv_set_transparency(kw, obj): 
     1205    color = kw.get('color', None) 
     1206    if (isinstance(color, np.ndarray) 
     1207            and color.shape[-1] == 4 
     1208            and (color[..., 3] != 1.0).any()): 
     1209        obj.material.transparent = True 
     1210        obj.material.side = "FrontSide" 
     1211 
     1212def ipv_axes(): 
     1213    """ 
     1214    Build a matplotlib style Axes interface for ipyvolume 
     1215    """ 
     1216    import ipyvolume as ipv 
     1217 
     1218    class Axes(object): 
     1219        """ 
     1220        Matplotlib Axes3D style interface to ipyvolume renderer. 
     1221        """ 
     1222        # pylint: disable=no-self-use,no-init 
     1223        # transparency can be achieved by setting the following: 
     1224        #    mesh.color = [r, g, b, alpha] 
     1225        #    mesh.material.transparent = True 
     1226        #    mesh.material.side = "FrontSide" 
     1227        # smooth(ish) rotation can be achieved by setting: 
     1228        #    slide.continuous_update = True 
     1229        #    figure.animation = 0. 
     1230        #    mesh.material.x = x 
     1231        #    mesh.material.y = y 
     1232        #    mesh.material.z = z 
     1233        # maybe need to synchronize update of x/y/z to avoid shimmy when moving 
     1234        def plot(self, x, y, z, **kw): 
     1235            """mpl style plot interface for ipyvolume""" 
     1236            _ipv_fix_color(kw) 
     1237            x, y, z = make_vec(x, y, z) 
     1238            ipv.plot(x, y, z, **kw) 
     1239        def plot_surface(self, x, y, z, **kw): 
     1240            """mpl style plot_surface interface for ipyvolume""" 
     1241            facecolors = kw.pop('facecolors', None) 
     1242            if facecolors is not None: 
     1243                kw['color'] = facecolors 
     1244            _ipv_fix_color(kw) 
     1245            x, y, z = make_vec(x, y, z) 
     1246            h = ipv.plot_surface(x, y, z, **kw) 
     1247            _ipv_set_transparency(kw, h) 
     1248            #h.material.side = "DoubleSide" 
     1249            return h 
     1250        def plot_trisurf(self, x, y, triangles=None, Z=None, **kw): 
     1251            """mpl style plot_trisurf interface for ipyvolume""" 
     1252            kw.pop('linewidth', None) 
     1253            _ipv_fix_color(kw) 
     1254            x, y, z = make_vec(x, y, Z) 
     1255            if triangles is not None: 
     1256                triangles = np.asarray(triangles) 
     1257            h = ipv.plot_trisurf(x, y, z, triangles=triangles, **kw) 
     1258            _ipv_set_transparency(kw, h) 
     1259            return h 
     1260        def scatter(self, x, y, z, **kw): 
     1261            """mpl style scatter interface for ipyvolume""" 
     1262            x, y, z = make_vec(x, y, z) 
     1263            map_colors(z, kw) 
     1264            marker = kw.get('marker', None) 
     1265            kw['marker'] = _IPV_MARKERS.get(marker, marker) 
     1266            h = ipv.scatter(x, y, z, **kw) 
     1267            _ipv_set_transparency(kw, h) 
     1268            return h 
     1269        def contourf(self, x, y, v, zdir='z', offset=0, levels=None, **kw): 
     1270            """mpl style contourf interface for ipyvolume""" 
     1271            # pylint: disable=unused-argument 
     1272            # Don't use contour for now (although we might want to later) 
     1273            self.pcolor(x, y, v, zdir='z', offset=offset, **kw) 
     1274        def pcolor(self, x, y, v, zdir='z', offset=0, **kw): 
     1275            """mpl style pcolor interface for ipyvolume""" 
     1276            # pylint: disable=unused-argument 
     1277            x, y, v = make_vec(x, y, v) 
     1278            image = make_image(v, kw) 
     1279            xmin, xmax = x.min(), x.max() 
     1280            ymin, ymax = y.min(), y.max() 
     1281            x = np.array([[xmin, xmax], [xmin, xmax]]) 
     1282            y = np.array([[ymin, ymin], [ymax, ymax]]) 
     1283            z = x*0 + offset 
     1284            u = np.array([[0., 1], [0, 1]]) 
     1285            v = np.array([[0., 0], [1, 1]]) 
     1286            h = ipv.plot_mesh(x, y, z, u=u, v=v, texture=image, wireframe=False) 
     1287            _ipv_set_transparency(kw, h) 
     1288            h.material.side = "DoubleSide" 
     1289            return h 
     1290        def text(self, *args, **kw): 
     1291            """mpl style text interface for ipyvolume""" 
     1292            pass 
     1293        def set_xlim(self, limits): 
     1294            """mpl style set_xlim interface for ipyvolume""" 
     1295            ipv.xlim(*limits) 
     1296        def set_ylim(self, limits): 
     1297            """mpl style set_ylim interface for ipyvolume""" 
     1298            ipv.ylim(*limits) 
     1299        def set_zlim(self, limits): 
     1300            """mpl style set_zlim interface for ipyvolume""" 
     1301            ipv.zlim(*limits) 
     1302        def set_axes_on(self): 
     1303            """mpl style set_axes_on interface for ipyvolume""" 
     1304            ipv.style.axis_on() 
     1305        def set_axis_off(self): 
     1306            """mpl style set_axes_off interface for ipyvolume""" 
     1307            ipv.style.axes_off() 
     1308    return Axes() 
     1309 
     1310def _ipv_plot(calculator, draw_shape, size, view, jitter, dist, mesh, projection): 
     1311    from IPython.display import display 
     1312    import ipywidgets as widgets 
     1313    import ipyvolume as ipv 
     1314 
     1315    axes = ipv_axes() 
     1316 
     1317    def _draw(view, jitter): 
     1318        camera = ipv.gcf().camera 
     1319        #print(ipv.gcf().__dict__.keys()) 
     1320        #print(dir(ipv.gcf())) 
     1321        ipv.figure(animation=0.)  # no animation when updating object mesh 
     1322 
     1323        # set small jitter as 0 if multiple pd dims 
     1324        dims = sum(v > 0 for v in jitter) 
     1325        limit = [0, 0.5, 5, 5][dims] 
     1326        jitter = [0 if v < limit else v for v in jitter] 
     1327 
     1328        ## Visualize as person on globe 
     1329        #draw_beam(axes, (0, 0)) 
     1330        #draw_sphere(axes, radius=0.5) 
     1331        #draw_person_on_sphere(axes, view, radius=0.5) 
     1332 
     1333        ## Move beam instead of shape 
     1334        #draw_beam(axes, view=(-view[0], -view[1])) 
     1335        #draw_jitter(axes, view=(0,0,0), jitter=(0,0,0)) 
     1336 
     1337        ## Move shape and draw scattering 
     1338        draw_beam(axes, (0, 0), steps=25) 
     1339        draw_jitter(axes, view, jitter, dist=dist, size=size, alpha=1.0, 
     1340                    draw_shape=draw_shape, projection=projection) 
     1341        draw_mesh(axes, view, jitter, dist=dist, n=mesh, radius=0.95, 
     1342                  projection=projection) 
     1343        draw_scattering(calculator, axes, view, jitter, dist=dist) 
     1344 
     1345        draw_axes(axes, origin=(-1, -1, -1.1)) 
     1346        ipv.style.box_off() 
     1347        ipv.style.axes_off() 
     1348        ipv.xyzlabel(" ", " ", " ") 
     1349 
     1350        ipv.gcf().camera = camera 
     1351        ipv.show() 
     1352 
     1353 
     1354    trange, prange = (-180., 180., 1.), (-180., 180., 1.) 
     1355    dtrange, dprange = (0., 180., 1.), (0., 180., 1.) 
     1356 
     1357    ## Super simple interfaca, but uses non-ascii variable namese 
     1358    # Ξ φ ψ Δξ Δφ Δψ 
     1359    #def update(**kw): 
     1360    #    view = kw['Ξ'], kw['φ'], kw['ψ'] 
     1361    #    jitter = kw['Δξ'], kw['Δφ'], kw['Δψ'] 
     1362    #    draw(view, jitter) 
     1363    #widgets.interact(update, Ξ=trange, φ=prange, ψ=prange, Δξ=dtrange, Δφ=dprange, Δψ=dprange) 
     1364 
     1365    def _update(theta, phi, psi, dtheta, dphi, dpsi): 
     1366        _draw(view=(theta, phi, psi), jitter=(dtheta, dphi, dpsi)) 
     1367 
     1368    def _slider(name, slice, init=0.): 
     1369        return widgets.FloatSlider( 
     1370            value=init, 
     1371            min=slice[0], 
     1372            max=slice[1], 
     1373            step=slice[2], 
     1374            description=name, 
     1375            disabled=False, 
     1376            #continuous_update=True, 
     1377            continuous_update=False, 
     1378            orientation='horizontal', 
     1379            readout=True, 
     1380            readout_format='.1f', 
     1381            ) 
     1382    theta = _slider(u'Ξ', trange, view[0]) 
     1383    phi = _slider(u'φ', prange, view[1]) 
     1384    psi = _slider(u'ψ', prange, view[2]) 
     1385    dtheta = _slider(u'Δξ', dtrange, jitter[0]) 
     1386    dphi = _slider(u'Δφ', dprange, jitter[1]) 
     1387    dpsi = _slider(u'Δψ', dprange, jitter[2]) 
     1388    fields = { 
     1389        'theta': theta, 'phi': phi, 'psi': psi, 
     1390        'dtheta': dtheta, 'dphi': dphi, 'dpsi': dpsi, 
     1391    } 
     1392    ui = widgets.HBox([ 
     1393        widgets.VBox([theta, phi, psi]), 
     1394        widgets.VBox([dtheta, dphi, dpsi]) 
     1395    ]) 
     1396 
     1397    out = widgets.interactive_output(_update, fields) 
     1398    display(ui, out) 
     1399 
     1400 
     1401_ENGINES = { 
     1402    "matplotlib": _mpl_plot, 
     1403    "mpl": _mpl_plot, 
     1404    #"plotly": _plotly_plot, 
     1405    "ipvolume": _ipv_plot, 
     1406    "ipv": _ipv_plot, 
     1407} 
     1408PLOT_ENGINE = _ENGINES["matplotlib"] 
     1409def set_plotter(name): 
     1410    """ 
     1411    Setting the plotting engine to matplotlib/ipyvolume or equivalently mpl/ipv. 
     1412    """ 
     1413    global PLOT_ENGINE 
     1414    PLOT_ENGINE = _ENGINES[name] 
     1415 
    8031416def main(): 
     1417    """ 
     1418    Command line interface to the jitter viewer. 
     1419    """ 
    8041420    parser = argparse.ArgumentParser( 
    8051421        description="Display jitter", 
     
    8111427    parser.add_argument('-s', '--size', type=str, default='10,40,100', 
    8121428                        help='a,b,c lengths') 
     1429    parser.add_argument('-v', '--view', type=str, default='0,0,0', 
     1430                        help='initial view angles') 
     1431    parser.add_argument('-j', '--jitter', type=str, default='0,0,0', 
     1432                        help='initial angular dispersion') 
    8131433    parser.add_argument('-d', '--distribution', choices=DISTRIBUTIONS, 
    8141434                        default=DISTRIBUTIONS[0], 
     
    8191439                        help='oriented shape') 
    8201440    opts = parser.parse_args() 
    821     size = tuple(int(v) for v in opts.size.split(',')) 
    822     run(opts.shape, size=size, 
     1441    size = tuple(float(v) for v in opts.size.split(',')) 
     1442    view = tuple(float(v) for v in opts.view.split(',')) 
     1443    jitter = tuple(float(v) for v in opts.jitter.split(',')) 
     1444    run(opts.shape, size=size, view=view, jitter=jitter, 
    8231445        mesh=opts.mesh, dist=opts.distribution, 
    8241446        projection=opts.projection) 
  • sasmodels/kernel.py

    rcd28947 ra34b811  
    2323# pylint: enable=unused-import 
    2424 
     25 
    2526class KernelModel(object): 
     27    """ 
     28    Model definition for the compute engine. 
     29    """ 
    2630    info = None  # type: ModelInfo 
    2731    dtype = None # type: np.dtype 
    2832    def make_kernel(self, q_vectors): 
    2933        # type: (List[np.ndarray]) -> "Kernel" 
     34        """ 
     35        Instantiate a kernel for evaluating the model at *q_vectors*. 
     36        """ 
    3037        raise NotImplementedError("need to implement make_kernel") 
    3138 
    3239    def release(self): 
    3340        # type: () -> None 
     41        """ 
     42        Free resources associated with the kernel. 
     43        """ 
    3444        pass 
    3545 
     46 
    3647class Kernel(object): 
    37     #: kernel dimension, either "1d" or "2d" 
     48    """ 
     49    Instantiated model for the compute engine, applied to a particular *q*. 
     50 
     51    Subclasses should define :meth:`_call_kernel` to evaluate the kernel over 
     52    its inputs. 
     53    """ 
     54    #: Kernel dimension, either "1d" or "2d". 
    3855    dim = None  # type: str 
     56    #: Model info. 
    3957    info = None  # type: ModelInfo 
    40     results = None # type: List[np.ndarray] 
     58    #: Numerical precision for the computation. 
    4159    dtype = None  # type: np.dtype 
     60    #: q values at which the kernel is to be evaluated 
     61    q_input = None  # type: Any 
     62    #: Place to hold result of :meth:`_call_kernel` for subclass. 
     63    result = None # type: np.ndarray 
    4264 
    4365    def Iq(self, call_details, values, cutoff, magnetic): 
     
    5577        built into the model, and do not need an additional scale. 
    5678        """ 
    57         _, F2, _, shell_volume, _ = self.Fq(call_details, values, cutoff, magnetic, 
    58                               effective_radius_type=0) 
     79        _, F2, _, shell_volume, _ = self.Fq(call_details, values, cutoff, 
     80                                            magnetic, radius_effective_mode=0) 
    5981        combined_scale = values[0]/shell_volume 
    6082        background = values[1] 
     
    6284    __call__ = Iq 
    6385 
    64     def Fq(self, call_details, values, cutoff, magnetic, effective_radius_type=0): 
     86    def Fq(self, call_details, values, cutoff, magnetic, 
     87           radius_effective_mode=0): 
    6588        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
    6689        r""" 
     
    79102        .. math:: 
    80103 
    81             I(q) = \text{scale} * P (1 + <F>^2/<F^2> (S - 1)) + \text{background} 
     104            I(q) = \text{scale} P (1 + <F>^2/<F^2> (S - 1)) + \text{background} 
    82105                 = \text{scale}/<V> (<F^2> + <F>^2 (S - 1)) + \text{background} 
    83106 
     
    88111        .. math:: 
    89112 
    90             P(q) = \frac{\sum w_k F^2(q, x_k) / \sum w_k}{\sum w_k V_k / \sum w_k} 
     113            P(q)=\frac{\sum w_k F^2(q, x_k) / \sum w_k}{\sum w_k V_k / \sum w_k} 
    91114 
    92115        The form factor itself is scaled by volume and contrast to compute the 
     
    120143        volume fraction of the particles.  The model can have several 
    121144        different ways to compute effective radius, with the 
    122         *effective_radius_type* parameter used to select amongst them.  The 
     145        *radius_effective_mode* parameter used to select amongst them.  The 
    123146        volume fraction of particles should be determined from the total 
    124147        volume fraction of the form, not just the shell volume fraction. 
     
    129152        hollow and solid shapes. 
    130153        """ 
    131         self._call_kernel(call_details, values, cutoff, magnetic, effective_radius_type) 
     154        self._call_kernel(call_details, values, cutoff, magnetic, 
     155                          radius_effective_mode) 
    132156        #print("returned",self.q_input.q, self.result) 
    133157        nout = 2 if self.info.have_Fq and self.dim == '1d' else 1 
     
    141165        form_volume = self.result[nout*self.q_input.nq + 1]/total_weight 
    142166        shell_volume = self.result[nout*self.q_input.nq + 2]/total_weight 
    143         effective_radius = self.result[nout*self.q_input.nq + 3]/total_weight 
     167        radius_effective = self.result[nout*self.q_input.nq + 3]/total_weight 
    144168        if shell_volume == 0.: 
    145169            shell_volume = 1. 
    146         F1 = self.result[1:nout*self.q_input.nq:nout]/total_weight if nout == 2 else None 
     170        F1 = (self.result[1:nout*self.q_input.nq:nout]/total_weight 
     171              if nout == 2 else None) 
    147172        F2 = self.result[0:nout*self.q_input.nq:nout]/total_weight 
    148         return F1, F2, effective_radius, shell_volume, form_volume/shell_volume 
     173        return F1, F2, radius_effective, shell_volume, form_volume/shell_volume 
    149174 
    150175    def release(self): 
    151176        # type: () -> None 
     177        """ 
     178        Free resources associated with the kernel instance. 
     179        """ 
    152180        pass 
    153181 
    154     def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
     182    def _call_kernel(self, call_details, values, cutoff, magnetic, 
     183                     radius_effective_mode): 
    155184        # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
    156185        """ 
  • sasmodels/kernel_iq.c

    r12f4c19 rde032da  
    2727//      parameters in the parameter table. 
    2828//  CALL_VOLUME(form, shell, table) : assign form and shell values 
    29 //  CALL_EFFECTIVE_RADIUS(type, table) : call the R_eff function 
     29//  CALL_RADIUS_EFFECTIVE(mode, table) : call the R_eff function 
    3030//  CALL_IQ(q, table) : call the Iq function for 1D calcs. 
    3131//  CALL_IQ_A(q, table) : call the Iq function with |q| for 2D data. 
     
    8585static void set_spin_weights(double in_spin, double out_spin, double weight[6]) 
    8686{ 
     87 
     88  double norm; 
    8789  in_spin = clip(in_spin, 0.0, 1.0); 
    8890  out_spin = clip(out_spin, 0.0, 1.0); 
     
    9496  // However, since the weights are applied to the final intensity and 
    9597  // are not interned inside the I(q) function, we want the full 
    96   // weight and not the square root.  Any function using 
    97   // set_spin_weights as part of calculating an amplitude will need to 
    98   // manually take that square root, but there is currently no such 
    99   // function. 
    100   weight[0] = (1.0-in_spin) * (1.0-out_spin); // dd 
    101   weight[1] = (1.0-in_spin) * out_spin;       // du 
    102   weight[2] = in_spin * (1.0-out_spin);       // ud 
    103   weight[3] = in_spin * out_spin;             // uu 
     98  // weight and not the square root.  Anyway no function will ever use 
     99  // set_spin_weights as part of calculating an amplitude, as the weights are 
     100  // related to polarisation efficiency of the instrument. The weights serve to 
     101  // construct various magnet scattering cross sections, which are linear combinations 
     102  // of the spin-resolved cross sections. The polarisation efficiency e_in and e_out 
     103  // are parameters ranging from 0.5 (unpolarised) beam to 1 (perfect optics). 
     104  // For in_spin or out_spin <0.5 one assumes a CS, where the spin is reversed/flipped 
     105  // with respect to the initial supermirror polariser. The actual polarisation efficiency 
     106  // in this case is however e_in/out = 1-in/out_spin. 
     107 
     108  if (out_spin < 0.5){norm=1-out_spin;} 
     109  else{norm=out_spin;} 
     110 
     111 
     112// The norm is needed to make sure that the scattering cross sections are 
     113//correctly weighted, such that the sum of spin-resolved measurements adds up to 
     114// the unpolarised or half-polarised scattering cross section. No intensity weighting 
     115// needed on the incoming polariser side (assuming that a user), has normalised 
     116// to the incoming flux with polariser in for SANSPOl and unpolarised beam, respectively. 
     117 
     118 
     119  weight[0] = (1.0-in_spin) * (1.0-out_spin) / norm; // dd 
     120  weight[1] = (1.0-in_spin) * out_spin / norm;       // du 
     121  weight[2] = in_spin * (1.0-out_spin) / norm;       // ud 
     122  weight[3] = in_spin * out_spin / norm;             // uu 
    104123  weight[4] = weight[1]; // du.imag 
    105124  weight[5] = weight[2]; // ud.imag 
     
    119138    switch (xs) { 
    120139      default: // keep compiler happy; condition ensures xs in [0,1,2,3] 
    121       case 0: // uu => sld - D M_perpx 
     140      case 0: // dd => sld - D M_perpx 
    122141          return sld - px*perp; 
    123       case 1: // ud.real => -D M_perpy 
     142      case 1: // du.real => -D M_perpy 
    124143          return py*perp; 
    125       case 2: // du.real => -D M_perpy 
     144      case 2: // ud.real => -D M_perpy 
    126145          return py*perp; 
    127       case 3: // dd => sld + D M_perpx 
     146      case 3: // uu => sld + D M_perpx 
    128147          return sld + px*perp; 
    129148    } 
     
    285304    pglobal double *result,       // nq+1 return values, again with padding 
    286305    const double cutoff,          // cutoff in the dispersity weight product 
    287     int32_t effective_radius_type // which effective radius to compute 
     306    int32_t radius_effective_mode // which effective radius to compute 
    288307    ) 
    289308{ 
     
    703722      weighted_form += weight * form; 
    704723      weighted_shell += weight * shell; 
    705       if (effective_radius_type != 0) { 
    706         weighted_radius += weight * CALL_EFFECTIVE_RADIUS(effective_radius_type, local_values.table); 
     724      if (radius_effective_mode != 0) { 
     725        weighted_radius += weight * CALL_RADIUS_EFFECTIVE(radius_effective_mode, local_values.table); 
    707726      } 
    708727      BUILD_ROTATION(); 
  • sasmodels/kernelcl.py

    r8064d5e ra34b811  
    5353from __future__ import print_function 
    5454 
     55import sys 
    5556import os 
    5657import warnings 
     
    6162    from time import perf_counter as clock 
    6263except ImportError: # CRUFT: python < 3.3 
    63     import sys 
    6464    if sys.platform.count("darwin") > 0: 
    6565        from time import time as clock 
     
    6969import numpy as np  # type: ignore 
    7070 
    71 # Attempt to setup opencl. This may fail if the pyopencl package is not 
     71# Attempt to setup OpenCL. This may fail if the pyopencl package is not 
    7272# installed or if it is installed but there are no devices available. 
    7373try: 
     
    7575    from pyopencl import mem_flags as mf 
    7676    from pyopencl.characterize import get_fast_inaccurate_build_options 
    77     # Ask OpenCL for the default context so that we know that one exists 
     77    # Ask OpenCL for the default context so that we know that one exists. 
    7878    cl.create_some_context(interactive=False) 
    7979    HAVE_OPENCL = True 
     
    9696# pylint: enable=unused-import 
    9797 
    98 # CRUFT: pyopencl < 2017.1  (as of June 2016 needs quotes around include path) 
     98 
     99# CRUFT: pyopencl < 2017.1 (as of June 2016 needs quotes around include path). 
    99100def quote_path(v): 
     101    # type: (str) -> str 
    100102    """ 
    101103    Quote the path if it is not already quoted. 
     
    107109    return '"'+v+'"' if v and ' ' in v and not v[0] in "\"'-" else v 
    108110 
     111 
    109112def fix_pyopencl_include(): 
     113    # type: (None) -> None 
    110114    """ 
    111115    Monkey patch pyopencl to allow spaces in include file path. 
    112116    """ 
    113     import pyopencl as cl 
    114     if hasattr(cl, '_DEFAULT_INCLUDE_OPTIONS'): 
    115         cl._DEFAULT_INCLUDE_OPTIONS = [quote_path(v) for v in cl._DEFAULT_INCLUDE_OPTIONS] 
     117    # pylint: disable=protected-access 
     118    import pyopencl 
     119    if hasattr(pyopencl, '_DEFAULT_INCLUDE_OPTIONS'): 
     120        pyopencl._DEFAULT_INCLUDE_OPTIONS = [ 
     121            quote_path(v) for v in pyopencl._DEFAULT_INCLUDE_OPTIONS 
     122            ] 
     123 
    116124 
    117125if HAVE_OPENCL: 
     
    126134MAX_LOOPS = 2048 
    127135 
    128  
    129136# Pragmas for enable OpenCL features.  Be sure to protect them so that they 
    130137# still compile even if OpenCL is not present. 
     
    141148""" 
    142149 
     150 
    143151def use_opencl(): 
     152    # type: () -> bool 
     153    """Return True if OpenCL is the default computational engine""" 
    144154    sas_opencl = os.environ.get("SAS_OPENCL", "OpenCL").lower() 
    145155    return HAVE_OPENCL and sas_opencl != "none" and not sas_opencl.startswith("cuda") 
    146156 
     157 
    147158ENV = None 
    148159def reset_environment(): 
     160    # type: () -> None 
    149161    """ 
    150162    Call to create a new OpenCL context, such as after a change to SAS_OPENCL. 
     
    152164    global ENV 
    153165    ENV = GpuEnvironment() if use_opencl() else None 
     166 
    154167 
    155168def environment(): 
     
    169182    return ENV 
    170183 
     184 
    171185def has_type(device, dtype): 
    172186    # type: (cl.Device, np.dtype) -> bool 
     
    179193        return "cl_khr_fp64" in device.extensions 
    180194    else: 
    181         # Not supporting F16 type since it isn't accurate enough 
     195        # Not supporting F16 type since it isn't accurate enough. 
    182196        return False 
     197 
    183198 
    184199def get_warp(kernel, queue): 
     
    190205        cl.kernel_work_group_info.PREFERRED_WORK_GROUP_SIZE_MULTIPLE, 
    191206        queue.device) 
     207 
    192208 
    193209def compile_model(context, source, dtype, fast=False): 
     
    211227        source_list.insert(0, _F64_PRAGMA) 
    212228 
    213     # Note: USE_SINCOS makes the intel cpu slower under opencl 
     229    # Note: USE_SINCOS makes the Intel CPU slower under OpenCL. 
    214230    if context.devices[0].type == cl.device_type.GPU: 
    215231        source_list.insert(0, "#define USE_SINCOS\n") 
     
    218234    source = "\n".join(source_list) 
    219235    program = cl.Program(context, source).build(options=options) 
     236 
    220237    #print("done with "+program) 
    221238    return program 
    222239 
    223240 
    224 # for now, this returns one device in the context 
    225 # TODO: create a context that contains all devices on all platforms 
     241# For now, this returns one device in the context. 
     242# TODO: Create a context that contains all devices on all platforms. 
    226243class GpuEnvironment(object): 
    227244    """ 
    228     GPU context, with possibly many devices, and one queue per device. 
    229  
    230     Because the environment can be reset during a live program (e.g., if the 
    231     user changes the active GPU device in the GUI), everything associated 
    232     with the device context must be cached in the environment and recreated 
    233     if the environment changes.  The *cache* attribute is a simple dictionary 
    234     which holds keys and references to objects, such as compiled kernels and 
    235     allocated buffers.  The running program should check in the cache for 
    236     long lived objects and create them if they are not there.  The program 
    237     should not hold onto cached objects, but instead only keep them active 
    238     for the duration of a function call.  When the environment is destroyed 
    239     then the *release* method for each active cache item is called before 
    240     the environment is freed.  This means that each cl buffer should be 
    241     in its own cache entry. 
     245    GPU context for OpenCL, with possibly many devices and one queue per device. 
    242246    """ 
    243247    def __init__(self): 
    244248        # type: () -> None 
    245         # find gpu context 
     249        # Find gpu context. 
    246250        context_list = _create_some_context() 
    247251 
     
    257261                self.context[dtype] = None 
    258262 
    259         # Build a queue for each context 
     263        # Build a queue for each context. 
    260264        self.queue = {} 
    261265        context = self.context[F32] 
     
    267271            self.queue[F64] = cl.CommandQueue(context, context.devices[0]) 
    268272 
    269         # Byte boundary for data alignment 
     273        ## Byte boundary for data alignment. 
    270274        #self.data_boundary = max(context.devices[0].min_data_type_align_size 
    271275        #                         for context in self.context.values()) 
    272276 
    273         # Cache for compiled programs, and for items in context 
     277        # Cache for compiled programs, and for items in context. 
    274278        self.compiled = {} 
    275         self.cache = {} 
    276279 
    277280    def has_type(self, dtype): 
     
    288291        """ 
    289292        # Note: PyOpenCL caches based on md5 hash of source, options and device 
    290         # so we don't really need to cache things for ourselves.  I'll do so 
    291         # anyway just to save some data munging time. 
     293        # but I'll do so as well just to save some data munging time. 
    292294        tag = generate.tag_source(source) 
    293295        key = "%s-%s-%s%s"%(name, dtype, tag, ("-fast" if fast else "")) 
    294         # Check timestamp on program 
     296        # Check timestamp on program. 
    295297        program, program_timestamp = self.compiled.get(key, (None, np.inf)) 
    296298        if program_timestamp < timestamp: 
     
    305307        return program 
    306308 
    307     def free_buffer(self, key): 
    308         if key in self.cache: 
    309             self.cache[key].release() 
    310             del self.cache[key] 
    311  
    312     def __del__(self): 
    313         for v in self.cache.values(): 
    314             release = getattr(v, 'release', lambda: None) 
    315             release() 
    316         self.cache = {} 
    317  
    318 _CURRENT_ID = 0 
    319 def unique_id(): 
    320     global _CURRENT_ID 
    321     _CURRENT_ID += 1 
    322     return _CURRENT_ID 
    323309 
    324310def _create_some_context(): 
     
    333319    which one (and not a CUDA device, or no GPU). 
    334320    """ 
    335     # Assume we do not get here if SAS_OPENCL is None or CUDA 
     321    # Assume we do not get here if SAS_OPENCL is None or CUDA. 
    336322    sas_opencl = os.environ.get('SAS_OPENCL', 'opencl') 
    337323    if sas_opencl.lower() != 'opencl': 
    338         # Setting PYOPENCL_CTX as a SAS_OPENCL to create cl context 
     324        # Setting PYOPENCL_CTX as a SAS_OPENCL to create cl context. 
    339325        os.environ["PYOPENCL_CTX"] = sas_opencl 
    340326 
     
    343329            return [cl.create_some_context(interactive=False)] 
    344330        except Exception as exc: 
     331            # TODO: Should warnings instead be put into logging.warn? 
    345332            warnings.warn(str(exc)) 
    346             warnings.warn("pyopencl.create_some_context() failed") 
    347             warnings.warn("the environment variable 'SAS_OPENCL' or 'PYOPENCL_CTX' might not be set correctly") 
     333            warnings.warn( 
     334                "pyopencl.create_some_context() failed.  The environment " 
     335                "variable 'SAS_OPENCL' or 'PYOPENCL_CTX' might not be set " 
     336                "correctly") 
    348337 
    349338    return _get_default_context() 
     339 
    350340 
    351341def _get_default_context(): 
     
    360350    # is running may increase throughput. 
    361351    # 
    362     # Macbook pro, base install: 
     352    # MacBook Pro, base install: 
    363353    #     {'Apple': [Intel CPU, NVIDIA GPU]} 
    364     # Macbook pro, base install: 
     354    # MacBook Pro, base install: 
    365355    #     {'Apple': [Intel CPU, Intel GPU]} 
    366     # 2 x nvidia 295 with Intel and NVIDIA opencl drivers installed 
     356    # 2 x NVIDIA 295 with Intel and NVIDIA opencl drivers install: 
    367357    #     {'Intel': [CPU], 'NVIDIA': [GPU, GPU, GPU, GPU]} 
    368358    gpu, cpu = None, None 
     
    387377            else: 
    388378                # System has cl.device_type.ACCELERATOR or cl.device_type.CUSTOM 
    389                 # Intel Phi for example registers as an accelerator 
     379                # Intel Phi for example registers as an accelerator. 
    390380                # Since the user installed a custom device on their system 
    391381                # and went through the pain of sorting out OpenCL drivers for 
     
    394384                gpu = device 
    395385 
    396     # order the devices by gpu then by cpu; when searching for an available 
     386    # Order the devices by gpu then by cpu; when searching for an available 
    397387    # device by data type they will be checked in this order, which means 
    398388    # that if the gpu supports double then the cpu will never be used (though 
     
    421411    that the compiler is allowed to take shortcuts. 
    422412    """ 
     413    info = None  # type: ModelInfo 
     414    source = ""  # type: str 
     415    dtype = None  # type: np.dtype 
     416    fast = False  # type: bool 
     417    _program = None  # type: cl.Program 
     418    _kernels = None  # type: Dict[str, cl.Kernel] 
     419 
    423420    def __init__(self, source, model_info, dtype=generate.F32, fast=False): 
    424421        # type: (Dict[str,str], ModelInfo, np.dtype, bool) -> None 
     
    427424        self.dtype = dtype 
    428425        self.fast = fast 
    429         self.timestamp = generate.ocl_timestamp(self.info) 
    430         self._cache_key = unique_id() 
    431426 
    432427    def __getstate__(self): 
     
    437432        # type: (Tuple[ModelInfo, str, np.dtype, bool]) -> None 
    438433        self.info, self.source, self.dtype, self.fast = state 
     434        self._program = self._kernels = None 
    439435 
    440436    def make_kernel(self, q_vectors): 
     
    442438        return GpuKernel(self, q_vectors) 
    443439 
    444     @property 
    445     def Iq(self): 
    446         return self._fetch_kernel('Iq') 
    447  
    448     def fetch_kernel(self, name): 
     440    def get_function(self, name): 
    449441        # type: (str) -> cl.Kernel 
    450442        """ 
     
    452444        does not already exist. 
    453445        """ 
    454         gpu = environment() 
    455         key = self._cache_key 
    456         if key not in gpu.cache: 
    457             program = gpu.compile_program( 
    458                 self.info.name, 
    459                 self.source['opencl'], 
    460                 self.dtype, 
    461                 self.fast, 
    462                 self.timestamp) 
    463             variants = ['Iq', 'Iqxy', 'Imagnetic'] 
    464             names = [generate.kernel_name(self.info, k) for k in variants] 
    465             kernels = [getattr(program, k) for k in names] 
    466             data = dict((k, v) for k, v in zip(variants, kernels)) 
    467             # keep a handle to program so GC doesn't collect 
    468             data['program'] = program 
    469             gpu.cache[key] = data 
    470         else: 
    471             data = gpu.cache[key] 
    472         return data[name] 
    473  
    474 # TODO: check that we don't need a destructor for buffers which go out of scope 
     446        if self._program is None: 
     447            self._prepare_program() 
     448        return self._kernels[name] 
     449 
     450    def _prepare_program(self): 
     451        # type: (str) -> None 
     452        env = environment() 
     453        timestamp = generate.ocl_timestamp(self.info) 
     454        program = env.compile_program( 
     455            self.info.name, 
     456            self.source['opencl'], 
     457            self.dtype, 
     458            self.fast, 
     459            timestamp) 
     460        variants = ['Iq', 'Iqxy', 'Imagnetic'] 
     461        names = [generate.kernel_name(self.info, k) for k in variants] 
     462        functions = [getattr(program, k) for k in names] 
     463        self._kernels = {k: v for k, v in zip(variants, functions)} 
     464        # Keep a handle to program so GC doesn't collect. 
     465        self._program = program 
     466 
     467 
     468# TODO: Check that we don't need a destructor for buffers which go out of scope. 
    475469class GpuInput(object): 
    476470    """ 
     
    494488    def __init__(self, q_vectors, dtype=generate.F32): 
    495489        # type: (List[np.ndarray], np.dtype) -> None 
    496         # TODO: do we ever need double precision q? 
     490        # TODO: Do we ever need double precision q? 
    497491        self.nq = q_vectors[0].size 
    498492        self.dtype = np.dtype(dtype) 
    499493        self.is_2d = (len(q_vectors) == 2) 
    500         # TODO: stretch input based on get_warp() 
    501         # not doing it now since warp depends on kernel, which is not known 
     494        # TODO: Stretch input based on get_warp(). 
     495        # Not doing it now since warp depends on kernel, which is not known 
    502496        # at this point, so instead using 32, which is good on the set of 
    503497        # architectures tested so far. 
     
    512506            self.q[:self.nq] = q_vectors[0] 
    513507        self.global_size = [self.q.shape[0]] 
    514         self._cache_key = unique_id() 
    515  
    516     @property 
    517     def q_b(self): 
    518         """Lazy creation of q buffer so it can survive context reset""" 
     508        #print("creating inputs of size", self.global_size) 
     509 
     510        # Transfer input value to GPU. 
    519511        env = environment() 
    520         key = self._cache_key 
    521         if key not in env.cache: 
    522             context = env.context[self.dtype] 
    523             #print("creating inputs of size", self.global_size) 
    524             buffer = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
    525                                hostbuf=self.q) 
    526             env.cache[key] = buffer 
    527         return env.cache[key] 
     512        context = env.context[self.dtype] 
     513        self.q_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
     514                             hostbuf=self.q) 
    528515 
    529516    def release(self): 
    530517        # type: () -> None 
    531518        """ 
    532         Free the buffer associated with the q value 
    533         """ 
    534         environment().free_buffer(id(self)) 
     519        Free the buffer associated with the q value. 
     520        """ 
     521        if self.q_b is not None: 
     522            self.q_b.release() 
     523            self.q_b = None 
    535524 
    536525    def __del__(self): 
     
    538527        self.release() 
    539528 
     529 
    540530class GpuKernel(Kernel): 
    541531    """ 
     
    544534    *model* is the GpuModel object to call 
    545535 
    546     The following attributes are defined: 
    547  
    548     *info* is the module information 
    549  
    550     *dtype* is the kernel precision 
    551  
    552     *dim* is '1d' or '2d' 
    553  
    554     *result* is a vector to contain the results of the call 
    555  
    556     The resulting call method takes the *pars*, a list of values for 
    557     the fixed parameters to the kernel, and *pd_pars*, a list of (value,weight) 
    558     vectors for the polydisperse parameters.  *cutoff* determines the 
    559     integration limits: any points with combined weight less than *cutoff* 
    560     will not be calculated. 
     536    The kernel is derived from :class:`Kernel`, providing the 
     537    :meth:`call_kernel` method to evaluate the kernel for a given set of 
     538    parameters.  Because of the need to move the q values to the GPU before 
     539    evaluation, the kernel is instantiated for a particular set of q vectors, 
     540    and can be called many times without transfering q each time. 
    561541 
    562542    Call :meth:`release` when done with the kernel instance. 
    563543    """ 
     544    #: SAS model information structure. 
     545    info = None  # type: ModelInfo 
     546    #: Kernel precision. 
     547    dtype = None  # type: np.dtype 
     548    #: Kernel dimensions (1d or 2d). 
     549    dim = ""  # type: str 
     550    #: Calculation results, updated after each call to :meth:`_call_kernel`. 
     551    result = None  # type: np.ndarray 
     552 
    564553    def __init__(self, model, q_vectors): 
    565         # type: (cl.Kernel, np.dtype, ModelInfo, List[np.ndarray]) -> None 
     554        # type: (GpuModel, List[np.ndarray]) -> None 
    566555        dtype = model.dtype 
    567556        self.q_input = GpuInput(q_vectors, dtype) 
    568557        self._model = model 
    569         # F16 isn't sufficient, so don't support it 
    570         self._as_dtype = np.float64 if dtype == generate.F64 else np.float32 
    571         self._cache_key = unique_id() 
    572  
    573         # attributes accessed from the outside 
     558 
     559        # Attributes accessed from the outside. 
    574560        self.dim = '2d' if self.q_input.is_2d else '1d' 
    575561        self.info = model.info 
    576         self.dtype = model.dtype 
    577  
    578         # holding place for the returned value 
     562        self.dtype = dtype 
     563 
     564        # Converter to translate input to target type. 
     565        self._as_dtype = np.float64 if dtype == generate.F64 else np.float32 
     566 
     567        # Holding place for the returned value. 
    579568        nout = 2 if self.info.have_Fq and self.dim == '1d' else 1 
    580         extra_q = 4  # total weight, form volume, shell volume and R_eff 
    581         self.result = np.empty(self.q_input.nq*nout+extra_q, dtype) 
    582  
    583     @property 
    584     def _result_b(self): 
    585         """Lazy creation of result buffer so it can survive context reset""" 
     569        extra_q = 4  # Total weight, form volume, shell volume and R_eff. 
     570        self.result = np.empty(self.q_input.nq*nout + extra_q, dtype) 
     571 
     572        # Allocate result value on GPU. 
    586573        env = environment() 
    587         key = self._cache_key 
    588         if key not in env.cache: 
    589             context = env.context[self.dtype] 
    590             width = ((self.result.size+31)//32)*32 * self.dtype.itemsize 
    591             buffer = cl.Buffer(context, mf.READ_WRITE, width) 
    592             env.cache[key] = buffer 
    593         return env.cache[key] 
    594  
    595     def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
    596         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
     574        context = env.context[self.dtype] 
     575        width = ((self.result.size+31)//32)*32 * self.dtype.itemsize 
     576        self._result_b = cl.Buffer(context, mf.READ_WRITE, width) 
     577 
     578    def _call_kernel(self, call_details, values, cutoff, magnetic, 
     579                     radius_effective_mode): 
     580        # type: (CallDetails, np.ndarray, float, bool, int) -> np.ndarray 
    597581        env = environment() 
    598582        queue = env.queue[self._model.dtype] 
    599583        context = queue.context 
    600584 
    601         # Arrange data transfer to/from card 
    602         q_b = self.q_input.q_b 
    603         result_b = self._result_b 
     585        # Arrange data transfer to card. 
    604586        details_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
    605587                              hostbuf=call_details.buffer) 
     
    607589                             hostbuf=values) 
    608590 
     591        # Setup kernel function and arguments. 
    609592        name = 'Iq' if self.dim == '1d' else 'Imagnetic' if magnetic else 'Iqxy' 
    610         kernel = self._model.fetch_kernel(name) 
     593        kernel = self._model.get_function(name) 
    611594        kernel_args = [ 
    612             np.uint32(self.q_input.nq), None, None, 
    613             details_b, values_b, q_b, result_b, 
    614             self._as_dtype(cutoff), 
    615             np.uint32(effective_radius_type), 
     595            np.uint32(self.q_input.nq),  # Number of inputs. 
     596            None,  # Placeholder for pd_start. 
     597            None,  # Placeholder for pd_stop. 
     598            details_b,  # Problem definition. 
     599            values_b,  # Parameter values. 
     600            self.q_input.q_b,  # Q values. 
     601            self._result_b,   # Result storage. 
     602            self._as_dtype(cutoff),  # Probability cutoff. 
     603            np.uint32(radius_effective_mode),  # R_eff mode. 
    616604        ] 
     605 
     606        # Call kernel and retrieve results. 
    617607        #print("Calling OpenCL") 
    618608        #call_details.show(values) 
    619         #Call kernel and retrieve results 
    620609        wait_for = None 
    621610        last_nap = clock() 
     
    628617                               *kernel_args, wait_for=wait_for)] 
    629618            if stop < call_details.num_eval: 
    630                 # Allow other processes to run 
     619                # Allow other processes to run. 
    631620                wait_for[0].wait() 
    632621                current_time = clock() 
     
    634623                    time.sleep(0.001) 
    635624                    last_nap = current_time 
    636         cl.enqueue_copy(queue, self.result, result_b, wait_for=wait_for) 
     625        cl.enqueue_copy(queue, self.result, self._result_b, wait_for=wait_for) 
    637626        #print("result", self.result) 
    638627 
    639         # Free buffers 
    640         for v in (details_b, values_b): 
    641             if v is not None: 
    642                 v.release() 
     628        # Free buffers. 
     629        details_b.release() 
     630        values_b.release() 
    643631 
    644632    def release(self): 
     
    647635        Release resources associated with the kernel. 
    648636        """ 
    649         environment().free_buffer(id(self)) 
    650637        self.q_input.release() 
     638        if self._result_b is not None: 
     639            self._result_b.release() 
     640            self._result_b = None 
    651641 
    652642    def __del__(self): 
  • sasmodels/kernelcuda.py

    rf872fd1 ra34b811  
    5959 
    6060import os 
    61 import warnings 
    6261import logging 
    6362import time 
    6463import re 
     64import atexit 
    6565 
    6666import numpy as np  # type: ignore 
    6767 
    6868 
    69 # Attempt to setup cuda. This may fail if the pycuda package is not 
     69# Attempt to setup CUDA. This may fail if the pycuda package is not 
    7070# installed or if it is installed but there are no devices available. 
    7171try: 
     
    107107MAX_LOOPS = 2048 
    108108 
     109 
    109110def use_cuda(): 
    110     env = os.environ.get("SAS_OPENCL", "").lower() 
    111     return HAVE_CUDA and (env == "" or env.startswith("cuda")) 
     111    # type: None -> bool 
     112    """Returns True if CUDA is the default compute engine.""" 
     113    sas_opencl = os.environ.get("SAS_OPENCL", "CUDA").lower() 
     114    return HAVE_CUDA and sas_opencl.startswith("cuda") 
     115 
    112116 
    113117ENV = None 
     
    122126    ENV = GpuEnvironment() if use_cuda() else None 
    123127 
     128 
    124129def environment(): 
    125130    # type: () -> "GpuEnvironment" 
     
    132137        if not HAVE_CUDA: 
    133138            raise RuntimeError("CUDA startup failed with ***" 
    134                             + CUDA_ERROR + "***; using C compiler instead") 
     139                               + CUDA_ERROR + "***; using C compiler instead") 
    135140        reset_environment() 
    136141        if ENV is None: 
     
    138143    return ENV 
    139144 
     145 
     146# PyTest is not freeing ENV, so make sure it gets freed. 
     147atexit.register(lambda: ENV.release() if ENV is not None else None) 
     148 
     149 
    140150def has_type(dtype): 
    141151    # type: (np.dtype) -> bool 
     
    143153    Return true if device supports the requested precision. 
    144154    """ 
    145     # Assume the nvidia card supports 32-bit and 64-bit floats. 
    146     # TODO: check if pycuda support F16 
     155    # Assume the NVIDIA card supports 32-bit and 64-bit floats. 
     156    # TODO: Check if pycuda support F16. 
    147157    return dtype in (generate.F32, generate.F64) 
    148158 
    149159 
    150160FUNCTION_PATTERN = re.compile(r"""^ 
    151   (?P<space>\s*)                   # initial space 
    152   (?P<qualifiers>^(?:\s*\b\w+\b\s*)+) # one or more qualifiers before function 
    153   (?P<function>\s*\b\w+\b\s*[(])      # function name plus open parens 
     161  (?P<space>\s*)                       # Initial space. 
     162  (?P<qualifiers>^(?:\s*\b\w+\b\s*)+)  # One or more qualifiers before function. 
     163  (?P<function>\s*\b\w+\b\s*[(])       # Function name plus open parens. 
    154164  """, re.VERBOSE|re.MULTILINE) 
    155165 
     
    158168  """, re.VERBOSE|re.MULTILINE) 
    159169 
     170 
    160171def _add_device_tag(match): 
    161172    # type: (None) -> str 
    162     # Note: should be re.Match, but that isn't a simple type 
     173    # Note: Should be re.Match, but that isn't a simple type. 
    163174    """ 
    164175    replace qualifiers with __device__ qualifiers if needed 
     
    173184        return "".join((space, "__device__ ", qualifiers, function)) 
    174185 
     186 
    175187def mark_device_functions(source): 
    176188    # type: (str) -> str 
     
    180192    return FUNCTION_PATTERN.sub(_add_device_tag, source) 
    181193 
     194 
    182195def show_device_functions(source): 
    183196    # type: (str) -> str 
     
    186199    """ 
    187200    for match in FUNCTION_PATTERN.finditer(source): 
    188         print(match.group('qualifiers').replace('\n',r'\n'), match.group('function'), '(') 
     201        print(match.group('qualifiers').replace('\n', r'\n'), 
     202              match.group('function'), '(') 
    189203    return source 
     204 
    190205 
    191206def compile_model(source, dtype, fast=False): 
     
    212227    #options = ['--verbose', '-E'] 
    213228    options = ['--use_fast_math'] if fast else None 
    214     program = SourceModule(source, no_extern_c=True, options=options) # include_dirs=[...] 
     229    program = SourceModule(source, no_extern_c=True, options=options) #, include_dirs=[...]) 
    215230 
    216231    #print("done with "+program) 
     
    218233 
    219234 
    220 # for now, this returns one device in the context 
    221 # TODO: create a context that contains all devices on all platforms 
     235# For now, this returns one device in the context. 
     236# TODO: Create a context that contains all devices on all platforms. 
    222237class GpuEnvironment(object): 
    223238    """ 
    224     GPU context, with possibly many devices, and one queue per device. 
     239    GPU context for CUDA. 
    225240    """ 
    226241    context = None # type: cuda.Context 
    227242    def __init__(self, devnum=None): 
    228243        # type: (int) -> None 
    229         # Byte boundary for data alignment 
    230         #self.data_boundary = max(d.min_data_type_align_size 
    231         #                         for d in self.context.devices) 
    232         self.compiled = {} 
    233244        env = os.environ.get("SAS_OPENCL", "").lower() 
    234245        if devnum is None and env.startswith("cuda:"): 
    235246            devnum = int(env[5:]) 
     247 
    236248        # Set the global context to the particular device number if one is 
    237249        # given, otherwise use the default context.  Perhaps this will be set 
     
    242254            self.context = make_default_context() 
    243255 
     256        ## Byte boundary for data alignment. 
     257        #self.data_boundary = max(d.min_data_type_align_size 
     258        #                         for d in self.context.devices) 
     259 
     260        # Cache for compiled programs, and for items in context. 
     261        self.compiled = {} 
     262 
    244263    def release(self): 
     264        """Free the CUDA device associated with this context.""" 
    245265        if self.context is not None: 
    246266            self.context.pop() 
     
    262282        Compile the program for the device in the given context. 
    263283        """ 
    264         # Note: PyOpenCL caches based on md5 hash of source, options and device 
    265         # so we don't really need to cache things for ourselves.  I'll do so 
    266         # anyway just to save some data munging time. 
     284        # Note: PyCuda (probably) caches but I'll do so as well just to 
     285        # save some data munging time. 
    267286        tag = generate.tag_source(source) 
    268287        key = "%s-%s-%s%s"%(name, dtype, tag, ("-fast" if fast else "")) 
    269         # Check timestamp on program 
     288        # Check timestamp on program. 
    270289        program, program_timestamp = self.compiled.get(key, (None, np.inf)) 
    271290        if program_timestamp < timestamp: 
     
    277296        return program 
    278297 
     298 
    279299class GpuModel(KernelModel): 
    280300    """ 
     
    292312    that the compiler is allowed to take shortcuts. 
    293313    """ 
    294     info = None # type: ModelInfo 
    295     source = "" # type: str 
    296     dtype = None # type: np.dtype 
    297     fast = False # type: bool 
    298     program = None # type: SourceModule 
    299     _kernels = None # type: List[cuda.Function] 
     314    info = None  # type: ModelInfo 
     315    source = ""  # type: str 
     316    dtype = None  # type: np.dtype 
     317    fast = False  # type: bool 
     318    _program = None # type: SourceModule 
     319    _kernels = None  # type: Dict[str, cuda.Function] 
    300320 
    301321    def __init__(self, source, model_info, dtype=generate.F32, fast=False): 
     
    305325        self.dtype = dtype 
    306326        self.fast = fast 
    307         self.program = None # delay program creation 
    308         self._kernels = None 
    309327 
    310328    def __getstate__(self): 
     
    315333        # type: (Tuple[ModelInfo, str, np.dtype, bool]) -> None 
    316334        self.info, self.source, self.dtype, self.fast = state 
    317         self.program = None 
     335        self._program = self._kernels = None 
    318336 
    319337    def make_kernel(self, q_vectors): 
    320338        # type: (List[np.ndarray]) -> "GpuKernel" 
    321         if self.program is None: 
    322             compile_program = environment().compile_program 
    323             timestamp = generate.ocl_timestamp(self.info) 
    324             self.program = compile_program( 
    325                 self.info.name, 
    326                 self.source['opencl'], 
    327                 self.dtype, 
    328                 self.fast, 
    329                 timestamp) 
    330             variants = ['Iq', 'Iqxy', 'Imagnetic'] 
    331             names = [generate.kernel_name(self.info, k) for k in variants] 
    332             kernels = [self.program.get_function(k) for k in names] 
    333             self._kernels = dict((k, v) for k, v in zip(variants, kernels)) 
    334         is_2d = len(q_vectors) == 2 
    335         if is_2d: 
    336             kernel = [self._kernels['Iqxy'], self._kernels['Imagnetic']] 
    337         else: 
    338             kernel = [self._kernels['Iq']]*2 
    339         return GpuKernel(kernel, self.dtype, self.info, q_vectors) 
    340  
    341     def release(self): 
    342         # type: () -> None 
    343         """ 
    344         Free the resources associated with the model. 
    345         """ 
    346         if self.program is not None: 
    347             self.program = None 
    348  
    349     def __del__(self): 
    350         # type: () -> None 
    351         self.release() 
    352  
    353 # TODO: check that we don't need a destructor for buffers which go out of scope 
     339        return GpuKernel(self, q_vectors) 
     340 
     341    def get_function(self, name): 
     342        # type: (str) -> cuda.Function 
     343        """ 
     344        Fetch the kernel from the environment by name, compiling it if it 
     345        does not already exist. 
     346        """ 
     347        if self._program is None: 
     348            self._prepare_program() 
     349        return self._kernels[name] 
     350 
     351    def _prepare_program(self): 
     352        # type: (str) -> None 
     353        env = environment() 
     354        timestamp = generate.ocl_timestamp(self.info) 
     355        program = env.compile_program( 
     356            self.info.name, 
     357            self.source['opencl'], 
     358            self.dtype, 
     359            self.fast, 
     360            timestamp) 
     361        variants = ['Iq', 'Iqxy', 'Imagnetic'] 
     362        names = [generate.kernel_name(self.info, k) for k in variants] 
     363        functions = [program.get_function(k) for k in names] 
     364        self._kernels = {k: v for k, v in zip(variants, functions)} 
     365        # Keep a handle to program so GC doesn't collect. 
     366        self._program = program 
     367 
     368 
     369# TODO: Check that we don't need a destructor for buffers which go out of scope. 
    354370class GpuInput(object): 
    355371    """ 
     
    373389    def __init__(self, q_vectors, dtype=generate.F32): 
    374390        # type: (List[np.ndarray], np.dtype) -> None 
    375         # TODO: do we ever need double precision q? 
     391        # TODO: Do we ever need double precision q? 
    376392        self.nq = q_vectors[0].size 
    377393        self.dtype = np.dtype(dtype) 
    378394        self.is_2d = (len(q_vectors) == 2) 
    379         # TODO: stretch input based on get_warp() 
    380         # not doing it now since warp depends on kernel, which is not known 
     395        # TODO: Stretch input based on get_warp(). 
     396        # Not doing it now since warp depends on kernel, which is not known 
    381397        # at this point, so instead using 32, which is good on the set of 
    382398        # architectures tested so far. 
    383399        if self.is_2d: 
    384             # Note: 16 rather than 15 because result is 1 longer than input. 
    385             width = ((self.nq+16)//16)*16 
     400            width = ((self.nq+15)//16)*16 
    386401            self.q = np.empty((width, 2), dtype=dtype) 
    387402            self.q[:self.nq, 0] = q_vectors[0] 
    388403            self.q[:self.nq, 1] = q_vectors[1] 
    389404        else: 
    390             # Note: 32 rather than 31 because result is 1 longer than input. 
    391             width = ((self.nq+32)//32)*32 
     405            width = ((self.nq+31)//32)*32 
    392406            self.q = np.empty(width, dtype=dtype) 
    393407            self.q[:self.nq] = q_vectors[0] 
    394408        self.global_size = [self.q.shape[0]] 
    395409        #print("creating inputs of size", self.global_size) 
     410 
     411        # Transfer input value to GPU. 
    396412        self.q_b = cuda.to_device(self.q) 
    397413 
     
    399415        # type: () -> None 
    400416        """ 
    401         Free the memory. 
     417        Free the buffer associated with the q value. 
    402418        """ 
    403419        if self.q_b is not None: 
     
    409425        self.release() 
    410426 
     427 
    411428class GpuKernel(Kernel): 
    412429    """ 
    413430    Callable SAS kernel. 
    414431 
    415     *kernel* is the GpuKernel object to call 
    416  
    417     *model_info* is the module information 
    418  
    419     *q_vectors* is the q vectors at which the kernel should be evaluated 
    420  
    421     *dtype* is the kernel precision 
    422  
    423     The resulting call method takes the *pars*, a list of values for 
    424     the fixed parameters to the kernel, and *pd_pars*, a list of (value,weight) 
    425     vectors for the polydisperse parameters.  *cutoff* determines the 
    426     integration limits: any points with combined weight less than *cutoff* 
    427     will not be calculated. 
     432    *model* is the GpuModel object to call 
     433 
     434    The kernel is derived from :class:`Kernel`, providing the 
     435    :meth:`call_kernel` method to evaluate the kernel for a given set of 
     436    parameters.  Because of the need to move the q values to the GPU before 
     437    evaluation, the kernel is instantiated for a particular set of q vectors, 
     438    and can be called many times without transfering q each time. 
    428439 
    429440    Call :meth:`release` when done with the kernel instance. 
    430441    """ 
    431     def __init__(self, kernel, dtype, model_info, q_vectors): 
    432         # type: (cl.Kernel, np.dtype, ModelInfo, List[np.ndarray]) -> None 
     442    #: SAS model information structure. 
     443    info = None  # type: ModelInfo 
     444    #: Kernel precision. 
     445    dtype = None  # type: np.dtype 
     446    #: Kernel dimensions (1d or 2d). 
     447    dim = ""  # type: str 
     448    #: Calculation results, updated after each call to :meth:`_call_kernel`. 
     449    result = None  # type: np.ndarray 
     450 
     451    def __init__(self, model, q_vectors): 
     452        # type: (GpuModel, List[np.ndarray]) -> None 
     453        dtype = model.dtype 
    433454        self.q_input = GpuInput(q_vectors, dtype) 
    434         self.kernel = kernel 
    435         # F16 isn't sufficient, so don't support it 
     455        self._model = model 
     456 
     457        # Attributes accessed from the outside. 
     458        self.dim = '2d' if self.q_input.is_2d else '1d' 
     459        self.info = model.info 
     460        self.dtype = dtype 
     461 
     462        # Converter to translate input to target type. 
    436463        self._as_dtype = np.float64 if dtype == generate.F64 else np.float32 
    437464 
    438         # attributes accessed from the outside 
    439         self.dim = '2d' if self.q_input.is_2d else '1d' 
    440         self.info = model_info 
    441         self.dtype = dtype 
    442  
    443         # holding place for the returned value 
     465        # Holding place for the returned value. 
    444466        nout = 2 if self.info.have_Fq and self.dim == '1d' else 1 
    445         extra_q = 4  # total weight, form volume, shell volume and R_eff 
    446         self.result = np.empty(self.q_input.nq*nout+extra_q, dtype) 
    447  
    448         # Inputs and outputs for each kernel call 
    449         # Note: res may be shorter than res_b if global_size != nq 
     467        extra_q = 4  # Total weight, form volume, shell volume and R_eff. 
     468        self.result = np.empty(self.q_input.nq*nout + extra_q, dtype) 
     469 
     470        # Allocate result value on GPU. 
    450471        width = ((self.result.size+31)//32)*32 * self.dtype.itemsize 
    451         self.result_b = cuda.mem_alloc(width) 
    452         self._need_release = [self.result_b] 
    453  
    454     def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
    455         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    456         # Arrange data transfer to card 
     472        self._result_b = cuda.mem_alloc(width) 
     473 
     474    def _call_kernel(self, call_details, values, cutoff, magnetic, 
     475                     radius_effective_mode): 
     476        # type: (CallDetails, np.ndarray, float, bool, int) -> np.ndarray 
     477 
     478        # Arrange data transfer to card. 
    457479        details_b = cuda.to_device(call_details.buffer) 
    458480        values_b = cuda.to_device(values) 
    459481 
    460         kernel = self.kernel[1 if magnetic else 0] 
    461         args = [ 
    462             np.uint32(self.q_input.nq), None, None, 
    463             details_b, values_b, self.q_input.q_b, self.result_b, 
    464             self._as_dtype(cutoff), 
    465             np.uint32(effective_radius_type), 
     482        # Setup kernel function and arguments. 
     483        name = 'Iq' if self.dim == '1d' else 'Imagnetic' if magnetic else 'Iqxy' 
     484        kernel = self._model.get_function(name) 
     485        kernel_args = [ 
     486            np.uint32(self.q_input.nq),  # Number of inputs. 
     487            None,  # Placeholder for pd_start. 
     488            None,  # Placeholder for pd_stop. 
     489            details_b,  # Problem definition. 
     490            values_b,  # Parameter values. 
     491            self.q_input.q_b,  # Q values. 
     492            self._result_b,   # Result storage. 
     493            self._as_dtype(cutoff),  # Probability cutoff. 
     494            np.uint32(radius_effective_mode),  # R_eff mode. 
    466495        ] 
    467496        grid = partition(self.q_input.nq) 
    468         #print("Calling OpenCL") 
     497 
     498        # Call kernel and retrieve results. 
     499        #print("Calling CUDA") 
    469500        #call_details.show(values) 
    470         # Call kernel and retrieve results 
    471501        last_nap = time.clock() 
    472502        step = 100000000//self.q_input.nq + 1 
     
    475505            stop = min(start + step, call_details.num_eval) 
    476506            #print("queuing",start,stop) 
    477             args[1:3] = [np.int32(start), np.int32(stop)] 
    478             kernel(*args, **grid) 
     507            kernel_args[1:3] = [np.int32(start), np.int32(stop)] 
     508            kernel(*kernel_args, **grid) 
    479509            if stop < call_details.num_eval: 
    480510                sync() 
    481                 # Allow other processes to run 
     511                # Allow other processes to run. 
    482512                current_time = time.clock() 
    483513                if current_time - last_nap > 0.5: 
     
    485515                    last_nap = current_time 
    486516        sync() 
    487         cuda.memcpy_dtoh(self.result, self.result_b) 
     517        cuda.memcpy_dtoh(self.result, self._result_b) 
    488518        #print("result", self.result) 
    489519 
     
    496526        Release resources associated with the kernel. 
    497527        """ 
    498         for p in self._need_release: 
    499             p.free() 
    500         self._need_release = [] 
     528        self.q_input.release() 
     529        if self._result_b is not None: 
     530            self._result_b.free() 
     531            self._result_b = None 
    501532 
    502533    def __del__(self): 
     
    512543    Note: Maybe context.synchronize() is sufficient. 
    513544    """ 
    514     #return # The following works in C++; don't know what pycuda is doing 
    515     # Create an event with which to synchronize 
     545    # Create an event with which to synchronize. 
    516546    done = cuda.Event() 
    517547 
     
    519549    done.record() 
    520550 
    521     #line added to not hog resources 
     551    # Make sure we don't hog resource while waiting to sync. 
    522552    while not done.query(): 
    523553        time.sleep(0.01) 
     
    525555    # Block until the GPU executes the kernel. 
    526556    done.synchronize() 
     557 
    527558    # Clean up the event; I don't think they can be reused. 
    528559    del done 
  • sasmodels/kerneldll.py

    re44432d ra34b811  
    100100# pylint: enable=unused-import 
    101101 
    102 # Compiler output is a byte stream that needs to be decode in python 3 
     102# Compiler output is a byte stream that needs to be decode in python 3. 
    103103decode = (lambda s: s) if sys.version_info[0] < 3 else (lambda s: s.decode('utf8')) 
    104104 
     
    115115        COMPILER = "tinycc" 
    116116    elif "VCINSTALLDIR" in os.environ: 
    117         # If vcvarsall.bat has been called, then VCINSTALLDIR is in the environment 
    118         # and we can use the MSVC compiler.  Otherwise, if tinycc is available 
    119         # the use it.  Otherwise, hope that mingw is available. 
     117        # If vcvarsall.bat has been called, then VCINSTALLDIR is in the 
     118        # environment and we can use the MSVC compiler.  Otherwise, if 
     119        # tinycc is available then use it.  Otherwise, hope that mingw 
     120        # is available. 
    120121        COMPILER = "msvc" 
    121122    else: 
     
    124125    COMPILER = "unix" 
    125126 
    126 ARCH = "" if ct.sizeof(ct.c_void_p) > 4 else "x86"  # 4 byte pointers on x86 
     127ARCH = "" if ct.sizeof(ct.c_void_p) > 4 else "x86"  # 4 byte pointers on x86. 
    127128if COMPILER == "unix": 
    128     # Generic unix compile 
    129     # On mac users will need the X code command line tools installed 
     129    # Generic unix compile. 
     130    # On Mac users will need the X code command line tools installed. 
    130131    #COMPILE = "gcc-mp-4.7 -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm -lgomp" 
    131132    CC = "cc -shared -fPIC -std=c99 -O2 -Wall".split() 
    132     # add openmp support if not running on a mac 
     133    # Add OpenMP support if not running on a Mac. 
    133134    if sys.platform != "darwin": 
    134         # OpenMP seems to be broken on gcc 5.4.0 (ubuntu 16.04.9) 
     135        # OpenMP seems to be broken on gcc 5.4.0 (ubuntu 16.04.9). 
    135136        # Shut it off for all unix until we can investigate. 
    136137        #CC.append("-fopenmp") 
     
    144145    # vcomp90.dll on the path.  One may be found here: 
    145146    #       C:/Windows/winsxs/x86_microsoft.vc90.openmp*/vcomp90.dll 
    146     # Copy this to the python directory and uncomment the OpenMP COMPILE 
    147     # TODO: remove intermediate OBJ file created in the directory 
    148     # TODO: maybe don't use randomized name for the c file 
    149     # TODO: maybe ask distutils to find MSVC 
     147    # Copy this to the python directory and uncomment the OpenMP COMPILE. 
     148    # TODO: Remove intermediate OBJ file created in the directory. 
     149    # TODO: Maybe don't use randomized name for the c file. 
     150    # TODO: Maybe ask distutils to find MSVC. 
    150151    CC = "cl /nologo /Ox /MD /W3 /GS- /DNDEBUG".split() 
    151152    if "SAS_OPENMP" in os.environ: 
     
    172173ALLOW_SINGLE_PRECISION_DLLS = True 
    173174 
    174 def compile(source, output): 
     175 
     176def compile_model(source, output): 
    175177    # type: (str, str) -> None 
    176178    """ 
     
    183185    logging.info(command_str) 
    184186    try: 
    185         # need shell=True on windows to keep console box from popping up 
     187        # Need shell=True on windows to keep console box from popping up. 
    186188        shell = (os.name == 'nt') 
    187189        subprocess.check_output(command, shell=shell, stderr=subprocess.STDOUT) 
     
    192194        raise RuntimeError("compile failed.  File is in %r"%source) 
    193195 
     196 
    194197def dll_name(model_info, dtype): 
    195198    # type: (ModelInfo, np.dtype) ->  str 
     
    202205    basename += ARCH + ".so" 
    203206 
    204     # Hack to find precompiled dlls 
     207    # Hack to find precompiled dlls. 
    205208    path = joinpath(generate.DATA_PATH, '..', 'compiled_models', basename) 
    206209    if os.path.exists(path): 
     
    242245        raise ValueError("16 bit floats not supported") 
    243246    if dtype == F32 and not ALLOW_SINGLE_PRECISION_DLLS: 
    244         dtype = F64  # Force 64-bit dll 
    245     # Note: dtype may be F128 for long double precision 
     247        dtype = F64  # Force 64-bit dll. 
     248    # Note: dtype may be F128 for long double precision. 
    246249 
    247250    dll = dll_path(model_info, dtype) 
     
    254257        need_recompile = dll_time < newest_source 
    255258    if need_recompile: 
    256         # Make sure the DLL path exists 
     259        # Make sure the DLL path exists. 
    257260        if not os.path.exists(SAS_DLL_PATH): 
    258261            os.makedirs(SAS_DLL_PATH) 
     
    262265        with os.fdopen(system_fd, "w") as file_handle: 
    263266            file_handle.write(source) 
    264         compile(source=filename, output=dll) 
    265         # comment the following to keep the generated c file 
    266         # Note: if there is a syntax error then compile raises an error 
     267        compile_model(source=filename, output=dll) 
     268        # Comment the following to keep the generated C file. 
     269        # Note: If there is a syntax error then compile raises an error 
    267270        # and the source file will not be deleted. 
    268271        os.unlink(filename) 
     
    303306        self.dllpath = dllpath 
    304307        self._dll = None  # type: ct.CDLL 
    305         self._kernels = None # type: List[Callable, Callable] 
     308        self._kernels = None  # type: List[Callable, Callable] 
    306309        self.dtype = np.dtype(dtype) 
    307310 
     
    338341        # type: (List[np.ndarray]) -> DllKernel 
    339342        q_input = PyInput(q_vectors, self.dtype) 
    340         # Note: pickle not supported for DllKernel 
     343        # Note: DLL is lazy loaded. 
    341344        if self._dll is None: 
    342345            self._load_dll() 
     
    358361        self._dll = None 
    359362 
     363 
    360364class DllKernel(Kernel): 
    361365    """ 
     
    379383    def __init__(self, kernel, model_info, q_input): 
    380384        # type: (Callable[[], np.ndarray], ModelInfo, PyInput) -> None 
    381         #,model_info,q_input) 
     385        dtype = q_input.dtype 
     386        self.q_input = q_input 
    382387        self.kernel = kernel 
     388 
     389        # Attributes accessed from the outside. 
     390        self.dim = '2d' if q_input.is_2d else '1d' 
    383391        self.info = model_info 
    384         self.q_input = q_input 
    385         self.dtype = q_input.dtype 
    386         self.dim = '2d' if q_input.is_2d else '1d' 
    387         # leave room for f1/f2 results in case we need to compute beta for 1d models 
     392        self.dtype = dtype 
     393 
     394        # Converter to translate input to target type. 
     395        self._as_dtype = (np.float32 if dtype == generate.F32 
     396                          else np.float64 if dtype == generate.F64 
     397                          else np.float128) 
     398 
     399        # Holding place for the returned value. 
    388400        nout = 2 if self.info.have_Fq else 1 
    389         # +4 for total weight, shell volume, effective radius, form volume 
    390         self.result = np.empty(q_input.nq*nout + 4, self.dtype) 
    391         self.real = (np.float32 if self.q_input.dtype == generate.F32 
    392                      else np.float64 if self.q_input.dtype == generate.F64 
    393                      else np.float128) 
    394  
    395     def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
    396         # type: (CallDetails, np.ndarray, np.ndarray, float, bool, int) -> np.ndarray 
     401        extra_q = 4  # Total weight, form volume, shell volume and R_eff. 
     402        self.result = np.empty(self.q_input.nq*nout + extra_q, dtype) 
     403 
     404    def _call_kernel(self, call_details, values, cutoff, magnetic, 
     405                     radius_effective_mode): 
     406        # type: (CallDetails, np.ndarray, float, bool, int) -> np.ndarray 
     407 
     408        # Setup kernel function and arguments. 
    397409        kernel = self.kernel[1 if magnetic else 0] 
    398         args = [ 
    399             self.q_input.nq, # nq 
    400             None, # pd_start 
    401             None, # pd_stop pd_stride[MAX_PD] 
    402             call_details.buffer.ctypes.data, # problem 
    403             values.ctypes.data,  # pars 
    404             self.q_input.q.ctypes.data, # q 
    405             self.result.ctypes.data,   # results 
    406             self.real(cutoff), # cutoff 
    407             effective_radius_type, # cutoff 
     410        kernel_args = [ 
     411            self.q_input.nq,  # Number of inputs. 
     412            None,  # Placeholder for pd_start. 
     413            None,  # Placeholder for pd_stop. 
     414            call_details.buffer.ctypes.data,  # Problem definition. 
     415            values.ctypes.data,  # Parameter values. 
     416            self.q_input.q.ctypes.data,  # Q values. 
     417            self.result.ctypes.data,   # Result storage. 
     418            self._as_dtype(cutoff),  # Probability cutoff. 
     419            radius_effective_mode,  # R_eff mode. 
    408420        ] 
     421 
     422        # Call kernel and retrieve results. 
    409423        #print("Calling DLL") 
    410424        #call_details.show(values) 
    411425        step = 100 
     426        # TODO: Do we need the explicit sleep like the OpenCL and CUDA loops? 
    412427        for start in range(0, call_details.num_eval, step): 
    413428            stop = min(start + step, call_details.num_eval) 
    414             args[1:3] = [start, stop] 
    415             kernel(*args) # type: ignore 
     429            kernel_args[1:3] = [start, stop] 
     430            kernel(*kernel_args) # type: ignore 
    416431 
    417432    def release(self): 
    418433        # type: () -> None 
    419434        """ 
    420         Release any resources associated with the kernel. 
     435        Release resources associated with the kernel. 
    421436        """ 
    422         self.q_input.release() 
     437        # TODO: OpenCL/CUDA allocate q_input in __init__ and free it in release. 
     438        # Should we be doing the same for DLL? 
     439        #self.q_input.release() 
     440        pass 
     441 
     442    def __del__(self): 
     443        # type: () -> None 
     444        self.release() 
  • sasmodels/kernelpy.py

    raa8c6e0 ra34b811  
    1616    from numpy import cbrt 
    1717except ImportError: 
    18     def cbrt(x): return x ** (1.0/3.0) 
     18    def cbrt(x): 
     19        """Return cubed root of x.""" 
     20        return x ** (1.0/3.0) 
    1921 
    2022from .generate import F64 
     
    3335logger = logging.getLogger(__name__) 
    3436 
     37 
    3538class PyModel(KernelModel): 
    3639    """ 
     
    3841    """ 
    3942    def __init__(self, model_info): 
    40         # Make sure Iq is available and vectorized 
     43        # Make sure Iq is available and vectorized. 
    4144        _create_default_functions(model_info) 
    4245        self.info = model_info 
    4346        self.dtype = np.dtype('d') 
    44         logger.info("make python model " + self.info.name) 
     47        logger.info("make python model %s", self.info.name) 
    4548 
    4649    def make_kernel(self, q_vectors): 
     50        """Instantiate the python kernel with input *q_vectors*""" 
    4751        q_input = PyInput(q_vectors, dtype=F64) 
    4852        return PyKernel(self.info, q_input) 
     
    5357        """ 
    5458        pass 
     59 
    5560 
    5661class PyInput(object): 
     
    9196        self.q = None 
    9297 
     98 
    9399class PyKernel(Kernel): 
    94100    """ 
     
    131137        parameter_vector = np.empty(len(partable.call_parameters)-2, 'd') 
    132138 
    133         # Create views into the array to hold the arguments 
     139        # Create views into the array to hold the arguments. 
    134140        offset = 0 
    135141        kernel_args, volume_args = [], [] 
     
    166172        volume = model_info.form_volume 
    167173        shell = model_info.shell_volume 
    168         radius = model_info.effective_radius 
     174        radius = model_info.radius_effective 
    169175        self._volume = ((lambda: (shell(*volume_args), volume(*volume_args))) if shell and volume 
    170176                        else (lambda: [volume(*volume_args)]*2) if volume 
     
    174180                        else (lambda mode: 1.0)) 
    175181 
    176  
    177  
    178     def _call_kernel(self, call_details, values, cutoff, magnetic, effective_radius_type): 
     182    def _call_kernel(self, call_details, values, cutoff, magnetic, radius_effective_mode): 
    179183        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    180184        if magnetic: 
     
    182186        #print("Calling python kernel") 
    183187        #call_details.show(values) 
    184         radius = ((lambda: 0.0) if effective_radius_type == 0 
    185                   else (lambda: self._radius(effective_radius_type))) 
     188        radius = ((lambda: 0.0) if radius_effective_mode == 0 
     189                  else (lambda: self._radius(radius_effective_mode))) 
    186190        self.result = _loops( 
    187191            self._parameter_vector, self._form, self._volume, radius, 
     
    195199        self.q_input.release() 
    196200        self.q_input = None 
     201 
    197202 
    198203def _loops(parameters,    # type: np.ndarray 
     
    254259        total = np.zeros(nq, 'd') 
    255260        for loop_index in range(call_details.num_eval): 
    256             # update polydispersity parameter values 
     261            # Update polydispersity parameter values. 
    257262            if p0_index == p0_length: 
    258263                pd_index = (loop_index//pd_stride)%pd_length 
     
    265270            p0_index += 1 
    266271            if weight > cutoff: 
    267                 # Call the scattering function 
     272                # Call the scattering function. 
    268273                # Assume that NaNs are only generated if the parameters are bad; 
    269274                # exclude all q for that NaN.  Even better would be to have an 
     
    273278                    continue 
    274279 
    275                 # update value and norm 
     280                # Update value and norm. 
    276281                total += weight * Iq 
    277282                weight_norm += weight 
     
    293298    any functions that are not already marked as vectorized. 
    294299    """ 
    295     # Note: must call create_vector_Iq before create_vector_Iqxy 
     300    # Note: Must call create_vector_Iq before create_vector_Iqxy. 
    296301    _create_vector_Iq(model_info) 
    297302    _create_vector_Iqxy(model_info) 
  • sasmodels/list_pars.py

    r2d81cfe rb297ba9  
    1616from .compare import columnize 
    1717 
    18 def find_pars(type=None): 
     18def find_pars(kind=None): 
    1919    """ 
    2020    Find all parameters in all models. 
     
    2626        model_info = load_model_info(name) 
    2727        for p in model_info.parameters.kernel_parameters: 
    28             if type is None or p.type == type: 
     28            if kind is None or p.type == kind: 
    2929                partable.setdefault(p.name, []) 
    3030                partable[p.name].append(name) 
    3131    return partable 
    3232 
    33 def list_pars(names_only=True, type=None): 
     33def list_pars(names_only=True, kind=None): 
    3434    """ 
    3535    Print all parameters in all models. 
     
    3838    occurs in. 
    3939    """ 
    40     partable = find_pars(type) 
     40    partable = find_pars(kind) 
    4141    if names_only: 
    4242        print(columnize(list(sorted(partable.keys())))) 
     
    5757        help="list models which use this argument") 
    5858    parser.add_argument( 
    59         'type', default="any", nargs='?', 
     59        'kind', default="any", nargs='?', 
    6060        metavar="volume|orientation|sld|none|any", 
    6161        choices=['volume', 'orientation', 'sld', None, 'any'], 
    62         type=lambda v: None if v == 'any' else '' if v == 'none' else v, 
    63         help="only list arguments of the given type") 
     62        type=lambda v: None if v == 'any' else None if v == 'none' else v, 
     63        help="only list arguments of the given kind") 
    6464    args = parser.parse_args() 
    6565 
    66     list_pars(names_only=not args.verbose, type=args.type) 
     66    list_pars(names_only=not args.verbose, kind=args.kind) 
    6767 
    6868if __name__ == "__main__": 
  • sasmodels/mixture.py

    r39a06c9 rb297ba9  
    120120 
    121121    def random(): 
     122        """Random set of model parameters for mixture model""" 
    122123        combined_pars = {} 
    123124        for k, part in enumerate(parts): 
     
    149150 
    150151class MixtureModel(KernelModel): 
     152    """ 
     153    Model definition for mixture of models. 
     154    """ 
    151155    def __init__(self, model_info, parts): 
    152156        # type: (ModelInfo, List[KernelModel]) -> None 
     
    165169        kernels = [part.make_kernel(q_vectors) for part in self.parts] 
    166170        return MixtureKernel(self.info, kernels) 
     171    make_kernel.__doc__ = KernelModel.make_kernel.__doc__ 
    167172 
    168173    def release(self): 
    169174        # type: () -> None 
    170         """ 
    171         Free resources associated with the model. 
    172         """ 
     175        """Free resources associated with the model.""" 
    173176        for part in self.parts: 
    174177            part.release() 
     178    release.__doc__ = KernelModel.release.__doc__ 
    175179 
    176180 
    177181class MixtureKernel(Kernel): 
     182    """ 
     183    Instantiated kernel for mixture of models. 
     184    """ 
    178185    def __init__(self, model_info, kernels): 
    179186        # type: (ModelInfo, List[Kernel]) -> None 
     
    185192        self.results = []  # type: List[np.ndarray] 
    186193 
    187     def __call__(self, call_details, values, cutoff, magnetic): 
     194    def Iq(self, call_details, values, cutoff, magnetic): 
    188195        # type: (CallDetails, np.ndarray, np.ndarry, float, bool) -> np.ndarray 
    189196        scale, background = values[0:2] 
     
    191198        # remember the parts for plotting later 
    192199        self.results = []  # type: List[np.ndarray] 
    193         parts = MixtureParts(self.info, self.kernels, call_details, values) 
     200        parts = _MixtureParts(self.info, self.kernels, call_details, values) 
    194201        for kernel, kernel_details, kernel_values in parts: 
    195202            #print("calling kernel", kernel.info.name) 
     
    208215        return scale*total + background 
    209216 
     217    Iq.__doc__ = Kernel.Iq.__doc__ 
     218    __call__ = Iq 
     219 
    210220    def release(self): 
    211221        # type: () -> None 
     222        """Free resources associated with the kernel.""" 
    212223        for k in self.kernels: 
    213224            k.release() 
    214225 
    215226 
    216 class MixtureParts(object): 
     227# Note: _MixtureParts doesn't implement iteration correctly, and only allows 
     228# a single iterator to be active at once.  It doesn't matter in this case 
     229# since _MixtureParts is only used in one place, but it is not clean style. 
     230class _MixtureParts(object): 
     231    """ 
     232    Mixture component iterator. 
     233    """ 
    217234    def __init__(self, model_info, kernels, call_details, values): 
    218235        # type: (ModelInfo, List[Kernel], CallDetails, np.ndarray) -> None 
     
    223240        self.values = values 
    224241        self.spin_index = model_info.parameters.npars + 2 
     242        # The following are redefined by __iter__, but set them here so that 
     243        # lint complains a little less. 
     244        self.part_num = -1 
     245        self.par_index = -1 
     246        self.mag_index = -1 
    225247        #call_details.show(values) 
    226248 
  • sasmodels/model_test.py

    r5024a56 r8795b6f  
    55Usage:: 
    66 
    7     python -m sasmodels.model_test [opencl|cuda|dll] model1 model2 ... 
    8  
    9     if model1 is 'all', then all except the remaining models will be tested 
     7    python -m sasmodels.model_test [opencl|cuda|dll|all] model1 model2 ... 
     8 
     9If model1 is 'all', then all except the remaining models will be tested. 
     10Subgroups are also possible, such as 'py', 'single' or '1d'.  See 
     11:func:`core.list_models` for details. 
    1012 
    1113Each model is tested using the default parameters at q=0.1, (qx, qy)=(0.1, 0.1), 
     
    4547from __future__ import print_function 
    4648 
     49import argparse 
    4750import sys 
    4851import unittest 
     
    5861import numpy as np  # type: ignore 
    5962 
    60 from . import core 
    6163from .core import list_models, load_model_info, build_model 
    6264from .direct_model import call_kernel, call_Fq 
     
    8991    suite = unittest.TestSuite() 
    9092 
    91     if models[0] in core.KINDS: 
     93    try: 
     94        # See if the first model parses as a model group 
     95        group = list_models(models[0]) 
    9296        skip = models[1:] 
    93         models = list_models(models[0]) 
    94     else: 
     97        models = group 
     98    except Exception: 
    9599        skip = [] 
    96100    for model_name in models: 
     
    132136        test_method_name = "test_%s_python" % model_info.id 
    133137        test = ModelTestCase(test_name, model_info, 
    134                                 test_method_name, 
    135                                 platform="dll",  # so that 
    136                                 dtype="double", 
    137                                 stash=stash) 
     138                             test_method_name, 
     139                             platform="dll",  # so that 
     140                             dtype="double", 
     141                             stash=stash) 
    138142        suite.addTest(test) 
    139143    else:   # kernel implemented in C 
     
    144148            test_method_name = "test_%s_dll" % model_info.id 
    145149            test = ModelTestCase(test_name, model_info, 
    146                                     test_method_name, 
    147                                     platform="dll", 
    148                                     dtype="double", 
    149                                     stash=stash) 
     150                                 test_method_name, 
     151                                 platform="dll", 
     152                                 dtype="double", 
     153                                 stash=stash) 
    150154            suite.addTest(test) 
    151155 
     
    159163            # presence of *single=False* in the model file. 
    160164            test = ModelTestCase(test_name, model_info, 
    161                                     test_method_name, 
    162                                     platform="ocl", dtype=None, 
    163                                     stash=stash) 
     165                                 test_method_name, 
     166                                 platform="ocl", dtype=None, 
     167                                 stash=stash) 
    164168            #print("defining", test_name) 
    165169            suite.addTest(test) 
     
    167171        # test using cuda if desired and available 
    168172        if 'cuda' in loaders and use_cuda(): 
    169             test_name = "%s-cuda"%model_name 
     173            test_name = "%s-cuda" % model_info.id 
    170174            test_method_name = "test_%s_cuda" % model_info.id 
    171175            # Using dtype=None so that the models that are only 
     
    174178            # presence of *single=False* in the model file. 
    175179            test = ModelTestCase(test_name, model_info, 
    176                                     test_method_name, 
    177                                     platform="cuda", dtype=None, 
    178                                     stash=stash) 
     180                                 test_method_name, 
     181                                 platform="cuda", dtype=None, 
     182                                 stash=stash) 
    179183            #print("defining", test_name) 
    180184            suite.addTest(test) 
     
    236240                    s_name = pars.pop('@S') 
    237241                    ps_test = [pars] + list(test[1:]) 
     242                    #print("PS TEST PARAMS!!!",ps_test) 
    238243                    # build the P@S model 
    239244                    s_info = load_model_info(s_name) 
     
    242247                                           platform=self.platform) 
    243248                    # run the tests 
     249                    #self.info = ps_model.info 
     250                    #print("SELF.INFO PARAMS!!!",[p.id for p in self.info.parameters.call_parameters]) 
     251                    #print("PS MODEL PARAMETERS:",[p.id for p in ps_model.info.parameters.call_parameters]) 
    244252                    results.append(self.run_one(ps_model, ps_test)) 
    245253 
     
    299307            """Run a single test case.""" 
    300308            user_pars, x, y = test[:3] 
    301             pars = expand_pars(self.info.parameters, user_pars) 
    302             invalid = invalid_pars(self.info.parameters, pars) 
     309            #print("PS MODEL PARAMETERS:",[p.id for p in model.info.parameters.call_parameters]) 
     310            pars = expand_pars(model.info.parameters, user_pars) 
     311            invalid = invalid_pars(model.info.parameters, pars) 
    303312            if invalid: 
    304313                raise ValueError("Unknown parameters in test: " + ", ".join(invalid)) 
     
    324333            else: 
    325334                y1 = y 
    326                 y2 = test[3] if not isinstance(test[3], list) else [test[3]] 
    327                 F1, F2, R_eff, volume, volume_ratio = call_Fq(kernel, pars) 
    328                 if F1 is not None:  # F1 is none for models with Iq instead of Fq 
    329                     self._check_vectors(x, y1, F1, 'F') 
    330                 self._check_vectors(x, y2, F2, 'F^2') 
     335                y2 = test[3] if isinstance(test[3], list) else [test[3]] 
     336                F, Fsq, R_eff, volume, volume_ratio = call_Fq(kernel, pars) 
     337                if F is not None:  # F is none for models with Iq instead of Fq 
     338                    self._check_vectors(x, y1, F, 'F') 
     339                self._check_vectors(x, y2, Fsq, 'F^2') 
    331340                self._check_scalar(test[4], R_eff, 'R_eff') 
    332341                self._check_scalar(test[5], volume, 'volume') 
    333342                self._check_scalar(test[6], volume_ratio, 'form:shell ratio') 
    334                 return F2 
     343                return Fsq 
    335344 
    336345        def _check_scalar(self, target, actual, name): 
    337             if target is None: 
    338                 # smoke test --- make sure it runs and produces a value 
    339                 self.assertTrue(not np.isnan(actual), 
    340                                 'invalid %s: %s' % (name, actual)) 
    341             elif np.isnan(target): 
    342                 # make sure nans match 
    343                 self.assertTrue(np.isnan(actual), 
    344                                 '%s: expected:%s; actual:%s' 
    345                                 % (name, target, actual)) 
    346             else: 
    347                 # is_near does not work for infinite values, so also test 
    348                 # for exact values. 
    349                 self.assertTrue(target == actual or is_near(target, actual, 5), 
    350                                 '%s: expected:%s; actual:%s' 
    351                                 % (name, target, actual)) 
     346            self.assertTrue(is_near(target, actual, 5), 
     347                            '%s: expected:%s; actual:%s' 
     348                            % (name, target, actual)) 
    352349 
    353350        def _check_vectors(self, x, target, actual, name='I'): 
     
    359356                             '%s(...) returned wrong length'%name) 
    360357            for xi, yi, actual_yi in zip(x, target, actual): 
    361                 if yi is None: 
    362                     # smoke test --- make sure it runs and produces a value 
    363                     self.assertTrue(not np.isnan(actual_yi), 
    364                                     'invalid %s(%s): %s' % (name, xi, actual_yi)) 
    365                 elif np.isnan(yi): 
    366                     # make sure nans match 
    367                     self.assertTrue(np.isnan(actual_yi), 
    368                                     '%s(%s): expected:%s; actual:%s' 
    369                                     % (name, xi, yi, actual_yi)) 
    370                 else: 
    371                     # is_near does not work for infinite values, so also test 
    372                     # for exact values. 
    373                     self.assertTrue(yi == actual_yi or is_near(yi, actual_yi, 5), 
    374                                     '%s(%s); expected:%s; actual:%s' 
    375                                     % (name, xi, yi, actual_yi)) 
     358                self.assertTrue(is_near(yi, actual_yi, 5), 
     359                                '%s(%s): expected:%s; actual:%s' 
     360                                % (name, xi, target, actual)) 
    376361 
    377362    return ModelTestCase 
     
    385370    invalid = [] 
    386371    for par in sorted(pars.keys()): 
    387         # special handling of R_eff mode, which is not a usual parameter 
     372        # Ignore the R_eff mode parameter when checking for valid parameters. 
     373        # It is an allowed parameter for a model even though it does not exist 
     374        # in the parameter table.  The call_Fq() function pops it from the 
     375        # parameter list and sends it directly to kernel.Fq(). 
    388376        if par == product.RADIUS_MODE_ID: 
    389377            continue 
     
    401389    """ 
    402390    Returns true if *actual* is within *digits* significant digits of *target*. 
    403     """ 
    404     import math 
    405     shift = 10**math.ceil(math.log10(abs(target))) 
    406     return abs(target-actual)/shift < 1.5*10**-digits 
     391 
     392    *taget* zero and inf should match *actual* zero and inf.  If you want to 
     393    accept eps for zero, choose a value such as 1e-10, which must match up to 
     394    +/- 1e-15 when *digits* is the default value of 5. 
     395 
     396    If *target* is None, then just make sure that *actual* is not NaN. 
     397 
     398    If *target* is NaN, make sure *actual* is NaN. 
     399    """ 
     400    if target is None: 
     401        # target is None => actual cannot be NaN 
     402        return not np.isnan(actual) 
     403    elif target == 0.: 
     404        # target is 0. => actual must be 0. 
     405        # Note: if small values are allowed, then use maybe test zero against eps instead? 
     406        return actual == 0. 
     407    elif np.isfinite(target): 
     408        shift = np.ceil(np.log10(abs(target))) 
     409        return abs(target-actual) < 1.5*10**(shift-digits) 
     410    elif target == actual: 
     411        # target is inf => actual must be inf of same sign 
     412        return True 
     413    else: 
     414        # target is NaN => actual must be NaN 
     415        return np.isnan(target) == np.isnan(actual) 
    407416 
    408417# CRUFT: old interface; should be deprecated and removed 
    409418def run_one(model_name): 
     419    # type: (str) -> str 
     420    """ 
     421    [Deprecated] Run the tests associated with *model_name*. 
     422 
     423    Use the following instead:: 
     424 
     425        succss, output = check_model(load_model_info(model_name)) 
     426    """ 
    410427    # msg = "use check_model(model_info) rather than run_one(model_name)" 
    411428    # warnings.warn(msg, category=DeprecationWarning, stacklevel=2) 
     
    416433        return output 
    417434 
    418     success, output = check_model(model_info) 
     435    _, output = check_model(model_info) 
    419436    return output 
    420437 
    421438def check_model(model_info): 
    422     # type: (ModelInfo) -> str 
     439    # type: (ModelInfo) -> Tuple[bool, str] 
    423440    """ 
    424441    Run the tests for a single model, capturing the output. 
     
    474491 
    475492 
    476 def main(*models): 
    477     # type: (*str) -> int 
    478     """ 
    479     Run tests given is models. 
    480  
    481     Returns 0 if success or 1 if any tests fail. 
    482     """ 
    483     try: 
    484         from xmlrunner import XMLTestRunner as TestRunner 
    485         test_args = {'output': 'logs'} 
    486     except ImportError: 
    487         from unittest import TextTestRunner as TestRunner 
    488         test_args = {} 
    489  
    490     if models and models[0] == '-v': 
    491         verbosity = 2 
    492         models = models[1:] 
    493     else: 
    494         verbosity = 1 
    495     if models and models[0] == 'opencl': 
    496         if not use_opencl(): 
    497             print("opencl is not available") 
    498             return 1 
    499         loaders = ['opencl'] 
    500         models = models[1:] 
    501     elif models and models[0] == 'cuda': 
    502         if not use_cuda(): 
    503             print("cuda is not available") 
    504             return 1 
    505         loaders = ['cuda'] 
    506         models = models[1:] 
    507     elif models and models[0] == 'dll': 
    508         # TODO: test if compiler is available? 
    509         loaders = ['dll'] 
    510         models = models[1:] 
    511     else: 
    512         loaders = ['dll'] 
    513         if use_opencl(): 
    514             loaders.append('opencl') 
    515         if use_cuda(): 
    516             loaders.append('cuda') 
    517     if not models: 
    518         print("""\ 
    519 usage: 
    520   python -m sasmodels.model_test [-v] [opencl|cuda|dll] model1 model2 ... 
    521  
    522 If -v is included on the command line, then use verbose output. 
    523  
    524 If no platform is specified, then models will be tested with dll, and 
    525 if available, OpenCL and CUDA; the compute target is ignored for pure python models. 
    526  
    527 If model1 is 'all', then all except the remaining models will be tested. 
    528  
    529 """) 
    530  
    531         return 1 
    532  
    533     runner = TestRunner(verbosity=verbosity, **test_args) 
    534     result = runner.run(make_suite(loaders, models)) 
    535     return 1 if result.failures or result.errors else 0 
    536  
    537  
    538493def model_tests(): 
    539494    # type: () -> Iterator[Callable[[], None]] 
     
    549504        loaders.append('cuda') 
    550505    tests = make_suite(loaders, ['all']) 
    551     def build_test(test): 
     506    def _build_test(test): 
    552507        # In order for nosetest to show the test name, wrap the test.run_all 
    553508        # instance in function that takes the test name as a parameter which 
     
    567522 
    568523    for test in tests: 
    569         yield build_test(test) 
     524        yield _build_test(test) 
     525 
     526 
     527def main(): 
     528    # type: (*str) -> int 
     529    """ 
     530    Run tests given is models. 
     531 
     532    Returns 0 if success or 1 if any tests fail. 
     533    """ 
     534    try: 
     535        from xmlrunner import XMLTestRunner as TestRunner 
     536        test_args = {'output': 'logs'} 
     537    except ImportError: 
     538        from unittest import TextTestRunner as TestRunner 
     539        test_args = {} 
     540 
     541    parser = argparse.ArgumentParser(description="Test SasModels Models") 
     542    parser.add_argument("-v", "--verbose", action="store_const", 
     543                        default=1, const=2, help="Use verbose output") 
     544    parser.add_argument("-e", "--engine", default="all", 
     545                        help="Engines on which to run the test.  " 
     546                        "Valid values are opencl, cuda, dll, and all. " 
     547                        "Defaults to all if no value is given") 
     548    parser.add_argument("models", nargs="*", 
     549                        help="The names of the models to be tested.  " 
     550                        "If the first model is 'all', then all but the listed " 
     551                        "models will be tested.  See core.list_models() for " 
     552                        "names of other groups, such as 'py' or 'single'.") 
     553    opts = parser.parse_args() 
     554 
     555    if opts.engine == "opencl": 
     556        if not use_opencl(): 
     557            print("opencl is not available") 
     558            return 1 
     559        loaders = ['opencl'] 
     560    elif opts.engine == "dll": 
     561        loaders = ["dll"] 
     562    elif opts.engine == "cuda": 
     563        if not use_cuda(): 
     564            print("cuda is not available") 
     565            return 1 
     566        loaders = ['cuda'] 
     567    elif opts.engine == "all": 
     568        loaders = ['dll'] 
     569        if use_opencl(): 
     570            loaders.append('opencl') 
     571        if use_cuda(): 
     572            loaders.append('cuda') 
     573    else: 
     574        print("unknown engine " + opts.engine) 
     575        return 1 
     576 
     577    runner = TestRunner(verbosity=opts.verbose, **test_args) 
     578    result = runner.run(make_suite(loaders, opts.models)) 
     579    return 1 if result.failures or result.errors else 0 
    570580 
    571581 
    572582if __name__ == "__main__": 
    573     sys.exit(main(*sys.argv[1:])) 
     583    sys.exit(main()) 
  • sasmodels/modelinfo.py

    rc1799d3 ra34b811  
    265265    other sld parameters. The volume parameters are used for calls 
    266266    to form_volume within the kernel (required for volume normalization), 
    267     to shell_volume (for hollow shapes), and to effective_radius (for 
     267    to shell_volume (for hollow shapes), and to radius_effective (for 
    268268    structure factor interactions) respectively. 
    269269 
     
    290290    example, might be used to set the value of a shape parameter. 
    291291 
    292     These values are set by :func:`make_parameter_table` and 
     292    Control parameters are used for variant models such as :ref:`rpa` which 
     293    have different cases with different parameters, as well as models 
     294    like :ref:`spherical_sld` with its user defined number of shells. 
     295    The control parameter should appear in the parameter table along with the 
     296    parameters it is is controlling.  For variant models, use *[CASES]* in 
     297    place of the parameter limits Within the parameter definition table, 
     298    with case names such as:: 
     299 
     300         CASES = ["diblock copolymer", "triblock copolymer", ...] 
     301 
     302    This should give *limits=[[case1, case2, ...]]*, but the model loader 
     303    translates it to *limits=[0, len(CASES)-1]*, and adds *choices=CASES* to 
     304    the :class:`Parameter` definition. Note that models can use a list of 
     305    cases as a parameter without it being a control parameter.  Either way, 
     306    the parameter is sent to the model evaluator as *float(choice_num)*, 
     307    where choices are numbered from 0. :meth:`ModelInfo.get_hidden_parameters` 
     308    will determine which parameers to display. 
     309 
     310    The class contructor should not be called directly, but instead the 
     311    parameter table is built using :func:`make_parameter_table` and 
    293312    :func:`parse_parameter` therein. 
    294313    """ 
     
    822841    info.structure_factor = getattr(kernel_module, 'structure_factor', False) 
    823842    # TODO: find Fq by inspection 
    824     info.effective_radius_type = getattr(kernel_module, 'effective_radius_type', None) 
     843    info.radius_effective_modes = getattr(kernel_module, 'radius_effective_modes', None) 
    825844    info.have_Fq = getattr(kernel_module, 'have_Fq', False) 
    826845    info.profile_axes = getattr(kernel_module, 'profile_axes', ['x', 'y']) 
     
    829848    info.source = getattr(kernel_module, 'source', []) 
    830849    info.c_code = getattr(kernel_module, 'c_code', None) 
    831     info.effective_radius = getattr(kernel_module, 'effective_radius', None) 
     850    info.radius_effective = getattr(kernel_module, 'radius_effective', None) 
    832851    # TODO: check the structure of the tests 
    833852    info.tests = getattr(kernel_module, 'tests', []) 
     
    845864    info.single = getattr(kernel_module, 'single', not callable(info.Iq)) 
    846865    info.random = getattr(kernel_module, 'random', None) 
    847  
    848     # multiplicity info 
    849     control_pars = [p.id for p in parameters.kernel_parameters if p.is_control] 
    850     default_control = control_pars[0] if control_pars else None 
    851     info.control = getattr(kernel_module, 'control', default_control) 
    852866    info.hidden = getattr(kernel_module, 'hidden', None) # type: ignore 
     867 
     868    # Set control flag for explicitly set parameters, e.g., in the RPA model. 
     869    control = getattr(kernel_module, 'control', None) 
     870    if control is not None: 
     871        parameters[control].is_control = True 
    853872 
    854873    if callable(info.Iq) and parameters.has_2d: 
     
    903922    #: *sphere*hardsphere* or *cylinder+sphere*. 
    904923    composition = None      # type: Optional[Tuple[str, List[ModelInfo]]] 
    905     #: Name of the control parameter for a variant model such as :ref:`rpa`. 
    906     #: The *control* parameter should appear in the parameter table, with 
    907     #: limits defined as *[CASES]*, for case names such as 
    908     #: *CASES = ["diblock copolymer", "triblock copolymer", ...]*. 
    909     #: This should give *limits=[[case1, case2, ...]]*, but the 
    910     #: model loader translates this to *limits=[0, len(CASES)-1]*, and adds 
    911     #: *choices=CASES* to the :class:`Parameter` definition. Note that 
    912     #: models can use a list of cases as a parameter without it being a 
    913     #: control parameter.  Either way, the parameter is sent to the model 
    914     #: evaluator as *float(choice_num)*, where choices are numbered from 0. 
    915     #: See also :attr:`hidden`. 
    916     control = None          # type: str 
    917924    #: Different variants require different parameters.  In order to show 
    918925    #: just the parameters needed for the variant selected by :attr:`control`, 
     
    954961    #: List of options for computing the effective radius of the shape, 
    955962    #: or None if the model is not usable as a form factor model. 
    956     effective_radius_type = None   # type: List[str] 
     963    radius_effective_modes = None   # type: List[str] 
    957964    #: List of C source files used to define the model.  The source files 
    958965    #: should define the *Iq* function, and possibly *Iqac* or *Iqabc* if the 
     
    978985    #: C code, either defined as a string, or in the sources. 
    979986    shell_volume = None      # type: Union[None, str, Callable[[np.ndarray], float]] 
     987    #: Computes the effective radius of the shape given the volume parameters. 
     988    #: Only needed for models defined in python that can be used for 
     989    #: monodisperse approximation for non-dilute solutions, P@S.  The first 
     990    #: argument is the integer effective radius mode, with default 0. 
     991    radius_effective = None  # type: Union[None, Callable[[int, np.ndarray], float]] 
    980992    #: Returns *I(q, a, b, ...)* for parameters *a*, *b*, etc. defined 
    981993    #: by the parameter table.  *Iq* can be defined as a python function, or 
     
    9891001    #: will return *I(q, a, b, ...)*.  Multiplicity parameters are sent as 
    9901002    #: pointers to doubles.  Constants in floating point expressions should 
    991     #: include the decimal point. See :mod:`generate` for more details. 
    992     Iq = None               # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
     1003    #: include the decimal point. See :mod:`generate` for more details. If 
     1004    #: *have_Fq* is True, then Iq should return an interleaved array of 
     1005    #: $[\sum F(q_1), \sum F^2(q_1), \ldots, \sum F(q_n), \sum F^2(q_n)]$. 
     1006    Iq = None               # type: Union[None, str, Callable[[...], np.ndarray]] 
     1007    #: Returns *I(qx, qy, a, b, ...)*.  The interface follows :attr:`Iq`. 
     1008    Iqxy = None             # type: Union[None, str, Callable[[...], np.ndarray]] 
    9931009    #: Returns *I(qab, qc, a, b, ...)*.  The interface follows :attr:`Iq`. 
    994     Iqac = None             # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
     1010    Iqac = None             # type: Union[None, str, Callable[[...], np.ndarray]] 
    9951011    #: Returns *I(qa, qb, qc, a, b, ...)*.  The interface follows :attr:`Iq`. 
    996     Iqabc = None            # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
     1012    Iqabc = None            # type: Union[None, str, Callable[[...], np.ndarray]] 
    9971013    #: Returns *I(qx, qy, a, b, ...)*.  The interface follows :attr:`Iq`. 
    998     Imagnetic = None        # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
     1014    Imagnetic = None        # type: Union[None, str, Callable[[...], np.ndarray]] 
    9991015    #: Returns a model profile curve *x, y*.  If *profile* is defined, this 
    10001016    #: curve will appear in response to the *Show* button in SasView.  Use 
  • sasmodels/models/_spherepy.py

    r304c775 ra34b811  
    3333to the output of the software provided by the NIST (Kline, 2006). 
    3434 
    35  
    3635References 
    3736---------- 
     
    4039John Wiley and Sons, New York, (1955) 
    4140 
     41Source 
     42------ 
     43 
     44`_spherepy.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/_spherepy.py>`_ 
     45`sphere.c <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/sphere.c>`_ 
     46 
     47Authorship and Verification 
     48---------------------------- 
     49 
     50* **Author: P Kienzle** 
     51* **Last Modified by:** 
    4252* **Last Reviewed by:** S King and P Parker **Date:** 2013/09/09 and 2014/01/06 
     53* **Source added by :** Steve King **Date:** March 25, 2019 
    4354""" 
    4455 
     
    6980 
    7081def form_volume(radius): 
     82    """Calculate volume for sphere""" 
    7183    return 1.333333333333333 * pi * radius ** 3 
    7284 
    73 def effective_radius(mode, radius): 
    74     return radius 
     85def radius_effective(mode, radius): 
     86    """Calculate R_eff for sphere""" 
     87    return radius if mode else 0. 
    7588 
    7689def Iq(q, sld, sld_solvent, radius): 
     90    """Calculate I(q) for sphere""" 
    7791    #print "q",q 
    7892    #print "sld,r",sld,sld_solvent,radius 
  • sasmodels/models/adsorbed_layer.py

    r2d81cfe r0507e09  
    4949   Layers*, *Macromol. Symp.*, 190 (2002) 33-42. 
    5050 
     51Source 
     52------ 
     53 
     54`adsorbed_layer.py <https://github.com/SasView/sasmodels/blob/master/sasmodels/models/adsorbed_layer.py>`_ 
     55 
    5156Authorship and Verification 
    5257---------------------------- 
     
    5560* **Last Modified by:** Paul Kienzle **Date:** April 14, 2016 
    5661* **Last Reviewed by:** Steve King **Date:** March 18, 2016 
     62* **Source added by :** Steve King **Date:** March 25, 2019 
    5763""" 
    5864 
     
    8793def Iq(q, second_moment, adsorbed_amount, density_shell, radius, 
    8894       volfraction, sld_shell, sld_solvent): 
     95    """Return I(q) for adsorbed layer model.""" 
    8996    with errstate(divide='ignore'): 
    9097        aa = ((sld_shell - sld_solvent)/density_shell * adsorbed_amount) / q