In [2]:
import sys
sys.path.append("..")
import splitwavepy as sw

import numpy as np
import matplotlib.pyplot as plt

In [13]:
x = np.random.rand(500)
y = np.random.rand(500)
degs = np.arange(-90,90,2)
slags = np.arange(0,50)
window = sw.Window(51,0)

In [51]:
def eigvalcov(data):
    """
    return sorted eigenvalues of covariance matrix
    lambda1 first, lambda2 second
    """
    return np.sort(np.linalg.eigvalsh(np.cov(data)))

def grideigval(x, y, degs, slags, window, **kwargs):
    """
    Grid search for splitting parameters applied to data.
    
    lags = 1-D array of sample shifts to search over, if None an attempt at finding sensible values is made
    degs = 1-D array of rotations to search over, if None an attempt at finding sensible values is made
    window = Window object (if None will guess an appropriate window)
    rcvcorr = receiver correction parameters in tuple (fast,lag) 
    srccorr = source correction parameters in tuple (fast,lag) 
    """
    
    # if 'pol' in kwargs:
    #     def silverchan(x,y):
    #         # measure energy on components
    #         ux, uy = rotate(ux,uy,-degs[0,ii]+pol)
    #         return transenergy(ux,uy)
    # else:
    #     def silverchan(x,y):
    #         return
                                 
    # grid of degs and lags to search over
    degs, lags = np.meshgrid(degs,slags)
    shape = degs.shape
    lam1 = np.zeros(shape)
    lam2 = np.zeros(shape)
#     cov = np.zeros(3,3,shape)
    
    # avoid using "dots" in loops for performance
    rotate = sw.core.rotate
    lag = sw.core.lag
    chop = sw.core.chop
    cov = np.cov
    vstack = np.vstack
    
#     covstack
    
    for ii in np.arange(shape[1]):
        tx, ty = rotate(x,y,degs[0,ii])
        for jj in np.arange(shape[0]):
            # remove splitting so use inverse operator (negative lag)
            ux, uy = lag(tx,ty,-lags[jj,ii])
            # chop to analysis window
            ux, uy = chop(ux,uy,window=window)
            # measure eigenvalues of covariance matrix
#             covstack = cov(vstack((ux,uy)))
#             lam2[jj,ii], lam1[jj,ii] = eigvalcov(np.vstack((ux,uy)))
            np.cov(np.vstack((ux,uy)))
#             pass
#     eigvals = np.linalg.eigvalsh(cov)
            
#     return degs,lags,lam1,lam2

In [52]:
%timeit grideigval(x,y,degs,slags,window)


1 loop, best of 3: 199 ms per loop

In [ ]: