The following are code examples for showing how to use . They are extracted from open source Python projects. You can vote up the examples you like or vote down the exmaples you don’t like. You can also save this page to your account.
Example 1
def encode(img_path, wm_path, res_path, alpha):
img = cv2.imread(img_path)
img_f = np.fft.fft2(img)
height, width, channel = np.shape(img)
watermark = cv2.imread(wm_path)
wm_height, wm_width = watermark.shape[0], watermark.shape[1]
x, y = range(height / 2), range(width)
random.seed(height + width)
random.shuffle(x)
random.shuffle(y)
tmp = np.zeros(img.shape)
for i in range(height / 2):
for j in range(width):
if x[i] < wm_height and y[j] < wm_width:
tmp[i][j] = watermark[x[i]][y[j]]
tmp[height - 1 - i][width - 1 - j] = tmp[i][j]
res_f = img_f + alpha * tmp
res = np.fft.ifft2(res_f)
res = np.real(res)
cv2.imwrite(res_path, res, [int(cv2.IMWRITE_JPEG_QUALITY), 100]) Example 2
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 Example 3
def quaternion_real(quaternion):
"""Return real part of quaternion.
>>> quaternion_real([3, 0, 1, 2])
3.0
"""
return float(quaternion[0]) Example 4
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 Example 5
def save_image(imager, grid_data, grid_norm, output_file):
"""Makes an image from gridded visibilities and saves it to a FITS file.
Args:
imager (oskar.Imager): Handle to configured imager.
grid_data (numpy.ndarray): Final visibility grid.
grid_norm (float): Grid normalisation to apply.
output_file (str): Name of output FITS file to write.
"""
# Make the image (take the FFT, normalise, and apply grid correction).
imager.finalise_plane(grid_data, grid_norm)
grid_data = numpy.real(grid_data)
# Trim the image if required.
border = (imager.plane_size - imager.image_size) // 2
if border > 0:
end = border + imager.image_size
grid_data = grid_data[border:end, border:end]
# Write the FITS file.
hdr = fits.header.Header()
fits.writeto(output_file, grid_data, hdr, clobber=True) Example 6
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 Example 7
def kaiser_bessel_ft(u, J, alpha, kb_m, d):
'''
Interpolation weight for given J/alpha/kb-m
'''
u = u * (1.0 + 0.0j)
import scipy.special
z = numpy.sqrt((2 * numpy.pi * (J / 2) * u) ** 2.0 - alpha ** 2.0)
nu = d / 2 + kb_m
y = ((2 * numpy.pi) ** (d / 2)) * ((J / 2) ** d) * (alpha ** kb_m) / \
scipy.special.iv(kb_m, alpha) * scipy.special.jv(nu, z) / (z ** nu)
y = numpy.real(y)
return y Example 8
def mdct(x, odd=True):
""" Calculate modified discrete cosine transform of input signal
Parameters
----------
X : array_like
The input signal
odd : boolean, optional
Switch to oddly stacked transform. Defaults to :code:`True`.
Returns
-------
out : array_like
The output signal
"""
return numpy.real(cmdct(x, odd=odd)) * numpy.sqrt(2) Example 9
def ift(self):
return MyImage(np.real(fft.ifft2(fft.fftshift(self.ft)))) Example 10
def ift(self):
self.imgifft = MyImage(np.real(fft.ifft2(fft.fftshift(self.imgfft)))) Example 11
def get_polar_t(self):
mag = self.get_magnitude()
sizeimg = np.real(self.imgfft).shape
pol = np.zeros(sizeimg)
for x in range(sizeimg[0]):
for y in range(sizeimg[1]):
my = y - sizeimg[1] / 2
mx = x - sizeimg[0] / 2
if mx != 0:
phi = np.arctan(my / float(mx))
else:
phi = 0
r = np.sqrt(mx**2 + my **2)
ix = map_range(phi, -np.pi, np.pi, sizeimg[0], 0)
iy = map_range(r, 0, sizeimg[0], 0, sizeimg[1])
if ix >= 0 and ix < sizeimg[0] and iy >= 0 and iy < sizeimg[1]:
pol[x][y] = mag.data[int(ix)][int(iy)]
pol = MyImage(pol)
pol.limit(1)
return pol Example 12
def operate(self, x):
"""
Apply the separable filter to the signal vector *x*.
"""
X = NP.fft.fftn(x, s=self.k_full)
if NP.isrealobj(self.h) and NP.isrealobj(x):
y = NP.real(NP.fft.ifftn(self.H * X))
else:
y = NP.fft.ifftn(self.H * X)
if self.mode == 'full' or self.mode == 'circ':
return y
elif self.mode == 'valid':
slice_list = []
for i in range(self.ndim):
if self.m[i]-1 == 0:
slice_list.append(slice(None, None, None))
else:
slice_list.append(slice(self.m[i]-1, -(self.m[i]-1), None))
return y[slice_list]
else:
assert(False) Example 13
def chain(mpas, astype=None):
"""Computes the tensor product of MPAs given in ``*args`` by adding more
sites to the array.
:param mpas: Iterable of MPAs in the order as they should appear in the
chain
:param astype: dtype of the returned MPA. If ``None``, use the type of the
first MPA.
:returns: MPA of length ``len(args[0]) + ... + len(args[-1])``
.. todo:: Make this canonicalization aware
.. todo:: Raise warning when casting complex to real dtype
"""
mpas = iter(mpas)
try:
first = next(mpas)
except StopIteration:
raise ValueError('Argument `mpas` is an empty list')
rest = (lt for mpa in mpas for lt in mpa.lt)
if astype is None:
astype = type(first)
return astype(it.chain(first.lt, rest)) Example 14
def fftDf( df , part = "abs") :
#Handle series or DataFrame
if type(df) == pd.Series :
df = pd.DataFrame(df)
ise = True
else :
ise = False
res = pd.DataFrame( index = np.fft.rfftfreq( df.index.size, d = dx( df ) ) )
for col in df.columns :
if part == "abs" :
res["FFT_"+col] = np.abs( np.fft.rfft(df[col]) ) / (0.5*df.index.size)
elif part == "real" :
res["FFT_"+col] = np.real( np.fft.rfft(df[col]) ) / (0.5*df.index.size)
elif part == "imag" :
res["FFT_"+col] = np.imag( np.fft.rfft(df[col]) ) / (0.5*df.index.size)
if ise :
return res.iloc[:,0]
else :
return res Example 15
def derivFFT(df, n=1 ) :
""" Deriv a signal trought FFT, warning, edge can be a bit noisy...
indexList : channel to derive
n : order of derivation
"""
deriv = []
for iSig in range(df.shape[1]) :
fft = np.fft.fft( df.values[:,iSig] ) #FFT
freq = np.fft.fftfreq( df.shape[0] , dx(df) )
from copy import deepcopy
fft0 = deepcopy(fft)
if n>0 :
fft *= (1j * 2*pi* freq[:])**n #Derivation in frequency domain
else :
fft[-n:] *= (1j * 2*pi* freq[-n:])**n
fft[0:-n] = 0.
tts = np.real(np.fft.ifft(fft))
tts -= tts[0]
deriv.append( tts ) #Inverse FFT
return pd.DataFrame( data = np.transpose(deriv), index = df.index , columns = [ "DerivFFT("+ x +")" for x in df.columns ] ) Example 16
def process(self, wave):
wave.check_mono()
if wave.sample_rate != self.sr:
raise Exception("Wrong sample rate")
n = int(np.ceil(2 * wave.num_frames / float(self.w_len)))
m = (n + 1) * self.w_len / 2
swindow = self.make_signal_window(n)
win_ratios = [self.window / swindow[t * self.w_len / 2 :
t * self.w_len / 2 + self.w_len]
for t in range(n)]
wave = wave.zero_pad(0, int(m - wave.num_frames))
wave = audio.Wave(signal.hilbert(wave), wave.sample_rate)
result = np.zeros((self.n_bins, n))
for b in range(self.n_bins):
w = self.widths[b]
wc = 1 / np.square(w + 1)
filter = self.filters[b]
band = fftfilt(filter, wave.zero_pad(0, int(2 * w))[:,0])
band = band[int(w) : int(w + m), np.newaxis]
for t in range(n):
frame = band[t * self.w_len / 2:
t * self.w_len / 2 + self.w_len,:] * win_ratios[t]
result[b, t] = wc * np.real(np.conj(np.dot(frame.conj().T, frame)))
return audio.Spectrogram(result, self.sr, self.w_len, self.w_len / 2) Example 17
def test_psi(adjcube):
"""Tests retrieval of the wave functions and eigenvalues.
"""
from pydft.bases.fourier import psi, O, H
cell = adjcube
V = QHO(cell)
W = W4(cell)
Ns = W.shape[1]
Psi, epsilon = psi(V, W, cell, forceR=False)
#Make sure that the eigenvalues are real.
assert np.sum(np.imag(epsilon)) < 1e-13
checkI = np.dot(Psi.conjugate().T, O(Psi, cell))
assert abs(np.sum(np.diag(checkI))-Ns) < 1e-13 # Should be the identity
assert np.abs(np.sum(checkI)-Ns) < 1e-13
checkD = np.dot(Psi.conjugate().T, H(V, Psi, cell))
diagsum = np.sum(np.diag(checkD))
assert np.abs(np.sum(checkD)-diagsum) < 1e-12 # Should be diagonal
# Should match the diagonal elements of previous matrix
assert np.allclose(np.diag(checkD), epsilon) Example 18
def psi(V, W, cell, forceR=True):
"""Calculates the normalized wave functions using the basis coefficients.
Args:
V (pydft.potential.Potential): describing the potential for the
particles.
W (numpy.ndarray): wave function sample points.
cell (pydft.geometry.Cell): describing the unit cell and sampling
points.
forceR (bool): forces the result to be real.
"""
WN = Y(W, cell)
mu = np.dot(WN.conjugate().T, H(V, WN, cell))
epsilon, D = np.linalg.eig(mu)
if forceR:
epsilon = np.real(epsilon)
return (np.dot(WN, D), epsilon) Example 19
def Idag(v=None, cell=None):
"""Computes the complex conjugate of the `I` operator for Fourier basis.
Args:
v (numpy.ndarray): if None, then return the matrix :math:`\mathbb{I^\dag}`,
else, return :math:`\mathbb{I^\dag}\cdot v`.
cell (pydft.geometry.Cell): that describes the unit cell and
sampling points for real and reciprocal space.
"""
#It turns out that for Fourier, the complex conjugate only differs by a -1
#on the i (symmetric in R and G), so that we can return J instead.
cell = get_cell(cell)
def ifft(X):
FB = np.fft.ifftn(np.reshape(X, cell.S, order='F'))
return np.reshape(FB, X.shape, order='F')
return matprod(ifft, v)*np.prod(cell.S) Example 20
def Jdag(v=None, cell=None):
"""Computes the complex conjugate of the `J` operator for Fourier basis.
Args:
v (numpy.ndarray): if None, then return the matrix :math:`\mathbb{J^\dag}`,
else, return :math:`\mathbb{J^\dag}\cdot v`.
cell (pydft.geometry.Cell): that describes the unit cell and
sampling points for real and reciprocal space.
"""
#It turns out that for Fourier, the complex conjugate only differs by a -1
#on the i (symmetric in R and G), so that we can return I instead.
cell = get_cell(cell)
def fft(X):
FF = np.fft.fftn(np.reshape(X, cell.S, order='F'))
return np.reshape(FF, X.shape, order='F')
return matprod(fft, v)/np.prod(cell.S) Example 21
def plotProfileXZplane(Ax,X,Z,Hx,Hz,Flag):
FS = 20
if Flag == 'Hp':
Ax.streamplot(X,Z,Hx,Hz,color='b',linewidth=3.5,arrowsize=2)
Ax.set_title('Primary Field',fontsize=FS+6)
elif Flag == 'Hs_real':
Ax.streamplot(X,Z,Hx,Hz,color='r',linewidth=3.5,arrowsize=2)
Ax.set_title('Secondary Field (real)',fontsize=FS+6)
elif Flag == 'Hs_imag':
Ax.streamplot(X,Z,Hx,Hz,color='r',linewidth=3.5,arrowsize=2)
Ax.set_title('Secondary Field (imaginary)',fontsize=FS+6)
Ax.set_xbound(np.min(X),np.max(X))
Ax.set_ybound(np.min(Z),np.max(Z))
Ax.set_xlabel('X [m]',fontsize=FS+2)
Ax.set_ylabel('Z [m]',fontsize=FS+2,labelpad=-10)
Ax.tick_params(labelsize=FS-2) Example 22
def decode(ori_path, img_path, res_path, alpha):
ori = cv2.imread(ori_path)
img = cv2.imread(img_path)
ori_f = np.fft.fft2(ori)
img_f = np.fft.fft2(img)
height, width = ori.shape[0], ori.shape[1]
watermark = (ori_f - img_f) / alpha
watermark = np.real(watermark)
res = np.zeros(watermark.shape)
random.seed(height + width)
x = range(height / 2)
y = range(width)
random.shuffle(x)
random.shuffle(y)
for i in range(height / 2):
for j in range(width):
res[x[i]][y[j]] = watermark[i][j]
cv2.imwrite(res_path, res, [int(cv2.IMWRITE_JPEG_QUALITY), 100]) Example 23
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 Example 24
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) Example 25
def __init__(self,jet,kernels,k,x,y,pt,subpixel):
self.jet = jet
self.kernels = kernels
self.k = k
self.x = x
self.y = y
re = np.real(jet)
im = np.imag(jet)
self.mag = np.sqrt(re*re + im*im)
self.phase = np.arctan2(re,im)
if subpixel:
d = np.array([[pt.X()-x],[pt.Y()-y]])
comp = np.dot(self.k,d)
self.phase -= comp.flatten()
self.jet = self.mag*np.exp(1.0j*self.phase) Example 26
def calculate_katz_centrality(graph):
"""
Compute the katz centrality for nodes.
"""
# if not graph.is_directed():
# raise nx.NetworkXError( \
# "katz_centrality() not defined for undirected graphs.")
print "\n\tCalculating Katz Centrality..."
print "\tWarning: This might take a long time larger pedigrees."
g = graph
A = nx.adjacency_matrix(g)
from scipy import linalg as LA
max_eign = float(np.real(max(LA.eigvals(A.todense()))))
print "\t-Max.Eigenvalue(A) ", round(max_eign, 3)
kt = nx.katz_centrality(g, tol=1.0e-4, alpha=1/max_eign-0.01, beta=1.0, max_iter=999999)
nx.set_node_attributes(g, 'katz', kt)
katz_sorted = sorted(kt.items(), key=itemgetter(1), reverse=True)
for key, value in katz_sorted[0:10]:
print "\t > ", key, round(value, 4)
return g, kt Example 27
def pad_cpu2gpu(x, sz, offset=(0,0), dtype='real'):
block_size = (16, 16 ,1)
grid_size = (int(np.ceil(np.float32(sz[1])/block_size[1])),
int(np.ceil(np.float32(sz[0])/block_size[0])))
sx = x.shape
if x.__class__ == np.ndarray:
x = np.array(x).astype(np.float32)
x_gpu = cua.to_gpu(x)
elif x.__class__ == cua.GPUArray:
x_gpu = x
if dtype == 'real':
mod = cu.module_from_buffer(cubin)
zeroPadKernel = mod.get_function("zeroPadKernel")
x_padded_gpu = cua.zeros(tuple((int(sz[0]),int(sz[1]))), np.float32)
zeroPadKernel(x_padded_gpu.gpudata, np.int32(sz[0]), np.int32(sz[1]),
x_gpu.gpudata, np.int32(sx[0]), np.int32(sx[1]),
np.int32(offset[0]), np.int32(offset[1]),
block=block_size, grid=grid_size)
elif dtype == 'complex':
mod = cu.module_from_buffer(cubin)
#mod = SourceModule(open('gputools.cu').read(), keep=True)
zeroPadComplexKernel = mod.get_function("zeroPadComplexKernel")
x_padded_gpu = cua.zeros(tuple((int(sz[0]),int(sz[1]))), np.complex64)
zeroPadComplexKernel(x_padded_gpu.gpudata, np.int32(sz[0]), np.int32(sz[1]),
x_gpu.gpudata, np.int32(sx[0]), np.int32(sx[1]),
np.int32(offset[0]), np.int32(offset[1]),
block=block_size, grid=grid_size)
return x_padded_gpu Example 28
def cnvinv_gradfun(self, z, sz, y_gpu, alpha=0., beta=0.):
"""
Computes gradient used for 'lbfgsb' mode of deconv method.
See deconv for details.
"""
if z.__class__ == np.ndarray:
z = np.array(np.reshape(z,sz)).astype(np.float32)
z_gpu = cua.to_gpu(z)
grad_gpu = self.cnvtp(self.res_gpu)
# Thikonov regularization
# alpha > 0: Thikonov on the gradient of z
if alpha > 0:
grad_gpu += alpha * self.lz_gpu
# beta > 0: Thikonov on z
if beta > 0:
grad_gpu += beta * z_gpu
grad = -np.real(grad_gpu.get())
grad = grad.flatten()
return grad.astype(np.float64) Example 29
def psfScale(D, wavelength, pixSize):
"""
Return the PSF scale appropriate for the required pixel size, wavelength and telescope diameter
The aperture is padded by this amount; resultant pix scale is lambda/D/psf_scale, so for instance full frame 256 pix
for 3.5 m at 532 nm is 256*5.32e-7/3.5/3 = 2.67 arcsec for psf_scale = 3
Args:
D (real): telescope diameter in m
wavelength (real): wavelength in Angstrom
pixSize (real): pixel size in arcsec
Returns:
real: psf scale
"""
DInCm = D * 100.0
wavelengthInCm = wavelength * 1e-8
return 206265.0 * wavelengthInCm / (DInCm * pixSize) Example 30
def fwd_filter(img, S):
img_w, img_h, ch = img.shape
F = pad2(img)
F.data[F.mask] = 0. # make sure its zero-filled!
# Forward transform
specF = np.fft.fft2(F.data.astype(float), axes=(0, 1))
specN = np.fft.fft2(1. - F.mask.astype(float), axes=(0, 1))
specS = np.fft.fft2(S[::-1, ::-1])
out = np.real(np.fft.ifft2(specF * specS[:, :, np.newaxis], axes=(0, 1)))
norm = np.real(np.fft.ifft2(specN * specS[:, :, np.newaxis], axes=(0, 1)))
eps = 1e-15
norm = np.maximum(norm, eps)
out /= norm
out = out[-img_w:, -img_h:]
out[img.mask] = 0.
return np.ma.MaskedArray(data=out, mask=img.mask) Example 31
def kernel_impute(img, S):
F = pad2(img)
F.data[F.mask] = 0. # make sure its zero-filled!
img_w, img_h, img_ch = img.shape
Q = S
specF = np.fft.fft2(F.data.astype(float), axes=(0, 1))
specN = np.fft.fft2(1. - F.mask.astype(float), axes=(0, 1))
specQ = np.fft.fft2(Q[::-1, ::-1])
numer = np.real(np.fft.ifft2(specF * specQ[:, :, np.newaxis], axes=(0, 1)))
denom = np.real(np.fft.ifft2(specN * specQ[:, :, np.newaxis], axes=(0, 1)))
eps = 1e-15
fill = numer/(denom+eps)
fill = fill[-img_w:, -img_h:]
image = img.data.copy()
# img = img.copy()
image[img.mask] = fill[img.mask]
mask = np.zeros_like(img.mask, dtype=bool)
return np.ma.MaskedArray(data=image, mask=mask) Example 32
def get_neg_log_post(Phi, sigma_J_list, ROI_list, G, MMT, q, Sigma_E, GL,
nu, V, prior_on = False):
eps = 1E-13
p = Phi.shape[0]
n_ROI = len(sigma_J_list)
Qu = Phi.dot(Phi.T)
G_Sigma_G = np.zeros(MMT.shape)
for i in range(n_ROI):
G_Sigma_G += sigma_J_list[i]**2 * np.dot(G[:,ROI_list[i]], G[:,ROI_list[i]].T)
cov = Sigma_E + G_Sigma_G + GL.dot(Qu).dot(GL.T)
inv_cov = np.linalg.inv(cov)
eigs = np.real(np.linalg.eigvals(cov)) + eps
log_det_cov = np.sum(np.log(eigs))
result = q*log_det_cov + np.trace(MMT.dot(inv_cov))
if prior_on:
inv_Q = np.linalg.inv(Qu)
#det_Q = np.linalg.det(Qu)
log_det_Q = np.sum(np.log(np.diag(Phi)**2))
result = result + np.float(nu+p+1)*log_det_Q+ np.trace(V.dot(inv_Q))
return result
#==============================================================================
# update both Qu and Sigma_J, gradient of Qu and Sigma J Example 33
def get_neg_log_post(Phi, sigma_J_list, ROI_list, G, MMT, q, Sigma_E, GL,
nu, V, prior_on = False):
eps = 1E-13
p = Phi.shape[0]
n_ROI = len(sigma_J_list)
Qu = Phi.dot(Phi.T)
G_Sigma_G = np.zeros(MMT.shape)
for i in range(n_ROI):
G_Sigma_G += sigma_J_list[i]**2 * np.dot(G[:,ROI_list[i]], G[:,ROI_list[i]].T)
cov = Sigma_E + G_Sigma_G + GL.dot(Qu).dot(GL.T)
inv_cov = np.linalg.inv(cov)
eigs = np.real(np.linalg.eigvals(cov)) + eps
log_det_cov = np.sum(np.log(eigs))
result = q*log_det_cov + np.trace(MMT.dot(inv_cov))
if prior_on:
inv_Q = np.linalg.inv(Qu)
#det_Q = np.linalg.det(Qu)
log_det_Q = np.sum(np.log(np.diag(Phi)**2))
result = result + np.float(nu+p+1)*log_det_Q+ np.trace(V.dot(inv_Q))
return result
#==============================================================================
# update both Qu and Sigma_J, gradient of Qu and Sigma J Example 34
def __init__(self, qubit_names, quad="real"):
super(PulseCalibration, self).__init__()
self.qubit_names = qubit_names if isinstance(qubit_names, list) else [qubit_names]
self.qubit = [QubitFactory(qubit_name) for qubit_name in qubit_names] if isinstance(qubit_names, list) else QubitFactory(qubit_names)
self.filename = 'None'
self.exp = None
self.axis_descriptor = None
self.cw_mode = False
self.saved_settings = config.load_meas_file(config.meas_file)
self.settings = deepcopy(self.saved_settings) #make a copy for used during calibration
self.quad = quad
if quad == "real":
self.quad_fun = np.real
elif quad == "imag":
self.quad_fun = np.imag
elif quad == "amp":
self.quad_fun = np.abs
elif quad == "phase":
self.quad_fun = np.angle
else:
raise ValueError('Quadrature to calibrate must be one of ("real", "imag", "amp", "phase").')
self.plot = self.init_plot() Example 35
def find_null_offset(xpts, powers, default=0.0):
"""Finds the offset corresponding to the minimum power using a fit to the measured data"""
def model(x, a, b, c):
return a*(x - b)**2 + c
powers = np.power(10, powers/10.)
min_idx = np.argmin(powers)
try:
fit = curve_fit(model, xpts, powers, p0=[1, xpts[min_idx], powers[min_idx]])
except RuntimeError:
logger.warning("Mixer null offset fit failed.")
return default, np.zeros(len(powers))
best_offset = np.real(fit[0][1])
best_offset = np.minimum(best_offset, xpts[-1])
best_offset = np.maximum(best_offset, xpts[0])
xpts_fine = np.linspace(xpts[0],xpts[-1],101)
fit_pts = np.array([np.real(model(x, *fit[0])) for x in xpts_fine])
if min(fit_pts)<0: fit_pts-=min(fit_pts)-1e-10 #prevent log of a negative number
return best_offset, xpts_fine, 10*np.log10(fit_pts) Example 36
def make_layout(self):
self.lay = QtWidgets.QHBoxLayout()
self.lay.setContentsMargins(0, 0, 0, 0)
self.real = FloatSpinBox(label=self.labeltext,
min=self.minimum,
max=self.maximum,
increment=self.singleStep,
log_increment=self.log_increment,
halflife_seconds=self.halflife_seconds,
decimals=self.decimals)
self.imag = FloatSpinBox(label=self.labeltext,
min=self.minimum,
max=self.maximum,
increment=self.singleStep,
log_increment=self.log_increment,
halflife_seconds=self.halflife_seconds,
decimals=self.decimals)
self.real.value_changed.connect(self.value_changed)
self.lay.addWidget(self.real)
self.label = QtWidgets.QLabel(" + j")
self.lay.addWidget(self.label)
self.imag.value_changed.connect(self.value_changed)
self.lay.addWidget(self.imag)
self.setLayout(self.lay)
self.setFocusPolicy(QtCore.Qt.ClickFocus) Example 37
def save_image(imager, grid_data, grid_norm, output_file):
"""Makes an image from gridded visibilities and saves it to a FITS file.
Args:
imager (oskar.Imager): Handle to configured imager.
grid_data (numpy.ndarray): Final visibility grid.
grid_norm (float): Grid normalisation to apply.
output_file (str): Name of output FITS file to write.
"""
# Make the image (take the FFT, normalise, and apply grid correction).
imager.finalise_plane(grid_data, grid_norm)
grid_data = numpy.real(grid_data)
# Trim the image if required.
border = (imager.plane_size - imager.image_size) // 2
if border > 0:
end = border + imager.image_size
grid_data = grid_data[border:end, border:end]
# Write the FITS file.
hdr = fits.header.Header()
fits.writeto(output_file, grid_data, hdr, clobber=True) Example 38
def show_image(im: Image, fig=None, title: str = '', pol=0, chan=0, cm='rainbow'):
""" Show an Image with coordinates using matplotlib
:param im:
:param fig:
:param title:
:return:
"""
import matplotlib.pyplot as plt
assert isinstance(im, Image)
if not fig:
fig = plt.figure()
plt.clf()
fig.add_subplot(111, projection=im.wcs.sub(['longitude', 'latitude']))
if len(im.data.shape) == 4:
plt.imshow(numpy.real(im.data[chan, pol, :, :]), origin='lower', cmap=cm)
elif len(im.data.shape) == 2:
plt.imshow(numpy.real(im.data[:, :]), origin='lower', cmap=cm)
plt.xlabel('RA---SIN')
plt.ylabel('DEC--SIN')
plt.title(title)
plt.colorbar()
return fig Example 39
def convolve_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: Twice convolved image [nscales, nscales, nx, ny]
"""
nscales, nx, ny = scalestack.shape
convolved_shape = [nscales, nscales, nx, ny]
convolved = numpy.zeros(convolved_shape)
ximg = numpy.fft.fftshift(numpy.fft.fft2(numpy.fft.fftshift(img)))
xscaleshape = [nscales, nx, ny]
xscale = numpy.zeros(xscaleshape, dtype='complex')
for s in range(nscales):
xscale[s] = numpy.fft.fftshift(numpy.fft.fft2(numpy.fft.fftshift(scalestack[s, ...])))
for s in range(nscales):
for p in range(nscales):
xmult = ximg * xscale[p] * numpy.conjugate(xscale[s])
convolved[s, p, ...] = numpy.real(numpy.fft.ifftshift(numpy.fft.ifft2(numpy.fft.ifftshift(xmult))))
return convolved Example 40
def plot_waveforms(waveforms, figTitle=''):
channels = waveforms.keys()
# plot
plots = []
for (ct, chan) in enumerate(channels):
fig = bk.figure(title=figTitle + repr(chan),
plot_width=800,
plot_height=350,
y_range=[-1.05, 1.05],
x_axis_label=u'Time (?s)')
fig.background_fill_color = config.plotBackground
if config.gridColor:
fig.xgrid.grid_line_color = config.gridColor
fig.ygrid.grid_line_color = config.gridColor
waveformToPlot = waveforms[chan]
xpts = np.linspace(0, len(waveformToPlot) / chan.phys_chan.sampling_rate
/ 1e-6, len(waveformToPlot))
fig.line(xpts, np.real(waveformToPlot), color='red')
fig.line(xpts, np.imag(waveformToPlot), color='blue')
plots.append(fig)
bk.show(column(*plots)) Example 41
def merge_waveform(n, chAB, chAm1, chAm2, chBm1, chBm2):
'''
Builds packed I and Q waveforms from the nth mini LL, merging in marker data.
'''
wfAB = np.array([], dtype=np.complex)
for entry in chAB['linkList'][n % len(chAB['linkList'])]:
if not entry.isTimeAmp:
wfAB = np.append(wfAB, chAB['wfLib'][entry.key])
else:
wfAB = np.append(wfAB, chAB['wfLib'][entry.key][0] *
np.ones(entry.length * entry.repeat))
wfAm1 = marker_waveform(chAm1['linkList'][n % len(chAm1['linkList'])],
chAm1['wfLib'])
wfAm2 = marker_waveform(chAm2['linkList'][n % len(chAm2['linkList'])],
chAm2['wfLib'])
wfBm1 = marker_waveform(chBm1['linkList'][n % len(chBm1['linkList'])],
chBm1['wfLib'])
wfBm2 = marker_waveform(chBm2['linkList'][n % len(chBm2['linkList'])],
chBm2['wfLib'])
wfA = pack_waveform(np.real(wfAB), wfAm1, wfAm2)
wfB = pack_waveform(np.imag(wfAB), wfBm1, wfBm2)
return wfA, wfB Example 42
def tune_everything(x0squared, C, T, gmin, gmax):
# First tune based on dynamic range
if C==0:
dr=gmax/gmin
mustar=((np.sqrt(dr)-1)/(np.sqrt(dr)+1))**2
alpha_star = (1+np.sqrt(mustar))**2/gmax
return alpha_star,mustar
dist_to_opt = x0squared
grad_var = C
max_curv = gmax
min_curv = gmin
const_fact = dist_to_opt * min_curv**2 / 2 / grad_var
coef = [-1, 3, -(3 + const_fact), 1]
roots = np.roots(coef)
roots = roots[np.real(roots) > 0]
roots = roots[np.real(roots) < 1]
root = roots[np.argmin(np.imag(roots) ) ]
assert root > 0 and root < 1 and np.absolute(root.imag) < 1e-6
dr = max_curv / min_curv
assert max_curv >= min_curv
mu = max( ( (np.sqrt(dr) - 1) / (np.sqrt(dr) + 1) )**2, root**2)
lr_min = (1 - np.sqrt(mu) )**2 / min_curv
lr_max = (1 + np.sqrt(mu) )**2 / max_curv
alpha_star = lr_min
mustar = mu
return alpha_star, mustar Example 43
def reflection_from_matrix(matrix):
"""Return mirror plane point and normal vector from reflection matrix.
>>> v0 = numpy.random.random(3) - 0.5
>>> v1 = numpy.random.random(3) - 0.5
>>> M0 = reflection_matrix(v0, v1)
>>> point, normal = reflection_from_matrix(M0)
>>> M1 = reflection_matrix(point, normal)
>>> is_same_transform(M0, M1)
True
"""
M = numpy.array(matrix, dtype=numpy.float64, copy=False)
# normal: unit eigenvector corresponding to eigenvalue -1
l, V = numpy.linalg.eig(M[:3, :3])
i = numpy.where(abs(numpy.real(l) + 1.0) < 1e-8)[0]
if not len(i):
raise ValueError("no unit eigenvector corresponding to eigenvalue -1")
normal = numpy.real(V[:, i[0]]).squeeze()
# point: any unit eigenvector corresponding to eigenvalue 1
l, V = numpy.linalg.eig(M)
i = numpy.where(abs(numpy.real(l) - 1.0) < 1e-8)[0]
if not len(i):
raise ValueError("no unit eigenvector corresponding to eigenvalue 1")
point = numpy.real(V[:, i[-1]]).squeeze()
point /= point[3]
return point, normal Example 44
def rotation_from_matrix(matrix):
"""Return rotation angle and axis from rotation matrix.
>>> angle = (random.random() - 0.5) * (2*math.pi)
>>> direc = numpy.random.random(3) - 0.5
>>> point = numpy.random.random(3) - 0.5
>>> R0 = rotation_matrix(angle, direc, point)
>>> angle, direc, point = rotation_from_matrix(R0)
>>> R1 = rotation_matrix(angle, direc, point)
>>> is_same_transform(R0, R1)
True
"""
R = numpy.array(matrix, dtype=numpy.float64, copy=False)
R33 = R[:3, :3]
# direction: unit eigenvector of R33 corresponding to eigenvalue of 1
l, W = numpy.linalg.eig(R33.T)
i = numpy.where(abs(numpy.real(l) - 1.0) < 1e-8)[0]
if not len(i):
raise ValueError("no unit eigenvector corresponding to eigenvalue 1")
direction = numpy.real(W[:, i[-1]]).squeeze()
# point: unit eigenvector of R33 corresponding to eigenvalue of 1
l, Q = numpy.linalg.eig(R)
i = numpy.where(abs(numpy.real(l) - 1.0) < 1e-8)[0]
if not len(i):
raise ValueError("no unit eigenvector corresponding to eigenvalue 1")
point = numpy.real(Q[:, i[-1]]).squeeze()
point /= point[3]
# rotation angle depending on direction
cosa = (numpy.trace(R33) - 1.0) / 2.0
if abs(direction[2]) > 1e-8:
sina = (R[1, 0] + (cosa-1.0)*direction[0]*direction[1]) / direction[2]
elif abs(direction[1]) > 1e-8:
sina = (R[0, 2] + (cosa-1.0)*direction[0]*direction[2]) / direction[1]
else:
sina = (R[2, 1] + (cosa-1.0)*direction[1]*direction[2]) / direction[0]
angle = math.atan2(sina, cosa)
return angle, direction, point Example 45
def scale_from_matrix(matrix):
"""Return scaling factor, origin and direction from scaling matrix.
>>> factor = random.random() * 10 - 5
>>> origin = numpy.random.random(3) - 0.5
>>> direct = numpy.random.random(3) - 0.5
>>> S0 = scale_matrix(factor, origin)
>>> factor, origin, direction = scale_from_matrix(S0)
>>> S1 = scale_matrix(factor, origin, direction)
>>> is_same_transform(S0, S1)
True
>>> S0 = scale_matrix(factor, origin, direct)
>>> factor, origin, direction = scale_from_matrix(S0)
>>> S1 = scale_matrix(factor, origin, direction)
>>> is_same_transform(S0, S1)
True
"""
M = numpy.array(matrix, dtype=numpy.float64, copy=False)
M33 = M[:3, :3]
factor = numpy.trace(M33) - 2.0
try:
# direction: unit eigenvector corresponding to eigenvalue factor
l, V = numpy.linalg.eig(M33)
i = numpy.where(abs(numpy.real(l) - factor) < 1e-8)[0][0]
direction = numpy.real(V[:, i]).squeeze()
direction /= vector_norm(direction)
except IndexError:
# uniform scaling
factor = (factor + 2.0) / 3.0
direction = None
# origin: any eigenvector corresponding to eigenvalue 1
l, V = numpy.linalg.eig(M)
i = numpy.where(abs(numpy.real(l) - 1.0) < 1e-8)[0]
if not len(i):
raise ValueError("no eigenvector corresponding to eigenvalue 1")
origin = numpy.real(V[:, i[-1]]).squeeze()
origin /= origin[3]
return factor, origin, direction Example 46
def reflection_from_matrix(matrix):
"""Return mirror plane point and normal vector from reflection matrix.
>>> v0 = numpy.random.random(3) - 0.5
>>> v1 = numpy.random.random(3) - 0.5
>>> M0 = reflection_matrix(v0, v1)
>>> point, normal = reflection_from_matrix(M0)
>>> M1 = reflection_matrix(point, normal)
>>> is_same_transform(M0, M1)
True
"""
M = numpy.array(matrix, dtype=numpy.float64, copy=False)
# normal: unit eigenvector corresponding to eigenvalue -1
w, V = numpy.linalg.eig(M[:3, :3])
i = numpy.where(abs(numpy.real(w) + 1.0) < 1e-8)[0]
if not len(i):
raise ValueError("no unit eigenvector corresponding to eigenvalue -1")
normal = numpy.real(V[:, i[0]]).squeeze()
# point: any unit eigenvector corresponding to eigenvalue 1
w, V = numpy.linalg.eig(M)
i = numpy.where(abs(numpy.real(w) - 1.0) < 1e-8)[0]
if not len(i):
raise ValueError("no unit eigenvector corresponding to eigenvalue 1")
point = numpy.real(V[:, i[-1]]).squeeze()
point /= point[3]
return point, normal Example 47
def rotation_from_matrix(matrix):
"""Return rotation angle and axis from rotation matrix.
>>> angle = (random.random() - 0.5) * (2*math.pi)
>>> direc = numpy.random.random(3) - 0.5
>>> point = numpy.random.random(3) - 0.5
>>> R0 = rotation_matrix(angle, direc, point)
>>> angle, direc, point = rotation_from_matrix(R0)
>>> R1 = rotation_matrix(angle, direc, point)
>>> is_same_transform(R0, R1)
True
"""
R = numpy.array(matrix, dtype=numpy.float64, copy=False)
R33 = R[:3, :3]
# direction: unit eigenvector of R33 corresponding to eigenvalue of 1
w, W = numpy.linalg.eig(R33.T)
i = numpy.where(abs(numpy.real(w) - 1.0) < 1e-8)[0]
if not len(i):
raise ValueError("no unit eigenvector corresponding to eigenvalue 1")
direction = numpy.real(W[:, i[-1]]).squeeze()
# point: unit eigenvector of R33 corresponding to eigenvalue of 1
w, Q = numpy.linalg.eig(R)
i = numpy.where(abs(numpy.real(w) - 1.0) < 1e-8)[0]
if not len(i):
raise ValueError("no unit eigenvector corresponding to eigenvalue 1")
point = numpy.real(Q[:, i[-1]]).squeeze()
point /= point[3]
# rotation angle depending on direction
cosa = (numpy.trace(R33) - 1.0) / 2.0
if abs(direction[2]) > 1e-8:
sina = (R[1, 0] + (cosa-1.0)*direction[0]*direction[1]) / direction[2]
elif abs(direction[1]) > 1e-8:
sina = (R[0, 2] + (cosa-1.0)*direction[0]*direction[2]) / direction[1]
else:
sina = (R[2, 1] + (cosa-1.0)*direction[1]*direction[2]) / direction[0]
angle = math.atan2(sina, cosa)
return angle, direction, point Example 48
def scale_from_matrix(matrix):
"""Return scaling factor, origin and direction from scaling matrix.
>>> factor = random.random() * 10 - 5
>>> origin = numpy.random.random(3) - 0.5
>>> direct = numpy.random.random(3) - 0.5
>>> S0 = scale_matrix(factor, origin)
>>> factor, origin, direction = scale_from_matrix(S0)
>>> S1 = scale_matrix(factor, origin, direction)
>>> is_same_transform(S0, S1)
True
>>> S0 = scale_matrix(factor, origin, direct)
>>> factor, origin, direction = scale_from_matrix(S0)
>>> S1 = scale_matrix(factor, origin, direction)
>>> is_same_transform(S0, S1)
True
"""
M = numpy.array(matrix, dtype=numpy.float64, copy=False)
M33 = M[:3, :3]
factor = numpy.trace(M33) - 2.0
try:
# direction: unit eigenvector corresponding to eigenvalue factor
w, V = numpy.linalg.eig(M33)
i = numpy.where(abs(numpy.real(w) - factor) < 1e-8)[0][0]
direction = numpy.real(V[:, i]).squeeze()
direction /= vector_norm(direction)
except IndexError:
# uniform scaling
factor = (factor + 2.0) / 3.0
direction = None
# origin: any eigenvector corresponding to eigenvalue 1
w, V = numpy.linalg.eig(M)
i = numpy.where(abs(numpy.real(w) - 1.0) < 1e-8)[0]
if not len(i):
raise ValueError("no eigenvector corresponding to eigenvalue 1")
origin = numpy.real(V[:, i[-1]]).squeeze()
origin /= origin[3]
return factor, origin, direction Example 49
def rotation_from_matrix(matrix):
"""Return rotation angle and axis from rotation matrix.
>>> angle = (random.random() - 0.5) * (2*math.pi)
>>> direc = numpy.random.random(3) - 0.5
>>> point = numpy.random.random(3) - 0.5
>>> R0 = rotation_matrix(angle, direc, point)
>>> angle, direc, point = rotation_from_matrix(R0)
>>> R1 = rotation_matrix(angle, direc, point)
>>> is_same_transform(R0, R1)
True
"""
R = numpy.array(matrix, dtype=numpy.float64, copy=False)
R33 = R[:3, :3]
# direction: unit eigenvector of R33 corresponding to eigenvalue of 1
w, W = numpy.linalg.eig(R33.T)
i = numpy.where(abs(numpy.real(w) - 1.0) < 1e-8)[0]
if not len(i):
raise ValueError("no unit eigenvector corresponding to eigenvalue 1")
direction = numpy.real(W[:, i[-1]]).squeeze()
# point: unit eigenvector of R33 corresponding to eigenvalue of 1
w, Q = numpy.linalg.eig(R)
i = numpy.where(abs(numpy.real(w) - 1.0) < 1e-8)[0]
if not len(i):
raise ValueError("no unit eigenvector corresponding to eigenvalue 1")
point = numpy.real(Q[:, i[-1]]).squeeze()
point /= point[3]
# rotation angle depending on direction
cosa = (numpy.trace(R33) - 1.0) / 2.0
if abs(direction[2]) > 1e-8:
sina = (R[1, 0] + (cosa-1.0)*direction[0]*direction[1]) / direction[2]
elif abs(direction[1]) > 1e-8:
sina = (R[0, 2] + (cosa-1.0)*direction[0]*direction[2]) / direction[1]
else:
sina = (R[2, 1] + (cosa-1.0)*direction[1]*direction[2]) / direction[0]
angle = math.atan2(sina, cosa)
return angle, direction, point
# Function to translate handshape coding to degrees of rotation, adduction, flexion Example 50
def rotation_from_matrix(matrix):
"""Return rotation angle and axis from rotation matrix.
>>> angle = (random.random() - 0.5) * (2*math.pi)
>>> direc = numpy.random.random(3) - 0.5
>>> point = numpy.random.random(3) - 0.5
>>> R0 = rotation_matrix(angle, direc, point)
>>> angle, direc, point = rotation_from_matrix(R0)
>>> R1 = rotation_matrix(angle, direc, point)
>>> is_same_transform(R0, R1)
True
"""
R = numpy.array(matrix, dtype=numpy.float64, copy=False)
R33 = R[:3, :3]
# direction: unit eigenvector of R33 corresponding to eigenvalue of 1
w, W = numpy.linalg.eig(R33.T)
i = numpy.where(abs(numpy.real(w) - 1.0) < 1e-8)[0]
if not len(i):
raise ValueError("no unit eigenvector corresponding to eigenvalue 1")
direction = numpy.real(W[:, i[-1]]).squeeze()
# point: unit eigenvector of R33 corresponding to eigenvalue of 1
w, Q = numpy.linalg.eig(R)
i = numpy.where(abs(numpy.real(w) - 1.0) < 1e-8)[0]
if not len(i):
raise ValueError("no unit eigenvector corresponding to eigenvalue 1")
point = numpy.real(Q[:, i[-1]]).squeeze()
point /= point[3]
# rotation angle depending on direction
cosa = (numpy.trace(R33) - 1.0) / 2.0
if abs(direction[2]) > 1e-8:
sina = (R[1, 0] + (cosa-1.0)*direction[0]*direction[1]) / direction[2]
elif abs(direction[1]) > 1e-8:
sina = (R[0, 2] + (cosa-1.0)*direction[0]*direction[2]) / direction[1]
else:
sina = (R[2, 1] + (cosa-1.0)*direction[1]*direction[2]) / direction[0]
angle = math.atan2(sina, cosa)
return angle, direction, point
# Function to translate handshape coding to degrees of rotation, adduction, flexion 