Source code for pyrocko.parstack

# http://pyrocko.org - GPLv3
#
# The Pyrocko Developers, 21st Century
# ---|P------/S----------~Lg----------

'''
Fast weight-delay-and-sum implementation for seismic array processing and
migration.
'''

import numpy as num

from . import parstack_ext

try:
    range = xrange
except NameError:
    pass


def parstack(arrays, offsets, shifts, weights, method,
             lengthout=-1,
             offsetout=0,
             result=None,
             nparallel=None,
             dtype=num.float64,
             impl='openmp'):

    if nparallel is None:
        import multiprocessing
        nparallel = multiprocessing.cpu_count()

    narrays = offsets.size
    assert len(arrays) == narrays
    nshifts = shifts.size // narrays
    assert shifts.shape == (nshifts, narrays)
    shifts = num.reshape(shifts, (nshifts*narrays))
    assert weights.shape == (nshifts, narrays)
    weights = num.reshape(weights, (nshifts*narrays))

    weights = weights.astype(dtype, copy=False)
    arrays = [arr.astype(dtype, copy=False) for arr in arrays]

    if impl == 'openmp':
        parstack_impl = parstack_ext.parstack
    elif impl == 'numpy':
        parstack_impl = parstack_numpy

    result, offset = parstack_impl(
        arrays, offsets, shifts, weights, method,
        lengthout, offsetout, result, nparallel)

    if method == 0:
        nsamps = result.size // nshifts
        result = result.reshape((nshifts, nsamps))

    return result, offset


def get_offset_and_length(arrays, offsets, shifts):
    narrays = offsets.size
    nshifts = shifts.size // narrays
    if shifts.ndim == 2:
        shifts = num.reshape(shifts, (nshifts*narrays))

    lengths = num.array([a.size for a in arrays], dtype=int)
    imin = offsets[0] + shifts[0]
    imax = imin + lengths[0]
    for iarray in range(len(arrays)):
        istarts = offsets[iarray] + shifts[iarray::narrays]
        iends = istarts + lengths[iarray]
        imin = min(imin, num.amin(istarts))
        imax = max(imax, num.amax(iends))

    return imin, imax - imin


def parstack_numpy(
        arrays,
        offsets,
        shifts,
        weights,
        method,
        lengthout,
        offsetout,
        result,
        nparallel):

    # nparallel is ignored here

    narrays = offsets.size

    lengths = num.array([a.size for a in arrays], dtype=int)
    if lengthout < 0:
        imin, nsamp = get_offset_and_length(arrays, offsets, shifts)
    else:
        nsamp = lengthout
        imin = offsetout

    nshifts = shifts.size // narrays
    result = num.zeros(nsamp*nshifts, dtype=float)

    for ishift in range(nshifts):
        for iarray in range(narrays):
            istart = offsets[iarray] + shifts[ishift*narrays + iarray]
            weight = weights[ishift*narrays + iarray]
            istart_r = ishift*nsamp + istart - imin
            jstart = max(0, imin - istart)
            jstop = max(jstart, min(nsamp - istart + imin, lengths[iarray]))
            result[istart_r+jstart:istart_r+jstop] += \
                arrays[iarray][jstart:jstop] * weight

    if method == 0:
        return result, imin
    elif method == 1:
        return num.amax(result.reshape((nshifts, nsamp)), axis=1), imin


[docs]def argmax(a, nparallel=1): ''' Same as numpys' argmax for 2 dimensional arrays along axis 0 but more memory efficient and twice as fast. ''' return parstack_ext.argmax(a, nparallel)