Changeset bc248f8 in sasmodels


Ignore:
Timestamp:
Mar 4, 2018 8:43:57 PM (5 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, 49d1f8b8
Parents:
1ee5ac6
Message:

add support for instrumental resolution to multiple scattering calculation and create 1D fitting example

Files:
1 added
1 edited

Legend:

Unmodified
Added
Removed
  • 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__": 
Note: See TracChangeset for help on using the changeset viewer.