Python numpy 模块，conjugate() 实例源码
我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用numpy.conjugate()。
def correlate(self, imgfft):
#Very much related to the convolution theorem, the cross-correlation
#theorem states that the Fourier transform of the cross-correlation of
#two functions is equal to the product of the individual Fourier
#transforms, where one of them has been complex conjugated:
if self.imgfft is not 0 or imgfft.imgfft is not 0:
imgcj = np.conjugate(self.imgfft)
imgft = imgfft.imgfft
prod = deepcopy(imgcj)
for x in range(imgcj.shape[0]):
for y in range(imgcj.shape[0]):
prod[x][y] = imgcj[x][y] * imgft[x][y]
cc = Corr( np.real(fft.ifft2(fft.fftshift(prod)))) # real image of the correlation
# adjust to center
cc.data = np.roll(cc.data, int(cc.data.shape[0] / 2), axis = 0)
cc.data = np.roll(cc.data, int(cc.data.shape[1] / 2), axis = 1)
else:
raise FFTnotInit()
return cc
def genSpectra(time,dipole,signal):
fw, frequency = pade(time,dipole)
fw_sig, frequency = pade(time,signal,alternate=True)
fw_re = np.real(fw)
fw_im = np.imag(fw)
fw_abs = fw_re**2 + fw_im**2
#spectra = (fw_re*17.32)/(np.pi*field*damp_const)
#spectra = (fw_re*17.32*514.220652)/(np.pi*field*damp_const)
#numerator = np.imag((fw*np.conjugate(fw_sig)))
numerator = np.imag(fw_abs*np.conjugate(fw_sig))
#numerator = np.abs((fw*np.conjugate(fw_sig)))
#numerator = np.abs(fw)
denominator = np.real(np.conjugate(fw_sig)*fw_sig)
#denominator = 1.0
spectra = ((4.0*27.21138602*2*frequency*np.pi*(numerator))/(3.0*137.036*denominator))
spectra *= 1.0/100.0
#plt.plot(frequency*27.2114,fourier)
#plt.show()
return frequency, spectra
def test_generic_methods(self):
# Tests some MaskedArray methods.
a = array([1, 3, 2])
assert_equal(a.any(), a._data.any())
assert_equal(a.all(), a._data.all())
assert_equal(a.argmax(), a._data.argmax())
assert_equal(a.argmin(), a._data.argmin())
assert_equal(a.choose(0, 1, 2, 3, 4), a._data.choose(0, 1, 2, 3, 4))
assert_equal(a.compress([1, 0, 1]), a._data.compress([1, 0, 1]))
assert_equal(a.conj(), a._data.conj())
assert_equal(a.conjugate(), a._data.conjugate())
m = array([[1, 2], [3, 4]])
assert_equal(m.diagonal(), m._data.diagonal())
assert_equal(a.sum(), a._data.sum())
assert_equal(a.take([1, 2]), a._data.take([1, 2]))
assert_equal(m.transpose(), m._data.transpose())
def test_testArrayMethods(self):
a = array([1, 3, 2])
self.assertTrue(eq(a.any(), a._data.any()))
self.assertTrue(eq(a.all(), a._data.all()))
self.assertTrue(eq(a.argmax(), a._data.argmax()))
self.assertTrue(eq(a.argmin(), a._data.argmin()))
self.assertTrue(eq(a.choose(0, 1, 2, 3, 4),
a._data.choose(0, 1, 2, 3, 4)))
self.assertTrue(eq(a.compress([1, 0, 1]), a._data.compress([1, 0, 1])))
self.assertTrue(eq(a.conj(), a._data.conj()))
self.assertTrue(eq(a.conjugate(), a._data.conjugate()))
m = array([[1, 2], [3, 4]])
self.assertTrue(eq(m.diagonal(), m._data.diagonal()))
self.assertTrue(eq(a.sum(), a._data.sum()))
self.assertTrue(eq(a.take([1, 2]), a._data.take([1, 2])))
self.assertTrue(eq(m.transpose(), m._data.transpose()))
def test_spherical_harmonic(scheme):
'''Assert the norm of the spherical harmonic
Y_1^1(phi, theta) = -1/2 sqrt(3/2/pi) * exp(i*phi) * sin(theta)
is indeed 1, i.e.,
int_0^2pi int_0^pi
Y_1^1(phi, theta) * conj(Y_1^1(phi, theta)) * sin(theta)
dphi dtheta = 1.
'''
def spherical_harmonic_11(azimuthal, polar):
# y00 = 1.0 / numpy.sqrt(4*numpy.pi)
y11 = -0.5 * numpy.sqrt(3.0/2.0/numpy.pi) \
* numpy.exp(1j*azimuthal) * numpy.sin(polar)
return y11 * numpy.conjugate(y11)
val = quadpy.sphere.integrate_spherical(
spherical_harmonic_11,
rule=scheme
)
assert abs(val - 1.0) < 1.0e-15
return
def __init__(self, n):
self.degree = n
self.points = -numpy.cos(numpy.pi * (numpy.arange(n) + 0.5) / n)
# n -= 1
N = numpy.arange(1, n, 2)
length = len(N)
m = n - length
K = numpy.arange(m)
v0 = numpy.concatenate([
2 * numpy.exp(1j*numpy.pi*K/n) / (1 - 4*K**2),
numpy.zeros(length+1)
])
v1 = v0[:-1] + numpy.conjugate(v0[:0:-1])
w = numpy.fft.ifft(v1)
assert max(w.imag) < 1.0e-15
self.weights = w.real
return
def CoreShellScatteringFunction(mCore,mShell,wavelength,dCore,dShell,minAngle=0, maxAngle=180, angularResolution=0.5, normed=False):
# http://pymiescatt.readthedocs.io/en/latest/forwardCS.html#CoreShellScatteringFunction
xCore = np.pi*dCore/wavelength
xShell = np.pi*dShell/wavelength
theta = np.linspace(minAngle,maxAngle,int((maxAngle-minAngle)/angularResolution))*np.pi/180
thetaSteps = len(theta)
SL = np.zeros(thetaSteps)
SR = np.zeros(thetaSteps)
SU = np.zeros(thetaSteps)
for j in range(thetaSteps):
u = np.cos(theta[j])
S1,S2 = CoreShellS1S2(mCore,mShell,xCore,xShell,u)
SL[j] = (np.sum((np.conjugate(S1)*S1))).real
SR[j] = (np.sum((np.conjugate(S2)*S2))).real
SU[j] = (SR[j]+SL[j])/2
if normed:
SL /= np.max(SL)
SR /= np.max(SR)
SU /= np.max(SU)
return theta,SL,SR,SU
def CoreShellScatteringFunction(mCore,mShell,wavelength,dCore,dShell,minAngle=0, maxAngle=180, angularResolution=0.5, normed=False):
# http://pymiescatt.readthedocs.io/en/latest/forwardCS.html#CoreShellScatteringFunction
xCore = np.pi*dCore/wavelength
xShell = np.pi*dShell/wavelength
theta = np.linspace(minAngle,maxAngle,int((maxAngle-minAngle)/angularResolution))*np.pi/180
thetaSteps = len(theta)
SL = np.zeros(thetaSteps)
SR = np.zeros(thetaSteps)
SU = np.zeros(thetaSteps)
for j in range(thetaSteps):
u = np.cos(theta[j])
S1,S2 = CoreShellS1S2(mCore,mShell,xCore,xShell,u)
SL[j] = (np.sum((np.conjugate(S1)*S1))).real
SR[j] = (np.sum((np.conjugate(S2)*S2))).real
SU[j] = (SR[j]+SL[j])/2
if normed:
SL /= np.max(SL)
SR /= np.max(SR)
SU /= np.max(SU)
return theta,SL,SR,SU
def gain_substitution_scalar(gain, x, xwt):
nants, nchan, nrec, _ = gain.shape
newgain = numpy.ones_like(gain, dtype='complex')
gwt = numpy.zeros_like(gain, dtype='float')
# We are going to work with Jones 2x2 matrix formalism so everything has to be
# converted to that format
x = x.reshape(nants, nants, nchan, nrec, nrec)
xwt = xwt.reshape(nants, nants, nchan, nrec, nrec)
for ant1 in range(nants):
for chan in range(nchan):
# Loop over e.g. 'RR', 'LL, or 'xx', 'YY' ignoring cross terms
top = numpy.sum(x[:, ant1, chan, 0, 0] *
gain[:, chan, 0, 0] * xwt[:, ant1, chan, 0, 0], axis=0)
bot = numpy.sum((gain[:, chan, 0, 0] * numpy.conjugate(gain[:, chan, 0, 0]) *
xwt[:, ant1, chan, 0, 0]).real, axis=0)
if bot > 0.0:
newgain[ant1, chan, 0, 0] = top / bot
gwt[ant1, chan, 0, 0] = bot
else:
newgain[ant1, chan, 0, 0] = 0.0
gwt[ant1, chan, 0, 0] = 0.0
return newgain, gwt
def convolve_scalestack(scalestack, img):
"""Convolve img by the specified scalestack, returning the resulting stack
:param scalestack: stack containing the scales
:param img: Image to be convolved
:return: stack
"""
convolved = numpy.zeros(scalestack.shape)
ximg = numpy.fft.fftshift(numpy.fft.fft2(numpy.fft.fftshift(img)))
nscales = scalestack.shape[0]
for iscale in range(nscales):
xscale = numpy.fft.fftshift(numpy.fft.fft2(numpy.fft.fftshift(scalestack[iscale, :, :])))
xmult = ximg * numpy.conjugate(xscale)
convolved[iscale, :, :] = numpy.real(numpy.fft.ifftshift(numpy.fft.ifft2(numpy.fft.ifftshift(xmult))))
return convolved
def test_generic_methods(self):
# Tests some MaskedArray methods.
a = array([1, 3, 2])
assert_equal(a.any(), a._data.any())
assert_equal(a.all(), a._data.all())
assert_equal(a.argmax(), a._data.argmax())
assert_equal(a.argmin(), a._data.argmin())
assert_equal(a.choose(0, 1, 2, 3, 4), a._data.choose(0, 1, 2, 3, 4))
assert_equal(a.compress([1, 0, 1]), a._data.compress([1, 0, 1]))
assert_equal(a.conj(), a._data.conj())
assert_equal(a.conjugate(), a._data.conjugate())
m = array([[1, 2], [3, 4]])
assert_equal(m.diagonal(), m._data.diagonal())
assert_equal(a.sum(), a._data.sum())
assert_equal(a.take([1, 2]), a._data.take([1, 2]))
assert_equal(m.transpose(), m._data.transpose())
def test_testArrayMethods(self):
a = array([1, 3, 2])
self.assertTrue(eq(a.any(), a._data.any()))
self.assertTrue(eq(a.all(), a._data.all()))
self.assertTrue(eq(a.argmax(), a._data.argmax()))
self.assertTrue(eq(a.argmin(), a._data.argmin()))
self.assertTrue(eq(a.choose(0, 1, 2, 3, 4),
a._data.choose(0, 1, 2, 3, 4)))
self.assertTrue(eq(a.compress([1, 0, 1]), a._data.compress([1, 0, 1])))
self.assertTrue(eq(a.conj(), a._data.conj()))
self.assertTrue(eq(a.conjugate(), a._data.conjugate()))
m = array([[1, 2], [3, 4]])
self.assertTrue(eq(m.diagonal(), m._data.diagonal()))
self.assertTrue(eq(a.sum(), a._data.sum()))
self.assertTrue(eq(a.take([1, 2]), a._data.take([1, 2])))
self.assertTrue(eq(m.transpose(), m._data.transpose()))
def test_generic_methods(self):
# Tests some MaskedArray methods.
a = array([1, 3, 2])
assert_equal(a.any(), a._data.any())
assert_equal(a.all(), a._data.all())
assert_equal(a.argmax(), a._data.argmax())
assert_equal(a.argmin(), a._data.argmin())
assert_equal(a.choose(0, 1, 2, 3, 4), a._data.choose(0, 1, 2, 3, 4))
assert_equal(a.compress([1, 0, 1]), a._data.compress([1, 0, 1]))
assert_equal(a.conj(), a._data.conj())
assert_equal(a.conjugate(), a._data.conjugate())
m = array([[1, 2], [3, 4]])
assert_equal(m.diagonal(), m._data.diagonal())
assert_equal(a.sum(), a._data.sum())
assert_equal(a.take([1, 2]), a._data.take([1, 2]))
assert_equal(m.transpose(), m._data.transpose())
def test_testArrayMethods(self):
a = array([1, 3, 2])
self.assertTrue(eq(a.any(), a._data.any()))
self.assertTrue(eq(a.all(), a._data.all()))
self.assertTrue(eq(a.argmax(), a._data.argmax()))
self.assertTrue(eq(a.argmin(), a._data.argmin()))
self.assertTrue(eq(a.choose(0, 1, 2, 3, 4),
a._data.choose(0, 1, 2, 3, 4)))
self.assertTrue(eq(a.compress([1, 0, 1]), a._data.compress([1, 0, 1])))
self.assertTrue(eq(a.conj(), a._data.conj()))
self.assertTrue(eq(a.conjugate(), a._data.conjugate()))
m = array([[1, 2], [3, 4]])
self.assertTrue(eq(m.diagonal(), m._data.diagonal()))
self.assertTrue(eq(a.sum(), a._data.sum()))
self.assertTrue(eq(a.take([1, 2]), a._data.take([1, 2])))
self.assertTrue(eq(m.transpose(), m._data.transpose()))
def cumsum_snr(freqs, detector, data, hpf, hxf, theta, phi, psi, zeroFreq=False, normalizeTemplate=True ):
"""
returns the cumulative sum of the snr as a function of frequency
does NOT maximize over the phase at coalescence
"""
template = detector.project(freqs, hpf, hxf, theta, phi, psi, zeroFreq=zeroFreq)
PSD = detector.PSD(freqs)
deltaF = freqs[1]-freqs[0]
ans = 2*np.cumsum(deltaF*np.conjugate(data)*template/PSD).real
if normalizeTemplate:
ans /= np.sum(deltaF*np.conjugate(template)*template/PSD).real**0.5
return ans
#------------------------
def __newlambdagammaC(self, theta, l):
""" Apply Eqns. 25-27 in Vidal 2003 to update lambda^C and gamma^C
(lambda and gamma of this qbit).
"""
gamma_ket = self.coefs[l+1].lam
gamma_bra = np.conjugate(gamma_ket)
Gamma_star = np.conjugate(self.coefs[l+1].gamma)
inputs = [Gamma_star, theta, gamma_bra, gamma_ket]
Gamma_star_idx = [1, -3, -2]
theta_idx = [-1, 1, -4, -5]
gamma_bra_idx = [-6]
gamma_ket_idx = [-7]
idx = [Gamma_star_idx, theta_idx, gamma_bra_idx, gamma_ket_idx]
contract_me = scon(inputs, idx)
svd_me = np.einsum('agibggg', contract_me)
evals, evecs = la.eigh(svd_me)
return evals, evecs
def paulidouble(i, j, tensor=True):
pauli_i = pauli(i)
pauli_j = pauli(j)
outer = np.zeros((4, 4), dtype=np.complex)
outer[:2, :2] = pauli_i
outer[2:, 2:] = pauli_j
if tensor:
outer.shape = (2, 2, 2, 2)
return outer
# def paulitwo_left(i):
# return np.kron(pauli(i), pauli(0))
# def paulitwo_right(i):
# return np.kron(pauli(0), pauli(i))
# def newrho_DK(Lket, theta_ij):
# Lbra = np.conjugate(L_before)
# theta_star = np.conjugate(theta_ij)
# in_bracket = scon(Lbra, Lket, theta_ij, theta_star,
# [1], [1], [2, 3, 1,
def test_generic_methods(self):
# Tests some MaskedArray methods.
a = array([1, 3, 2])
assert_equal(a.any(), a._data.any())
assert_equal(a.all(), a._data.all())
assert_equal(a.argmax(), a._data.argmax())
assert_equal(a.argmin(), a._data.argmin())
assert_equal(a.choose(0, 1, 2, 3, 4), a._data.choose(0, 1, 2, 3, 4))
assert_equal(a.compress([1, 0, 1]), a._data.compress([1, 0, 1]))
assert_equal(a.conj(), a._data.conj())
assert_equal(a.conjugate(), a._data.conjugate())
m = array([[1, 2], [3, 4]])
assert_equal(m.diagonal(), m._data.diagonal())
assert_equal(a.sum(), a._data.sum())
assert_equal(a.take([1, 2]), a._data.take([1, 2]))
assert_equal(m.transpose(), m._data.transpose())
def test_testArrayMethods(self):
a = array([1, 3, 2])
self.assertTrue(eq(a.any(), a._data.any()))
self.assertTrue(eq(a.all(), a._data.all()))
self.assertTrue(eq(a.argmax(), a._data.argmax()))
self.assertTrue(eq(a.argmin(), a._data.argmin()))
self.assertTrue(eq(a.choose(0, 1, 2, 3, 4),
a._data.choose(0, 1, 2, 3, 4)))
self.assertTrue(eq(a.compress([1, 0, 1]), a._data.compress([1, 0, 1])))
self.assertTrue(eq(a.conj(), a._data.conj()))
self.assertTrue(eq(a.conjugate(), a._data.conjugate()))
m = array([[1, 2], [3, 4]])
self.assertTrue(eq(m.diagonal(), m._data.diagonal()))
self.assertTrue(eq(a.sum(), a._data.sum()))
self.assertTrue(eq(a.take([1, 2]), a._data.take([1, 2])))
self.assertTrue(eq(m.transpose(), m._data.transpose()))
def qisunitary(self,op):
mat = op[1]
(r,c) = mat.shape
if r != c:
return False
invmat = np.conjugate(np.transpose(mat))
pmat = np.asarray(mat * invmat)
for i in range(r):
for j in range(c):
if i != j:
if np.absolute(pmat[i][j]) > self.maxerr:
return False
else:
if np.absolute(pmat[i][j]-1.0) > self.maxerr:
return False
return True
def __init__(self, gain):
self.gain = gain * np.pi # pi term puts angle output to pi range
# components
self.conjugate = Conjugate()
self.complex_mult = ComplexMultiply()
self.angle = Angle()
self.out = Sfix()
# specify component delay
self._delay = self.conjugate._delay + \
self.complex_mult._delay + \
self.angle._delay + 1
# constants
self.gain_sfix = Const(Sfix(self.gain, 3, -14))
def test_generic_methods(self):
# Tests some MaskedArray methods.
a = array([1, 3, 2])
assert_equal(a.any(), a._data.any())
assert_equal(a.all(), a._data.all())
assert_equal(a.argmax(), a._data.argmax())
assert_equal(a.argmin(), a._data.argmin())
assert_equal(a.choose(0, 1, 2, 3, 4), a._data.choose(0, 1, 2, 3, 4))
assert_equal(a.compress([1, 0, 1]), a._data.compress([1, 0, 1]))
assert_equal(a.conj(), a._data.conj())
assert_equal(a.conjugate(), a._data.conjugate())
m = array([[1, 2], [3, 4]])
assert_equal(m.diagonal(), m._data.diagonal())
assert_equal(a.sum(), a._data.sum())
assert_equal(a.take([1, 2]), a._data.take([1, 2]))
assert_equal(m.transpose(), m._data.transpose())
def test_testArrayMethods(self):
a = array([1, 3, 2])
self.assertTrue(eq(a.any(), a._data.any()))
self.assertTrue(eq(a.all(), a._data.all()))
self.assertTrue(eq(a.argmax(), a._data.argmax()))
self.assertTrue(eq(a.argmin(), a._data.argmin()))
self.assertTrue(eq(a.choose(0, 1, 2, 3, 4),
a._data.choose(0, 1, 2, 3, 4)))
self.assertTrue(eq(a.compress([1, 0, 1]), a._data.compress([1, 0, 1])))
self.assertTrue(eq(a.conj(), a._data.conj()))
self.assertTrue(eq(a.conjugate(), a._data.conjugate()))
m = array([[1, 2], [3, 4]])
self.assertTrue(eq(m.diagonal(), m._data.diagonal()))
self.assertTrue(eq(a.sum(), a._data.sum()))
self.assertTrue(eq(a.take([1, 2]), a._data.take([1, 2])))
self.assertTrue(eq(m.transpose(), m._data.transpose()))
def test_poly(self):
assert_array_almost_equal(np.poly([3, -np.sqrt(2), np.sqrt(2)]),
[1, -3, -2, 6])
# From matlab docs
A = [[1, 2, 3], [4, 5, 6], [7, 8, 0]]
assert_array_almost_equal(np.poly(A), [1, -6, -72, -27])
# Should produce real output for perfect conjugates
assert_(np.isrealobj(np.poly([+1.082j, +2.613j, -2.613j, -1.082j])))
assert_(np.isrealobj(np.poly([0+1j, -0+-1j, 1+2j,
1-2j, 1.+3.5j, 1-3.5j])))
assert_(np.isrealobj(np.poly([1j, -1j, 1+2j, 1-2j, 1+3j, 1-3.j])))
assert_(np.isrealobj(np.poly([1j, -1j, 1+2j, 1-2j])))
assert_(np.isrealobj(np.poly([1j, -1j, 2j, -2j])))
assert_(np.isrealobj(np.poly([1j, -1j])))
assert_(np.isrealobj(np.poly([1, -1])))
assert_(np.iscomplexobj(np.poly([1j, -1.0000001j])))
np.random.seed(42)
a = np.random.randn(100) + 1j*np.random.randn(100)
assert_(np.isrealobj(np.poly(np.concatenate((a, np.conjugate(a))))))
def test_generic_methods(self):
# Tests some MaskedArray methods.
a = array([1, 3, 2])
assert_equal(a.any(), a._data.any())
assert_equal(a.all(), a._data.all())
assert_equal(a.argmax(), a._data.argmax())
assert_equal(a.argmin(), a._data.argmin())
assert_equal(a.choose(0, 1, 2, 3, 4), a._data.choose(0, 1, 2, 3, 4))
assert_equal(a.compress([1, 0, 1]), a._data.compress([1, 0, 1]))
assert_equal(a.conj(), a._data.conj())
assert_equal(a.conjugate(), a._data.conjugate())
m = array([[1, 2], [3, 4]])
assert_equal(m.diagonal(), m._data.diagonal())
assert_equal(a.sum(), a._data.sum())
assert_equal(a.take([1, 2]), a._data.take([1, 2]))
assert_equal(m.transpose(), m._data.transpose())
def test_testArrayMethods(self):
a = array([1, 3, 2])
self.assertTrue(eq(a.any(), a._data.any()))
self.assertTrue(eq(a.all(), a._data.all()))
self.assertTrue(eq(a.argmax(), a._data.argmax()))
self.assertTrue(eq(a.argmin(), a._data.argmin()))
self.assertTrue(eq(a.choose(0, 1, 2, 3, 4),
a._data.choose(0, 1, 2, 3, 4)))
self.assertTrue(eq(a.compress([1, 0, 1]), a._data.compress([1, 0, 1])))
self.assertTrue(eq(a.conj(), a._data.conj()))
self.assertTrue(eq(a.conjugate(), a._data.conjugate()))
m = array([[1, 2], [3, 4]])
self.assertTrue(eq(m.diagonal(), m._data.diagonal()))
self.assertTrue(eq(a.sum(), a._data.sum()))
self.assertTrue(eq(a.take([1, 2]), a._data.take([1, 2])))
self.assertTrue(eq(m.transpose(), m._data.transpose()))
def get_pmf_xi(self):
"""Return the values of the variable ``xi``.
The components ``xi`` make up the probability mass function, i.e.
:math:`\\xi(k) = pmf(k) = Pr(X = k)`.
"""
chi = np.empty(self.number_trials + 1, dtype=complex)
chi[0] = 1
half_number_trials = int(
self.number_trials / 2 + self.number_trials % 2)
# set first half of chis:
chi[1:half_number_trials + 1] = self.get_chi(
np.arange(1, half_number_trials + 1))
# set second half of chis:
chi[half_number_trials + 1:self.number_trials + 1] = np.conjugate(
chi[1:self.number_trials - half_number_trials + 1] [::-1])
chi /= self.number_trials + 1
xi = np.fft.fft(chi)
if self.check_xi_are_real(xi):
xi = xi.real
else:
raise TypeError("pmf / xi values have to be real.")
xi += np.finfo(type(xi[0])).eps
return xi
def test_generic_methods(self):
# Tests some MaskedArray methods.
a = array([1, 3, 2])
assert_equal(a.any(), a._data.any())
assert_equal(a.all(), a._data.all())
assert_equal(a.argmax(), a._data.argmax())
assert_equal(a.argmin(), a._data.argmin())
assert_equal(a.choose(0, 1, 2, 3, 4), a._data.choose(0, 1, 2, 3, 4))
assert_equal(a.compress([1, 0, 1]), a._data.compress([1, 0, 1]))
assert_equal(a.conj(), a._data.conj())
assert_equal(a.conjugate(), a._data.conjugate())
m = array([[1, 2], [3, 4]])
assert_equal(m.diagonal(), m._data.diagonal())
assert_equal(a.sum(), a._data.sum())
assert_equal(a.take([1, 2]), a._data.take([1, 2]))
assert_equal(m.transpose(), m._data.transpose())
def test_testArrayMethods(self):
a = array([1, 3, 2])
self.assertTrue(eq(a.any(), a._data.any()))
self.assertTrue(eq(a.all(), a._data.all()))
self.assertTrue(eq(a.argmax(), a._data.argmax()))
self.assertTrue(eq(a.argmin(), a._data.argmin()))
self.assertTrue(eq(a.choose(0, 1, 2, 3, 4),
a._data.choose(0, 1, 2, 3, 4)))
self.assertTrue(eq(a.compress([1, 0, 1]), a._data.compress([1, 0, 1])))
self.assertTrue(eq(a.conj(), a._data.conj()))
self.assertTrue(eq(a.conjugate(), a._data.conjugate()))
m = array([[1, 2], [3, 4]])
self.assertTrue(eq(m.diagonal(), m._data.diagonal()))
self.assertTrue(eq(a.sum(), a._data.sum()))
self.assertTrue(eq(a.take([1, 2]), a._data.take([1, 2])))
self.assertTrue(eq(m.transpose(), m._data.transpose()))
def _solve_lens_equation(self, complex_value):
"""
Solve the lens equation for the given point (in complex coordinates).
"""
complex_conjugate = np.conjugate(complex_value)
return complex_value - (1. / (1. + self.q)) * (
(1./complex_conjugate) + (self.q / (complex_conjugate - self.s)))
def _jacobian_determinant_ok_WM95(self, source_x, source_y):
"""determinants of lens equation Jacobian for verified roots"""
roots_ok_bar = np.conjugate(self._polynomial_roots_ok_WM95(
source_x=source_x, source_y=source_y))
# Variable X_bar is conjugate of variable X.
denominator_1 = self._position_z1_WM95 - roots_ok_bar
add_1 = self.mass_1 / denominator_1**2
denominator_2 = self._position_z2_WM95 - roots_ok_bar
add_2 = self.mass_2 / denominator_2**2
derivative = add_1 + add_2
# Can we use Utils.complex_fsum()-like function here?
return 1.-derivative*np.conjugate(derivative)
# Can we use Utils.complex_fsum()-like function here?
def adj(self,x):
"""Returns Hermitian adjoint of a matrix"""
assert x.shape[0] == x.shape[1]
return np.conjugate(x).T
def test_conjugate(self):
a = np.array([1-1j, 1+1j, 23+23.0j])
ac = a.conj()
assert_equal(a.real, ac.real)
assert_equal(a.imag, -ac.imag)
assert_equal(ac, a.conjugate())
assert_equal(ac, np.conjugate(a))
a = np.array([1-1j, 1+1j, 23+23.0j], 'F')
ac = a.conj()
assert_equal(a.real, ac.real)
assert_equal(a.imag, -ac.imag)
assert_equal(ac, a.conjugate())
assert_equal(ac, np.conjugate(a))
a = np.array([1, 2, 3])
ac = a.conj()
assert_equal(a, ac)
assert_equal(ac, a.conjugate())
assert_equal(ac, np.conjugate(a))
a = np.array([1.0, 2.0, 3.0])
ac = a.conj()
assert_equal(a, ac)
assert_equal(ac, a.conjugate())
assert_equal(ac, np.conjugate(a))
a = np.array([1-1j, 1+1j, 1, 2.0], object)
ac = a.conj()
assert_equal(ac, [k.conjugate() for k in a])
assert_equal(ac, a.conjugate())
assert_equal(ac, np.conjugate(a))
a = np.array([1-1j, 1, 2.0, 'f'], object)
assert_raises(AttributeError, lambda: a.conj())
assert_raises(AttributeError, lambda: a.conjugate())
def test_var_values(self):
for mat in [self.rmat, self.cmat, self.omat]:
for axis in [0, 1, None]:
msqr = _mean(mat * mat.conj(), axis=axis)
mean = _mean(mat, axis=axis)
tgt = msqr - mean * mean.conjugate()
res = _var(mat, axis=axis)
assert_almost_equal(res, tgt)
def test_oddfeatures_1(self):
# Test of other odd features
x = arange(20)
x = x.reshape(4, 5)
x.flat[5] = 12
assert_(x[1, 0] == 12)
z = x + 10j * x
assert_equal(z.real, x)
assert_equal(z.imag, 10 * x)
assert_equal((z * conjugate(z)).real, 101 * x * x)
z.imag[...] = 0.0
x = arange(10)
x[3] = masked
assert_(str(x[3]) == str(masked))
c = x >= 8
assert_(count(where(c, masked, masked)) == 0)
assert_(shape(where(c, masked, masked)) == c.shape)
z = masked_where(c, x)
assert_(z.dtype is x.dtype)
assert_(z[3] is masked)
assert_(z[4] is not masked)
assert_(z[7] is not masked)
assert_(z[8] is masked)
assert_(z[9] is masked)
assert_equal(x, z)
def test_basic_ufuncs(self):
# Test various functions such as sin, cos.
(x, y, a10, m1, m2, xm, ym, z, zm, xf) = self.d
assert_equal(np.cos(x), cos(xm))
assert_equal(np.cosh(x), cosh(xm))
assert_equal(np.sin(x), sin(xm))
assert_equal(np.sinh(x), sinh(xm))
assert_equal(np.tan(x), tan(xm))
assert_equal(np.tanh(x), tanh(xm))
assert_equal(np.sqrt(abs(x)), sqrt(xm))
assert_equal(np.log(abs(x)), log(xm))
assert_equal(np.log10(abs(x)), log10(xm))
assert_equal(np.exp(x), exp(xm))
assert_equal(np.arcsin(z), arcsin(zm))
assert_equal(np.arccos(z), arccos(zm))
assert_equal(np.arctan(z), arctan(zm))
assert_equal(np.arctan2(x, y), arctan2(xm, ym))
assert_equal(np.absolute(x), absolute(xm))
assert_equal(np.angle(x + 1j*y), angle(xm + 1j*ym))
assert_equal(np.angle(x + 1j*y, deg=True), angle(xm + 1j*ym, deg=True))
assert_equal(np.equal(x, y), equal(xm, ym))
assert_equal(np.not_equal(x, y), not_equal(xm, ym))
assert_equal(np.less(x, y), less(xm, ym))
assert_equal(np.greater(x, y), greater(xm, ym))
assert_equal(np.less_equal(x, y), less_equal(xm, ym))
assert_equal(np.greater_equal(x, y), greater_equal(xm, ym))
assert_equal(np.conjugate(x), conjugate(xm))
def test_testUfuncs1(self):
# Test various functions such as sin, cos.
(x, y, a10, m1, m2, xm, ym, z, zm, xf, s) = self.d
self.assertTrue(eq(np.cos(x), cos(xm)))
self.assertTrue(eq(np.cosh(x), cosh(xm)))
self.assertTrue(eq(np.sin(x), sin(xm)))
self.assertTrue(eq(np.sinh(x), sinh(xm)))
self.assertTrue(eq(np.tan(x), tan(xm)))
self.assertTrue(eq(np.tanh(x), tanh(xm)))
with np.errstate(divide='ignore', invalid='ignore'):
self.assertTrue(eq(np.sqrt(abs(x)), sqrt(xm)))
self.assertTrue(eq(np.log(abs(x)), log(xm)))
self.assertTrue(eq(np.log10(abs(x)), log10(xm)))
self.assertTrue(eq(np.exp(x), exp(xm)))
self.assertTrue(eq(np.arcsin(z), arcsin(zm)))
self.assertTrue(eq(np.arccos(z), arccos(zm)))
self.assertTrue(eq(np.arctan(z), arctan(zm)))
self.assertTrue(eq(np.arctan2(x, y), arctan2(xm, ym)))
self.assertTrue(eq(np.absolute(x), absolute(xm)))
self.assertTrue(eq(np.equal(x, y), equal(xm, ym)))
self.assertTrue(eq(np.not_equal(x, y), not_equal(xm, ym)))
self.assertTrue(eq(np.less(x, y), less(xm, ym)))
self.assertTrue(eq(np.greater(x, y), greater(xm, ym)))
self.assertTrue(eq(np.less_equal(x, y), less_equal(xm, ym)))
self.assertTrue(eq(np.greater_equal(x, y), greater_equal(xm, ym)))
self.assertTrue(eq(np.conjugate(x), conjugate(xm)))
self.assertTrue(eq(np.concatenate((x, y)), concatenate((xm, ym))))
self.assertTrue(eq(np.concatenate((x, y)), concatenate((x, y))))
self.assertTrue(eq(np.concatenate((x, y)), concatenate((xm, y))))
self.assertTrue(eq(np.concatenate((x, y, x)), concatenate((x, ym, x))))
def test_testUfuncRegression(self):
f_invalid_ignore = [
'sqrt', 'arctanh', 'arcsin', 'arccos',
'arccosh', 'arctanh', 'log', 'log10', 'divide',
'true_divide', 'floor_divide', 'remainder', 'fmod']
for f in ['sqrt', 'log', 'log10', 'exp', 'conjugate',
'sin', 'cos', 'tan',
'arcsin', 'arccos', 'arctan',
'sinh', 'cosh', 'tanh',
'arcsinh',
'arccosh',
'arctanh',
'absolute', 'fabs', 'negative',
'floor', 'ceil',
'logical_not',
'add', 'subtract', 'multiply',
'divide', 'true_divide', 'floor_divide',
'remainder', 'fmod', 'hypot', 'arctan2',
'equal', 'not_equal', 'less_equal', 'greater_equal',
'less', 'greater',
'logical_and', 'logical_or', 'logical_xor']:
try:
uf = getattr(umath, f)
except AttributeError:
uf = getattr(fromnumeric, f)
mf = getattr(np.ma, f)
args = self.d[:uf.nin]
with np.errstate():
if f in f_invalid_ignore:
np.seterr(invalid='ignore')
if f in ['arctanh', 'log', 'log10']:
np.seterr(divide='ignore')
ur = uf(*args)
mr = mf(*args)
self.assertTrue(eq(ur.filled(0), mr.filled(0), f))
self.assertTrue(eqmask(ur.mask, mr.mask))
def power_process(data, sfreq, toffset, modulus, integration, log_scale, zscale, title):
""" Break power by modulus and display each block. Integration here acts
as a pure average on the power level data.
"""
if modulus:
block = 0
block_size = integration * modulus
block_toffset = toffset
while block < len(data) / block_size:
vblock = data[block * block_size:block * block_size + modulus]
pblock = (vblock * numpy.conjugate(vblock)).real
# complete integration
for idx in range(1, integration):
vblock = data[block * block_size + idx *
modulus:block * block_size + idx * modulus + modulus]
pblock += (vblock * numpy.conjugate(vblock)).real
pblock /= integration
yield power_plot(pblock, sfreq, block_toffset, log_scale, zscale, title)
block += 1
block_toffset += block_size / sfreq
else:
pdata = (data * numpy.conjugate(data)).real
yield power_plot(pdata, sfreq, toffset, log_scale, zscale, title)
def MatrixElements(m,wavelength,diameter,mu):
# http://pymiescatt.readthedocs.io/en/latest/forward.html#MatrixElements
x = np.pi*diameter/wavelength
# B&H eqs. 4.77, where mu=cos(theta)
S1,S2 = MieS1S2(m,x,mu)
S11 = 0.5*(np.abs(S2)**2+np.abs(S1)**2)
S12 = 0.5*(np.abs(S2)**2-np.abs(S1)**2)
S33 = 0.5*(np.conjugate(S2)*S1+S2*np.conjugate(S1))
S34 = 0.5j*(S1*np.conjugate(S2)-S2*np.conjugate(S1))
return S11, S12, S33, S34
def CoreShellS1S2(mCore,mShell,xCore,xShell,mu):
# http://pymiescatt.readthedocs.io/en/latest/forwardCS.html#CoreShellS1S2
nmax = np.round(2+xShell+4*(xShell**(1/3)))
an,bn = CoreShell_ab(mCore,mShell,xCore,xShell)
pin,taun = MiePiTau(mu,nmax)
n = np.arange(1,int(nmax)+1)
n2 = (2*n+1)/(n*(n+1))
pin *= n2
taun *= n2
S1=np.sum(an*np.conjugate(pin))+np.sum(bn*np.conjugate(taun))
S2=np.sum(an*np.conjugate(taun))+np.sum(bn*np.conjugate(pin))
return S1,S2
def CoreShellMatrixElements(mCore,mShell,xCore,xShell,mu):
# http://pymiescatt.readthedocs.io/en/latest/forwardCS.html#CoreShellMatrixElements
S1,S2 = CoreShellS1S2(mCore,mShell,xCore,xShell,mu)
S11 = 0.5*(np.abs(S2)**2+np.abs(S1)**2)
S12 = 0.5*(np.abs(S2)**2-np.abs(S1)**2)
S33 = 0.5*(np.conjugate(S2)*S1+S2*np.conjugate(S1))
S34 = 0.5j*(S1*np.conjugate(S2)-S2*np.conjugate(S1))
return S11, S12, S33, S34
def MatrixElements(m,wavelength,diameter,mu):
# http://pymiescatt.readthedocs.io/en/latest/forward.html#MatrixElements
x = np.pi*diameter/wavelength
# B&H eqs. 4.77, where mu=cos(theta)
S1,S2 = MieS1S2(m,x,mu)
S11 = 0.5*(np.abs(S2)**2+np.abs(S1)**2)
S12 = 0.5*(np.abs(S2)**2-np.abs(S1)**2)
S33 = 0.5*(np.conjugate(S2)*S1+S2*np.conjugate(S1))
S34 = 0.5j*(S1*np.conjugate(S2)-S2*np.conjugate(S1))
return S11, S12, S33, S34
def CoreShellS1S2(mCore,mShell,xCore,xShell,mu):
# http://pymiescatt.readthedocs.io/en/latest/forwardCS.html#CoreShellS1S2
nmax = np.round(2+xShell+4*(xShell**(1/3)))
an,bn = CoreShell_ab(mCore,mShell,xCore,xShell)
pin,taun = MiePiTau(mu,nmax)
n = np.arange(1,int(nmax)+1)
n2 = (2*n+1)/(n*(n+1))
pin *= n2
taun *= n2
S1=np.sum(an*np.conjugate(pin))+np.sum(bn*np.conjugate(taun))
S2=np.sum(an*np.conjugate(taun))+np.sum(bn*np.conjugate(pin))
return S1,S2
def CoreShellMatrixElements(mCore,mShell,xCore,xShell,mu):
# http://pymiescatt.readthedocs.io/en/latest/forwardCS.html#CoreShellMatrixElements
S1,S2 = CoreShellS1S2(mCore,mShell,xCore,xShell,mu)
S11 = 0.5*(np.abs(S2)**2+np.abs(S1)**2)
S12 = 0.5*(np.abs(S2)**2-np.abs(S1)**2)
S33 = 0.5*(np.conjugate(S2)*S1+S2*np.conjugate(S1))
S34 = 0.5j*(S1*np.conjugate(S2)-S2*np.conjugate(S1))
return S11, S12, S33, S34
def sum_visibility(vis: Visibility, direction: SkyCoord) -> numpy.array:
""" Direct Fourier summation in a given direction
:param vis: Visibility to be summed
:param direction: Direction of summation
:return: flux[nch,npol], weight[nch,pol]
"""
# TODO: Convert to Visibility or remove?
assert isinstance(vis, Visibility) or isinstance(vis, BlockVisibility), vis
svis = copy_visibility(vis)
l, m, n = skycoord_to_lmn(direction, svis.phasecentre)
phasor = numpy.conjugate(simulate_point(svis.uvw, l, m))
# Need to put correct mapping here
_, frequency = get_frequency_map(svis, None)
frequency = list(frequency)
nchan = max(frequency) + 1
npol = svis.polarisation_frame.npol
flux = numpy.zeros([nchan, npol])
weight = numpy.zeros([nchan, npol])
coords = svis.vis, svis.weight, phasor, list(frequency)
for v, wt, p, ic in zip(*coords):
for pol in range(npol):
flux[ic, pol] += numpy.real(wt[pol] * v[pol] * p)
weight[ic, pol] += wt[pol]
flux[weight > 0.0] = flux[weight > 0.0] / weight[weight > 0.0]
flux[weight <= 0.0] = 0.0
return flux, weight
def gain_substitution_vector(gain, x, xwt):
nants, nchan, nrec, _ = gain.shape
newgain = numpy.ones_like(gain, dtype='complex')
if nrec > 0:
newgain[..., 0, 1] = 0.0
newgain[..., 1, 0] = 0.0
gwt = numpy.zeros_like(gain, dtype='float')
# We are going to work with Jones 2x2 matrix formalism so everything has to be
# converted to that format
x = x.reshape(nants, nants, nchan, nrec, nrec)
xwt = xwt.reshape(nants, nants, nchan, nrec, nrec)
if nrec > 0:
gain[..., 0, 1] = 0.0
gain[..., 1, 0] = 0.0
for ant1 in range(nants):
for chan in range(nchan):
# Loop over e.g. 'RR', 'LL, or 'xx', 'YY' ignoring cross terms
for rec in range(nrec):
top = numpy.sum(x[:, ant1, chan, rec, rec] *
gain[:, chan, rec, rec] * xwt[:, ant1, chan, rec, rec], axis=0)
bot = numpy.sum((gain[:, chan, rec, rec] * numpy.conjugate(gain[:, chan, rec, rec]) *
xwt[:, ant1, chan, rec, rec]).real, axis=0)
if bot > 0.0:
newgain[ant1, chan, rec, rec] = top / bot
gwt[ant1, chan, rec, rec] = bot
else:
newgain[ant1, chan, rec, rec] = 0.0
gwt[ant1, chan, rec, rec] = 0.0
return newgain, gwt
def gain_substitution_matrix(gain, x, xwt):
nants, nchan, nrec, _ = gain.shape
newgain = numpy.ones_like(gain, dtype='complex')
gwt = numpy.zeros_like(gain, dtype='float')
# We are going to work with Jones 2x2 matrix formalism so everything has to be
# converted to that format
x = x.reshape(nants, nants, nchan, nrec, nrec)
xwt = xwt.reshape(nants, nants, nchan, nrec, nrec)
# Write these loops out explicitly. Derivation of these vector equations is tedious but they are
# structurally identical to the scalar case with the following changes
# Vis -> 2x2 coherency vector, g-> 2x2 Jones matrix, *-> matmul, conjugate->Hermitean transpose (.H)
for ant1 in range(nants):
for chan in range(nchan):
top = 0.0
bot = 0.0
for ant2 in range(nants):
if ant1 != ant2:
xmat = x[ant2, ant1, chan]
xwtmat = xwt[ant2, ant1, chan]
g2 = gain[ant2, chan]
top += xmat * xwtmat * g2
bot += numpy.conjugate(g2) * xwtmat * g2
newgain[ant1, chan][bot > 0.0] = top[bot > 0.0] / bot[bot > 0.0]
newgain[ant1, chan][bot <= 0.0] = 0.0
gwt[ant1, chan] = bot.real
return newgain, gwt
def solution_residual_scalar(gain, x, xwt):
"""Calculate residual across all baselines of gain for point source equivalent visibilities
:param gain: gain [nant, ...]
:param x: Point source equivalent visibility [nant, ...]
:param xwt: Point source equivalent weight [nant, ...]
:return: residual[...]
"""
nants, nchan, nrec, _ = gain.shape
x = x.reshape(nants, nants, nchan, nrec, nrec)
xwt = xwt.reshape(nants, nants, nchan, nrec, nrec)
residual = numpy.zeros([nchan, nrec, nrec])
sumwt = numpy.zeros([nchan, nrec, nrec])
for ant1 in range(nants):
for ant2 in range(nants):
for chan in range(nchan):
error = x[ant2, ant1, chan, 0, 0] - \
gain[ant1, chan, 0, 0] * numpy.conjugate(gain[ant2, chan, 0, 0])
residual += (error * xwt[ant2, ant1, chan, 0, 0] * numpy.conjugate(error)).real
sumwt += xwt[ant2, ant1, chan, 0, 0]
residual[sumwt > 0.0] = numpy.sqrt(residual[sumwt > 0.0] / sumwt[sumwt > 0.0])
residual[sumwt <= 0.0] = 0.0
return residual
def solution_residual_vector(gain, x, xwt):
"""Calculate residual across all baselines of gain for point source equivalent visibilities
Vector case i.e. off-diagonals of gains are zero
:param gain: gain [nant, ...]
:param x: Point source equivalent visibility [nant, ...]
:param xwt: Point source equivalent weight [nant, ...]
:return: residual[...]
"""
nants, nchan, nrec, _ = gain.shape
x = x.reshape(nants, nants, nchan, nrec, nrec)
x[..., 1, 0] = 0.0
x[..., 0, 1] = 0.0
xwt = xwt.reshape(nants, nants, nchan, nrec, nrec)
xwt[..., 1, 0] = 0.0
xwt[..., 0, 1] = 0.0
residual = numpy.zeros([nchan, nrec, nrec])
sumwt = numpy.zeros([nchan, nrec, nrec])
for ant1 in range(nants):
for ant2 in range(nants):
for chan in range(nchan):
for rec in range(nrec):
error = x[ant2, ant1, chan, rec, rec] - \
gain[ant1, chan, rec, rec] * numpy.conjugate(gain[ant2, chan, rec, rec])
residual += (error * xwt[ant2, ant1, chan, rec, rec] * numpy.conjugate(error)).real
sumwt += xwt[ant2, ant1, chan, rec, rec]
residual[sumwt > 0.0] = numpy.sqrt(residual[sumwt > 0.0] / sumwt[sumwt > 0.0])
residual[sumwt <= 0.0] = 0.0
return residual