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
python类conj()的实例源码
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)
