Changeset b3703f5 in sasmodels


Ignore:
Timestamp:
Mar 26, 2018 3:50:35 PM (7 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:
c11d09f, c462169
Parents:
802c412
Message:

lint reduction

Location:
sasmodels
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/guyou.py

    rd86f0fc rb3703f5  
    195195        plt.plot(x, y, 'g') 
    196196 
    197     for long in range(-limit, limit+1, step): 
    198         x, y = guyou(scale*long, scale*long_line) 
     197    for longitude in range(-limit, limit+1, step): 
     198        x, y = guyou(scale*longitude, scale*long_line) 
    199199        plt.plot(x, y, 'b') 
    200200    #plt.xlabel('longitude') 
    201201    plt.ylabel('latitude') 
     202    plt.title('forward transform') 
    202203 
    203204    plt.subplot(212) 
    … …  
    211212    plt.xlabel('longitude') 
    212213    plt.ylabel('latitude') 
     214    plt.title('inverse transform') 
    213215 
    214216def main(): 
     217    """Show the Guyou transformation""" 
    215218    plot_grid() 
    216219    import matplotlib.pyplot as plt 

  • sasmodels/jitter.py

    r802c412 rb3703f5  
    9999    _draw_crystal(axes, size, view, jitter, atoms=atoms) 
    100100 
    101 def _draw_crystal(axes, size, view, jitter, steps=None, alpha=1, atoms=None): 
     101def _draw_crystal(axes, size, view, jitter, atoms=None): 
    102102    atoms, size = np.asarray(atoms, 'd').T, np.asarray(size, 'd') 
    103103    x, y, z = atoms*size[:, None] 
    … …  
    298298    # TODO: try Kent distribution instead of a gaussian warped by projection 
    299299 
    300     t = np.linspace(-1, 1, n) 
    301     weights = np.ones_like(t) 
     300    dist_x = np.linspace(-1, 1, n) 
     301    weights = np.ones_like(dist_x) 
    302302    if dist == 'gaussian': 
    303         t *= 3 
    304         weights = exp(-0.5*t**2) 
     303        dist_x *= 3 
     304        weights = exp(-0.5*dist_x**2) 
    305305    elif dist == 'rectangle': 
    306306        # Note: uses sasmodels ridiculous definition of rectangle width 
    307         t *= sqrt(3) 
     307        dist_x *= sqrt(3) 
    308308    elif dist == 'uniform': 
    309309        pass 
    … …  
    392392    z = np.matrix([[0], [0], [radius]]) 
    393393    points = np.hstack([_rotate(theta_i, phi_j)*z 
    394                         for theta_i in dtheta*t 
    395                         for phi_j in dphi*t]) 
    396     w = np.array([_weight(theta_i, phi_j, w_i, w_j) 
    397                   for w_i, theta_i in zip(weights, dtheta*t) 
    398                   for w_j, phi_j in zip(weights, dphi*t)]) 
    399     #print(max(w), min(w), min(w[w>0])) 
    400     points = points[:, w > 0] 
    401     w = w[w > 0] 
    402     w /= max(w) 
     394                        for theta_i in dtheta*dist_x 
     395                        for phi_j in dphi*dist_x]) 
     396    dist_w = np.array([_weight(theta_i, phi_j, w_i, w_j) 
     397                       for w_i, theta_i in zip(weights, dtheta*dist_x) 
     398                       for w_j, phi_j in zip(weights, dphi*dist_x)]) 
     399    #print(max(dist_w), min(dist_w), min(dist_w[dist_w > 0])) 
     400    points = points[:, dist_w > 0] 
     401    dist_w = dist_w[dist_w > 0] 
     402    dist_w /= max(dist_w) 
    403403 
    404404    # rotate relative to beam 
    … …  
    407407    x, y, z = [np.array(v).flatten() for v in points] 
    408408    #plt.figure(2); plt.clf(); plt.hist(z, bins=np.linspace(-1, 1, 51)) 
    409     axes.scatter(x, y, z, c=w, marker='o', vmin=0., vmax=1.) 
     409    axes.scatter(x, y, z, c=dist_w, marker='o', vmin=0., vmax=1.) 
    410410 
    411411def draw_labels(axes, view, jitter, text): 
    … …  
    430430    """Construct a matrix to rotate points about *x* by *angle* degrees.""" 
    431431    angle = radians(angle) 
    432     R = [[1, 0, 0], 
    433          [0, +cos(angle), -sin(angle)], 
    434          [0, +sin(angle), +cos(angle)]] 
    435     return np.matrix(R) 
     432    rot = [[1, 0, 0], 
     433           [0, +cos(angle), -sin(angle)], 
     434           [0, +sin(angle), +cos(angle)]] 
     435    return np.matrix(rot) 
    436436 
    437437def Ry(angle): 
    438438    """Construct a matrix to rotate points about *y* by *angle* degrees.""" 
    439439    angle = radians(angle) 
    440     R = [[+cos(angle), 0, +sin(angle)], 
    441          [0, 1, 0], 
    442          [-sin(angle), 0, +cos(angle)]] 
    443     return np.matrix(R) 
     440    rot = [[+cos(angle), 0, +sin(angle)], 
     441           [0, 1, 0], 
     442           [-sin(angle), 0, +cos(angle)]] 
     443    return np.matrix(rot) 
    444444 
    445445def Rz(angle): 
    446446    """Construct a matrix to rotate points about *z* by *angle* degrees.""" 
    447447    angle = radians(angle) 
    448     R = [[+cos(angle), -sin(angle), 0], 
    449          [+sin(angle), +cos(angle), 0], 
    450          [0, 0, 1]] 
    451     return np.matrix(R) 
     448    rot = [[+cos(angle), -sin(angle), 0], 
     449           [+sin(angle), +cos(angle), 0], 
     450           [0, 0, 1]] 
     451    return np.matrix(rot) 
    452452 
    453453def transform_xyz(view, jitter, x, y, z): 
    … …  
    564564    elif 1: 
    565565        axes.contourf(qx/qx.max(), qy/qy.max(), Iqxy, zdir='z', offset=-1.1, 
    566                     levels=np.linspace(vmin, vmax, 24)) 
     566                      levels=np.linspace(vmin, vmax, 24)) 
    567567    else: 
    568568        axes.pcolormesh(qx, qy, Iqxy) 
    … …  
    583583    for details. 
    584584    """ 
    585     from sasmodels.core import load_model_info, build_model 
     585    from sasmodels.core import load_model_info, build_model as build_sasmodel 
    586586    from sasmodels.data import empty_data2D 
    587587    from sasmodels.direct_model import DirectModel 
    588588 
    589589    model_info = load_model_info(model_name) 
    590     model = build_model(model_info) #, dtype='double!') 
     590    model = build_sasmodel(model_info) #, dtype='double!') 
    591591    q = np.linspace(-qmax, qmax, n) 
    592592    data = empty_data2D(q, q) 

  • sasmodels/multiscat.py

    rd86f0fc rb3703f5  
    7373import argparse 
    7474import time 
    75 import os.path 
    7675 
    7776import numpy as np 
    … …  
    8180from sasmodels import core 
    8281from sasmodels import compare 
    83 from sasmodels import resolution2d 
    8482from sasmodels.resolution import Resolution, bin_edges 
    85 from sasmodels.data import empty_data1D, empty_data2D, plot_data 
    8683from sasmodels.direct_model import call_kernel 
    8784import sasmodels.kernelcl 
    … …  
    106103USE_FAST = True  # OpenCL faster, less accurate math 
    107104 
    108 class NumpyCalculator: 
     105class ICalculator: 
     106    """ 
     107    Multiple scattering calculator 
     108    """ 
     109    def fft(self, Iq): 
     110        """ 
     111        Compute the forward FFT for an image, real -> complex. 
     112        """ 
     113        raise NotImplementedError() 
     114 
     115    def ifft(self, Iq): 
     116        """ 
     117        Compute the inverse FFT for an image, complex -> complex. 
     118        """ 
     119        raise NotImplementedError() 
     120 
     121    def mulitple_scattering(self, Iq): 
     122        r""" 
     123        Compute multiple scattering for I(q) given scattering probability p. 
     124 
     125        Given a probability p of scattering with the thickness, the expected 
     126        number of scattering events, $\lambda$ is $-\log(1 - p)$, giving a 
     127        Poisson weighted sum of single, double, triple, etc. scattering patterns. 
     128        The number of patterns used is based on coverage (default 99%). 
     129        """ 
     130        raise NotImplementedError() 
     131 
     132class NumpyCalculator(ICalculator): 
     133    """ 
     134    Multiple scattering calculator using numpy fft. 
     135    """ 
    109136    def __init__(self, dims=None, dtype=PRECISION): 
    110137        self.dtype = dtype 
    111138        self.complex_dtype = np.dtype('F') if dtype == np.dtype('f') else np.dtype('D') 
    112         pass 
    113139 
    114140    def fft(self, Iq): 
    … …  
    127153 
    128154    def multiple_scattering(self, Iq, p, coverage=0.99): 
    129         r""" 
    130         Compute multiple scattering for I(q) given scattering probability p. 
    131  
    132         Given a probability p of scattering with the thickness, the expected 
    133         number of scattering events, $\lambda$ is $-\log(1 - p)$, giving a 
    134         Poisson weighted sum of single, double, triple, etc. scattering patterns. 
    135         The number of patterns used is based on coverage (default 99%). 
    136         """ 
    137155        #t0 = time.time() 
    138156        coeffs = scattering_coeffs(p, coverage) 
    … …  
    140158        scale = np.sum(Iq) 
    141159        frame = _forward_shift(Iq/scale, dtype=self.dtype) 
    142         F = np.fft.fft2(frame) 
    143         F_convolved = F * np.polyval(poly, F) 
    144         frame = np.fft.ifft2(F_convolved) 
     160        fourier_frame = np.fft.fft2(frame) 
     161        convolved = fourier_frame * np.polyval(poly, fourier_frame) 
     162        frame = np.fft.ifft2(convolved) 
    145163        result = scale * _inverse_shift(frame.real, dtype=self.dtype) 
    146164        #print("numpy multiscat time", time.time()-t0) 
    … …  
    173191""" 
    174192 
    175 class OpenclCalculator(NumpyCalculator): 
     193class OpenclCalculator(ICalculator): 
     194    """ 
     195    Multiple scattering calculator using OpenCL via pyfft. 
     196    """ 
    176197    polyval1f = None 
    177198    polyval1d = None 
    … …  
    180201        context = env.get_context(dtype) 
    181202        if dtype == np.dtype('f'): 
    182             if self.polyval1f is None: 
     203            if OpenclCalculator.polyval1f is None: 
    183204                program = sasmodels.kernelcl.compile_model( 
    184205                    context, POLYVAL1_KERNEL, dtype, fast=USE_FAST) 
    … …  
    187208            self.dtype = dtype 
    188209            self.complex_dtype = np.dtype('F') 
    189             self.polyval1 = self.polyval1f 
     210            self.polyval1 = OpenclCalculator.polyval1f 
    190211        else: 
    191             if self.polyval1d is None: 
     212            if OpenclCalculator.polyval1d is None: 
    192213                program = sasmodels.kernelcl.compile_model( 
    193214                    context, POLYVAL1_KERNEL, dtype, fast=False) 
    … …  
    196217            self.dtype = dtype 
    197218            self.complex_type = np.dtype('D') 
    198             self.polyval1 = self.polyval1d 
     219            self.polyval1 = OpenclCalculator.polyval1d 
    199220        self.queue = env.get_queue(dtype) 
    200221        self.plan = pyfft.cl.Plan(dims, queue=self.queue) 
    … …  
    229250        gpu_poly = cl_array.to_device(self.queue, poly) 
    230251        self.plan.execute(gpu_data.data) 
    231         degree, n = poly.shape[0], frame.shape[0]*frame.shape[1] 
     252        degree, data_size= poly.shape[0], frame.shape[0]*frame.shape[1] 
    232253        self.polyval1( 
    233             self.queue, [n], None, 
    234             np.int32(degree), gpu_poly.data, np.int32(n), gpu_data.data) 
     254            self.queue, [data_size], None, 
     255            np.int32(degree), gpu_poly.data, np.int32(data_size), gpu_data.data) 
    235256        self.plan.execute(gpu_data.data, inverse=True) 
    236257        frame = gpu_data.get() 
    … …  
    251272    """ 
    252273    if transform is None: 
    253         nx, ny = Iq.shape 
    254         transform = Calculator(dims=(nx*2, ny*2), dtype=dtype) 
     274        n_x, n_y = Iq.shape 
     275        transform = Calculator(dims=(n_x*2, n_y*2), dtype=dtype) 
    255276    scale = np.sum(Iq) 
    256277    frame = _forward_shift(Iq/scale, dtype=dtype) 
    … …  
    528549def parse_pars(model, opts): 
    529550    # type: (ModelInfo, argparse.Namespace) -> Dict[str, float] 
     551    """ 
     552    Parse par=val arguments from the command line. 
     553    """ 
    530554 
    531555    seed = np.random.randint(1000000) if opts.random and opts.seed < 0 else opts.seed 
    … …  
    541565        'is2d': opts.is2d, 
    542566    } 
    543     pars, pars2 = compare.parse_pars(compare_opts) 
     567    # Note: sascomp allows comparison on a pair of models, so ignore the second. 
     568    pars, _ = compare.parse_pars(compare_opts) 
    544569    return pars 
    545570 
    … …  
    550575        formatter_class=argparse.ArgumentDefaultsHelpFormatter, 
    551576        ) 
    552     parser.add_argument('-p', '--probability', type=float, default=0.1, help="scattering probability") 
    553     parser.add_argument('-n', '--nq', type=int, default=1024, help='number of mesh points') 
    554     parser.add_argument('-q', '--qmax', type=float, default=0.5, help='max q') 
    555     parser.add_argument('-w', '--window', type=float, default=2.0, help='q calc = q max * window') 
    556     parser.add_argument('-2', '--2d', dest='is2d', action='store_true', help='oriented sample') 
    557     parser.add_argument('-s', '--seed', default=-1, help='random pars with given seed') 
    558     parser.add_argument('-r', '--random', action='store_true', help='random pars with random seed') 
    559     parser.add_argument('-o', '--outfile', type=str, default="", help='random pars with random seed') 
    560     parser.add_argument('model', type=str, help='sas model name such as cylinder') 
    561     parser.add_argument('pars', type=str, nargs='*', help='model parameters such as radius=30') 
     577    parser.add_argument('-p', '--probability', type=float, default=0.1, 
     578                        help="scattering probability") 
     579    parser.add_argument('-n', '--nq', type=int, default=1024, 
     580                        help='number of mesh points') 
     581    parser.add_argument('-q', '--qmax', type=float, default=0.5, 
     582                        help='max q') 
     583    parser.add_argument('-w', '--window', type=float, default=2.0, 
     584                        help='q calc = q max * window') 
     585    parser.add_argument('-2', '--2d', dest='is2d', action='store_true', 
     586                        help='oriented sample') 
     587    parser.add_argument('-s', '--seed', default=-1, 
     588                        help='random pars with given seed') 
     589    parser.add_argument('-r', '--random', action='store_true', 
     590                        help='random pars with random seed') 
     591    parser.add_argument('-o', '--outfile', type=str, default="", 
     592                        help='random pars with random seed') 
     593    parser.add_argument('model', type=str, 
     594                        help='sas model name such as cylinder') 
     595    parser.add_argument('pars', type=str, nargs='*', 
     596                        help='model parameters such as radius=30') 
    562597    opts = parser.parse_args() 
    563598    assert opts.nq%2 == 0, "require even # points" 
    … …  
    607642            plotxy((res._q_steps, res._q_steps), res.Iqxy+background) 
    608643            pylab.title("total scattering for p=%g" % probability) 
     644            if res.resolution is not None: 
     645                pylab.figure() 
     646                plotxy((res._q_steps, res._q_steps), result) 
     647                pylab.title("total scattering with resolution") 
    609648    else: 
    610649        q = res._q 
    … …  
    624663            # Plot 1D pattern for partial scattering 
    625664            pylab.loglog(q, res.Iq+background, label="total for p=%g"%probability) 
     665            if res.resolution is not None: 
     666                pylab.loglog(q, result, label="total with dQ") 
    626667            #new_annulus = annular_average(res._radius, res.Iqxy, res._edges) 
    627668            #pylab.loglog(q, new_annulus+background, label="new total for p=%g"%probability) 
Note: See TracChangeset for help on using the changeset viewer.