Changeset d86f0fc in sasmodels for sasmodels/jitter.py


Ignore:
Timestamp:
Mar 26, 2018 3:52:24 PM (6 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
802c412
Parents:
f4ae8c4
Message:

lint reduction

File:
1 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/jitter.py

    r199bd07 rd86f0fc  
    77""" 
    88from __future__ import division, print_function 
    9  
    10 import sys, os 
    11 sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.realpath(__file__)))) 
    129 
    1310import argparse 
    … …  
    5047def draw_ellipsoid(ax, size, view, jitter, steps=25, alpha=1): 
    5148    """Draw an ellipsoid.""" 
    52     a,b,c = size 
     49    a, b, c = size 
    5350    u = np.linspace(0, 2 * np.pi, steps) 
    5451    v = np.linspace(0, np.pi, steps) 
    … …  
    6158 
    6259    draw_labels(ax, view, jitter, [ 
    63          ('c+', [ 0, 0, c], [ 1, 0, 0]), 
    64          ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
    65          ('a+', [ a, 0, 0], [ 0, 0, 1]), 
    66          ('a-', [-a, 0, 0], [ 0, 0,-1]), 
    67          ('b+', [ 0, b, 0], [-1, 0, 0]), 
    68          ('b-', [ 0,-b, 0], [-1, 0, 0]), 
     60        ('c+', [+0, +0, +c], [+1, +0, +0]), 
     61        ('c-', [+0, +0, -c], [+0, +0, -1]), 
     62        ('a+', [+a, +0, +0], [+0, +0, +1]), 
     63        ('a-', [-a, +0, +0], [+0, +0, -1]), 
     64        ('b+', [+0, +b, +0], [-1, +0, +0]), 
     65        ('b-', [+0, -b, +0], [-1, +0, +0]), 
    6966    ]) 
    7067 
    7168def draw_sc(ax, size, view, jitter, steps=None, alpha=1): 
     69    """Draw points for simple cubic paracrystal""" 
    7270    atoms = _build_sc() 
    7371    _draw_crystal(ax, size, view, jitter, atoms=atoms) 
    7472 
    7573def draw_fcc(ax, size, view, jitter, steps=None, alpha=1): 
     74    """Draw points for face-centered cubic paracrystal""" 
    7675    # Build the simple cubic crystal 
    7776    atoms = _build_sc() 
    … …  
    8887 
    8988def draw_bcc(ax, size, view, jitter, steps=None, alpha=1): 
     89    """Draw points for body-centered cubic paracrystal""" 
    9090    # Build the simple cubic crystal 
    9191    atoms = _build_sc() 
    … …  
    127127    """Draw a parallelepiped.""" 
    128128    a, b, c = size 
    129     x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1]) 
    130     y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1]) 
    131     z = c*np.array([ 1, 1, 1, 1,-1,-1,-1,-1]) 
     129    x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1]) 
     130    y = b*np.array([+1, +1, -1, -1, +1, +1, -1, -1]) 
     131    z = c*np.array([+1, +1, +1, +1, -1, -1, -1, -1]) 
    132132    tri = np.array([ 
    133133        # counter clockwise triangles 
    134134        # z: up/down, x: right/left, y: front/back 
    135         [0,1,2], [3,2,1], # top face 
    136         [6,5,4], [5,6,7], # bottom face 
    137         [0,2,6], [6,4,0], # right face 
    138         [1,5,7], [7,3,1], # left face 
    139         [2,3,6], [7,6,3], # front face 
    140         [4,1,0], [5,1,4], # back face 
     135        [0, 1, 2], [3, 2, 1], # top face 
     136        [6, 5, 4], [5, 6, 7], # bottom face 
     137        [0, 2, 6], [6, 4, 0], # right face 
     138        [1, 5, 7], [7, 3, 1], # left face 
     139        [2, 3, 6], [7, 6, 3], # front face 
     140        [4, 1, 0], [5, 1, 4], # back face 
    141141    ]) 
    142142 
    … …  
    149149    # rotate that face. 
    150150    if 1: 
    151         x = a*np.array([ 1,-1, 1,-1, 1,-1, 1,-1]) 
    152         y = b*np.array([ 1, 1,-1,-1, 1, 1,-1,-1]) 
    153         z = c*np.array([ 1, 1, 1, 1,-1,-1,-1,-1]) 
     151        x = a*np.array([+1, -1, +1, -1, +1, -1, +1, -1]) 
     152        y = b*np.array([+1, +1, -1, -1, +1, +1, -1, -1]) 
     153        z = c*np.array([+1, +1, +1, +1, -1, -1, -1, -1]) 
    154154        x, y, z = transform_xyz(view, jitter, x, y, abs(z)+0.001) 
    155         ax.plot_trisurf(x, y, triangles=tri, Z=z, color=[1,0.6,0.6], alpha=alpha) 
     155        ax.plot_trisurf(x, y, triangles=tri, Z=z, color=[1, 0.6, 0.6], alpha=alpha) 
    156156 
    157157    draw_labels(ax, view, jitter, [ 
    158          ('c+', [ 0, 0, c], [ 1, 0, 0]), 
    159          ('c-', [ 0, 0,-c], [ 0, 0,-1]), 
    160          ('a+', [ a, 0, 0], [ 0, 0, 1]), 
    161          ('a-', [-a, 0, 0], [ 0, 0,-1]), 
    162          ('b+', [ 0, b, 0], [-1, 0, 0]), 
    163          ('b-', [ 0,-b, 0], [-1, 0, 0]), 
     158        ('c+', [+0, +0, +c], [+1, +0, +0]), 
     159        ('c-', [+0, +0, -c], [+0, +0, -1]), 
     160        ('a+', [+a, +0, +0], [+0, +0, +1]), 
     161        ('a-', [-a, +0, +0], [+0, +0, -1]), 
     162        ('b+', [+0, +b, +0], [-1, +0, +0]), 
     163        ('b-', [+0, -b, +0], [-1, +0, +0]), 
    164164    ]) 
    165165 
    … …  
    187187    cloud = [ 
    188188        [-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], 
     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], 
    215215    ] 
    216216    dtheta, dphi, dpsi = jitter 
    … …  
    228228    for v in 'xyz': 
    229229        a, b, c = size 
    230         lim = np.sqrt(a**2+b**2+c**2) 
     230        lim = np.sqrt(a**2 + b**2 + c**2) 
    231231        getattr(ax, 'set_'+v+'lim')([-lim, lim]) 
    232232        getattr(ax, v+'axis').label.set_text(v) 
    … …  
    297297        Should allow free movement in theta, but phi is distorted. 
    298298    """ 
    299     theta, phi, psi = view 
    300     dtheta, dphi, dpsi = jitter 
    301  
    302299    t = np.linspace(-1, 1, n) 
    303300    weights = np.ones_like(t) 
    … …  
    314311 
    315312    if projection == 'equirectangular':  #define PROJECTION 1 
    316         def rotate(theta_i, phi_j): 
     313        def _rotate(theta_i, phi_j): 
    317314            return Rx(phi_j)*Ry(theta_i) 
    318         def weight(theta_i, phi_j, wi, wj): 
     315        def _weight(theta_i, phi_j, wi, wj): 
    319316            return wi*wj*abs(cos(radians(theta_i))) 
    320317    elif projection == 'sinusoidal':  #define PROJECTION 2 
    321         def rotate(theta_i, phi_j): 
     318        def _rotate(theta_i, phi_j): 
    322319            latitude = theta_i 
    323320            scale = cos(radians(latitude)) 
    … …  
    325322            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
    326323            return Rx(longitude)*Ry(latitude) 
    327         def weight(theta_i, phi_j, wi, wj): 
     324        def _weight(theta_i, phi_j, wi, wj): 
    328325            latitude = theta_i 
    329326            scale = cos(radians(latitude)) 
    … …  
    331328            return w*wi*wj 
    332329    elif projection == 'guyou':  #define PROJECTION 3  (eventually?) 
    333         def rotate(theta_i, phi_j): 
     330        def _rotate(theta_i, phi_j): 
    334331            from guyou import guyou_invert 
    335332            #latitude, longitude = guyou_invert([theta_i], [phi_j]) 
    336333            longitude, latitude = guyou_invert([phi_j], [theta_i]) 
    337334            return Rx(longitude[0])*Ry(latitude[0]) 
    338         def weight(theta_i, phi_j, wi, wj): 
     335        def _weight(theta_i, phi_j, wi, wj): 
    339336            return wi*wj 
    340337    elif projection == 'azimuthal_equidistance':  # Note: Rz Ry, not Rx Ry 
    341         def rotate(theta_i, phi_j): 
     338        def _rotate(theta_i, phi_j): 
    342339            latitude = sqrt(theta_i**2 + phi_j**2) 
    343340            longitude = degrees(np.arctan2(phi_j, theta_i)) 
    344341            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
    345342            return Rz(longitude)*Ry(latitude) 
    346         def weight(theta_i, phi_j, wi, wj): 
     343        def _weight(theta_i, phi_j, wi, wj): 
    347344            # Weighting for each point comes from the integral: 
    348345            #     \int\int I(q, lat, log) sin(lat) dlat dlog 
    … …  
    377374            return w*wi*wj if latitude < 180 else 0 
    378375    elif projection == 'azimuthal_equal_area': 
    379         def rotate(theta_i, phi_j): 
     376        def _rotate(theta_i, phi_j): 
    380377            R = min(1, sqrt(theta_i**2 + phi_j**2)/180) 
    381378            latitude = 180-degrees(2*np.arccos(R)) 
    … …  
    383380            #print("(%+7.2f, %+7.2f) => (%+7.2f, %+7.2f)"%(theta_i, phi_j, latitude, longitude)) 
    384381            return Rz(longitude)*Ry(latitude) 
    385         def weight(theta_i, phi_j, wi, wj): 
     382        def _weight(theta_i, phi_j, wi, wj): 
    386383            latitude = sqrt(theta_i**2 + phi_j**2) 
    387384            w = sin(radians(latitude))/latitude if latitude != 0 else 1 
    … …  
    391388 
    392389    # mesh in theta, phi formed by rotating z 
     390    dtheta, dphi, dpsi = jitter 
    393391    z = np.matrix([[0], [0], [radius]]) 
    394     points = np.hstack([rotate(theta_i, phi_j)*z 
     392    points = np.hstack([_rotate(theta_i, phi_j)*z 
    395393                        for theta_i in dtheta*t 
    396394                        for phi_j in dphi*t]) 
    397395    # select just the active points (i.e., those with phi < 180 
    398     w = np.array([weight(theta_i, phi_j, wi, wj) 
     396    w = np.array([_weight(theta_i, phi_j, wi, wj) 
    399397                  for wi, theta_i in zip(weights, dtheta*t) 
    400398                  for wj, phi_j in zip(weights, dphi*t)]) 
    401399    #print(max(w), min(w), min(w[w>0])) 
    402     points = points[:, w>0] 
    403     w = w[w>0] 
     400    points = points[:, w > 0] 
     401    w = w[w > 0] 
    404402    w /= max(w) 
    405403 
    … …  
    469467    x, y, z = [np.asarray(v) for v in (x, y, z)] 
    470468    shape = x.shape 
    471     points = np.matrix([x.flatten(),y.flatten(),z.flatten()]) 
     469    points = np.matrix([x.flatten(), y.flatten(), z.flatten()]) 
    472470    points = apply_jitter(jitter, points) 
    473471    points = orient_relative_to_beam(view, points) 
    … …  
    543541    theta, phi, psi = view 
    544542    theta_pd, phi_pd, psi_pd = [scale*v for v in jitter] 
    545     theta_pd_n, phi_pd_n, psi_pd_n = PD_N_TABLE[(theta_pd>0, phi_pd>0, psi_pd>0)] 
     543    theta_pd_n, phi_pd_n, psi_pd_n = PD_N_TABLE[(theta_pd > 0, phi_pd > 0, psi_pd > 0)] 
    546544    ## increase pd_n for testing jitter integration rather than simple viz 
    547545    #theta_pd_n, phi_pd_n, psi_pd_n = [5*v for v in (theta_pd_n, phi_pd_n, psi_pd_n)] 
    … …  
    571569    if 0: 
    572570        level = np.asarray(255*(Iqxy - vmin)/(vmax - vmin), 'i') 
    573         level[level<0] = 0 
     571        level[level < 0] = 0 
    574572        colors = plt.get_cmap()(level) 
    575573        ax.plot_surface(qx, qy, -1.1, rstride=1, cstride=1, facecolors=colors) 
    … …  
    618616    return calculator 
    619617 
    620 def select_calculator(model_name, n=150, size=(10,40,100)): 
     618def select_calculator(model_name, n=150, size=(10, 40, 100)): 
    621619    """ 
    622620    Create a model calculator for the given shape. 
    … …  
    641639        radius = 0.5*c 
    642640        calculator = build_model('sc_paracrystal', n=n, dnn=dnn, 
    643                                   d_factor=d_factor, radius=(1-d_factor)*radius, 
    644                                   background=0) 
     641                                 d_factor=d_factor, radius=(1-d_factor)*radius, 
     642                                 background=0) 
    645643    elif model_name == 'fcc_paracrystal': 
    646644        a = b = c 
    … …  
    650648        radius = sqrt(2)/4 * c 
    651649        calculator = build_model('fcc_paracrystal', n=n, dnn=dnn, 
    652                                   d_factor=d_factor, radius=(1-d_factor)*radius, 
    653                                   background=0) 
     650                                 d_factor=d_factor, radius=(1-d_factor)*radius, 
     651                                 background=0) 
    654652    elif model_name == 'bcc_paracrystal': 
    655653        a = b = c 
    … …  
    659657        radius = sqrt(3)/2 * c 
    660658        calculator = build_model('bcc_paracrystal', n=n, dnn=dnn, 
    661                                   d_factor=d_factor, radius=(1-d_factor)*radius, 
    662                                   background=0) 
     659                                 d_factor=d_factor, radius=(1-d_factor)*radius, 
     660                                 background=0) 
    663661    elif model_name == 'cylinder': 
    664662        calculator = build_model('cylinder', n=n, qmax=0.3, radius=b, length=c) 
    … …  
    685683    'cylinder', 
    686684    'fcc_paracrystal', 'bcc_paracrystal', 'sc_paracrystal', 
    687  ] 
     685] 
    688686 
    689687DRAW_SHAPES = { 
    … …  
    761759 
    762760    ## add control widgets to plot 
    763     axtheta  = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor) 
     761    axtheta = plt.axes([0.1, 0.15, 0.45, 0.04], axisbg=axcolor) 
    764762    axphi = plt.axes([0.1, 0.1, 0.45, 0.04], axisbg=axcolor) 
    765763    axpsi = plt.axes([0.1, 0.05, 0.45, 0.04], axisbg=axcolor) 
    … …  
    768766    spsi = Slider(axpsi, 'Psi', -180, 180, valinit=psi) 
    769767 
    770     axdtheta  = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor) 
     768    axdtheta = plt.axes([0.75, 0.15, 0.15, 0.04], axisbg=axcolor) 
    771769    axdphi = plt.axes([0.75, 0.1, 0.15, 0.04], axisbg=axcolor) 
    772     axdpsi= plt.axes([0.75, 0.05, 0.15, 0.04], axisbg=axcolor) 
     770    axdpsi = plt.axes([0.75, 0.05, 0.15, 0.04], axisbg=axcolor) 
    773771    # Note: using ridiculous definition of rectangle distribution, whose width 
    774772    # in sasmodels is sqrt(3) times the given width.  Divide by sqrt(3) to keep 
    … …  
    797795 
    798796    ## bind control widgets to view updater 
    799     stheta.on_changed(lambda v: update(v,'theta')) 
     797    stheta.on_changed(lambda v: update(v, 'theta')) 
    800798    sphi.on_changed(lambda v: update(v, 'phi')) 
    801799    spsi.on_changed(lambda v: update(v, 'psi')) 
    … …  
    815813        formatter_class=argparse.ArgumentDefaultsHelpFormatter, 
    816814        ) 
    817     parser.add_argument('-p', '--projection', choices=PROJECTIONS, default=PROJECTIONS[0], help='coordinate projection') 
    818     parser.add_argument('-s', '--size', type=str, default='10,40,100', help='a,b,c lengths') 
    819     parser.add_argument('-d', '--distribution', choices=DISTRIBUTIONS, default=DISTRIBUTIONS[0], help='jitter distribution') 
    820     parser.add_argument('-m', '--mesh', type=int, default=30, help='#points in theta-phi mesh') 
    821     parser.add_argument('shape', choices=SHAPES, nargs='?', default=SHAPES[0], help='oriented shape') 
     815    parser.add_argument('-p', '--projection', choices=PROJECTIONS, 
     816                        default=PROJECTIONS[0], 
     817                        help='coordinate projection') 
     818    parser.add_argument('-s', '--size', type=str, default='10,40,100', 
     819                        help='a,b,c lengths') 
     820    parser.add_argument('-d', '--distribution', choices=DISTRIBUTIONS, 
     821                        default=DISTRIBUTIONS[0], 
     822                        help='jitter distribution') 
     823    parser.add_argument('-m', '--mesh', type=int, default=30, 
     824                        help='#points in theta-phi mesh') 
     825    parser.add_argument('shape', choices=SHAPES, nargs='?', default=SHAPES[0], 
     826                        help='oriented shape') 
    822827    opts = parser.parse_args() 
    823828    size = tuple(int(v) for v in opts.size.split(',')) 
Note: See TracChangeset for help on using the changeset viewer.