Changeset b510b92 in sasmodels


Ignore:
Timestamp:
Mar 4, 2018 10:53:03 PM (4 years ago)
Author:
butler
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
367886f
Parents:
5bc373b (diff), bc248f8 (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-896

Files:
1 added
8 edited

Legend:

Unmodified
Added
Removed
  • example/fit.py

    rf4b36fa r1ee5ac6  
    179179 
    180180else: 
    181     print "No parameters for %s"%name 
     181    print("No parameters for %s"%name) 
    182182    sys.exit(1) 
    183183 
     
    192192else: 
    193193   problem = FitProblem(M) 
     194 
     195if __name__ == "__main__": 
     196   problem.plot() 
     197   import pylab; pylab.show() 
  • explore/multiscat.py

    re309e23 rbc248f8  
    6161The usual pinhole or slit resolution calculation can performed from these 
    6262calculated values. 
     63 
     64By default single precision OpenCL is used for the calculation.  Set the 
     65environment variable *SAS_OPENCL=none* to use double precision numpy FFT 
     66instead.  The OpenCL versions is about 10x faster on an elderly Mac with 
     67Intel HD 4000 graphics.  The single precision numerical artifacts don't 
     68seem to seriously impact overall accuracy, though they look pretty bad. 
    6369""" 
    6470 
     
    6773import argparse 
    6874import time 
     75import os.path 
    6976 
    7077import numpy as np 
     78from numpy import pi 
    7179from scipy.special import gamma 
    7280 
     
    7785from sasmodels.data import empty_data1D, empty_data2D, plot_data 
    7886from sasmodels.direct_model import call_kernel 
    79  
    80  
    81 class MultipleScattering(Resolution): 
    82     def __init__(self, qmax, nq, probability, is2d, window=2, power=0): 
    83         self.qmax = qmax 
    84         self.nq = nq 
    85         self.power = power 
    86         self.probability = probability 
    87         self.is2d = is2d 
    88         self.window = window 
    89  
    90         q_range = qmax * window 
    91         q = np.linspace(-q_range, q_range, nq) 
    92         qx, qy = np.meshgrid(q, q) 
    93  
    94         if is2d: 
    95             q_calc = [qx.flatten(), qy.flatten()] 
     87import sasmodels.kernelcl 
     88 
     89# TODO: select fast and accurate fft library 
     90# clFFT: https://github.com/clMathLibraries/clFFT (AMD's OpenCL) 
     91# - gpyfft: https://github.com/geggo/gpyfft (wraps clFFT) 
     92# - arrayfire: https://github.com/arrayfire (wraps clFFT and much more; +cuda) 
     93# pyFFT: https://github.com/fjarri-attic/pyfft (based on Apple's OpenCL; +cuda) 
     94# - Reikna: https://github.com/fjarri/reikna (evolved from pyfft) 
     95# genFFT: https://software.intel.com/en-us/articles/genFFT (Intel's OpenCL) 
     96# VexCL: https://vexcl.readthedocs.io/en/latest/ (c++ library) 
     97# TODO: switch to fftw when opencl is not available 
     98 
     99try: 
     100    import pyfft.cl 
     101    import pyopencl.array as cl_array 
     102    HAVE_OPENCL = sasmodels.kernelcl.use_opencl() 
     103except ImportError: 
     104    HAVE_OPENCL = False 
     105PRECISION = np.dtype('f' if HAVE_OPENCL else 'd')  # 'f' or 'd' 
     106USE_FAST = True  # OpenCL faster, less accurate math 
     107 
     108class NumpyCalculator: 
     109    def __init__(self, dims=None, dtype=PRECISION): 
     110        self.dtype = dtype 
     111        self.complex_dtype = np.dtype('F') if dtype == np.dtype('f') else np.dtype('D') 
     112        pass 
     113 
     114    def fft(self, Iq): 
     115        #t0 = time.time() 
     116        Iq = np.asarray(Iq, self.dtype) 
     117        result = np.fft.fft2(Iq) 
     118        #print("fft time", time.time()-t0) 
     119        return result 
     120 
     121    def ifft(self, fourier_frame): 
     122        #t0 = time.time() 
     123        fourier_frame = np.asarray(fourier_frame, self.complex_dtype) 
     124        result = np.fft.ifft2(fourier_frame) 
     125        #print("ifft time", time.time()-t0) 
     126        return result 
     127 
     128    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        """ 
     137        #t0 = time.time() 
     138        coeffs = scattering_coeffs(p, coverage) 
     139        poly = np.asarray(coeffs[::-1], dtype=self.dtype) 
     140        scale = np.sum(Iq) 
     141        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) 
     145        result = scale * _inverse_shift(frame.real, dtype=self.dtype) 
     146        #print("multiscat time", time.time()-t0) 
     147        return result 
     148 
     149# polyval1(c, x) computes (...((c0 x + c1) x + c2) x ... + cn) x 
     150# where c is an array of length *degree* and x is an array of 
     151# complex values (type double2) of length *n*. 2-D arrays can of 
     152# course be treated as 1-D arrays of length *nx* X *ny*. 
     153# When compiling with sasmodels.kernelcl.compile_model the double precision 
     154# types are converted to single precision as needed.  See the code in 
     155# sasmodels.generate._convert_type for details. 
     156POLYVAL1_KERNEL = """ 
     157kernel void polyval1( 
     158    const int degree, 
     159    global const double *coeff, 
     160    const int n, 
     161    global double2 *array) 
     162{ 
     163    int index = get_global_id(0); 
     164    if (index < n) { 
     165        const double2 x = array[index]; 
     166        double2 total = coeff[0]; 
     167        for (int k=1; k < degree; k++) { 
     168            total = fma(total, x, coeff[k]); 
     169        } 
     170        array[index] = total * x; 
     171    } 
     172} 
     173""" 
     174 
     175class OpenclCalculator(NumpyCalculator): 
     176    polyval1f = None 
     177    polyval1d = None 
     178    def __init__(self, dims, dtype=PRECISION): 
     179        env = sasmodels.kernelcl.environment() 
     180        context = env.get_context(dtype) 
     181        if dtype == np.dtype('f'): 
     182            if self.polyval1f is None: 
     183                program = sasmodels.kernelcl.compile_model( 
     184                    context, POLYVAL1_KERNEL, dtype, fast=USE_FAST) 
     185                # Assume context is always the same for a given dtype 
     186                OpenclCalculator.polyval1f = program.polyval1 
     187            self.dtype = dtype 
     188            self.complex_dtype = np.dtype('F') 
     189            self.polyval1 = self.polyval1f 
    96190        else: 
    97             q_range_corners = np.sqrt(2.) * q_range 
    98             nq_corners = int(np.sqrt(2.) * nq/2) 
    99             q_corners = np.linspace(0, q_range_corners, nq_corners+1)[1:] 
    100             q_calc = [q_corners] 
    101             self._qxy = np.sqrt(qx**2 + qy**2) 
    102             self._edges = bin_edges(q_corners) 
    103             self._norm = np.histogram(self._qxy, bins=self._edges)[0] 
    104  
    105         self.q_calc = q_calc 
    106         self.q_range = q_range 
    107  
    108     def apply(self, theory, background=0.0, outfile="", plot=False): 
     191            if self.polyval1d is None: 
     192                program = sasmodels.kernelcl.compile_model( 
     193                    context, POLYVAL1_KERNEL, dtype, fast=False) 
     194                # Assume context is always the same for a given dtype 
     195                OpenclCalculator.polyval1d = program.polyval1 
     196            self.dtype = dtype 
     197            self.complex_type = np.dtype('D') 
     198            self.polyval1 = self.polyval1d 
     199        self.queue = env.get_queue(dtype) 
     200        self.plan = pyfft.cl.Plan(dims, queue=self.queue) 
     201 
     202    def fft(self, Iq): 
     203        # forward transform 
    109204        #t0 = time.time() 
    110         if self.is2d: 
    111             Iq_calc = theory 
    112         else: 
    113             q_corners = self.q_calc[0] 
    114             Iq_calc = np.interp(self._qxy, q_corners, theory) 
    115         Iq_calc = Iq_calc.reshape(self.nq, self.nq) 
    116         #plotxy(Iq_calc); import pylab; pylab.figure() 
    117         if self.power > 0: 
    118             Iqxy = scattering_power(Iq_calc, self.power) 
    119             powers = [] 
    120         else: 
    121             Iqxy, powers = multiple_scattering( 
    122                 Iq_calc, self.probability, 
    123                 coverage=0.99, 
    124                 return_powers=(outfile!="") or plot, 
    125                 ) 
    126         #print("multiple scattering calc time", time.time()-t0) 
    127  
    128         #plotxy(Iqxy); import pylab; pylab.figure() 
    129         if self.is2d: 
    130             if outfile and powers: 
    131                 data = np.vstack([Ipower.flatten() for Ipower in powers]).T 
    132                 np.savetxt(outfile + "_powers.txt", data) 
    133             if outfile: 
    134                 data = np.vstack(Iq_calc).T 
    135                 np.savetxt(outfile + ".txt", data) 
    136             if plot: 
    137                 import pylab 
    138                 plotxy(Iq_calc) 
    139                 pylab.title("single scattering") 
    140                 pylab.figure() 
    141  
    142             return Iqxy + background 
    143         else: 
    144             # circular average, no anti-aliasing 
    145             Iq = np.histogram(self._qxy, bins=self._edges, weights=Iqxy)[0]/self._norm 
    146             if powers: 
    147                 Iq_powers = [np.histogram(self._qxy, bins=self._edges, weights=Ipower)[0]/self._norm 
    148                              for Ipower in powers]#  or powers[:5] for the first five 
    149  
    150             if outfile: 
    151                 data = np.vstack([q_corners, theory, Iq]).T 
    152                 np.savetxt(outfile + ".txt", data) 
    153             if outfile and powers: 
    154                 # circular average, no anti-aliasing for individual powers 
    155                 data = np.vstack([q_corners] + Iq_powers).T 
    156                 np.savetxt(outfile + "_powers.txt", data) 
    157             if plot: 
    158                 import pylab 
    159                 plotxy(Iqxy) 
    160                 pylab.title("multiple scattering") 
    161                 pylab.figure() 
    162             if plot and powers: 
    163                 import pylab 
    164                 L = -np.log(1-self.probability) 
    165                 pylab.loglog(q_corners, Iq, label="total %g"%self.probability) 
    166                 for n, Ipower in enumerate(Iq_powers): 
    167                     k = n+1 
    168                     w = L**(k-1)/gamma(k+1) 
    169                     pylab.loglog(q_corners, w*Ipower, label="scattering**%d"%k) 
    170                 pylab.legend() 
    171                 pylab.figure() 
    172  
    173             return q_corners, Iq + background 
    174  
    175 def scattering_power(Iq, n): 
     205        data = np.asarray(Iq, self.complex_dtype) 
     206        gpu_data = cl_array.to_device(self.queue, data) 
     207        self.plan.execute(gpu_data.data) 
     208        result = gpu_data.get() 
     209        #print("fft time", time.time()-t0) 
     210        return result 
     211 
     212    def ifft(self, fourier_frame): 
     213        # inverse transform 
     214        #t0 = time.time() 
     215        data = np.asarray(fourier_frame, self.complex_dtype) 
     216        gpu_data = cl_array.to_device(self.queue, data) 
     217        self.plan.execute(gpu_data.data, inverse=True) 
     218        result = gpu_data.get() 
     219        #print("ifft time", time.time()-t0) 
     220        return result 
     221 
     222    def multiple_scattering(self, Iq, p, coverage=0.99): 
     223        #t0 = time.time() 
     224        coeffs = scattering_coeffs(p, coverage) 
     225        scale = np.sum(Iq) 
     226        poly = np.asarray(coeffs[::-1], self.dtype) 
     227        frame = _forward_shift(Iq/scale, dtype=self.complex_dtype) 
     228        gpu_data = cl_array.to_device(self.queue, frame) 
     229        gpu_poly = cl_array.to_device(self.queue, poly) 
     230        self.plan.execute(gpu_data.data) 
     231        degree, n = poly.shape[0], frame.shape[0]*frame.shape[1] 
     232        self.polyval1( 
     233            self.queue, [n], None, 
     234            np.int32(degree), gpu_poly.data, np.int32(n), gpu_data.data) 
     235        self.plan.execute(gpu_data.data, inverse=True) 
     236        frame = gpu_data.get() 
     237        #result = scale * _inverse_shift(frame.real, dtype=self.dtype) 
     238        result = scale * _inverse_shift(frame.real, dtype=self.dtype) 
     239        #print("multiscat time", time.time()-t0) 
     240        return result 
     241 
     242Calculator = OpenclCalculator if HAVE_OPENCL else NumpyCalculator 
     243 
     244def scattering_powers(Iq, n, dtype='f', transform=None): 
     245    r""" 
     246    Calculate the scattering powers up to n. 
     247 
     248    This includes 1 even though it should just be Iq itself 
     249 
     250    The frames are unweighted; to weight scale by $\lambda^k e^{-\lambda}/k!$. 
    176251    """ 
    177     Calculate the nth scattering power as a distribution.  To get the 
    178     weighted contribution, scale by $\lambda^k e^{-\lambda}/k!$. 
    179     """ 
     252    if transform is None: 
     253        nx, ny = Iq.shape 
     254        transform = Calculator(dims=(nx*2, ny*2), dtype=dtype) 
    180255    scale = np.sum(Iq) 
    181     F = _forward_fft(Iq/scale) 
    182     result = _inverse_fft(F**n) 
    183     return result 
    184  
    185 def multiple_scattering(Iq, p, coverage=0.99, return_powers=False): 
    186     """ 
    187     Compute multiple scattering for I(q) given scattering probability p. 
    188  
    189     Given a probability p of scattering with the thickness, the expected 
    190     number of scattering events, $\lambda$ is $-\log(1 - p)$, giving a 
    191     Poisson weighted sum of single, double, triple, etc. scattering patterns. 
    192     The number of patterns used is based on coverage (default 99%). 
    193     """ 
     256    frame = _forward_shift(Iq/scale, dtype=dtype) 
     257    F = transform.fft(frame) 
     258    powers = [scale * _inverse_shift(transform.ifft(F**(k+1)).real, dtype=dtype) 
     259              for k in range(n)] 
     260    return powers 
     261 
     262def scattering_coeffs(p, coverage=0.99): 
    194263    L = -np.log(1-p) 
    195264    num_scatter = truncated_poisson_invcdf(coverage, L) 
    196  
    197     # Compute multiple scattering via convolution. 
    198     scale = np.sum(Iq) 
    199     F = _forward_fft(Iq/scale) 
    200     coeffs = [L**(k-1)/gamma(k+1) for k in range(1, num_scatter+1)] 
    201     multiple_scattering = F * np.polyval(coeffs[::-1], F) 
    202     result = scale * _inverse_fft(multiple_scattering) 
    203  
    204     if return_powers: 
    205         powers = [scale * _inverse_fft(F**k) for k in range(1, num_scatter+1)] 
    206     else: 
    207         powers = [] 
    208  
    209     return result, powers 
     265    coeffs = [L**k/gamma(k+2) for k in range(num_scatter)] 
     266    return coeffs 
    210267 
    211268def truncated_poisson_invcdf(p, L): 
    212     """ 
     269    r""" 
    213270    Return smallest k such that cdf(k; L) > p from the truncated Poisson 
    214271    probability excluding k=0 
     
    224281    return k 
    225282 
    226 def _forward_fft(Iq): 
     283def _forward_shift(Iq, dtype=PRECISION): 
    227284    # Prepare padded array and forward transform 
    228285    nq = Iq.shape[0] 
    229286    half_nq = nq//2 
    230     frame = np.zeros((2*nq, 2*nq)) 
     287    frame = np.zeros((2*nq, 2*nq), dtype=dtype) 
    231288    frame[:half_nq, :half_nq] = Iq[half_nq:, half_nq:] 
    232289    frame[-half_nq:, :half_nq] = Iq[:half_nq, half_nq:] 
    233290    frame[:half_nq, -half_nq:] = Iq[half_nq:, :half_nq] 
    234291    frame[-half_nq:, -half_nq:] = Iq[:half_nq, :half_nq] 
    235     fourier_frame = np.fft.fft2(frame) 
    236     return fourier_frame 
    237  
    238 def _inverse_fft(fourier_frame): 
     292    return frame 
     293 
     294def _inverse_shift(frame, dtype=PRECISION): 
    239295    # Invert the transform and recover the transformed data 
    240     nq = fourier_frame.shape[0]//2 
     296    nq = frame.shape[0]//2 
    241297    half_nq = nq//2 
    242     frame = np.fft.ifft2(fourier_frame).real 
    243     Iq = np.empty((nq, nq)) 
     298    Iq = np.empty((nq, nq), dtype=dtype) 
    244299    Iq[half_nq:, half_nq:] = frame[:half_nq, :half_nq] 
    245300    Iq[:half_nq, half_nq:] = frame[-half_nq:, :half_nq] 
     
    247302    Iq[:half_nq, :half_nq] = frame[-half_nq:, -half_nq:] 
    248303    return Iq 
     304 
     305 
     306class MultipleScattering(Resolution): 
     307    def __init__(self, qmin=None, qmax=None, nq=None, window=2, 
     308                 probability=None, coverage=0.99, 
     309                 is2d=False, resolution=None, 
     310                 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        """ 
     335        # Infer qmin, qmax from instrument resolution calculator, if present 
     336        if resolution is not None: 
     337            is2d = hasattr(resolution, 'qx_data') 
     338            if is2d: 
     339                # 2D data 
     340                if qmax is None: 
     341                    qx_calc, qy_calc = resolution.q_calc 
     342                    qmax = np.sqrt(np.max(qx_calc**2 + qy_calc**2)) 
     343                if qmin is None and nq is None: 
     344                    qx, qy = resolution.data.x_bins, resolution.data.y_bins 
     345                    if qx and qy: 
     346                        dx = (np.max(qx) - np.min(qx)) / len(qx) 
     347                        dy = (np.max(qy) - np.min(qy)) / len(qy) 
     348                    else: 
     349                        qx, qy = resolution.data.qx_data, resolution.data.qy_data 
     350                        steps = np.sqrt(len(qx)) 
     351                        dx = (np.max(qx) - np.min(qx)) / steps 
     352                        dy = (np.max(qy) - np.min(qy)) / steps 
     353                    qmin = min(dx, dy) 
     354            else: 
     355                # 1D data 
     356                if qmax is None: 
     357                    qmax = np.max(resolution.q_calc) 
     358                if qmin is None and nq is None: 
     359                    qmin = np.min(np.abs(resolution.q_calc)) 
     360 
     361        # estimate nq from qmin, qmax if not given explicitly 
     362        q_range = qmax * window 
     363        if nq is None: 
     364            nq = 2**np.ceil(np.log2(q_range/qmin)) 
     365        nq = int(nq) 
     366        # Compute available qmin based on nq 
     367        qmin = 2*q_range / nq 
     368        #print(nq) 
     369 
     370        # remember input parameters 
     371        self.qmax = qmax 
     372        self.qmin = qmin 
     373        self.nq = nq 
     374        self.probability = probability 
     375        self.coverage = coverage 
     376        self.is2d = is2d 
     377        self.window = window 
     378        self.resolution = resolution 
     379 
     380        # Determine the q values to calculate 
     381        q = np.linspace(-q_range, q_range, nq) 
     382        qx, qy = np.meshgrid(q, q) 
     383        if is2d: 
     384            q_calc = (qx.flatten(), qy.flatten()) 
     385        else: 
     386            # For 1-D patterns, compute q from the center to the corners and 
     387            # interpolate from there into the individual pixels.  Given that 
     388            # nq represents the points in [-qmax*windows, qmax*window], 
     389            # then using 2*sqrt(2)*nq/2 will oversample the points in q by 
     390            # a factor of two relative to the pixels. 
     391            q_range_to_corner = np.sqrt(2.) * q_range 
     392            nq_to_corner = 10*int(np.ceil(np.sqrt(2.) * nq)) 
     393            q_to_corner = np.linspace(0, q_range_to_corner, nq_to_corner+1)[1:] 
     394            q_calc = (q_to_corner,) 
     395            # Remember the q radii of the calculated points 
     396            self._radius = np.sqrt(qx**2 + qy**2) 
     397            #self._q = q_to_corner 
     398        self._q_steps = q 
     399        self.q_calc = q_calc 
     400 
     401        # TODO: use cleaner data representation than that from sasview 
     402        # Resolution function forwards underlying q data (for plotting, etc?) 
     403        if is2d: 
     404            if resolution is not None: 
     405                # forward resolution function info to multiscattering 
     406                self.qx_data = resolution.qx_data 
     407                self.qy_data = resolution.qy_data 
     408            else: 
     409                # no underlying resolution function, but make it look like there is 
     410                self.qx_data, self.qy_data = q_calc 
     411        else: 
     412            # 1-D radial profile is determined by the q values we need to 
     413            # compute, either for the calculated q values for the resolution 
     414            # function (if any) or for the raw q values desired 
     415            self._q = np.linspace(qmin, qmax, nq//(2*window)) 
     416            self._edges = bin_edges(self._q) 
     417            self._norm, _ = np.histogram(self._radius, bins=self._edges) 
     418            if resolution is not None: 
     419                self.q = resolution.q 
     420            else: 
     421                # no underlying resolution function, but make it look like there is 
     422                self.q = self._q 
     423 
     424        # Prepare the multiple scattering calculator (either numpy or OpenCL) 
     425        self.transform = Calculator((2*nq, 2*nq), dtype=dtype) 
     426 
     427    def apply(self, theory): 
     428        if self.is2d: 
     429            Iq_calc = theory 
     430        else: 
     431            Iq_calc = np.interp(self._radius, self.q_calc[0], theory) 
     432        Iq_calc = Iq_calc.reshape(self.nq, self.nq) 
     433 
     434        probability = self.probability() if callable(self.probability) else self.probability 
     435        coverage = self.coverage 
     436        #t0 = time.time() 
     437        Iqxy = self.transform.multiple_scattering(Iq_calc, probability, coverage) 
     438        #print("multiple scattering calc time", time.time()-t0) 
     439 
     440        # remember the intermediate result in case we want to see it later 
     441        self.Iqxy = Iqxy 
     442 
     443        if self.is2d: 
     444            if self.resolution is not None: 
     445                Iqxy = self.resolution.apply(Iqxy) 
     446            return Iqxy 
     447        else: 
     448            # remember the intermediate result in case we want to see it later 
     449            Iq = self.radial_profile(Iqxy) 
     450            self.Iq = Iq 
     451            if self.resolution is not None: 
     452                q = self._q 
     453                Iq_res = np.interp(np.abs(self.resolution.q_calc), q, self.Iq) 
     454                _ = """ 
     455                k = 6 
     456                print("q theory", q[:k]) 
     457                print("Iq theory", theory[:k]) 
     458                print("interp NaN?", np.any(np.isnan(Iq_calc))) 
     459                print("convolved NaN?", np.any(np.isnan(Iqxy))) 
     460                print("Iq intgrated", self.Iq[:k]) 
     461                print("radius", self._radius[self.nq/2,self.nq/2:self.nq/2+k]) 
     462                print("q edges", self._edges[:k+1]) 
     463                print("Iq norm", self._norm[:k]) 
     464                print("res q", self.resolution.q_calc[:k]) 
     465                print("Iq res", Iq_res[:k]) 
     466                #print(Iq) 
     467                #print(Iq_res) 
     468                """ 
     469                Iq = self.resolution.apply(Iq_res) 
     470            return Iq 
     471 
     472    def radial_profile(self, Iqxy): 
     473        # circular average, no anti-aliasing 
     474        Iq = np.histogram(self._radius, bins=self._edges, weights=Iqxy)[0]/self._norm 
     475        return Iq 
     476 
     477 
     478def annular_average(qxy, Iqxy, qbins): 
     479    """ 
     480    Compute annular average of points at 
     481    """ 
     482    qxy, Iqxy = qxy.flatten(), Iqxy.flatten() 
     483    index = np.argsort(qxy) 
     484    qxy, Iqxy = qxy[index], Iqxy[index] 
     485    print(qxy.shape, Iqxy.shape, index.shape, qbins.shape) 
     486    #values = rebin(np.vstack((0., qxy)), Iqxy, qbins) 
     487    integral = np.cumsum(Iqxy) 
     488    Io = np.diff(np.interp(qbins, qxy, integral, left=0.)) 
     489    # normalize by area of annulus 
     490    # TODO: doesn't properly account for box. 
     491    # Need to chop off chords that are greater than the box edges 
     492    # (left, right, top and bottom), then add back in the corners 
     493    # chopped off by both. https://en.wikipedia.org/wiki/Circular_segment 
     494    norms = np.diff(pi*qbins**2) 
     495    return Io/norms 
     496 
     497def rebin(x, I, xo): 
     498    """ 
     499    Rebin from edges *x*, bins I into edges *xo*. 
     500 
     501    *x* and *xo* should be monotonically increasing. 
     502 
     503    If *x* has duplicate values, then all corresponding values at I(x) will 
     504    be effectively summed into the same bin.  If *xo* has duplicate values, 
     505    the first bin will contain the entire contents and subsequent bins will 
     506    contain zeros. 
     507    """ 
     508    integral = np.cumsum(I) 
     509    Io = np.diff(np.interp(xo, x[1:], integral, left=0.)) 
     510    return Io 
     511 
    249512 
    250513def parse_pars(model, opts): 
     
    259522        'magnetic': False, 
    260523        'values': opts.pars, 
    261         'show_pars': True, 
     524        #'show_pars': True, 
     525        'show_pars': False, 
    262526        'is2d': opts.is2d, 
    263527    } 
     
    271535        formatter_class=argparse.ArgumentDefaultsHelpFormatter, 
    272536        ) 
    273     parser.add_argument('-k', '--power', type=int, default=0, help="show pattern for nth scattering") 
    274537    parser.add_argument('-p', '--probability', type=float, default=0.1, help="scattering probability") 
    275538    parser.add_argument('-n', '--nq', type=int, default=1024, help='number of mesh points') 
     
    287550    model = core.load_model(opts.model) 
    288551    pars = parse_pars(model, opts) 
    289     res = MultipleScattering(opts.qmax, opts.nq, opts.probability, opts.is2d, 
    290                              window=opts.window, power=opts.power) 
     552    res = MultipleScattering(qmax=opts.qmax, nq=opts.nq, window=opts.window, 
     553                             probability=opts.probability, is2d=opts.is2d) 
    291554    kernel = model.make_kernel(res.q_calc) 
    292555    #print(pars) 
    293556    bg = pars.get('background', 0.0) 
    294557    pars['background'] = 0.0 
    295     Iq_calc = call_kernel(kernel, pars) 
    296     Iq = res.apply(Iq_calc, background=bg, outfile=opts.outfile, plot=True) 
    297     plotxy(Iq) 
     558    theory = call_kernel(kernel, pars) 
     559    Iq = res.apply(theory) + bg 
     560    plot_and_save_powers(res, theory, Iq, outfile=opts.outfile, background=bg) 
     561 
     562def plot_and_save_powers(res, theory, result, plot=True, outfile="", background=0.): 
    298563    import pylab 
    299     if opts.power > 0: 
    300         pylab.title('scattering power %d'%opts.power) 
     564    probability, coverage = res.probability, res.coverage 
     565    weights = scattering_coeffs(probability, coverage) 
     566 
     567    # cribbed from MultipleScattering.apply 
     568    if res.is2d: 
     569        Iq_calc = theory 
    301570    else: 
    302         pylab.title('multiple scattering with fraction %g'%opts.probability) 
     571        Iq_calc = np.interp(res._radius, res.q_calc[0], theory) 
     572    Iq_calc = Iq_calc.reshape(res.nq, res.nq) 
     573 
     574    # Compute the scattering powers for 1, 2, ... n scattering events 
     575    powers = scattering_powers(Iq_calc, len(weights)) 
     576 
     577    #plotxy(Iqxy); import pylab; pylab.figure() 
     578    if res.is2d: 
     579        if outfile: 
     580            data = np.vstack([Ipower.flatten() for Ipower in powers]).T 
     581            np.savetxt(outfile + "_powers.txt", data) 
     582            data = np.vstack(Iq_calc).T 
     583            np.savetxt(outfile + ".txt", data) 
     584        if plot: 
     585            plotxy((res._q_steps, res._q_steps), Iq_calc) 
     586            pylab.title("single scattering F") 
     587            for k, v in enumerate(powers[1:]): 
     588                pylab.figure() 
     589                plotxy((res._q_steps, res._q_steps), v+background) 
     590                pylab.title("multiple scattering F^%d" % (k+2)) 
     591            pylab.figure() 
     592            plotxy((res._q_steps, res._q_steps), res.Iqxy+background) 
     593            pylab.title("total scattering for p=%g" % probability) 
     594    else: 
     595        q = res._q 
     596        Iq_powers = [res.radial_profile(Iqxy) for Iqxy in powers] 
     597        if outfile: 
     598            data = np.vstack([q, theory, res.Iq]).T 
     599            np.savetxt(outfile + ".txt", data) 
     600            # circular average, no anti-aliasing for individual powers 
     601            data = np.vstack([q] + Iq_powers).T 
     602            np.savetxt(outfile + "_powers.txt", data) 
     603        if plot: 
     604            # Plot 2D pattern for total scattering 
     605            plotxy((res._q_steps, res._q_steps), res.Iqxy+background) 
     606            pylab.title("total scattering for p=%g" % probability) 
     607            pylab.figure() 
     608 
     609            # Plot 1D pattern for partial scattering 
     610            pylab.loglog(q, res.Iq+background, label="total for p=%g"%probability) 
     611            #new_annulus = annular_average(res._radius, res.Iqxy, res._edges) 
     612            #pylab.loglog(q, new_annulus+background, label="new total for p=%g"%probability) 
     613            for n, (w, Ipower) in enumerate(zip(weights, Iq_powers)): 
     614                pylab.loglog(q, w*Ipower+background, label="scattering^%d"%(n+1)) 
     615            pylab.legend() 
     616            pylab.title('total scattering for p=%g' % probability) 
    303617    pylab.show() 
    304618 
    305 def plotxy(Iq): 
     619def plotxy(q, Iq): 
    306620    import pylab 
    307     if isinstance(Iq, tuple): 
    308         q, Iq = Iq 
    309         pylab.loglog(q, Iq) 
     621    # q is a tuple of (q,) or (qx, qy) 
     622    if len(q) == 1: 
     623        pylab.loglog(q[0], Iq) 
    310624    else: 
    311         data = Iq+0. 
    312         data[Iq <= 0] = np.min(Iq[Iq>0])/2 
     625        data = Iq.copy() 
     626        data[Iq <= 0] = np.min(Iq[Iq > 0])/2 
    313627        pylab.imshow(np.log10(data)) 
    314     #import pylab; pylab.show() 
    315628 
    316629if __name__ == "__main__": 
  • sasmodels/data.py

    r2d81cfe rf549e37  
    706706    else: 
    707707        vmin_scaled, vmax_scaled = vmin, vmax 
    708     plt.imshow(plottable.reshape(len(data.x_bins), len(data.y_bins)), 
     708    nx, ny = len(data.x_bins), len(data.y_bins) 
     709    x_bins, y_bins, image = _build_matrix(data, plottable) 
     710    plt.imshow(image, 
    709711               interpolation='nearest', aspect=1, origin='lower', 
    710712               extent=[xmin, xmax, ymin, ymax], 
     
    714716    return vmin, vmax 
    715717 
     718 
     719# === The following is modified from sas.sasgui.plottools.PlotPanel 
     720def _build_matrix(self, plottable): 
     721    """ 
     722    Build a matrix for 2d plot from a vector 
     723    Returns a matrix (image) with ~ square binning 
     724    Requirement: need 1d array formats of 
     725    self.data, self.qx_data, and self.qy_data 
     726    where each one corresponds to z, x, or y axis values 
     727 
     728    """ 
     729    # No qx or qy given in a vector format 
     730    if self.qx_data is None or self.qy_data is None \ 
     731            or self.qx_data.ndim != 1 or self.qy_data.ndim != 1: 
     732        return self.x_bins, self.y_bins, plottable 
     733 
     734    # maximum # of loops to fillup_pixels 
     735    # otherwise, loop could never stop depending on data 
     736    max_loop = 1 
     737    # get the x and y_bin arrays. 
     738    x_bins, y_bins = _get_bins(self) 
     739    # set zero to None 
     740 
     741    #Note: Can not use scipy.interpolate.Rbf: 
     742    # 'cause too many data points (>10000)<=JHC. 
     743    # 1d array to use for weighting the data point averaging 
     744    #when they fall into a same bin. 
     745    weights_data = np.ones([self.data.size]) 
     746    # get histogram of ones w/len(data); this will provide 
     747    #the weights of data on each bins 
     748    weights, xedges, yedges = np.histogram2d(x=self.qy_data, 
     749                                             y=self.qx_data, 
     750                                             bins=[y_bins, x_bins], 
     751                                             weights=weights_data) 
     752    # get histogram of data, all points into a bin in a way of summing 
     753    image, xedges, yedges = np.histogram2d(x=self.qy_data, 
     754                                           y=self.qx_data, 
     755                                           bins=[y_bins, x_bins], 
     756                                           weights=plottable) 
     757    # Now, normalize the image by weights only for weights>1: 
     758    # If weight == 1, there is only one data point in the bin so 
     759    # that no normalization is required. 
     760    image[weights > 1] = image[weights > 1] / weights[weights > 1] 
     761    # Set image bins w/o a data point (weight==0) as None (was set to zero 
     762    # by histogram2d.) 
     763    image[weights == 0] = None 
     764 
     765    # Fill empty bins with 8 nearest neighbors only when at least 
     766    #one None point exists 
     767    loop = 0 
     768 
     769    # do while loop until all vacant bins are filled up up 
     770    #to loop = max_loop 
     771    while (weights == 0).any(): 
     772        if loop >= max_loop:  # this protects never-ending loop 
     773            break 
     774        image = _fillup_pixels(self, image=image, weights=weights) 
     775        loop += 1 
     776 
     777    return x_bins, y_bins, image 
     778 
     779def _get_bins(self): 
     780    """ 
     781    get bins 
     782    set x_bins and y_bins into self, 1d arrays of the index with 
     783    ~ square binning 
     784    Requirement: need 1d array formats of 
     785    self.qx_data, and self.qy_data 
     786    where each one corresponds to  x, or y axis values 
     787    """ 
     788    # find max and min values of qx and qy 
     789    xmax = self.qx_data.max() 
     790    xmin = self.qx_data.min() 
     791    ymax = self.qy_data.max() 
     792    ymin = self.qy_data.min() 
     793 
     794    # calculate the range of qx and qy: this way, it is a little 
     795    # more independent 
     796    x_size = xmax - xmin 
     797    y_size = ymax - ymin 
     798 
     799    # estimate the # of pixels on each axes 
     800    npix_y = int(np.floor(np.sqrt(len(self.qy_data)))) 
     801    npix_x = int(np.floor(len(self.qy_data) / npix_y)) 
     802 
     803    # bin size: x- & y-directions 
     804    xstep = x_size / (npix_x - 1) 
     805    ystep = y_size / (npix_y - 1) 
     806 
     807    # max and min taking account of the bin sizes 
     808    xmax = xmax + xstep / 2.0 
     809    xmin = xmin - xstep / 2.0 
     810    ymax = ymax + ystep / 2.0 
     811    ymin = ymin - ystep / 2.0 
     812 
     813    # store x and y bin centers in q space 
     814    x_bins = np.linspace(xmin, xmax, npix_x) 
     815    y_bins = np.linspace(ymin, ymax, npix_y) 
     816 
     817    return x_bins, y_bins 
     818 
     819def _fillup_pixels(self, image=None, weights=None): 
     820    """ 
     821    Fill z values of the empty cells of 2d image matrix 
     822    with the average over up-to next nearest neighbor points 
     823 
     824    :param image: (2d matrix with some zi = None) 
     825 
     826    :return: image (2d array ) 
     827 
     828    :TODO: Find better way to do for-loop below 
     829 
     830    """ 
     831    # No image matrix given 
     832    if image is None or np.ndim(image) != 2 \ 
     833            or np.isfinite(image).all() \ 
     834            or weights is None: 
     835        return image 
     836    # Get bin size in y and x directions 
     837    len_y = len(image) 
     838    len_x = len(image[1]) 
     839    temp_image = np.zeros([len_y, len_x]) 
     840    weit = np.zeros([len_y, len_x]) 
     841    # do for-loop for all pixels 
     842    for n_y in range(len(image)): 
     843        for n_x in range(len(image[1])): 
     844            # find only null pixels 
     845            if weights[n_y][n_x] > 0 or np.isfinite(image[n_y][n_x]): 
     846                continue 
     847            else: 
     848                # find 4 nearest neighbors 
     849                # check where or not it is at the corner 
     850                if n_y != 0 and np.isfinite(image[n_y - 1][n_x]): 
     851                    temp_image[n_y][n_x] += image[n_y - 1][n_x] 
     852                    weit[n_y][n_x] += 1 
     853                if n_x != 0 and np.isfinite(image[n_y][n_x - 1]): 
     854                    temp_image[n_y][n_x] += image[n_y][n_x - 1] 
     855                    weit[n_y][n_x] += 1 
     856                if n_y != len_y - 1 and np.isfinite(image[n_y + 1][n_x]): 
     857                    temp_image[n_y][n_x] += image[n_y + 1][n_x] 
     858                    weit[n_y][n_x] += 1 
     859                if n_x != len_x - 1 and np.isfinite(image[n_y][n_x + 1]): 
     860                    temp_image[n_y][n_x] += image[n_y][n_x + 1] 
     861                    weit[n_y][n_x] += 1 
     862                # go 4 next nearest neighbors when no non-zero 
     863                # neighbor exists 
     864                if n_y != 0 and n_x != 0 and \ 
     865                        np.isfinite(image[n_y - 1][n_x - 1]): 
     866                    temp_image[n_y][n_x] += image[n_y - 1][n_x - 1] 
     867                    weit[n_y][n_x] += 1 
     868                if n_y != len_y - 1 and n_x != 0 and \ 
     869                        np.isfinite(image[n_y + 1][n_x - 1]): 
     870                    temp_image[n_y][n_x] += image[n_y + 1][n_x - 1] 
     871                    weit[n_y][n_x] += 1 
     872                if n_y != len_y and n_x != len_x - 1 and \ 
     873                        np.isfinite(image[n_y - 1][n_x + 1]): 
     874                    temp_image[n_y][n_x] += image[n_y - 1][n_x + 1] 
     875                    weit[n_y][n_x] += 1 
     876                if n_y != len_y - 1 and n_x != len_x - 1 and \ 
     877                        np.isfinite(image[n_y + 1][n_x + 1]): 
     878                    temp_image[n_y][n_x] += image[n_y + 1][n_x + 1] 
     879                    weit[n_y][n_x] += 1 
     880 
     881    # get it normalized 
     882    ind = (weit > 0) 
     883    image[ind] = temp_image[ind] / weit[ind] 
     884 
     885    return image 
     886 
     887 
    716888def demo(): 
    717889    # type: () -> None 
  • sasmodels/direct_model.py

    r2d81cfe rd18d6dd  
    345345            self._kernel = self._model.make_kernel(self._kernel_inputs) 
    346346 
     347        # Need to pull background out of resolution for multiple scattering 
     348        background = pars.get('background', 0.) 
     349        pars = pars.copy() 
     350        pars['background'] = 0. 
     351 
    347352        Iq_calc = call_kernel(self._kernel, pars, cutoff=cutoff) 
    348353        # Storing the calculated Iq values so that they can be plotted. 
     
    357362                np.reshape(Iq_calc, (self.resolution.ny, self.resolution.nx)) 
    358363            ) 
    359         return result 
     364        return result + background 
    360365 
    361366 
  • sasmodels/generate.py

    r108e70e r6cbdcd4  
    214214 
    215215EXTERNAL_DIR = 'sasmodels-data' 
    216 DATA_PATH = get_data_path(EXTERNAL_DIR, 'kernel_template.c') 
     216DATA_PATH = get_data_path(EXTERNAL_DIR, 'kernel_iq.c') 
    217217MODEL_PATH = joinpath(DATA_PATH, 'models') 
    218218 
     
    389389    # TODO: fails DRY; templates appear two places. 
    390390    model_templates = [joinpath(DATA_PATH, filename) 
    391                        for filename in ('kernel_header.c', 'kernel_iq.cl')] 
     391                       for filename in ('kernel_header.c', 'kernel_iq.c')] 
    392392    source_files = (model_sources(model_info) 
    393393                    + model_templates 
  • sasmodels/kernelcl.py

    rb4272a2 r6cbdcd4  
    208208    Build a model to run on the gpu. 
    209209 
    210     Returns the compiled program and its type.  The returned type will 
    211     be float32 even if the desired type is float64 if any of the 
    212     devices in the context do not support the cl_khr_fp64 extension. 
     210    Returns the compiled program and its type. 
     211 
     212    Raises an error if the desired precision is not available. 
    213213    """ 
    214214    dtype = np.dtype(dtype) 
  • sasmodels/models/core_shell_parallelepiped.py

    r97be877 r5bc373b  
    2020$A < B < C$. 
    2121 
    22 .. image:: img/core_shell_parallelepiped_geometry.jpg 
     22.. figure:: img/parallelepiped_geometry.jpg 
     23 
     24   Core of the core shell Parallelepiped with the corresponding definition 
     25   of sides. 
     26 
    2327 
    2428There are rectangular "slabs" of thickness $t_A$ that add to the $A$ dimension 
     
    2630$(=t_C)$ faces. The projection in the $AB$ plane is then 
    2731 
    28 .. image:: img/core_shell_parallelepiped_projection.jpg 
     32.. figure:: img/core_shell_parallelepiped_projection.jpg 
     33 
     34   AB cut through the core-shell parllelipiped showing the cross secion of 
     35   four of the six shell slabs 
    2936 
    3037The volume of the solid is 
     
    4956.. math:: 
    5057 
    51     F(Q) 
     58    F(q) 
    5259    &= (\rho_\text{core}-\rho_\text{solvent}) 
    5360       S(Q_A, A) S(Q_B, B) S(Q_C, C) \\ 
    5461    &+ (\rho_\text{A}-\rho_\text{solvent}) 
    55         \left[S(Q_A, A+2t_A) - S(Q_A, Q)\right] S(Q_B, B) S(Q_C, C) \\ 
     62        \left[S(Q_A, A+2t_A) - S(Q_A, A)\right] S(Q_B, B) S(Q_C, C) \\ 
    5663    &+ (\rho_\text{B}-\rho_\text{solvent}) 
    5764        S(Q_A, A) \left[S(Q_B, B+2t_B) - S(Q_B, B)\right] S(Q_C, C) \\ 
     
    6976.. math:: 
    7077 
    71     Q_A &= \sin\alpha \sin\beta \\ 
    72     Q_B &= \sin\alpha \cos\beta \\ 
    73     Q_C &= \cos\alpha 
     78    Q_A &= q \sin\alpha \sin\beta \\ 
     79    Q_B &= q \sin\alpha \cos\beta \\ 
     80    Q_C &= q \cos\alpha 
    7481 
    7582 
     
    95102and length $(C+2t_C)$ values, after appropriately sorting the three dimensions 
    96103to give an oblate or prolate particle, to give an effective radius, 
    97 for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
     104for $S(q)$ when $P(q) * S(q)$ is applied. 
    98105 
    99106For 2d data the orientation of the particle is required, described using 
  • sasmodels/models/parallelepiped.py

    ref07e95 r5bc373b  
    1010 
    1111 This model calculates the scattering from a rectangular parallelepiped 
    12  (\:numref:`parallelepiped-image`\). 
     12 (:numref:`parallelepiped-image`). 
    1313 If you need to apply polydispersity, see also :ref:`rectangular-prism`. 
    1414 
Note: See TracChangeset for help on using the changeset viewer.