我们从Python开源项目中,提取了以下50个代码示例,用于说明如何使用numpy.real()。
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])
def quaternion_real(quaternion): """Return real part of quaternion. >>> quaternion_real([3, 0, 1, 2]) 3.0 """ return float(quaternion[0])
def _ncc_c(x, y): """ >>> _ncc_c([1,2,3,4], [1,2,3,4]) array([ 0.13333333, 0.36666667, 0.66666667, 1. , 0.66666667, 0.36666667, 0.13333333]) >>> _ncc_c([1,1,1], [1,1,1]) array([ 0.33333333, 0.66666667, 1. , 0.66666667, 0.33333333]) >>> _ncc_c([1,2,3], [-1,-1,-1]) array([-0.15430335, -0.46291005, -0.9258201 , -0.77151675, -0.46291005]) """ den = np.array(norm(x) * norm(y)) den[den == 0] = np.Inf x_len = len(x) fft_size = 1<<(2*x_len-1).bit_length() cc = ifft(fft(x, fft_size) * np.conj(fft(y, fft_size))) cc = np.concatenate((cc[-(x_len-1):], cc[:x_len])) return np.real(cc) / den
def nufft_scale1(N, K, alpha, beta, Nmid): ''' calculate image space scaling factor ''' # import types # if alpha is types.ComplexType: alpha = numpy.real(alpha) # print('complex alpha may not work, but I just let it as') L = len(alpha) - 1 if L > 0: sn = numpy.zeros((N, 1)) n = numpy.arange(0, N).reshape((N, 1), order='F') i_gam_n_n0 = 1j * (2 * numpy.pi / K) * (n - Nmid) * beta for l1 in range(-L, L + 1): alf = alpha[abs(l1)] if l1 < 0: alf = numpy.conj(alf) sn = sn + alf * numpy.exp(i_gam_n_n0 * l1) else: sn = numpy.dot(alpha, numpy.ones((N, 1), dtype=numpy.float32)) return sn
def 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
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)
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
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 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)
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))
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
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 ] )
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)
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)
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)
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)
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)
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)
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])
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 fft_convolve(X,Y, inv = 0): XF = np.fft.rfft2(X) YF = np.fft.rfft2(Y) # YF0 = np.copy(YF) # YF.imag = 0 # XF.imag = 0 if inv == 1: # plt.imshow(np.real(YF)); plt.colorbar(); plt.show() YF = np.conj(YF) SF = XF*YF S = np.fft.irfft2(SF) n1,n2 = np.shape(S) S = np.roll(S,-n1/2+1,axis = 0) S = np.roll(S,-n2/2+1,axis = 1) return np.real(S)
def __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)
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
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)
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)
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)
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)
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
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()
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)
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)
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)
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
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
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))
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
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
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
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
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
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
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
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
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
def nufft_alpha_kb_fit(N, J, K): ''' find out parameters alpha and beta of scaling factor st['sn'] Note, when J = 1 , alpha is hardwired as [1,0,0...] (uniform scaling factor) ''' beta = 1 Nmid = (N - 1.0) / 2.0 if N > 40: L = 13 else: L = numpy.ceil(N / 3) nlist = numpy.arange(0, N) * 1.0 - Nmid (kb_a, kb_m) = kaiser_bessel('string', J, 'best', 0, K / N) if J > 1: sn_kaiser = 1 / kaiser_bessel_ft(nlist / K, J, kb_a, kb_m, 1.0) elif J == 1: # The case when samples are on regular grids sn_kaiser = numpy.ones((1, N), dtype=dtype) gam = 2 * numpy.pi / K X_ant = beta * gam * nlist.reshape((N, 1), order='F') X_post = numpy.arange(0, L + 1) X_post = X_post.reshape((1, L + 1), order='F') X = numpy.dot(X_ant, X_post) # [N,L] X = numpy.cos(X) sn_kaiser = sn_kaiser.reshape((N, 1), order='F').conj() X = numpy.array(X, dtype=dtype) sn_kaiser = numpy.array(sn_kaiser, dtype=dtype) coef = numpy.linalg.lstsq(numpy.nan_to_num(X), numpy.nan_to_num(sn_kaiser))[0] # (X \ sn_kaiser.H); alphas = coef if J > 1: alphas[0] = alphas[0] alphas[1:] = alphas[1:] / 2.0 elif J == 1: # cases on grids alphas[0] = 1.0 alphas[1:] = 0.0 alphas = numpy.real(alphas) return (alphas, beta)
def init_bh_tsne(samples, workdir, no_dims=DEFAULT_NO_DIMS, initial_dims=INITIAL_DIMENSIONS, perplexity=DEFAULT_PERPLEXITY, theta=DEFAULT_THETA, randseed=EMPTY_SEED, verbose=False, use_pca=DEFAULT_USE_PCA, max_iter=DEFAULT_MAX_ITERATIONS): if use_pca: samples = samples - np.mean(samples, axis=0) cov_x = np.dot(np.transpose(samples), samples) [eig_val, eig_vec] = np.linalg.eig(cov_x) # sorting the eigen-values in the descending order eig_vec = eig_vec[:, eig_val.argsort()[::-1]] if initial_dims > len(eig_vec): initial_dims = len(eig_vec) # truncating the eigen-vectors matrix to keep the most important vectors eig_vec = np.real(eig_vec[:, :initial_dims]) samples = np.dot(samples, eig_vec) # Assume that the dimensionality of the first sample is representative for # the whole batch sample_dim = len(samples[0]) sample_count = len(samples) # Note: The binary format used by bh_tsne is roughly the same as for # vanilla tsne with open(path_join(workdir, 'data.dat'), 'wb') as data_file: # Write the bh_tsne header data_file.write(pack('iiddii', sample_count, sample_dim, theta, perplexity, no_dims, max_iter)) # Then write the data for sample in samples: data_file.write(pack('{}d'.format(len(sample)), *sample)) # Write random seed if specified if randseed != EMPTY_SEED: data_file.write(pack('i', randseed))