def transform(self, X, y):
"""transform function"""
XMat = np.array(X)
yMat = np.array(y)
if XMat.shape[0] != yMat.shape[0]:
yMat = yMat.T
assert XMat.shape[0] == yMat.shape[0]
XMat -= XMat.mean(axis=0)
Sw, Sb = calc_Sw_Sb(XMat, yMat)
evals, evecs = eig(Sw, Sb)
np.ascontiguousarray(evals)
np.ascontiguousarray(evecs)
idx = np.argsort(evals)
idx = idx[::-1]
evecs = evecs[:, idx]
self.W = evecs[:, :self.n_components]
X_transformed = np.dot(XMat, self.W)
return X_transformed
python类eig()的实例源码
def eigenbasis(se, nb):
# generate number sector
ns1 = se.model.numbersector(nb)
# get the size of the basis
ns1size = ns1.basis.len # length of the number sector basis
# G1i = range(ns1size) # our Greens function?
# self energy
# sigma = self.sigma(nb, phi)
# Effective Hamiltonian
H1n = ns1.hamiltonian
# Complete diagonalization
E1, psi1r = linalg.eig(H1n.toarray(), left=False)
psi1l = np.conj(np.linalg.inv(psi1r)).T
# psi1l = np.conj(psi1r).T
# check for dark states (throw a warning if one shows up)
# if (nb > 0):
# Setup.check_for_dark_states(nb, E1)
return E1, psi1l, psi1r
def contruct_ellipse_parallel(pars):
Coor,cm,A_i,Vr,dims,dist,max_size,min_size,d=pars
dist_cm = coo_matrix(np.hstack([Coor[c].reshape(-1, 1) - cm[k]
for k, c in enumerate(['x', 'y', 'z'][:len(dims)])]))
Vr.append(dist_cm.T * spdiags(A_i.toarray().squeeze(),
0, d, d) * dist_cm / A_i.sum(axis=0))
if np.sum(np.isnan(Vr)) > 0:
raise Exception('You cannot pass empty (all zeros) components!')
D, V = eig(Vr[-1])
dkk = [np.min((max_size**2, np.max((min_size**2, dd.real)))) for dd in D]
# search indexes for each component
return np.sqrt(np.sum([(dist_cm * V[:, k])**2 / dkk[k] for k in range(len(dkk))], 0)) <= dist
#%% threshold_components
def get_gev_vector(target_psd_matrix, noise_psd_matrix):
"""
Returns the GEV beamforming vector.
:param target_psd_matrix: Target PSD matrix
with shape (bins, sensors, sensors)
:param noise_psd_matrix: Noise PSD matrix
with shape (bins, sensors, sensors)
:return: Set of beamforming vectors with shape (bins, sensors)
"""
bins, sensors, _ = target_psd_matrix.shape
beamforming_vector = np.empty((bins, sensors), dtype=np.complex)
for f in range(bins):
try:
eigenvals, eigenvecs = eigh(target_psd_matrix[f, :, :],
noise_psd_matrix[f, :, :])
except np.linalg.LinAlgError:
eigenvals, eigenvecs = eig(target_psd_matrix[f, :, :],
noise_psd_matrix[f, :, :])
beamforming_vector[f, :] = eigenvecs[:, np.argmax(eigenvals)]
return beamforming_vector
def get_gevd_vals_vecs(target_psd_matrix, noise_psd_matrix):
"""
Returns the eigenvalues and eigenvectors of GEVD
:param target_psd_matrix: Target PSD matrix
with shape (bins, sensors, sensors)
:param noise_psd_matrix: Noise PSD matrix
with shape (bins, sensors, sensors)
:return: Set of eigen values with shape (bins, sensors)
eigenvectors (bins, sensors, sensors)
"""
bins, sensors, _ = target_psd_matrix.shape
eigen_values = np.empty((bins, sensors), dtype=np.complex)
eigen_vectors = np.empty((bins, sensors, sensors), dtype=np.complex)
for f in range(bins):
try:
eigenvals, eigenvecs = eigh(target_psd_matrix[f, :, :],
noise_psd_matrix[f, :, :])
except np.linalg.LinAlgError:
eigenvals, eigenvecs = eig(target_psd_matrix[f, :, :],
noise_psd_matrix[f, :, :])
# values in increasing order
eigen_values[f,:] = eigenvals
eigen_vectors[f, :] = eigenvecs
return eigen_values, eigen_vectors
def get_gev_vector(target_psd_matrix, noise_psd_matrix):
"""
Returns the GEV beamforming vector.
:param target_psd_matrix: Target PSD matrix
with shape (bins, sensors, sensors)
:param noise_psd_matrix: Noise PSD matrix
with shape (bins, sensors, sensors)
:return: Set of beamforming vectors with shape (bins, sensors)
"""
bins, sensors, _ = target_psd_matrix.shape
beamforming_vector = np.empty((bins, sensors), dtype=np.complex)
for f in range(bins):
try:
eigenvals, eigenvecs = eigh(target_psd_matrix[f, :, :],
noise_psd_matrix[f, :, :])
except np.linalg.LinAlgError:
eigenvals, eigenvecs = eig(target_psd_matrix[f, :, :],
noise_psd_matrix[f, :, :])
beamforming_vector[f, :] = eigenvecs[:, np.argmax(eigenvals)]
return beamforming_vector
def get_dressed_info(H0):
# assign index of the dressed state according to the overall with bare state
w_c, v_c = la.eig(H0)
dressed_id=[]
for ii in range(len(v_c)):
index = np.argmax(np.abs(v_c[:, ii]))
if index not in dressed_id:
dressed_id.append(index)
else:
temp = (np.abs(v_c[:, ii])).tolist()
while index in dressed_id:
temp[index] = 0
index = np.argmax(temp)
dressed_id.append(index)
return w_c, v_c, dressed_id
def get_gev_vector(target_psd_matrix, noise_psd_matrix):
"""
Returns the GEV beamforming vector.
:param target_psd_matrix: Target PSD matrix
with shape (bins, sensors, sensors)
:param noise_psd_matrix: Noise PSD matrix
with shape (bins, sensors, sensors)
:return: Set of beamforming vectors with shape (bins, sensors)
"""
bins, sensors, _ = target_psd_matrix.shape
beamforming_vector = np.empty((bins, sensors), dtype=np.complex)
for f in range(bins):
try:
eigenvals, eigenvecs = eigh(target_psd_matrix[f, :, :],
noise_psd_matrix[f, :, :])
except np.linalg.LinAlgError:
eigenvals, eigenvecs = eig(target_psd_matrix[f, :, :],
noise_psd_matrix[f, :, :])
beamforming_vector[f, :] = eigenvecs[:, np.argmax(eigenvals)]
return beamforming_vector
def gen_A(self, V, n_classes, n_dims, return_S_b=False):
""" A = [B][inv(? ** .5)][Q^T] and assumes same number of data
in each class v. """
B = np.random.randint(-100, 100, (n_dims, n_dims)).astype(float)
big_V = np.matmul(V.T, V) # V is now a scatter matrix.
vals, vecs = eig(big_V)
A = B / np.sqrt(vals.real)
A = np.matmul(A, vecs.T)
D = np.matmul(np.matmul(vecs.T, big_V), vecs)
assert np.allclose(D, np.diag(vals))
if return_S_b is True:
S_b = 1 /n_classes * np.matmul(np.matmul(A, big_V), A.T)
x = np.matmul(A, V.T).T
S_b_empirical = 1 / n_classes * np.matmul(x.T, x)
assert np.allclose(S_b, S_b_empirical)
return A, S_b
else:
return A
def isNormal(A, method = 'definition'):
# use Schur inequality to determine whether it's normal
if method == 'Schur':
# initialize eigenValue
eigenValue = la.eig(A)[0]
if abs(np.sum(eigenValue**2) - la.norm(A, 'fro')**2) < 0.00001:
return True
else:
return False
# use definition
else:
if abs((A.conjugate().T.dot(A) - A.dot(A.conjugate().T)).all()) < 0.00001:
return True
else:
return False
def get_esystem(basis,traj_edges,test_set=None,delay=1):
"""
"""
if test_set is None:
test_set = basis
# Calculate Generator, Stiffness matrix
L = get_generator(basis,traj_edges,test_set=test_set,delay=delay,dt_eff=1)
# L = get_transop(basis,traj_edges,test_set=test_set,delay=delay)
S = get_stiffness_mat(basis,traj_edges,test_set=test_set,delay=delay)
# Calculate, sort eigensystem
evals, evecs_l, evecs_r = spl.eig(L,b=S,left=True,right=True,overwrite_a=False,overwrite_b=False)
idx = evals.argsort()[::-1]
evals = evals[idx]
evecs_l = evecs_l[:,idx]
evecs_r = evecs_r[:,idx]
# Expand eigenvectors into real space.
expanded_evecs_l = np.dot(test_set,evecs_l)
expanded_evecs_r = np.dot(basis,evecs_r)
return evals, expanded_evecs_l, expanded_evecs_r
def FLQTQEB(engine,app):
'''
This method calculates the Floquet quasi-energy bands.
'''
if app.path is None:
result=zeros((2,engine.nmatrix+1))
result[:,0]=array(xrange(2))
result[0,1:]=angle(eig(engine.evolution(ts=app.ts.mesh('t')))[0])/app.ts.volume('t')
result[1,1:]=result[0,1:]
else:
rank,mesh=app.path.rank(0),app.path.mesh(0)
result=zeros((rank,engine.nmatrix+1))
result[:,0]=mesh if mesh.ndim==1 else array(xrange(rank))
for i,paras in app.path('+'):
result[i,1:]=angle(eig(engine.evolution(ts=app.ts.mesh('t'),**paras))[0])/app.ts.volume('t')
name='%s_%s'%(engine,app.name)
if app.savedata: savetxt('%s/%s.dat'%(engine.dout,name),result)
if app.plot: app.figure('L',result,'%s/%s'%(engine.dout,name))
if app.returndata: return result
def transform(self, X, y):
"""transform function"""
XMat = np.array(X)
yMat = np.array(y)
if XMat.shape[0] != yMat.shape[0]:
yMat = yMat.T
assert XMat.shape[0] == yMat.shape[0]
XMat -= XMat.mean(axis=0)
Sw, Sb = calc_Sw_Sb(XMat, yMat)
if self.method == 'svd':
U, S, V = np.linalg.svd(Sw)
S = np.diag(S)
Sw_inversed = V * np.linalg.pinv(S) * U.T
A = Sw_inversed * Sb
elif self.method == 'auto':
A = np.linalg.pinv(Sw) * Sb
eigval, eigvec = np.linalg.eig(A)
eigval = eigval[0:self.n_components]
eigvec = eigvec[:, 0:self.n_components]
X_transformed = np.dot(XMat, eigvec)
self.W = eigvec
return X_transformed
def control_systems(request):
ct_sys, ref = request.param
Ac, Bc, Cc = ct_sys.data
Dc = np.zeros((Cc.shape[0], 1))
Q = np.eye(Ac.shape[0])
R = np.eye(Bc.shape[1] if len(Bc.shape) > 1 else 1)
Sc = linalg.solve_continuous_are(Ac, Bc.reshape(-1, 1), Q, R,)
Kc = linalg.solve(R, Bc.T @ Sc).reshape(1, -1)
ct_ctr = LTISystem(Kc)
evals = np.sort(np.abs(
linalg.eig(Ac, left=False, right=False, check_finite=False)
))
dT = 1/(2*evals[-1])
Tsim = (8/np.min(evals[~np.isclose(evals, 0)])
if np.sum(np.isclose(evals[np.nonzero(evals)], 0)) > 0
else 8
)
dt_data = signal.cont2discrete((Ac, Bc.reshape(-1, 1), Cc, Dc), dT)
Ad, Bd, Cd, Dd = dt_data[:-1]
Sd = linalg.solve_discrete_are(Ad, Bd.reshape(-1, 1), Q, R,)
Kd = linalg.solve(Bd.T @ Sd @ Bd + R, Bd.T @ Sd @ Ad)
dt_sys = LTISystem(Ad, Bd, dt=dT)
dt_sys.initial_condition = ct_sys.initial_condition
dt_ctr = LTISystem(Kd, dt=dT)
yield ct_sys, ct_ctr, dt_sys, dt_ctr, ref, Tsim
def eigenbasis(self, nb, phi=0):
"""
Calculates the generalized eigen-energies along with
the left and right eigen-basis.
Parameters
----------
nb : int
Number of bosons
phi : float
Phase factor for the relevant photonic state
"""
phi = 0 if self.local else phi
ckey = '{}-{}'.format(nb, phi)
if ckey not in self._cache['eigen']:
# generate number sector
ns1 = self.model.numbersector(nb)
# get the size of the basis
ns1size = ns1.basis.len # length of the number sector basis
# G1i = xrange(ns1size) # our Greens function?
# self energy
sigma = self.sigma(nb, phi)
# Effective Hamiltonian
H1n = ns1.hamiltonian + sigma
# Complete diagonalization
E1, psi1r = linalg.eig(H1n.toarray(), left=False)
psi1l = np.conj(np.linalg.inv(psi1r)).T
# check for dark states (throw a warning if one shows up)
# if (nb > 0):
# Setup.check_for_dark_states(nb, E1)
self._cache['eigen'][ckey] = (E1, psi1l, psi1r)
return self._cache['eigen'][ckey]
def fitEllipse (x, y):
x = x [:, np.newaxis]
y = y [:, np.newaxis]
D = np.hstack(( x*x , x*y , y*y , x, y, np.ones_like(x)))
S = np.dot (D.T ,D)
C = np.zeros ([6 , 6])
C[0, 2] = C[2, 0] = 2
C[1, 1] = -1
E ,V = eig( np.dot(inv(S),C ))
n = np.argmax(np.abs(E))
a = V[:, n]
if a[0] < 0 : a = -a
return a
def Aint(qL, qR, d):
""" Returns the Osher-Solomon jump matrix for A, in the dth direction
NB: eig function should be replaced with analytical function, if known
"""
ret = zeros(n, dtype=complex128)
?q = qR - qL
for i in range(N+1):
q = qL + nodes[i] * ?q
J = jacobian(q, d)
?, R = eig(J, overwrite_a=1, check_finite=0)
b = solve(R, ?q, check_finite=0)
ret += weights[i] * dot(R, abs(?)*b)
return ret.real
def PCA(self,data):
"""
returns: data transformed in 2 dims/columns + regenerated original data
pass in: data as 2D NumPy array
"""
dims_rescaled_data=2
m, n = data.shape
data -= data.mean(axis=0)
# calculate the covariance matrix
R = np.cov(data, rowvar=False)
# calculate eigenvectors & eigenvalues of the covariance matrix
# use 'eigh' rather than 'eig' since R is symmetric,
# the performance gain is substantial
evals, evecs = LA.eig(R)
# sort eigenvalue in decreasing order
idx = np.argsort(evals)[::-1]
evecs = evecs[:,idx]
# sort eigenvectors according to same index
evals = evals[idx]
# select the first n eigenvectors (n is desired dimension
# of rescaled data array, or dims_rescaled_data)
evecs = evecs[:, :dims_rescaled_data]
# carry out the transformation on the data using eigenvectors
# and return the re-scaled data, eigenvalues, and eigenvectors
return np.dot(evecs.T, data.T).T, evals, evecs
def plot_eig(model, maps, t, condition = 0, color = [0, 0, 1]):
import scipy.linalg as la
w = get_weights(model)
wrec = w[1]
bin_maps = maps.astype('int')
bin_t = bin_maps[condition, t, :].reshape(np.shape(bin_maps)[2], 1)
bin_mask = bin_t.dot(bin_t.T)
eva = la.eig(wrec * bin_mask)
plt.scatter(eva[0].real, eva[0].imag, 12, color)
plt.xlabel('real')
plt.ylabel('imaginary')
# Dave's visualization function
def compute_scores(self,Y):
if Y.shape[1]==1:
Y=self.get_Y_from_y(Y)
P,N=Y.shape
R=Y*Y.H/N
S,U=lg.eig(R)
eigenvalues=np.sort(np.real(S))[::-1]
scores=np.zeros(len(self.L_vect))
for index,L in enumerate(self.L_vect):
if self.corrected==1:
nb_free_parameters=L*(2*P-1-L) #see [WIL94]
else:
nb_free_parameters=L*(2*P-L) #see [WAX85]
num=np.prod(eigenvalues[L:])
den=np.mean(eigenvalues[L:])**(P-L)
penalty_term=penalty_term_IC(nb_free_parameters,N,method=self.method)
scores[index]=-2*N*np.log(num/den)+penalty_term
self.scores=scores
def hsvd(sys):
"""Compute Hankel singular values of a linear system.
Parameters
----------
sys : :data:`linear_system_like`
Linear system representation.
Returns
-------
``(len(sys),) np.array``
Hankel singular values.
See Also
--------
:class:`.Hankel`
:func:`.balanced_transformation`
References
----------
.. [#] Glover, Keith, and Jonathan R. Partington. "Bounds on the
achievable accuracy in model reduction." Modelling, robustness and
sensitivity reduction in control systems. Springer Berlin
Heidelberg, 1987. 95-118.
.. [#] https://www.mathworks.com/help/control/ref/hsvd.html
"""
sys = LinearSystem(sys)
R, O = control_gram(sys), observe_gram(sys)
# sort needed to be consistent across different versions
return np.sort(np.sqrt(abs(eig(np.dot(O, R))[0])))[::-1]
def isSimple(A):
'''
ask is matrix A is a simple matrix
'''
# check if A is a squre matrix
if A.shape[1] != A.shape[0]:
return False
eigenValues, eigenVectors = la.eig(A)
while eigenValues.shape[0] != 0:
#dictValues.update({eigenValues[0]: 1})
index = np.argwhere(abs(eigenValues - eigenValues[0]) < 0.00001)
algebraicMulti = len(index)
geometricMulti = eigenVectors[:, index].shape[1]
if algebraicMulti != geometricMulti:
return False
#dictValues.update({eigenValues[0]: len})
eigenValues = np.delete(eigenValues, index)
# stack another spaces of eigenvalue and eigenvector
return True
# compute the spectrum decomposition of matrix A
def spectrum_decomposition(A):
'''
compute the spectrum decomposition of matrix A
'''
# check if A is a simple matrix
if isSimple(A) != True:
raise Exception('non-simple matrix cannot be spectrum-decomposed')
eigenValues, eigenVectors = la.eig(A)
invVectors = la.inv(eigenVectors)
eigenVectors_row = eigenVectors.shape[0]
eigenVectors_column = eigenVectors.shape[1]
invVectors_row = invVectors.shape[0]
invVectors_column = invVectors.shape[1]
# an array of all distinct eigen values with their associated algebraic multiplicity
# keys: eigen values
# values: algebraic multiplicity
# {eigenValue: algebraicMulti}
dictValues = {}
while eigenValues.shape[0] != 0:
index = np.argwhere(abs(eigenValues - eigenValues[0]) < 0.00001)
spectrum = eigenVectors[:, index].reshape((eigenVectors_row, len(index))).dot(invVectors[index, :].reshape((len(index), invVectors_column)) )
dictValues.update({eigenValues[0]: spectrum})
eigenValues = np.delete(eigenValues, index)
eigenVectors = np.delete(eigenVectors, index, 1)
invVectors = np.delete(invVectors, index, 0)
return dictValues
def verSchur(A):
'''
a function that verifies Schur inequality
test whether the sum of all square eigenValue is less than or equal to Frobenius norm of the matrix
'''
# initialize eigenValue
eigenValue = la.eig(A)[0]
if np.sum(eigenValue**2) <= la.norm(A, 'fro')**2:
return True
else:
return False
# approximation of the upper bound of the absolute value of each eigenValue
def isSimple(A):
# check if A is a squre matrix
if A.shape[1] != A.shape[0]:
return False
eigenValues, eigenVectors = la.eig(A)
while (eigenValues.shape[0] != 0):
#dictValues.update({eigenValues[0]: 1})
index = np.argwhere(abs(eigenValues - eigenValues[0]) < 0.00001)
algebraicMulti = len(index)
geometricMulti = eigenVectors[:, index].shape[1]
if algebraicMulti != geometricMulti:
return False
#dictValues.update({eigenValues[0]: len})
eigenValues = np.delete(eigenValues, index)
# stack another spaces of eigenvalue and eigenvector
return True
def calcAvec(new, dQ, W, lambda_, active_set, M, positive):
# TODO: comment
r, c = np.nonzero(active_set)
# [r,c] = find(active_set);
Mm = -M.take(r, axis=0).take(r, axis=1)
Mm = (Mm + Mm.T) / 2
#% verify that there is no numerical instability
if len(Mm) > 1:
# print Mm.shape
eigMm, _ = scipy.linalg.eig(Mm)
eigMm = np.real(eigMm)
# check_here
else:
eigMm = Mm
if any(eigMm < 0):
np.min(eigMm)
#%error('The matrix Mm has negative eigenvalues')
flag = 1
b = np.sign(W)
if new >= 0:
b[new] = np.sign(dQ[new])
b = b[active_set == 1]
if len(Mm) > 1:
avec = np.linalg.solve(Mm, b)
else:
avec = b / Mm
if positive:
if new >= 0:
in_ = np.sum(active_set[:new])
if avec[in_] < 0:
# new;
#%error('new component of a is negative')
flag = 1
one_vec = np.ones(W.shape)
dQa = np.zeros(W.shape)
for j in range(len(r)):
dQa = dQa + np.expand_dims(avec[j] * M[:, r[j]], axis=1)
gamma_plus = (lambda_ - dQ) / (one_vec + dQa)
gamma_minus = (lambda_ + dQ) / (one_vec - dQa)
return avec, gamma_plus, gamma_minus
def eigensystem(matrix):
"""
Given an array-like 'matrix', returns:
- An array of eigenvalues
- An array of eigenvectors
- A rotation matrix that rotates the eigenbasis
into the original basis
Example:
mat = [[1,2,3],[2,4,5],[3,5,6]]
evals, evecs, rot = eigensystem(mat)
evals
array([ 11.34481428+0.j, -0.51572947+0.j, 0.17091519+0.j]
evecs
array([[-0.32798528, -0.59100905, -0.73697623],
[-0.73697623, -0.32798528, 0.59100905],
[ 0.59100905, -0.73697623, 0.32798528]])
rot
array([[-0.32798528, -0.73697623, 0.59100905],
[-0.59100905, -0.32798528, -0.73697623],
[-0.73697623, 0.59100905, 0.32798528]]))
This allows you to translate between original and eigenbases:
If [v1, v2, v3] are the components of a vector in eigenbasis e1, e2, e3
Then:
rot.dot([v1,v2,v3]) = [vx,vy,vz]
Will give the components in the original basis x, y, z
If [wx, wy, wz] are the components of a vector in original basis z, y, z
Then:
inv(rot).dot([wx,wy,wz]) = [w1,w2,w3]
Will give the components in the eigenbasis e1, e2, e3
inv(rot).dot(mat).dot(rot)
array([[evals[0], 0, 0]
[0, evals[1], 0]
[0, 0, evals[2]]])
Note: For symmetric input 'matrix', inv(rot) == evecs
"""
evals, emat = eig(matrix)
return evals, np.transpose(emat), emat
def solveCSE(en, rot, mu, R, VT):
n, m, oo = VT.shape
V = np.transpose(VT)
# find channels that are open, as defined by E' > 0
edash = en - np.diag(VT[:, :, -1])
openchann = edash > 0
nopen = edash[openchann].size
mx = matching_point(en, rot, V, R, mu)
if mx < oo-5:
out = leastsq(eigen, (en, ), args=(rot, mx, V, R, mu), xtol=0.01)
en = float(out[0])
# solve CSE according to Johnson renormalized Numerov method
WI = WImat(en, rot, V, R, mu)
RI = RImat(WI, mx)
wf = []
if nopen > 0:
oc = 0
for j, ed in enumerate(edash):
if ed > 0:
f = fmat(j, RI, WI, mx)
wf.append(wavefunction(WI, oc, f))
oc += 1
else:
f = fmat(0, RI, WI, mx)
wf.append(wavefunction(WI, nopen, f))
wf = np.array(wf)
wf = np.transpose(wf)
if nopen == 0:
wf = normalize(wf, R)
else:
K, AI, B = amplitude(wf, R, edash, mu) # shape = nopen x nopen
# K = BA-1 = U tan xi UT
eig, U = linalg.eig(K)
# form A^-1 U cos xi exp(i xi) UT
I = np.identity(nopen, dtype=complex)
xi = np.arctan(eig)*I
cosxi = np.cos(xi)*I
expxi = np.exp((0+1j)*xi)*I
expxiUT = expxi @ np.transpose(U)
cosxiexpxiUT = cosxi @ expxiUT
UcosxiexpxiUT = U @ cosxiexpxiUT
Norm = AI @ UcosxiexpxiUT # == (cu, su) complex
# complex wavefunction array oo x n x nopen
wf = wf @ Norm
return wf, en
def processing(self):
'''
Perform the maf transformation
'''
# covariance of original bands
sigma = numpy.ma.cov(self.variates_stack.T, allow_masked=True) # TODO: Jenkins can't run this lines with allow_masked =True, whyyyy??
# covariance for horizontal and vertical shifts
sigmadh = numpy.ma.cov(self.H.T, allow_masked=True)
sigmadv = numpy.ma.cov(self.V.T, allow_masked=True)
sigma = numpy.cov(self.variates_stack.T)
# covariance for horizontal and vertical shifts
sigmadh = numpy.cov(self.H.T)
sigmadv = numpy.cov(self.V.T)
# simple pooling of shifts
sigmad = 0.5 * (numpy.array(sigmadh) + numpy.array(sigmadv))
# evalues, vec1 = scipy.linalg.eig(sigmad, sigma)
evalues, vec1 = linalg.eig(sigmad, sigma)
# Sort eigen values from smallest to largest and apply this order to
# eigen vector matrix
sort_index = evalues.argsort()
evalues = evalues[sort_index]
vec1 = vec1[:, sort_index]
# autocorrelation
# ac= 1-0.5*vec1
HH = 1 / numpy.std(self.variates_stack, 0, ddof=1)
diagvariates = numpy.diag(HH)
invstderrmaf = numpy.diag((1 / numpy.sqrt(numpy.diag(vec1.T * sigma * vec1))))
HHH = numpy.zeros((self.number_of_maf_variates), float)
for b in range(self.number_of_maf_variates):
# logger.info("Calculating component %d of MAF transformation" % b)
HHH[b] = cmp(numpy.sum((diagvariates * sigma * vec1 * invstderrmaf)[b]), 0)
sgn = numpy.diag(HHH) # assure positiviy
v = numpy.dot(vec1, sgn)
N = numpy.shape(self.variates_stack)[0]
X = self.variates_stack - numpy.tile(numpy.mean(self.variates_stack, 0), (N, 1))
# scale v to give MAFs with unit variance
aux1 = numpy.dot(numpy.dot(v.T, sigma), v) # dispersion of MAFs
aux2 = 1 / numpy.sqrt(numpy.diag(aux1))
aux3 = numpy.tile(aux2.T, (self.number_of_maf_variates, 1))
v = v * aux3 # now dispersion is unit matrix
self.mafs = numpy.dot(X, v)
