Python numpy 模块，conj() 实例源码
我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用numpy.conj()。
def _ncc_c(x, y):
"""
>>> _ncc_c([1,2,3,4], [1,2,3,4])
array([ 0.13333333, 0.36666667, 0.66666667, 1. , 0.66666667,
0.36666667, 0.13333333])
>>> _ncc_c([1,1,1], [1,1,1])
array([ 0.33333333, 0.66666667, 1. , 0.66666667, 0.33333333])
>>> _ncc_c([1,2,3], [-1,-1,-1])
array([-0.15430335, -0.46291005, -0.9258201 , -0.77151675, -0.46291005])
"""
den = np.array(norm(x) * norm(y))
den[den == 0] = np.Inf
x_len = len(x)
fft_size = 1<<(2*x_len-1).bit_length()
cc = ifft(fft(x, fft_size) * np.conj(fft(y, fft_size)))
cc = np.concatenate((cc[-(x_len-1):], cc[:x_len]))
return np.real(cc) / den
def branch_power(self):
"""
"branch_power" computes the branch power
(passive sign convention)
:return:
* self.pb: real power for '.op' and complex power for '.ac'
"""
# check is branch voltages are available
if self.vb is None:
self.branch_voltage()
# check is branch current are available
if self.ib is None:
self.branch_current()
if self.analysis[0].lower() == '.ac':
self.pb = self.vb * np.conj(self.ib)
else:
self.pb = self.vb * self.ib
def mandelbrot(h, w, maxit):
"""
Returns an image of the Mandelbrot fractal of size (h,w).
"""
y, x = np.ogrid[-1.4:1.4:h * 1j, -2:0.8:w * 1j]
c = x + y * 1j
z = c
divtime = maxit + np.zeros(z.shape, dtype=int)
for i in range(maxit):
z = z ** 2 + c
diverge = z * np.conj(z) > 2 ** 2
div_now = diverge & (divtime == maxit)
divtime[div_now] = i + 100
z[diverge] = 2
logger.debug("Updating divtime")
recorder.record('divtime', divtime)
return divtime
def nufft_scale1(N, K, alpha, beta, Nmid):
'''
calculate image space scaling factor
'''
# import types
# if alpha is types.ComplexType:
alpha = numpy.real(alpha)
# print('complex alpha may not work, but I just let it as')
L = len(alpha) - 1
if L > 0:
sn = numpy.zeros((N, 1))
n = numpy.arange(0, N).reshape((N, 1), order='F')
i_gam_n_n0 = 1j * (2 * numpy.pi / K) * (n - Nmid) * beta
for l1 in range(-L, L + 1):
alf = alpha[abs(l1)]
if l1 < 0:
alf = numpy.conj(alf)
sn = sn + alf * numpy.exp(i_gam_n_n0 * l1)
else:
sn = numpy.dot(alpha, numpy.ones((N, 1), dtype=numpy.float32))
return sn
def nufft_T(N, J, K, alpha, beta):
'''
equation (29) and (26)Fessler's paper
create the overlapping matrix CSSC (diagonal dominent matrix)
of J points
and then find out the pseudo-inverse of CSSC '''
# import scipy.linalg
L = numpy.size(alpha) - 1
# print('L = ', L, 'J = ',J, 'a b', alpha,beta )
cssc = numpy.zeros((J, J))
[j1, j2] = numpy.mgrid[1:J + 1, 1:J + 1]
overlapping_mat = j2 - j1
for l1 in range(-L, L + 1):
for l2 in range(-L, L + 1):
alf1 = alpha[abs(l1)]
# if l1 < 0: alf1 = numpy.conj(alf1)
alf2 = alpha[abs(l2)]
# if l2 < 0: alf2 = numpy.conj(alf2)
tmp = overlapping_mat + beta * (l1 - l2)
tmp = dirichlet(1.0 * tmp / (1.0 * K / N))
cssc = cssc + alf1 * numpy.conj(alf2) * tmp
return mat_inv(cssc)
def precompute(self):
# CSR_W = cuda_cffi.cusparse.CSR.to_CSR(self.st['W_gpu'],diag_type=True)
# Dia_W_cpu = scipy.sparse.dia_matrix( (self.st['M'], self.st['M']),dtype=dtype)
# Dia_W_cpu = scipy.sparse.dia_matrix( ( self.st['W'], 0 ), shape=(self.st['M'], self.st['M']) )
# Dia_W_cpu = scipy.sparse.diags(self.st['W'], format="csr", dtype=dtype)
# CSR_W = cuda_cffi.cusparse.CSR.to_CSR(Dia_W_cpu)
self.st['pHp_gpu'] = self.CSRH.gemm(self.CSR)
self.st['pHp']=self.st['pHp_gpu'].get()
print('untrimmed',self.st['pHp'].nnz)
self.truncate_selfadjoint(1e-5)
print('trimmed', self.st['pHp'].nnz)
self.st['pHp_gpu'] = cuda_cffi.cusparse.CSR.to_CSR(self.st['pHp'])
# self.st['pHWp_gpu'] = self.CSR.conj().gemm(CSR_W,transA=cuda_cffi.cusparse.CUSPARSE_OPERATION_TRANSPOSE)
# self.st['pHWp_gpu'] = self.st['pHWp_gpu'].gemm(self.CSR, transA=cuda_cffi.cusparse.CUSPARSE_OPERATION_NON_TRANSPOSE)
def nufft_scale1(N, K, alpha, beta, Nmid):
'''
Calculate image space scaling factor
'''
# import types
# if alpha is types.ComplexType:
alpha = numpy.real(alpha)
# print('complex alpha may not work, but I just let it as')
L = len(alpha) - 1
if L > 0:
sn = numpy.zeros((N, 1))
n = numpy.arange(0, N).reshape((N, 1), order='F')
i_gam_n_n0 = 1j * (2 * numpy.pi / K) * (n - Nmid) * beta
for l1 in range(-L, L + 1):
alf = alpha[abs(l1)]
if l1 < 0:
alf = numpy.conj(alf)
sn = sn + alf * numpy.exp(i_gam_n_n0 * l1)
else:
sn = numpy.dot(alpha, numpy.ones((N, 1)))
return sn
def nufft_T(N, J, K, alpha, beta):
'''
The Equation (29) and (26) in Fessler and Sutton 2003.
Create the overlapping matrix CSSC (diagonal dominent matrix)
of J points and find out the pseudo-inverse of CSSC '''
# import scipy.linalg
L = numpy.size(alpha) - 1
# print('L = ', L, 'J = ',J, 'a b', alpha,beta )
cssc = numpy.zeros((J, J))
[j1, j2] = numpy.mgrid[1:J + 1, 1:J + 1]
overlapping_mat = j2 - j1
for l1 in range(-L, L + 1):
for l2 in range(-L, L + 1):
alf1 = alpha[abs(l1)]
# if l1 < 0: alf1 = numpy.conj(alf1)
alf2 = alpha[abs(l2)]
# if l2 < 0: alf2 = numpy.conj(alf2)
tmp = overlapping_mat + beta * (l1 - l2)
tmp = dirichlet(1.0 * tmp / (1.0 * K / N))
cssc = cssc + alf1 * alf2 * tmp
return mat_inv(cssc)
def _evaluate_windows(self, window_a, window_b):
"""
Calculate the FFT of both windows, correlate and transform back.
In order to decrease the error a mean subtraction is performed.
To compensate for the indexing during the FFT a FFT Shift is performed.
:param window_a: interrogation window
:param window_b: search window
:returns: correlation window
"""
fft_a = self._fa_fft(window_a - np.mean(window_a))
fft_b = self._fb_fft(window_b - np.mean(window_b))
fft_corr = fft_a*np.conj(fft_b)
inv_fft = self._ift_fft(fft_corr)
return np.fft.fftshift(inv_fft)
def _random_op(sites, ldim, hermitian=False, normalized=False, randstate=None,
dtype=np.complex_):
"""Returns a random operator of shape (ldim,ldim) * sites with local
dimension `ldim` living on `sites` sites in global form.
:param sites: Number of local sites
:param ldim: Local ldimension
:param hermitian: Return only the hermitian part (default False)
:param normalized: Normalize to Frobenius norm=1 (default False)
:param randstate: numpy.random.RandomState instance or None
:returns: numpy.ndarray of shape (ldim,ldim) * sites
>>> A = _random_op(3, 2); A.shape
(2, 2, 2, 2, 2, 2)
"""
op = _randfuncs[dtype]((ldim**sites,) * 2, randstate=randstate)
if hermitian:
op += np.transpose(op).conj()
if normalized:
op /= np.linalg.norm(op)
return op.reshape((ldim,) * 2 * sites)
def ccorr(a, b):
"""
Circular correlation of vectors
Computes the circular correlation of two vectors a and b via their
fast fourier transforms
a \ast b = \mathcal{F}^{-1}(\overline{\mathcal{F}(a)} \odot \mathcal{F}(b))
Parameter
---------
a: real valued array (shape N)
b: real valued array (shape N)
Returns
-------
c: real valued array (shape N), representing the circular
correlation of a and b
"""
return ifft(np.conj(fft(a)) * fft(b)).real
def fft_convolve(X,Y, inv = 0):
XF = np.fft.rfft2(X)
YF = np.fft.rfft2(Y)
# YF0 = np.copy(YF)
# YF.imag = 0
# XF.imag = 0
if inv == 1:
# plt.imshow(np.real(YF)); plt.colorbar(); plt.show()
YF = np.conj(YF)
SF = XF*YF
S = np.fft.irfft2(SF)
n1,n2 = np.shape(S)
S = np.roll(S,-n1/2+1,axis = 0)
S = np.roll(S,-n2/2+1,axis = 1)
return np.real(S)
def _invert(self):
""" Calculate the streamfunction given the potential vorticity.
The algorithm is:
1) Calculate wave potential vorticity
2) Invert for wave, pw, and vortex stremfunctions, pv.
3) Calculate geostrophic stremfunction, p = pv+pw.
"""
# the wavy PV
self.phich = self.fft(np.conj(self.phi))
self.phi2 = np.abs(self.phi)**2
self.jacph = self.jacobian_phic_phi()
self.gphi2h = -self.wv2*self.fft(self.phi2)
self.qwh = (0.5*self.gphi2h + 1j*self.jacph)/self.f/2.
self.qwh *= self.filtr
# invert for psi
self.ph = -self.wv2i*(self.qh-self.qwh)
# calculate in physical space
self.p = self.ifft(self.ph).real
def fst_delay_snd(fst, snd, samp_rate):
# Verify argument shape.
s1, s2 = fst.shape, snd.shape
if len(s1) != 1 or len(s2) != 1 or s1[0] != s2[0]:
raise Exception("Argument shape invalid, in 'fst_delay_snd' function")
length = s1[0]
half_len = int(length / 2)
Xfst = numpy.fft.fft(fst)
Xsnd_star = numpy.conj(numpy.fft.fft(snd))
Xall = numpy.zeros(length, dtype=numpy.complex64)
for i in range(0, length):
if Xsnd_star[i] == 0 or Xfst[i] == 0:
Xall[i] = 0
else:
Xall[i] = (Xsnd_star[i] * Xfst[i]) / abs(Xsnd_star[i]) / abs(Xfst[i])
R = numpy.fft.ifft(Xall)
max_pos = numpy.argmax(R)
if max_pos > half_len:
return -(length - 1 - max_pos) / samp_rate
else:
return max_pos / samp_rate
def fst_delay_snd(fst, snd, samp_rate):
# Verify argument shape.
s1, s2 = fst.shape, snd.shape
if len(s1) != 1 or len(s2) != 1 or s1[0] != s2[0]:
raise Exception("Argument shape invalid, in 'fst_delay_snd' function")
length = s1[0]
half_len = int(length / 2)
Xfst = numpy.fft.fft(fst)
Xsnd = numpy.fft.fft(snd)
Xsnd_star = numpy.conj(Xsnd)
Xall = numpy.zeros(length, dtype=numpy.complex64)
for i in range(0, length):
Xall[i] = (Xsnd_star[i] * Xfst[i]) / abs(Xsnd_star[i]) / abs(Xfst[i])
R = numpy.fft.ifft(Xall)
max_pos = numpy.argmax(R)
if max_pos > half_len:
delta = -(length - 1 - max_pos) / samp_rate
else:
delta = max_pos / samp_rate
return delta, R, Xfst, Xsnd
def test_probabilities():
p = 1
n_qubits = 2
# known set of angles for barbell
angles = [1.96348709, 4.71241069]
wf = np.array([-1.17642098e-05 - 1j*7.67538040e-06,
-7.67563580e-06 - 1j*7.07106781e-01,
-7.67563580e-06 - 1j*7.07106781e-01,
-1.17642098e-05 - 1j*7.67538040e-06])
fakeQVM = Mock(spec=qvm_module.QVMConnection())
fakeQVM.wavefunction = Mock(return_value=(Wavefunction(wf)))
inst = QAOA(fakeQVM, n_qubits, steps=p,
rand_seed=42)
true_probs = np.zeros_like(wf)
for xx in range(wf.shape[0]):
true_probs[xx] = np.conj(wf[xx]) * wf[xx]
probs = inst.probabilities(angles)
assert isinstance(probs, np.ndarray)
prob_true = np.zeros((2**inst.n_qubits, 1))
prob_true[1] = 0.5
prob_true[2] = 0.5
assert np.isclose(probs, prob_true).all()
def probabilities(self, angles):
"""
Computes the probability of each state given a particular set of angles.
:param angles: [list] A concatenated list of angles [betas]+[gammas]
:return: [list] The probabilities of each outcome given those angles.
"""
if isinstance(angles, list):
angles = np.array(angles)
assert angles.shape[0] == 2 * self.steps, "angles must be 2 * steps"
param_prog = self.get_parameterized_program()
prog = param_prog(angles)
wf = self.qvm.wavefunction(prog)
wf = wf.amplitudes.reshape((-1, 1))
probs = np.zeros_like(wf)
for xx in range(2 ** self.n_qubits):
probs[xx] = np.conj(wf[xx]) * wf[xx]
return probs
def correlate_periodic(a, v=None):
"""Cross-correlation of two 1-dimensional periodic sequences.
a and v must be sequences with the same length. If v is not specified, it is
assumed to be the same as a (i.e. the function computes auto-correlation).
:param a: input sequence #1
:param v: input sequence #2
:returns: discrete periodic cross-correlation of a and v
"""
a_fft = _np.fft.fft(_np.asarray(a))
if v is None:
v_cfft = a_fft.conj()
else:
v_cfft = _np.fft.fft(_np.asarray(v)).conj()
x = _np.fft.ifft(a_fft * v_cfft)
if _np.isrealobj(a) and (v is None or _np.isrealobj(v)):
x = x.real
return x
def eval_gradf(self):
""" Compute gradient in Fourier domain """
# Compute X D - S
self.Ryf = self.eval_Rf(self.Yf)
# Map to spatial domain to multiply by mask
Ry = sl.irfftn(self.Ryf, self.cri.Nv, self.cri.axisN)
# Multiply by mask
WRy = (self.W**2) * Ry
# Map back to frequency domain
WRyf = sl.rfftn(WRy, self.cri.Nv, self.cri.axisN)
gradf = sl.inner(np.conj(self.Zf), WRyf, axis=self.cri.axisK)
# Multiple channel signal, single channel dictionary
if self.cri.C > 1 and self.cri.Cd == 1:
gradf = np.sum(gradf, axis=self.cri.axisC, keepdims=True)
return gradf
def setcoef(self, Z):
"""Set coefficient array."""
# If the dictionary has a single channel but the input (and
# therefore also the coefficient map array) has multiple
# channels, the channel index and multiple image index have
# the same behaviour in the dictionary update equation: the
# simplest way to handle this is to just reshape so that the
# channels also appear on the multiple image index.
if self.cri.Cd == 1 and self.cri.C > 1:
Z = Z.reshape(self.cri.Nv + (1,) + (self.cri.Cx*self.cri.K,) +
(self.cri.M,))
self.Z = np.asarray(Z, dtype=self.dtype)
self.Zf = sl.rfftn(self.Z, self.cri.Nv, self.cri.axisN)
# Compute X^H S
self.ZSf = sl.inner(np.conj(self.Zf), self.Sf, self.cri.axisK)
def setcoef(self, Z):
"""Set coefficient array."""
# If the dictionary has a single channel but the input (and
# therefore also the coefficient map array) has multiple
# channels, the channel index and multiple image index have
# the same behaviour in the dictionary update equation: the
# simplest way to handle this is to just reshape so that the
# channels also appear on the multiple image index.
if self.cri.Cd == 1 and self.cri.C > 1:
Z = Z.reshape(self.cri.Nv + (1,) + (self.cri.Cx*self.cri.K,) +
(self.cri.M,))
self.Z = np.asarray(Z, dtype=self.dtype)
self.Zf = sl.rfftn(self.Z, self.cri.Nv, self.cri.axisN)
# Compute X^H S
self.ZSf = np.conj(self.Zf) * self.Sf
def xstep(self):
r"""Minimise Augmented Lagrangian with respect to
:math:`\mathbf{x}`.
"""
self.cgit = None
self.YU[:] = self.Y - self.U
self.block_sep0(self.YU)[:] += self.S
YUf = sl.rfftn(self.YU, None, self.cri.axisN)
b = sl.inner(np.conj(self.Zf), self.block_sep0(YUf),
axis=self.cri.axisK) + self.block_sep1(YUf)
self.Xf[:], cgit = sl.solvemdbi_cg(self.Zf, 1.0, b,
self.cri.axisM, self.cri.axisK,
self.opt['CG', 'StopTol'],
self.opt['CG', 'MaxIter'], self.Xf)
self.cgit = cgit
self.X = sl.irfftn(self.Xf, self.cri.Nv, self.cri.axisN)
self.xstep_check(b)
def xstep(self):
"""The xstep of the baseline consensus class from which this
class is derived is re-used to implement the xstep of the
modified algorithm by replacing ``self.ZSf``, which is constant
in the baseline algorithm, with a quantity derived from the
additional variables ``self.Y1`` and ``self.U1``. It is also
necessary to set the penalty parameter to unity for the duration
of the x step.
"""
self.YU1[:] = self.Y1 - self.U1
self.ZSf = np.conj(self.Zf) * (self.Sf + sl.rfftn(
self.YU1, None, self.cri.axisN))
rho = self.rho
self.rho = 1.0
super(ConvCnstrMODMaskDcpl_Consensus, self).xstep()
self.rho = rho
def xstep(self):
r"""Minimise Augmented Lagrangian with respect to
:math:`\mathbf{x}`.
"""
b = self.AHSf + self.rho*np.sum(np.conj(self.Gf)*
sl.rfftn(self.Y-self.U, axes=self.axes),
axis=self.Y.ndim-1)
self.Xf = b / (self.AHAf + self.rho*self.GHGf)
self.X = sl.irfftn(self.Xf, None, axes=self.axes)
if self.opt['LinSolveCheck']:
ax = (self.AHAf + self.rho*self.GHGf)*self.Xf
self.xrrs = sl.rrs(ax, b)
else:
self.xrrs = None
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 _find_start_symbol(iq_data):
'''
Correlate to find symbol boundaries
'''
corr_length = 2*(SYMBOL_LENGTH)
corr = np.empty(corr_length)
for k in range(corr_length):
leading = iq_data[k:k+GUARD_LENGTH]
trailing = iq_data[k+USEFUL_LENGTH:k+SYMBOL_LENGTH]
corr[k] = np.abs(np.dot(leading, np.conj(trailing)))
first_symbol = np.argmax(corr)%(SYMBOL_LENGTH)
return first_symbol
def get_power_spectral_density_matrix(observation, mask=None):
"""
Calculates the weighted power spectral density matrix.
This does not yet work with more than one target mask.
:param observation: Complex observations with shape (bins, sensors, frames)
:param mask: Masks with shape (bins, frames) or (bins, 1, frames)
:return: PSD matrix with shape (bins, sensors, sensors)
"""
bins, sensors, frames = observation.shape
if mask is None:
mask = np.ones((bins, frames))
if mask.ndim == 2:
mask = mask[:, np.newaxis, :]
normalization = np.maximum(np.sum(mask, axis=-1, keepdims=True), 1e-6)
psd = np.einsum('...dt,...et->...de', mask * observation,
observation.conj())
psd /= normalization
return psd
def get_mvdr_vector(atf_vector, noise_psd_matrix):
"""
Returns the MVDR beamforming vector.
:param atf_vector: Acoustic transfer function vector
with shape (..., bins, sensors)
:param noise_psd_matrix: Noise PSD matrix
with shape (bins, sensors, sensors)
:return: Set of beamforming vectors with shape (..., bins, sensors)
"""
while atf_vector.ndim > noise_psd_matrix.ndim - 1:
noise_psd_matrix = np.expand_dims(noise_psd_matrix, axis=0)
# Make sure matrix is hermitian
noise_psd_matrix = 0.5 * (
noise_psd_matrix + np.conj(noise_psd_matrix.swapaxes(-1, -2)))
numerator = solve(noise_psd_matrix, atf_vector)
denominator = np.einsum('...d,...d->...', atf_vector.conj(), numerator)
beamforming_vector = numerator / np.expand_dims(denominator, axis=-1)
return beamforming_vector
def apply_sdw_mwf(mix, target_psd_matrix, noise_psd_matrix, mu=1, corr=None):
"""
Apply speech distortion weighted MWF: h = Tpsd * e1 / (Tpsd + mu*Npsd)
:param mix: the signal complex FFT
:param target_psd_matrix (bins, sensors, sensors)
:param noise_psd_matrix
:param mu: the lagrange factor
:return
"""
bins, sensors, frames = mix.shape
ref_vector = np.zeros((sensors,1), dtype=np.float)
if corr is None:
ref_ch = 0
else: # choose the channel with highest correlation with the others
corr=corr.tolist()
while len(corr) > sensors:
corr.remove(np.min(corr))
ref_ch=np.argmax(corr)
ref_vector[ref_ch,0]=1
mwf_vector = solve(target_psd_matrix + mu*noise_psd_matrix, target_psd_matrix[:,:,ref_ch])
return np.einsum('...a,...at->...t', mwf_vector.conj(), mix)
def test_blas_cgemm(backend, m, n, k, alpha, beta, forward):
b = backend()
y = indigo.util.rand64c(m,n)
M = indigo.util.rand64c(m,k)
x = indigo.util.rand64c(k,n)
if not forward:
x, y = y, x
M_exp = np.conj(M.T)
else:
M_exp = M
y_exp = alpha * M_exp.dot(x) + beta * y
y_d = b.copy_array(y)
M_d = b.copy_array(M)
x_d = b.copy_array(x)
b.cgemm(y_d, M_d, x_d, alpha, beta, forward=forward)
y_act = y_d.to_host()
np.testing.assert_allclose(y_exp, y_act, atol=1e-3)
def _lagged_coherence_1freq(x, f, Fs, N_cycles=3, f_step=1):
"""Calculate lagged coherence of x at frequency f using the hanning-taper FFT method"""
Nsamp = int(np.ceil(N_cycles*Fs / f))
# For each N-cycle chunk, calculate phase
chunks = _nonoverlapping_chunks(x,Nsamp)
C = len(chunks)
hann_window = signal.hanning(Nsamp)
fourier_f = np.fft.fftfreq(Nsamp,1/float(Fs))
fourier_f_idx = _arg_closest_value(fourier_f,f)
fourier_coefsoi = np.zeros(C,dtype=complex)
for i2, c in enumerate(chunks):
fourier_coef = np.fft.fft(c*hann_window)
fourier_coefsoi[i2] = fourier_coef[fourier_f_idx]
lcs_num = 0
for i2 in range(C-1):
lcs_num += fourier_coefsoi[i2]*np.conj(fourier_coefsoi[i2+1])
lcs_denom = np.sqrt(np.sum(np.abs(fourier_coefsoi[:-1])**2)*np.sum(np.abs(fourier_coefsoi[1:])**2))
return np.abs(lcs_num/lcs_denom)
def _create_rotational_weights_for_elements(self, kpoint, transformation_matrix, vectors):
"""
Parameters
----------
kpoint : 1d array
Reciprocal space point in fractional coordinates for PC.
vectors : (..., natoms_p * ndims, nbands) array
Vectors for SC after translational projection.
"""
projected_vectors = self._rotational_projector.project_vectors(
vectors, kpoint, transformation_matrix)
nirreps, natoms_p, nelms, tmp, nbands = projected_vectors.shape
shape = (nirreps, natoms_p, nelms, natoms_p, nelms, nbands)
weights = np.zeros(shape, dtype=complex)
for i in range(nirreps):
for j in range(nbands):
weights[i, ..., j] = np.inner(
np.conj(projected_vectors[i, ..., j]), projected_vectors[i, ..., j])
return weights, projected_vectors
def get_power_spectral_density_matrix(observation, mask=None):
"""
Calculates the weighted power spectral density matrix.
This does not yet work with more than one target mask.
:param observation: Complex observations with shape (bins, sensors, frames)
:param mask: Masks with shape (bins, frames) or (bins, 1, frames)
:return: PSD matrix with shape (bins, sensors, sensors)
"""
bins, sensors, frames = observation.shape
if mask is None:
mask = np.ones((bins, frames))
if mask.ndim == 2:
mask = mask[:, np.newaxis, :]
normalization = np.maximum(np.sum(mask, axis=-1, keepdims=True), 1e-6)
psd = np.einsum('...dt,...et->...de', mask * observation,
observation.conj())
psd /= normalization
return psd
def get_mvdr_vector(atf_vector, noise_psd_matrix):
"""
Returns the MVDR beamforming vector.
:param atf_vector: Acoustic transfer function vector
with shape (..., bins, sensors)
:param noise_psd_matrix: Noise PSD matrix
with shape (bins, sensors, sensors)
:return: Set of beamforming vectors with shape (..., bins, sensors)
"""
while atf_vector.ndim > noise_psd_matrix.ndim - 1:
noise_psd_matrix = np.expand_dims(noise_psd_matrix, axis=0)
# Make sure matrix is hermitian
noise_psd_matrix = 0.5 * (
noise_psd_matrix + np.conj(noise_psd_matrix.swapaxes(-1, -2)))
numerator = solve(noise_psd_matrix, atf_vector)
denominator = np.einsum('...d,...d->...', atf_vector.conj(), numerator)
beamforming_vector = numerator / np.expand_dims(denominator, axis=-1)
return beamforming_vector
def autocorrelation(self):
"Autocorrelation as a function of time"
if self.__autocorrelation is not None:
return self.__autocorrelationTimeSeries, self.__autocorrelation
negT = -np.flipud(self.timeSeries[1:])
autocorrelationTime = np.hstack((negT, self.timeSeries))
self.__autocorrelationTimeSeries = autocorrelationTime
initialWF = self[0]
ACF = []
for WF in self:
ACF.append(WF.overlap(initialWF))
ACF = np.array(ACF)
negACF = np.conj(np.flipud(ACF[1:]))
totalACF = np.hstack((negACF, ACF))
self.__autocorrelation = totalACF
return self.__autocorrelationTimeSeries, self.__autocorrelation
def swap_Nq(fft_y, fu, fft_x, N):
f = fu[:, 0].copy()
fft_x[0] = f[0].real
fft_x[1:N//2] = 0.5*(f[1:N//2] + np.conj(f[:N//2:-1]))
fft_x[N//2] = f[N//2].real
fu[:N//2+1, 0] = fft_x[:N//2+1]
fu[N//2+1:, 0] = np.conj(fft_x[(N//2-1):0:-1])
fft_y[0] = f[0].imag
fft_y[1:N//2] = -0.5*1j*(f[1:N//2] - np.conj(f[:N//2:-1]))
fft_y[N//2] = f[N//2].imag
fft_y[N//2+1:] = np.conj(fft_y[(N//2-1):0:-1])
return fft_y
def nufft_alpha_kb_fit(N, J, K):
'''
find out parameters alpha and beta
of scaling factor st['sn']
Note, when J = 1 , alpha is hardwired as [1,0,0...]
(uniform scaling factor)
'''
beta = 1
Nmid = (N - 1.0) / 2.0
if N > 40:
L = 13
else:
L = numpy.ceil(N / 3)
nlist = numpy.arange(0, N) * 1.0 - Nmid
(kb_a, kb_m) = kaiser_bessel('string', J, 'best', 0, K / N)
if J > 1:
sn_kaiser = 1 / kaiser_bessel_ft(nlist / K, J, kb_a, kb_m, 1.0)
elif J == 1: # The case when samples are on regular grids
sn_kaiser = numpy.ones((1, N), dtype=dtype)
gam = 2 * numpy.pi / K
X_ant = beta * gam * nlist.reshape((N, 1), order='F')
X_post = numpy.arange(0, L + 1)
X_post = X_post.reshape((1, L + 1), order='F')
X = numpy.dot(X_ant, X_post) # [N,L]
X = numpy.cos(X)
sn_kaiser = sn_kaiser.reshape((N, 1), order='F').conj()
X = numpy.array(X, dtype=dtype)
sn_kaiser = numpy.array(sn_kaiser, dtype=dtype)
coef = numpy.linalg.lstsq(numpy.nan_to_num(X), numpy.nan_to_num(sn_kaiser))[0] # (X \ sn_kaiser.H);
alphas = coef
if J > 1:
alphas[0] = alphas[0]
alphas[1:] = alphas[1:] / 2.0
elif J == 1: # cases on grids
alphas[0] = 1.0
alphas[1:] = 0.0
alphas = numpy.real(alphas)
return (alphas, beta)
def iterate_l1(L, alpha, arg, beta, K, N, rr):
oversample_ratio = (1.0 * K / N)
for l1 in range(-L, L + 1):
alf = alpha[abs(l1)] * 1.0
if l1 < 0:
alf = numpy.conj(alf)
# r1 = numpy.sinc(1.0*(arg+1.0*l1*beta)/(1.0*K/N))
input_array = (arg + 1.0 * l1 * beta) / oversample_ratio
r1 = dirichlet(input_array.astype(numpy.float32))
rr = iterate_sum(rr, alf, r1)
return rr
def nufft_r(om, N, J, K, alpha, beta):
'''
equation (30) of Fessler's paper
'''
def iterate_sum(rr, alf, r1):
rr = rr + alf * r1
return rr
def iterate_l1(L, alpha, arg, beta, K, N, rr):
oversample_ratio = (1.0 * K / N)
import time
t0=time.time()
for l1 in range(-L, L + 1):
alf = alpha[abs(l1)] * 1.0
# if l1 < 0:
# alf = numpy.conj(alf)
# r1 = numpy.sinc(1.0*(arg+1.0*l1*beta)/(1.0*K/N))
input_array = (arg + 1.0 * l1 * beta) / oversample_ratio
r1 = dirichlet(input_array)
rr = iterate_sum(rr, alf, r1)
return rr
M = numpy.size(om) # 1D size
gam = 2.0 * numpy.pi / (K * 1.0)
nufft_offset0 = nufft_offset(om, J, K) # om/gam - nufft_offset , [M,1]
dk = 1.0 * om / gam - nufft_offset0 # om/gam - nufft_offset , [M,1]
arg = outer_sum(-numpy.arange(1, J + 1) * 1.0, dk)
L = numpy.size(alpha) - 1
# print('alpha',alpha)
rr = numpy.zeros((J, M), dtype=numpy.float32)
rr = iterate_l1(L, alpha, arg, beta, K, N, rr)
return (rr, arg)
def series_nfft(series,
oversample=4):
"""
note that output period units are [days] (so is frequency)
"""
M = len(series)
if not is_power_of_2(M):
raise ValueError('series length {} is not a power of 2'.format(len(series)))
N = M
if N % 2 == 1:
# number of frequency samples must be even
N += 1
N *= oversample
# re-grid time the interval [-1/2, 1/2) as required by nfft
time = series.index.astype(NP.int64) / 1e9
time -= time[0]
b = -0.5
a = (M - 1) / (M * time[-1])
x = a * time + b
# setup for nfft computation
plan = NFFT(N, M)
plan.x = x
plan.f = series.values
plan.precompute()
# compute nfft (note that conjugation is necessary because of the
# difference in transform sign convention)
x_nfft = NP.conj(plan.adjoint())
# calculate frequencies and periods
dt = ((series.index[-1] - series.index[0]) / M).total_seconds() / DAYS_TO_SECONDS
f_range = NP.fft.fftshift(NP.fft.fftfreq(N, dt))
T_range = 1 / f_range
return x_nfft, f_range, T_range
def spatFT(data, position_grid, order_max=10, spherical_harmonic_bases=None):
''' Spatial Fourier Transform
Parameters
----------
data : array_like
Data to be transformed, with signals in rows and frequency bins in columns
order_max : int, optional
Maximum transform order (Default: 10)
position_grid : array_like or io.SphericalGrid
Azimuths/Colatitudes/Gridweights of spatial sampling points
Returns
-------
Pnm : array_like
Spatial Fourier Coefficients with nm coeffs in rows and FFT bins in columns
'''
data = _np.atleast_2d(data)
number_of_signals, FFTLength = data.shape
position_grid = SphericalGrid(*position_grid)
# Re-generate spherical harmonic bases if they were not provided or their order is too small
if (spherical_harmonic_bases is None or
spherical_harmonic_bases.shape[0] < number_of_signals or
spherical_harmonic_bases.shape[1] < (order_max + 1) ** 2):
spherical_harmonic_bases = sph_harm_all(order_max, position_grid.azimuth, position_grid.colatitude)
spherical_harmonic_bases = (_np.conj(spherical_harmonic_bases).T * (4 * _np.pi * position_grid.weight))
return _np.inner(spherical_harmonic_bases, data.T)
def _random_state(sites, ldim, randstate=None):
"""Returns a random positive semidefinite operator of shape (ldim, ldim) *
sites normalized to Tr rho = 1, i.e. a mixed state with local dimension
`ldim` living on `sites` sites. Note that the returned state is positive
semidefinite only when interpreted in global form (see
:func:`tools.global_to_local`)
:param sites: Number of local sites
:param ldim: Local ldimension
:param randstate: numpy.random.RandomState instance or None
:returns: numpy.ndarray of shape (ldim, ldim) * sites
>>> from numpy.linalg import eigvalsh
>>> rho = _random_state(3, 2).reshape((2**3, 2**3))
>>> all(eigvalsh(rho) >= 0)
True
>>> np.abs(np.trace(rho) - 1) < 1e-6
True
"""
shape = (ldim**sites, ldim**sites)
mat = _zrandn(shape, randstate=randstate)
rho = np.conj(mat.T).dot(mat)
rho /= np.trace(rho)
return rho.reshape((ldim,) * 2 * sites)
####################################
# Factory functions for MPArrays #
####################################
def random_mpo(sites, ldim, rank, randstate=None, hermitian=False,
normalized=True, force_rank=False):
"""Returns an hermitian MPO with randomly choosen local tensors
:param sites: Number of sites
:param ldim: Local dimension
:param rank: Rank
:param randstate: numpy.random.RandomState instance or None
:param hermitian: Is the operator supposed to be hermitian
:param normalized: Operator should have unit norm
:param force_rank: If True, the rank is exaclty `rank`.
Otherwise, it might be reduced if we reach the maximum sensible rank.
:returns: randomly choosen matrix product operator
>>> mpo = random_mpo(4, 2, 10, force_rank=True)
>>> mpo.ranks, mpo.shape
((10, 10, 10), ((2, 2), (2, 2), (2, 2), (2, 2)))
>>> mpo.canonical_form
(0, 4)
"""
mpo = random_mpa(sites, (ldim,) * 2, rank, randstate=randstate,
force_rank=force_rank, dtype=np.complex_)
if hermitian:
# make mpa Herimitan in place, without increasing rank:
ltens = (l + l.swapaxes(1, 2).conj() for l in mpo.lt)
mpo = mp.MPArray(ltens)
if normalized:
# we do this with a copy to ensure the returned state is not
# normalized
mpo /= mp.norm(mpo.copy())
return mpo
def random_mpdo(sites, ldim, rank, randstate=np.random):
"""Returns a randomly choosen matrix product density operator (i.e.
positive semidefinite matrix product operator with trace 1).
:param sites: Number of sites
:param ldim: Local dimension
:param rank: Rank
:param randstate: numpy.random.RandomState instance
:returns: randomly choosen classicaly correlated matrix product density op.
>>> rho = random_mpdo(4, 2, 4)
>>> rho.ranks, rho.shape
((4, 4, 4), ((2, 2), (2, 2), (2, 2), (2, 2)))
>>> rho.canonical_form
(0, 4)
"""
# generate density matrix as a mixture of `rank` pure product states
psis = [random_mps(sites, ldim, 1, randstate=randstate) for _ in range(rank)]
weights = (lambda x: x / np.sum(x))(randstate.rand(rank))
rho = mp.sumup(mpsmpo.mps_to_mpo(psi) * weight
for weight, psi in zip(weights, psis))
# Scramble the local tensors
for n, rank in enumerate(rho.ranks):
unitary = _unitary_haar(rank, randstate)
rho.lt[n] = matdot(rho.lt[n], unitary)
rho.lt[n + 1] = matdot(np.transpose(unitary).conj(), rho.lt[n + 1])
rho /= mp.trace(rho)
return rho
def test_conjugations(nr_sites, local_dim, _, rgen, dtype):
op = factory._random_op(nr_sites, local_dim, randstate=rgen, dtype=dtype)
mpo = mp.MPArray.from_array(op, 2)
assert_array_almost_equal(np.conj(op), mpo.conj().to_array())
assert mpo.conj().dtype == dtype
mpo.canonicalize()
mpo_c = mpo.conj()
assert_correct_normalization(mpo_c)
def test_norm(nr_sites, local_dim, rank, dtype, rgen):
mp_psi = factory.random_mpa(nr_sites, local_dim, rank, randstate=rgen,
dtype=dtype)
psi = mp_psi.to_array()
assert_almost_equal(mp.inner(mp_psi, mp_psi), mp.norm(mp_psi)**2)
assert_almost_equal(np.sum(psi.conj() * psi), mp.norm(mp_psi)**2)
def SLshearrecadjoint2D(X, shearletSystem):
"""
Adjoint of (pseudo-)inverse of 2D data.
Note that this is also the (pseudo-)inverse of the adjoint.
Usage:
coeffs = SLshearrecadjoint2D(X, shearletSystem)
Input:
X : 2D data.
shearletSystem: Structure containing a shearlet system. This
should be the same system as the one
previously used for decomposition.
Output:
coeffs: X x Y x N array of shearlet coefficients.
"""
# skipping useGPU stuff...
# STUFF
Xfreq = np.divide(fftlib.fftshift(fftlib.fft2(fftlib.ifftshift(X))), shearletSystem["dualFrameWeights"])
coeffs = np.zeros(shearletSystem["shearlets"].shape, dtype=complex)
for j in range(shearletSystem["nShearlets"]):
coeffs[:,:,j] = fftlib.fftshift(fftlib.ifft2(fftlib.ifftshift(Xfreq*np.conj(shearletSystem["shearlets"][:,:,j]))))
return np.real(coeffs).astype(X.dtype)
#
##############################################################################
def test_simple_conjugate(self):
ref = np.conj(np.sqrt(np.complex(1, 1)))
def f(z):
return np.sqrt(np.conj(z))
yield check_complex_value, f, 1, 1, ref.real, ref.imag, False
#def test_branch_cut(self):
# _check_branch_cut(f, -1, 0, 1, -1)
def test_cabs_inf_nan(self):
x, y = [], []
# cabs(+-nan + nani) returns nan
x.append(np.nan)
y.append(np.nan)
yield check_real_value, np.abs, np.nan, np.nan, np.nan
x.append(np.nan)
y.append(-np.nan)
yield check_real_value, np.abs, -np.nan, np.nan, np.nan
# According to C99 standard, if exactly one of the real/part is inf and
# the other nan, then cabs should return inf
x.append(np.inf)
y.append(np.nan)
yield check_real_value, np.abs, np.inf, np.nan, np.inf
x.append(-np.inf)
y.append(np.nan)
yield check_real_value, np.abs, -np.inf, np.nan, np.inf
# cabs(conj(z)) == conj(cabs(z)) (= cabs(z))
def f(a):
return np.abs(np.conj(a))
def g(a, b):
return np.abs(np.complex(a, b))
xa = np.array(x, dtype=np.complex)
for i in range(len(xa)):
ref = g(x[i], y[i])
yield check_real_value, f, x[i], y[i], ref
def cor_fft(x1,x2,sigma):
dist11 = np.sum(np.square(x1))
dist22 = np.sum(np.square(x2))
if len(x1.shape)==2:
c = np.fft.ifft2((np.conj(fft(x1))*fft(x2)))
else:
c = np.fft.ifft2(np.sum(np.conj(fft(x1))*fft(x2),2))
dist= dist11-2*c+dist22
cor = np.exp(-1*dist/(sigma**2*x1.size))
cor = np.real(cor)
return cor