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)
In [ ]: