Changeset d86f0fc in sasmodels


Ignore:
Timestamp:
Mar 26, 2018 1: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

Location:
sasmodels
Files:
12 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/compare.py

    r3221de0 rd86f0fc  
    718718    model = core.build_model(model_info, dtype=dtype, platform="ocl") 
    719719    calculator = DirectModel(data, model, cutoff=cutoff) 
    720     engine_type = calculator._model.__class__.__name__.replace('Model','').upper() 
     720    engine_type = calculator._model.__class__.__name__.replace('Model', '').upper() 
    721721    bits = calculator._model.dtype.itemsize*8 
    722722    precision = "fast" if getattr(calculator._model, 'fast', False) else str(bits) 
    … …  
    14911491            vmin, vmax = limits 
    14921492            self.limits = vmax*1e-7, 1.3*vmax 
    1493             import pylab; pylab.clf() 
     1493            import pylab 
     1494            pylab.clf() 
    14941495            plot_models(self.opts, result, limits=self.limits) 
    14951496 

  • sasmodels/data.py

    rf549e37 rd86f0fc  
    706706    else: 
    707707        vmin_scaled, vmax_scaled = vmin, vmax 
    708     nx, ny = len(data.x_bins), len(data.y_bins) 
     708    #nx, ny = len(data.x_bins), len(data.y_bins) 
    709709    x_bins, y_bins, image = _build_matrix(data, plottable) 
    710710    plt.imshow(image, 
    … …  
    772772        if loop >= max_loop:  # this protects never-ending loop 
    773773            break 
    774         image = _fillup_pixels(self, image=image, weights=weights) 
     774        image = _fillup_pixels(image=image, weights=weights) 
    775775        loop += 1 
    776776 
    … …  
    817817    return x_bins, y_bins 
    818818 
    819 def _fillup_pixels(self, image=None, weights=None): 
     819def _fillup_pixels(image=None, weights=None): 
    820820    """ 
    821821    Fill z values of the empty cells of 2d image matrix 

  • sasmodels/generate.py

    ra0243f7 rd86f0fc  
    290290    import loops. 
    291291    """ 
    292     if (info.source and any(lib.startswith('lib/gauss') for lib in info.source)): 
    293         import os.path 
     292    if info.source and any(lib.startswith('lib/gauss') for lib in info.source): 
    294293        from .gengauss import gengauss 
    295         path = os.path.join(MODEL_PATH, "lib", "gauss%d.c"%n) 
    296         if not os.path.exists(path): 
     294        path = joinpath(MODEL_PATH, "lib", "gauss%d.c"%n) 
     295        if not exists(path): 
    297296            gengauss(n, path) 
    298297        info.source = ["lib/gauss%d.c"%n if lib.startswith('lib/gauss') 
    299                         else lib for lib in info.source] 
     298                       else lib for lib in info.source] 
    300299 
    301300def format_units(units): 
    … …  
    321320                     for w, h in zip(column_widths, PARTABLE_HEADERS)] 
    322321 
    323     sep = " ".join("="*w for w in column_widths) 
     322    underbar = " ".join("="*w for w in column_widths) 
    324323    lines = [ 
    325         sep, 
     324        underbar, 
    326325        " ".join("%-*s" % (w, h) 
    327326                 for w, h in zip(column_widths, PARTABLE_HEADERS)), 
    328         sep, 
     327        underbar, 
    329328        ] 
    330329    for p in pars: 
    … …  
    335334            "%*g" % (column_widths[3], p.default), 
    336335            ])) 
    337     lines.append(sep) 
     336    lines.append(underbar) 
    338337    return "\n".join(lines) 
    339338 
    … …  
    613612    """ 
    614613    spaces = " "*depth 
    615     sep = "\n" + spaces 
    616     return spaces + sep.join(s.split("\n")) 
     614    interline_separator = "\n" + spaces 
     615    return spaces + interline_separator.join(s.split("\n")) 
    617616 
    618617 
    … …  
    620619def load_template(filename): 
    621620    # type: (str) -> str 
     621    """ 
     622    Load template file from sasmodels resource directory. 
     623    """ 
    622624    path = joinpath(DATA_PATH, filename) 
    623625    mtime = getmtime(path) 

  • sasmodels/guyou.py

    r0d5a655 rd86f0fc  
    3131# 
    3232# 2017-11-01 Paul Kienzle 
    33 # * converted to python, using degrees rather than radians 
     33# * converted to python, with degrees rather than radians 
     34""" 
     35Convert between latitude-longitude and Guyou map coordinates. 
     36""" 
     37 
    3438from __future__ import division, print_function 
    3539 
    3640import numpy as np 
    37 from numpy import sqrt, pi, tan, cos, sin, log, exp, arctan2 as atan2, sign, radians, degrees 
     41from numpy import sqrt, pi, tan, cos, sin, sign, radians, degrees 
     42from numpy import sinh, arctan as atan 
     43 
     44# scipy version of special functions 
     45from scipy.special import ellipj as ellipticJ, ellipkinc as ellipticF 
    3846 
    3947_ = """ 
    … …  
    5967""" 
    6068 
    61 # scipy version of special functions 
    62 from scipy.special import ellipj as ellipticJ, ellipkinc as ellipticF 
    63 from numpy import sinh, sign, arctan as atan 
    64  
    6569def ellipticJi(u, v, m): 
    6670    scalar = np.isscalar(u) and np.isscalar(v) and np.isscalar(m) 
    6771    u, v, m = np.broadcast_arrays(u, v, m) 
    68     result = np.empty_like([u,u,u], 'D') 
    69     real = v==0 
    70     imag = u==0 
     72    result = np.empty_like([u, u, u], 'D') 
     73    real = (v == 0) 
     74    imag = (u == 0) 
    7175    mixed = ~(real|imag) 
    7276    result[:, real] = _ellipticJi_real(u[real], m[real]) 
    7377    result[:, imag] = _ellipticJi_imag(v[imag], m[imag]) 
    7478    result[:, mixed] = _ellipticJi(u[mixed], v[mixed], m[mixed]) 
    75     return result[0,:] if scalar else result 
     79    return result[0, :] if scalar else result 
    7680 
    7781def _ellipticJi_real(u, m): 
    … …  
    104108        result = np.empty_like(phi, 'D') 
    105109        index = (phi == 0) 
    106         result[index] = ellipticF(atan(sinh(abs(phi[index]))), 1-m[index]) * sign(psi[index]) 
     110        result[index] = ellipticF(atan(sinh(abs(phi[index]))), 
     111                                  1-m[index]) * sign(psi[index]) 
    107112        result[~index] = ellipticFi(phi[~index], psi[~index], m[~index]) 
    108113        return result.reshape(1)[0] if scalar else result 
    … …  
    117122    cotlambda2 = (-b + sqrt(b * b - 4 * c)) / 2 
    118123    re = ellipticF(atan(1 / sqrt(cotlambda2)), m) * sign(phi) 
    119     im = ellipticF(atan(sqrt(np.maximum(0,(cotlambda2 / cotphi2 - 1) / m))), 1 - m) * sign(psi) 
     124    im = ellipticF(atan(sqrt(np.maximum(0, (cotlambda2 / cotphi2 - 1) / m))), 
     125                   1 - m) * sign(psi) 
    120126    return re + 1j*im 
    121127 
    122 sqrt2 = sqrt(2) 
     128SQRT2 = sqrt(2) 
    123129 
    124130# [PAK] renamed k_ => cos_u, k => sin_u, k*k => sinsq_u to avoid k,K confusion 
    … …  
    127133# K = 3.165103454447431823666270142140819753058976299237578486994... 
    128134def guyou(lam, phi): 
     135    """Transform from (latitude, longitude) to point (x, y)""" 
    129136    # [PAK] wrap into [-pi/2, pi/2] radians 
    130137    x, y = np.asarray(lam), np.asarray(phi) 
    … …  
    135142 
    136143    # Compute constant K 
    137     cos_u = (sqrt2 - 1) / (sqrt2 + 1) 
     144    cos_u = (SQRT2 - 1) / (SQRT2 + 1) 
    138145    sinsq_u = 1 - cos_u**2 
    139146    K = ellipticF(pi/2, sinsq_u) 
    … …  
    144151    at = atan(r * (cos(lam) - 1j*sin(lam))) 
    145152    t = ellipticFi(at.real, at.imag, sinsq_u) 
    146     x, y = (-t.imag, sign(phi + (phi==0))*(0.5 * K - t.real)) 
     153    x, y = (-t.imag, sign(phi + (phi == 0))*(0.5 * K - t.real)) 
    147154 
    148155    # [PAK] convert to degrees, and return to original tile 
    … …  
    150157 
    151158def guyou_invert(x, y): 
     159    """Transform from point (x, y) on plot to (latitude, longitude)""" 
    152160    # [PAK] wrap into [-pi/2, pi/2] radians 
    153161    x, y = np.asarray(x), np.asarray(y) 
    … …  
    158166 
    159167    # compute constant K 
    160     cos_u = (sqrt2 - 1) / (sqrt2 + 1) 
     168    cos_u = (SQRT2 - 1) / (SQRT2 + 1) 
    161169    sinsq_u = 1 - cos_u**2 
    162170    K = ellipticF(pi/2, sinsq_u) 
    … …  
    174182 
    175183def plot_grid(): 
     184    """Plot the latitude-longitude grid for Guyou transform""" 
    176185    import matplotlib.pyplot as plt 
    177186    from numpy import linspace 
    … …  
    203212    plt.ylabel('latitude') 
    204213 
     214def main(): 
     215    plot_grid() 
     216    import matplotlib.pyplot as plt 
     217    plt.show() 
     218 
    205219if __name__ == "__main__": 
    206     plot_grid() 
    207     import matplotlib.pyplot as plt; plt.show() 
    208  
     220    main() 
    209221 
    210222_ = """ 

  • 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(',')) 

  • sasmodels/kernel_iq.c

    raadec17 rd86f0fc  
    173173// Apply the rotation matrix returned from qac_rotation to the point (qx,qy), 
    174174// returning R*[qx,qy]' = [qa,qc]' 
    175 static double 
     175static void 
    176176qac_apply( 
    177177    QACRotation *rotation, 
    … …  
    246246// Apply the rotation matrix returned from qabc_rotation to the point (qx,qy), 
    247247// returning R*[qx,qy]' = [qa,qb,qc]' 
    248 static double 
     248static void 
    249249qabc_apply( 
    250250    QABCRotation *rotation, 

  • sasmodels/kernelcl.py

    r6cbdcd4 rd86f0fc  
    151151        if not HAVE_OPENCL: 
    152152            raise RuntimeError("OpenCL startup failed with ***" 
    153                             + OPENCL_ERROR + "***; using C compiler instead") 
     153                               + OPENCL_ERROR + "***; using C compiler instead") 
    154154        reset_environment() 
    155155        if ENV is None: 

  • sasmodels/models/core_shell_cylinder.c

    r108e70e rd86f0fc  
    4848 
    4949 
    50 double Iqac(double qab, double qc, 
     50static double 
     51Iqac(double qab, double qc, 
    5152    double core_sld, 
    5253    double shell_sld, 

  • sasmodels/models/hollow_rectangular_prism.c

    r108e70e rd86f0fc  
    1 double form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness); 
    2 double Iq(double q, double sld, double solvent_sld, double length_a, 
    3           double b2a_ratio, double c2a_ratio, double thickness); 
    4  
    5 double form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
     1static double 
     2form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) 
    63{ 
    74    double length_b = length_a * b2a_ratio; 
    … …  
    1613} 
    1714 
    18 double Iq(double q, 
     15static double 
     16Iq(double q, 
    1917    double sld, 
    2018    double solvent_sld, 
    … …  
    8583} 
    8684 
    87 double Iqabc(double qa, double qb, double qc, 
     85static double 
     86Iqabc(double qa, double qb, double qc, 
    8887    double sld, 
    8988    double solvent_sld, 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.c

    r74768cb rd86f0fc  
    1 double form_volume(double length_a, double b2a_ratio, double c2a_ratio); 
    2 double Iq(double q, double sld, double solvent_sld, double length_a, 
    3           double b2a_ratio, double c2a_ratio); 
    4  
    5 double form_volume(double length_a, double b2a_ratio, double c2a_ratio) 
     1static double 
     2form_volume(double length_a, double b2a_ratio, double c2a_ratio) 
    63{ 
    74    double length_b = length_a * b2a_ratio; 
    … …  
    118} 
    129 
    13 double Iq(double q, 
     10static double 
     11Iq(double q, 
    1412    double sld, 
    1513    double solvent_sld, 

  • sasmodels/models/rectangular_prism.c

    r108e70e rd86f0fc  
    1 double form_volume(double length_a, double b2a_ratio, double c2a_ratio); 
    2 double Iq(double q, double sld, double solvent_sld, double length_a, 
    3           double b2a_ratio, double c2a_ratio); 
    4  
    5 double form_volume(double length_a, double b2a_ratio, double c2a_ratio) 
     1static double 
     2form_volume(double length_a, double b2a_ratio, double c2a_ratio) 
    63{ 
    74    return length_a * (length_a*b2a_ratio) * (length_a*c2a_ratio); 
    85} 
    96 
    10 double Iq(double q, 
     7static double 
     8Iq(double q, 
    119    double sld, 
    1210    double solvent_sld, 
    … …  
    7169 
    7270 
    73 double Iqabc(double qa, double qb, double qc, 
     71static double 
     72Iqabc(double qa, double qb, double qc, 
    7473    double sld, 
    7574    double solvent_sld, 

  • sasmodels/multiscat.py

    r49d1f8b8 rd86f0fc  
    261261 
    262262def scattering_coeffs(p, coverage=0.99): 
     263    r""" 
     264    Return the coefficients of the scattering powers for transmission 
     265    probability *p*.  This is just the corresponding values for the 
     266    Poisson distribution for $\lambda = -\ln(1-p)$ such that 
     267    $\sum_{k = 0 \ldots n} P(k; \lambda)$ is larger than *coverage*. 
     268    """ 
    263269    L = -np.log(1-p) 
    264270    num_scatter = truncated_poisson_invcdf(coverage, L) 
    … …  
    266272    return coeffs 
    267273 
    268 def truncated_poisson_invcdf(p, L): 
     274def truncated_poisson_invcdf(coverage, L): 
    269275    r""" 
    270     Return smallest k such that cdf(k; L) > p from the truncated Poisson 
     276    Return smallest k such that cdf(k; L) > coverage from the truncated Poisson 
    271277    probability excluding k=0 
    272278    """ 
    … …  
    275281    pmf = -np.exp(-L) / np.expm1(-L) 
    276282    k = 0 
    277     while cdf < p: 
     283    while cdf < coverage: 
    278284        k += 1 
    279285        pmf *= L/k 
    … …  
    305311 
    306312class MultipleScattering(Resolution): 
     313    r""" 
     314    Compute multiple scattering using Fourier convolution. 
     315 
     316    The fourier steps are determined by *qmax*, the maximum $q$ value 
     317    desired, *nq* the number of $q$ steps and *window*, the amount 
     318    of padding around the circular convolution.  The $q$ spacing 
     319    will be $\Delta q = 2 q_\mathrm{max} w / n_q$.  If *nq* is not 
     320    given it will use $n_q = 2^k$ such that $\Delta q < q_\mathrm{min}$. 
     321 
     322    *probability* is related to the expected number of scattering 
     323    events in the sample $\lambda$ as $p = 1 = e^{-\lambda}$.  As a 
     324    hack to allow probability to be a fitted parameter, the "value" 
     325    can be a function that takes no parameters and returns the current 
     326    value of the probability.  *coverage* determines how many scattering 
     327    steps to consider.  The default is 0.99, which sets $n$ such that 
     328    $1 \ldots n$ covers 99% of the Poisson probability mass function. 
     329 
     330    *is2d* is True then 2D scattering is used, otherwise it accepts 
     331    and returns 1D scattering. 
     332 
     333    *resolution* is the resolution function to apply after multiple 
     334    scattering.  If present, then the resolution $q$ vectors will provide 
     335    default values for *qmin*, *qmax* and *nq*. 
     336    """ 
    307337    def __init__(self, qmin=None, qmax=None, nq=None, window=2, 
    308338                 probability=None, coverage=0.99, 
    309339                 is2d=False, resolution=None, 
    310340                 dtype=PRECISION): 
    311         r""" 
    312         Compute multiple scattering using Fourier convolution. 
    313  
    314         The fourier steps are determined by *qmax*, the maximum $q$ value 
    315         desired, *nq* the number of $q$ steps and *window*, the amount 
    316         of padding around the circular convolution.  The $q$ spacing 
    317         will be $\Delta q = 2 q_\mathrm{max} w / n_q$.  If *nq* is not 
    318         given it will use $n_q = 2^k$ such that $\Delta q < q_\mathrm{min}$. 
    319  
    320         *probability* is related to the expected number of scattering 
    321         events in the sample $\lambda$ as $p = 1 = e^{-\lambda}$.  As a 
    322         hack to allow probability to be a fitted parameter, the "value" 
    323         can be a function that takes no parameters and returns the current 
    324         value of the probability.  *coverage* determines how many scattering 
    325         steps to consider.  The default is 0.99, which sets $n$ such that 
    326         $1 \ldots n$ covers 99% of the Poisson probability mass function. 
    327  
    328         *is2d* is True then 2D scattering is used, otherwise it accepts 
    329         and returns 1D scattering. 
    330  
    331         *resolution* is the resolution function to apply after multiple 
    332         scattering.  If present, then the resolution $q$ vectors will provide 
    333         default values for *qmin*, *qmax* and *nq*. 
    334         """ 
    335341        # Infer qmin, qmax from instrument resolution calculator, if present 
    336342        if resolution is not None: 
    … …  
    424430        # Prepare the multiple scattering calculator (either numpy or OpenCL) 
    425431        self.transform = Calculator((2*nq, 2*nq), dtype=dtype) 
     432 
     433        # Iq and Iqxy will be set during apply 
     434        self.Iq = None # type: np.ndarray 
     435        self.Iqxy = None # type: np.ndarray 
    426436 
    427437    def apply(self, theory): 
    … …  
    471481 
    472482    def radial_profile(self, Iqxy): 
     483        """ 
     484        Compute that radial profile for the given Iqxy grid.  The grid should 
     485        be defined as for 
     486        """ 
    473487        # circular average, no anti-aliasing 
    474488        Iq = np.histogram(self._radius, bins=self._edges, weights=Iqxy)[0]/self._norm 
    … …  
    478492def annular_average(qxy, Iqxy, qbins): 
    479493    """ 
    480     Compute annular average of points at 
     494    Compute annular average of points in *Iqxy* at *qbins*.  The $q_x$, $q_y$ 
     495    coordinates for *Iqxy* are given in *qxy*. 
    481496    """ 
    482497    qxy, Iqxy = qxy.flatten(), Iqxy.flatten() 
Note: See TracChangeset for help on using the changeset viewer.