python类real()的实例源码

encode.py 文件源码 项目:blind-watermark 作者: linyacool 项目源码 文件源码 阅读 36 收藏 0 点赞 0 评论 0
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])
transformations.py 文件源码 项目:Neural-Networks-for-Inverse-Kinematics 作者: paramrajpura 项目源码 文件源码 阅读 47 收藏 0 点赞 0 评论 0
def quaternion_real(quaternion):
    """Return real part of quaternion.

    >>> quaternion_real([3, 0, 1, 2])
    3.0

    """
    return float(quaternion[0])
kshape.py 文件源码 项目:rca-evaluation 作者: sieve-microservices 项目源码 文件源码 阅读 30 收藏 0 点赞 0 评论 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
pynufft_gpu.py 文件源码 项目:pynufft 作者: jyhmiinlin 项目源码 文件源码 阅读 36 收藏 0 点赞 0 评论 0
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
helper.py 文件源码 项目:pynufft 作者: jyhmiinlin 项目源码 文件源码 阅读 26 收藏 0 点赞 0 评论 0
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
transforms.py 文件源码 项目:mdct 作者: nils-werner 项目源码 文件源码 阅读 37 收藏 0 点赞 0 评论 0
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)
ImageFFT_class.py 文件源码 项目:DenoiseAverage 作者: Pella86 项目源码 文件源码 阅读 25 收藏 0 点赞 0 评论 0
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
ImageFFT_class.py 文件源码 项目:DenoiseAverage 作者: Pella86 项目源码 文件源码 阅读 40 收藏 0 点赞 0 评论 0
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
sepfilter.py 文件源码 项目:pyrsss 作者: butala 项目源码 文件源码 阅读 37 收藏 0 点赞 0 评论 0
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)
mparray.py 文件源码 项目:mpnum 作者: dseuss 项目源码 文件源码 阅读 47 收藏 0 点赞 0 评论 0
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))
TimeSignals.py 文件源码 项目:droppy 作者: BV-DR 项目源码 文件源码 阅读 49 收藏 0 点赞 0 评论 0
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
TimeSignals.py 文件源码 项目:droppy 作者: BV-DR 项目源码 文件源码 阅读 52 收藏 0 点赞 0 评论 0
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 ]  )
qerbt.py 文件源码 项目:untwist 作者: IoSR-Surrey 项目源码 文件源码 阅读 48 收藏 0 点赞 0 评论 0
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)
test_fourier.py 文件源码 项目:dft 作者: rosenbrockc 项目源码 文件源码 阅读 44 收藏 0 点赞 0 评论 0
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)
fourier.py 文件源码 项目:dft 作者: rosenbrockc 项目源码 文件源码 阅读 27 收藏 0 点赞 0 评论 0
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)
fourier.py 文件源码 项目:dft 作者: rosenbrockc 项目源码 文件源码 阅读 27 收藏 0 点赞 0 评论 0
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)
fourier.py 文件源码 项目:dft 作者: rosenbrockc 项目源码 文件源码 阅读 30 收藏 0 点赞 0 评论 0
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)
InductionSphereFEM.py 文件源码 项目:em_examples 作者: geoscixyz 项目源码 文件源码 阅读 37 收藏 0 点赞 0 评论 0
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)
decode.py 文件源码 项目:blind-watermark 作者: linyacool 项目源码 文件源码 阅读 34 收藏 0 点赞 0 评论 0
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])
spectrum.py 文件源码 项目:McMurchie-Davidson 作者: jjgoings 项目源码 文件源码 阅读 28 收藏 0 点赞 0 评论 0
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
wave_transform.py 文件源码 项目:semantic-segmentation 作者: albertbuchard 项目源码 文件源码 阅读 38 收藏 0 点赞 0 评论 0
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)
GaborJets.py 文件源码 项目:nimo 作者: wolfram2012 项目源码 文件源码 阅读 37 收藏 0 点赞 0 评论 0
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)
pedWorks.py 文件源码 项目:PedWorks 作者: BrnCPrz 项目源码 文件源码 阅读 27 收藏 0 点赞 0 评论 0
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
olaGPU.py 文件源码 项目:slitSpectrographBlind 作者: aasensio 项目源码 文件源码 阅读 27 收藏 0 点赞 0 评论 0
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)
wavefront.py 文件源码 项目:slitSpectrographBlind 作者: aasensio 项目源码 文件源码 阅读 33 收藏 0 点赞 0 评论 0
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)
filtering.py 文件源码 项目:uncover-ml 作者: GeoscienceAustralia 项目源码 文件源码 阅读 53 收藏 0 点赞 0 评论 0
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)
filtering.py 文件源码 项目:uncover-ml 作者: GeoscienceAustralia 项目源码 文件源码 阅读 40 收藏 0 点赞 0 评论 0
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)
ROI_cov_Qu_SigmaJ.py 文件源码 项目:MEEG_connectivity 作者: YingYang 项目源码 文件源码 阅读 35 收藏 0 点赞 0 评论 0
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
ROI_cov_Qu_SigmaJ.py 文件源码 项目:MEEG_connectivity 作者: YingYang 项目源码 文件源码 阅读 37 收藏 0 点赞 0 评论 0
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
pulse_calibration.py 文件源码 项目:Auspex 作者: BBN-Q 项目源码 文件源码 阅读 27 收藏 0 点赞 0 评论 0
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()


问题


面经


文章

微信
公众号

扫码关注公众号