我们从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 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 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 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_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 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 _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 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
