Python numpy.real函数代码示例

 def invert(self, estimate):
     """Invert the estimate to produce slopes.
     
     Parameters
     ----------
     estimate : array_like
         Phase estimate to invert.
     
     Returns
     -------
     xs : array_like
         Estimate of the x slopes.
     ys : array_like
         Estimate of the y slopes.
     
     
     """
     if self.manage_tt:
         estimate, ttx, tty = remove_tiptilt(self.ap, estimate)
     
     est_ft = fftpack.fftn(estimate) / 2.0
     
     xs_ft = self.gx * est_ft
     ys_ft = self.gy * est_ft
     
     xs = np.real(fftpack.ifftn(xs_ft))
     ys = np.real(fftpack.ifftn(ys_ft))
     
     if self.manage_tt and not self.suppress_tt:
         xs += ttx
         ys += tty
     
     return (xs, ys)

	def SaveData (self, fname, verbose = True):

		if (verbose):
			print ("  Saving measurement to %s ... " % fname)

		start = time.clock()

		f = h5py.File(fname, 'w')

		f['data_r'] = np.squeeze(np.real(self.data).transpose())
		f['data_i'] = np.squeeze(np.imag(self.data).transpose())

		if (self.noise != 0):
			f['noise_r'] = np.squeeze(np.real(self.noise).transpose())
			f['noise_i'] = np.squeeze(np.imag(self.noise).transpose())
		
		if (self.acs != 0):
			f['acs_r'] = np.squeeze(np.real(self.acs).transpose())
			f['acs_i'] = np.squeeze(np.imag(self.acs).transpose())
		
		if (self.sync.any() != 0):
                        f['sync']  = self.sync.transpose()

		f.close()

		if (verbose):
			print '    ... saved in %(time).1f s.\n' % {"time": time.clock()-start}

		return

def FFT_Correlation(x,y):
    """
    FFT-based correlation, much faster than numpy autocorr.
    x and y are row-based vectors of arbitrary lengths.
    This is a vectorized implementation of O(N*log(N)) flops.
    """

    lengthx = x.shape[0]
    lengthy = y.shape[0]

    x = np.reshape(x,(1,lengthx))
    y = np.reshape(y,(1,lengthy))

    length = np.array([lengthx, lengthy]).min()
    
    x = x[:length]
    y = y[:length]
    
    fftx = fft(x, 2 * length - 1, axis=1) #pad with zeros
    ffty = fft(y, 2 * length - 1, axis=1)

    corr_xy = fft.ifft(fftx * np.conjugate(ffty), axis=1)
    corr_xy = np.real(fft.fftshift(corr_xy, axes=1)) #should be no imaginary part

    corr_yx = fft.ifft(ffty * np.conjugate(fftx), axis=1)
    corr_yx = np.real(fft.fftshift(corr_yx, axes=1))

    corr = 0.5 * (corr_xy[:,length:] + corr_yx[:,length:]) / range(1,length)[::-1]
    return np.reshape(corr,corr.shape[1])

def update_all_baseline_plots(i, fig, crawler, lines, norm_cross=False, forward=True):

    if forward:
        try:
            crawler.forward()
        except EOFError as err:
            print err
            raw_input("End of File. Press enter to quit.")
            sys.exit()

    burst = crawler

    for k in range(len(BASELINES)):
        if k < 4:
            #autos
            lines[k].set_data(FREQS, 10*np.log10(burst.autos[BASELINES[k]]))
            #overlays
            lines[-(k+1)].set_data(FREQS,10*np.log10(burst.autos[BASELINES[k]]))

        elif norm_cross:

			norm_val = np.array(burst.cross[BASELINES[k]])/np.sqrt(np.array(burst.autos[BASELINES[k][0]*2])*np.array(burst.autos[BASELINES[k][1]*2]))
			lines[k]['real'].set_data(FREQS, np.real(norm_val))
			lines[k]['imag'].set_data(FREQS, np.imag(norm_val))
        else:
			lines[k].set_data(FREQS, 10*np.log10(np.abs(np.real(burst.cross[BASELINES[k]]))))



    return lines

    def _exportDataToText(self, file):
        """Saves textual data to a file

        """
        Nt = self.data.shape[0]
        N = self.data.shape[1]
        # all elements [real] + (all elements - diagonal) [imaginary] + time
        Na = N + 1 + N*(N-1) 

        out = numpy.zeros((Nt, Na), dtype=numpy.float64)   
        
        for i in range(Nt):
            #print("%%%%")
            # time
            out[i,0] = self.TimeAxis.data[i]
            #print(0)
            # populations
            for j in range(N):
               out[i,j+1] = numpy.real(self.data[i,j,j])
               #print(j+1)
            # coherences
            l = 0
            for j in range(N):
                for k in range(j+1,N):
                    out[i,N+1+l] = numpy.real(self.data[i,j,k])
                    #print(N+1+l)
                    l += 1
                    out[i,N+1+l] = numpy.imag(self.data[i,j,k])
                    #print(N+1+l)
                    l += 1
                    
        numpy.savetxt(file, out)

    def eta_fil(self, x, V_app, apprx=(0, 0, 0, 0)):
        m_eff = self.m_r * const.electron_mass

        mpmath.mp.dps = 20
        x0 = Symbol('x0')  # eta_fil
        x1 = Symbol('x1')  # eta_ac
        x2 = Symbol('x2')  # eta_hop
        x3 = Symbol('x3')  # V_tunnel

        f0 = const.Boltzmann * self.T / (1 - self.alpha) / const.elementary_charge / self.z * \
             ln(self.A_fil/self.A_ac*(exp(- self.alpha * const.elementary_charge * self.z / const.Boltzmann / self.T * x0) - 1) + 1) - x1# eta_ac = f(eta_fil) x1 = f(x0)
        f1 = x*2*const.Boltzmann*self.T/self.a/self.z/const.elementary_charge*\
             asinh(self.j_0et/self.j_0hop*(exp(- self.alpha * const.elementary_charge * self.z / const.Boltzmann / self.T * x0) - 1)) - x2# eta_hop = f(eta_fil)
        f2 = x1 - x0 + x2 - x3

        f3 = -V_app + ((self.C * 3 * sqrt(2 * m_eff * ((4+x3/2)*const.elementary_charge)) / 2 / x * (const.elementary_charge / const.Planck)**2 * \
             exp(- 4 * const.pi * x / const.Planck * sqrt(2 * m_eff * ((4+x3/2)*const.elementary_charge))) * self.A_fil*x3)
                       + (self.j_0et*self.A_fil*(exp(-self.alpha*const.elementary_charge*self.z*x0/const.Boltzmann/self.T) - 1))) * (self.R_el + self.R_S + self.rho_fil*(self.L - x) / self.A_fil) \
             + x3

        eta_fil, eta_ac, eta_hop, V_tunnel = nsolve((f0, f1, f2, f3), [x0, x1, x2, x3], apprx)
        eta_fil = np.real(np.complex128(eta_fil))
        eta_ac = np.real(np.complex128(eta_ac))
        eta_hop = np.real(np.complex128(eta_hop))
        V_tunnel = np.real(np.complex128(V_tunnel))
        current = ((self.C * 3 * sqrt(2 * m_eff * ((4+V_tunnel)*const.elementary_charge)) / 2 / x * (const.elementary_charge / const.Planck)**2 * \
            exp(- 4 * const.pi * x / const.Planck * sqrt(2 * m_eff * ((4+V_tunnel)*const.elementary_charge))) * self.A_fil*V_tunnel)
                       + (self.j_0et*self.A_fil*(exp(-self.alpha*const.elementary_charge*self.z*eta_fil/const.Boltzmann/self.T) - 1)))
        print(eta_fil, eta_ac, eta_hop, V_tunnel)
        # print(eta_ac - eta_fil + eta_hop - V_tunnel)
        return eta_fil, eta_ac, eta_hop, V_tunnel, current

def window_fn_matrix(Q,N,num_remov=None,save_tag=None,lms=None):
    Q = n.matrix(Q); N = n.matrix(N)
    Ninv = uf.pseudo_inverse(N,num_remov=None) # XXX want to remove dynamically
    #print Ninv 
    info = n.dot(Q.H,n.dot(Ninv,Q))
    M = uf.pseudo_inverse(info,num_remov=num_remov)
    W = n.dot(M,info)

    if save_tag!=None:
        foo = W[0,:]
        foo = n.real(n.array(foo))
        foo.shape = (foo.shape[1]),
        print foo.shape
        p.scatter(lms[:,0],foo,c=lms[:,1],cmap=mpl.cm.PiYG,s=50)
        p.xlabel('l (color is m)')
        p.ylabel('W_0,lm')
        p.title('First Row of Window Function Matrix')
        p.colorbar()
        p.savefig('{0}/{1}_W.pdf'.format(fig_loc,save_tag))
        p.clf()

        print 'W ',W.shape
        p.imshow(n.real(W))
        p.title('Window Function Matrix')
        p.colorbar()
        p.savefig('{0}/{1}_W_im.pdf'.format(fig_loc,save_tag))
        p.clf()


    return W

def _test_random(test_case,inarr,outarr,tol):
    tc = test_case
    # Test that applying a transform and its inverse to reasonably long, random
    # input gives back the (appropriately scaled) input. We must allow for numerical
    # error, and it seems more reliable to check using normwise error (than elementwise).
    #
    # First test IFFT(FFT(random))
    # The numpy randn(n) provides an array of n numbers drawn from standard normal
    if dtype(inarr).kind == 'c':
        inarr._data[:] = randn(len(inarr)) +1j*randn(len(inarr))
        # If we're going to do a HC2R transform we must worry about DC/Nyquist imaginary
        if dtype(outarr).kind == 'f':
            inarr._data[0] = real(inarr[0])
            if (len(outarr)%2)==0:
                inarr._data[len(inarr)-1] = real(inarr[len(inarr)-1])
    else:
        inarr._data[:] = randn(len(inarr))
    incopy = type(inarr)(inarr)
    outarr.clear()
    # An FFT followed by IFFT gives Array scaled by len(Array), but for
    # Time/FrequencySeries there should be no scaling.
    if type(inarr) == pycbc.types.Array:
        incopy *= len(inarr)
    with tc.context:
        pycbc.fft.fft(inarr, outarr)
        pycbc.fft.ifft(outarr, inarr)
        emsg="IFFT(FFT(random)) did not reproduce original array to within tolerance {0}".format(tol)
        if isinstance(incopy,ts) or isinstance(incopy,fs):
            tc.assertTrue(incopy.almost_equal_norm(inarr,tol=tol,dtol=tol),
                          msg=emsg)
        else:
            tc.assertTrue(incopy.almost_equal_norm(inarr,tol=tol),
                          msg=emsg)
    # Perform arithmetic on outarr and inarr to pull them off of the GPU:
    outarr *= 1.0
    inarr *= 1.0
    # Now the same for FFT(IFFT(random))
    if dtype(outarr).kind == 'c':
        outarr._data[:] = randn(len(outarr))+1j*randn(len(outarr))
        # If we're going to do a HC2R transform we must worry about DC/Nyquist imaginary
        if dtype(inarr).kind == 'f':
            outarr._data[0] = real(outarr[0])
            if (len(inarr)%2)==0:
                outarr._data[len(outarr)-1] = real(outarr[len(outarr)-1])
    else:
        outarr._data[:] = randn(len(outarr))
    inarr.clear()
    outcopy = type(outarr)(outarr)
    if type(outarr) == pycbc.types.Array:
        outcopy *= len(inarr)
    with tc.context:
        pycbc.fft.ifft(outarr, inarr)
        pycbc.fft.fft(inarr, outarr)
        emsg="FFT(IFFT(random)) did not reproduce original array to within tolerance {0}".format(tol)
        if isinstance(outcopy,ts) or isinstance(outcopy,fs):
            tc.assertTrue(outcopy.almost_equal_norm(outarr,tol=tol,dtol=tol),
                          msg=emsg)
        else:
            tc.assertTrue(outcopy.almost_equal_norm(outarr,tol=tol),
                          msg=emsg)

    def icd(self, G1, G2):
        """Incomplete Cholesky decomposition
        """
        
        # remove mean. avoid standard calculation N0 = eye(N)-1/N*ones(N);
        G1 = G1 - numpy.array(numpy.mean(G1, 0), ndmin=2, copy=False)
        G2 = G2 - numpy.array(numpy.mean(G2, 0), ndmin=2, copy=False)

        R, D = self.method(G1, G2, self.reg)
        
        #solve generalized eigenvalues problem
        betas, alphas = scipy.linalg.eig(R,D)
        ind = numpy.argsort(numpy.real(betas))
        max_ind = ind[-1]
        alpha = alphas[:, max_ind]
        alpha = alpha/numpy.linalg.norm(alpha)
        beta = numpy.real(betas[max_ind])
        
        N1 = G1.shape[1]
        alpha1 = alpha[:N1]
        alpha2 = alpha[N1:]
        
        y1 = dot(G1, alpha1)
        y2 = dot(G2, alpha2)

        self.alpha1 = alpha1
        self.alpha2 = alpha2
        
        return (y1, y2, beta)                

def closest_point_on_the_parabola(a, px, py):
    thingy1 = 2. * a * py
    thingy2 = np.sqrt(-3 + 0j)
    thingy3 = 2 ** (1 / 3.)
    thingy4 = 2 ** (2 / 3.)
    thingy = (-108. * a ** 4 * px + np.sqrt(11664. * a ** 8 * px ** 2 - 864. * a ** 6 * (-1 + thingy1) ** 3 + 0j)) ** (
    1 / 3.)
    Aone = (thingy3 * (-1. + thingy1))
    Atwo = thingy
    Athree = thingy
    Afour = (6. * thingy3 * a ** 2)
    Bone = ((1. + thingy2) * (-1. + thingy1))
    Btwo = (thingy4 * thingy)
    Bthree = ((1. - thingy2) * thingy)
    Bfour = (12. * thingy3 * a ** 2)
    Cone = (1. - thingy2) * (-1 + thingy1)
    Ctwo = thingy4 * thingy
    Cthree = (1. + thingy2) * thingy
    Cfour = 12. * thingy3 * a ** 2

    A = -np.real(Aone / Atwo + Athree / Afour)
    B = np.real(Bone / Btwo + Bthree / Bfour)
    C = np.real(Cone / Ctwo + Cthree / Cfour)

    solns = [A, B, C]
    solns_temp = []
    for soln in solns:
        solns_temp.append(np.abs(soln - px))

    val, idx = min((val, idx) for (idx, val) in enumerate(solns_temp))
    return solns[idx]

    def simulate(self, steps, time, collect_dynamic = False):
        """!
        @brief Performs static simulation of oscillatory network.
        
        @param[in] steps (uint): Number simulation steps.
        @param[in] time (double): Time of simulation.
        @param[in] collect_dynamic (bool): If True - returns whole dynamic of oscillatory network, otherwise returns only last values of dynamics.
        
        @return (list) Dynamic of oscillatory network. If argument 'collect_dynamic' is True, than return dynamic for the whole simulation time,
                 otherwise returns only last values (last step of simulation) of output dynamic.
        
        @see simulate()
        @see simulate_dynamic()
        
        """
        
        dynamic_amplitude, dynamic_time = ([], []) if collect_dynamic is False else ([self.__amplitude], [0]);
        
        step = time / steps;
        int_step = step / 10.0;
        
        for t in numpy.arange(step, time + step, step):
            self.__amplitude = self.__calculate(t, step, int_step);
            
            if collect_dynamic is True:
                dynamic_amplitude.append([ numpy.real(amplitude)[0] for amplitude in self.__amplitude ]);
                dynamic_time.append(t);
        
        if collect_dynamic is False:
            dynamic_amplitude.append([ numpy.real(amplitude)[0] for amplitude in self.__amplitude ]);
            dynamic_time.append(time);

        output_sync_dynamic = fsync_dynamic(dynamic_amplitude, dynamic_time);
        return output_sync_dynamic;

def instantaneous_frequency(data, fs, fk):
    """
    Instantaneous frequency of a signal.

    Computes the instantaneous frequency of the given data which can be
    windowed or not. The instantaneous frequency is determined by the time
    derivative of the analytic signal of the input data.

    :type data: :class:`~numpy.ndarray`
    :param data: Data to determine instantaneous frequency of.
    :param fs: Sampling frequency.
    :param fk: Coefficients for calculating time derivatives
        (calculated via central difference).
    :return: **omega[, domega]** - Instantaneous frequency of input data, Time
        derivative of instantaneous frequency (windowed only).
    """
    x = envelope(data)
    if len(x[0].shape) > 1:
        omega = np.zeros(x[0].shape[0], dtype=np.float64)
        i = 0
        for row in x[0]:
            f = np.real(row)
            h = np.imag(row)
            # faster alternative to calculate f_add
            f_add = np.hstack(([f[0]] * (np.size(fk) // 2), f, [f[np.size(f) - 1]] * (np.size(fk) // 2)))
            fd = signal.lfilter(fk, 1, f_add)
            # correct start and end values of time derivative
            fd = fd[np.size(fk) - 1 : np.size(fd)]
            # faster alternative to calculate h_add
            h_add = np.hstack(([h[0]] * (np.size(fk) // 2), h, [h[np.size(h) - 1]] * (np.size(fk) // 2)))
            hd = signal.lfilter(fk, 1, h_add)
            # correct start and end values of time derivative
            hd = hd[np.size(fk) - 1 : np.size(hd)]
            omega_win = abs(((f * hd - fd * h) / (f * f + h * h)) * fs / 2 / np.pi)
            omega[i] = np.median(omega_win)
            i = i + 1
        # faster alternative to calculate omega_add
        omega_add = np.hstack(
            ([omega[0]] * (np.size(fk) // 2), omega, [omega[np.size(omega) - 1]] * (np.size(fk) // 2))
        )
        domega = signal.lfilter(fk, 1, omega_add)
        # correct start and end values of time derivative
        domega = domega[np.size(fk) - 1 : np.size(domega)]
        return omega, domega
    else:
        omega = np.zeros(np.size(x[0]), dtype=np.float64)
        f = np.real(x[0])
        h = np.imag(x[0])
        # faster alternative to calculate f_add
        f_add = np.hstack(([f[0]] * (np.size(fk) // 2), f, [f[np.size(f) - 1]] * (np.size(fk) // 2)))
        fd = signal.lfilter(fk, 1, f_add)
        # correct start and end values of time derivative
        fd = fd[np.size(fk) - 1 : np.size(fd)]
        # faster alternative to calculate h_add
        h_add = np.hstack(([h[0]] * (np.size(fk) // 2), h, [h[np.size(h) - 1]] * (np.size(fk) // 2)))
        hd = signal.lfilter(fk, 1, h_add)
        # correct start and end values of time derivative
        hd = hd[np.size(fk) - 1 : np.size(hd)]
        omega = abs(((f * hd - fd * h) / (f * f + h * h)) * fs / 2 / np.pi)
        return omega

def getArgumentVariable(_ComplexVariable):

	#Debug
	'''
	print('l 31 Numscipier')
	print('_ComplexVariable is ')
	print(_ComplexVariable)
	print('')
	'''

	#import
	import numpy as np

	#return
	return 2.*np.arctan(
	np.imag(_ComplexVariable)/(
	        np.sqrt(
	            np.imag(
	                _ComplexVariable
	            )**2+np.real(
	                _ComplexVariable
	            )**2)+np.real(
	                _ComplexVariable
	            )
	    )
	);

def test_morlet():
    """Test morlet with and without zero mean"""
    Wz = morlet(1000, [10], 2., zero_mean=True)
    W = morlet(1000, [10], 2., zero_mean=False)

    assert_true(np.abs(np.mean(np.real(Wz[0]))) < 1e-5)
    assert_true(np.abs(np.mean(np.real(W[0]))) > 1e-3)

def _select_function(sort, typ):
    if typ in ['F','D']:
        if callable(sort):
            # assume the user knows what they're doing
            sfunction = sort
        elif sort == 'lhp':
            sfunction = lambda x,y: (np.real(x/y) < 0.0)
        elif sort == 'rhp':
            sfunction = lambda x,y: (np.real(x/y) >= 0.0)
        elif sort == 'iuc':
            sfunction = lambda x,y: (abs(x/y) <= 1.0)
        elif sort == 'ouc':
            sfunction = lambda x,y: (abs(x/y) > 1.0)
        else:
            raise ValueError("sort parameter must be None, a callable, or "
                "one of ('lhp','rhp','iuc','ouc')")
    elif typ in ['f','d']:
        if callable(sort):
            # assume the user knows what they're doing
            sfunction = sort
        elif sort == 'lhp':
            sfunction = lambda x,y,z: (np.real((x+y*1j)/z) < 0.0)
        elif sort == 'rhp':
            sfunction = lambda x,y,z: (np.real((x+y*1j)/z) >= 0.0)
        elif sort == 'iuc':
            sfunction = lambda x,y,z: (abs((x+y*1j)/z) <= 1.0)
        elif sort == 'ouc':
            sfunction = lambda x,y,z: (abs((x+y*1j)/z) > 1.0)
        else:
            raise ValueError("sort parameter must be None, a callable, or "
                "one of ('lhp','rhp','iuc','ouc')")
    else:  # to avoid an error later
        raise ValueError("dtype %s not understood" % typ)
    return sfunction