Python pyfftw.n_byte_align_empty函数代码示例

def dfm3(input,angles,Npad):
    # input: aligned data
    # angles: projection angles
    # N_pad: size of padded projection.

    input = np.double(input)
    (Nx,Ny,Nproj) = input.shape
    angles = np.double(angles)
    cen = np.floor(Ny/2.0)
    cen_pad = np.floor(Npad/2.0)
    pad_pre = np.ceil((Npad-Ny)/2.0); pad_post = np.floor((Npad-Ny)/2.0)
    
    # Initialization
    Nz = Ny
    w = np.zeros((Nx,Ny,np.int(Nz/2+1))) #store weighting factors
    v = pyfftw.n_byte_align_empty((Nx,Ny,Nz/2+1),16,dtype='complex128')
    v = np.zeros(v.shape) + 1j*np.zeros(v.shape)
    recon = pyfftw.n_byte_align_empty((Nx,Ny,Nz),16,dtype='float64')
    recon_fftw_object = pyfftw.FFTW(v,recon,direction='FFTW_BACKWARD',axes=(0,1,2))

    p = pyfftw.n_byte_align_empty((Nx,Npad),16,dtype='float64')
    pF = pyfftw.n_byte_align_empty((Nx,Npad/2+1),16,dtype='complex128')
    p_fftw_object = pyfftw.FFTW(p,pF,axes=(0,1))

    dk = np.double(Ny)/np.double(Npad)

    for a in range(0, Nproj):
        #print angles[a]
        ang = angles[a]*np.pi/180
        projection = input[:,:,a] #2D projection image
        p = np.lib.pad(projection,((0,0),(pad_pre,pad_post)),'constant',constant_values=(0,0)) #pad zeros
        p = np.fft.ifftshift(p) 
        p_fftw_object.update_arrays(p,pF)
        p_fftw_object()

        probjection_f = pF.copy()
        if ang<0:
            probjection_f = np.conj(pF.copy())
            probjection_f[1:,:] = np.flipud(probjection_f[1:,:])
            ang = np.pi+ang

        # Bilinear extrapolation
        for i in range(0, np.int(np.ceil(Npad/2))+1):
            ky = i*dk;  #kz = 0;
            ky_new = np.cos(ang)*ky #new coord. after rotation
            kz_new = np.sin(ang)*ky 
            sy = abs(np.floor(ky_new) - ky_new) #calculate weights
            sz = abs(np.floor(kz_new) - kz_new)
            for b in range(1,5): #bilinear extrapolation
                pz,py,weight = bilinear(kz_new,ky_new,sz,sy,Ny,b)
                if (py>=0 and py<Ny and pz>=0 and pz<Nz/2+1):
                    w[:,py,pz] = w[:,py,pz] + weight
                    v[:,py,pz] = v[:,py,pz] + weight * probjection_f[:,i]

    v[w!=0] = v[w!=0]/w[w!=0]
    recon_F = v.copy()
    recon_fftw_object.update_arrays(v,recon)
    recon_fftw_object()
    recon = np.fft.fftshift(recon)
    return recon.astype(np.float32)

    def test_flags(self):
        '''Test to see if the flags are correct
        '''
        fft = FFTW(self.input_array, self.output_array)
        self.assertEqual(fft.flags, ('FFTW_MEASURE',))

        fft = FFTW(self.input_array, self.output_array, 
                flags=('FFTW_DESTROY_INPUT', 'FFTW_UNALIGNED'))
        self.assertEqual(fft.flags, ('FFTW_DESTROY_INPUT', 'FFTW_UNALIGNED'))

        # Test an implicit flag
        _input_array = n_byte_align_empty(256, 16, dtype='complex64')
        _output_array = n_byte_align_empty(256, 16, dtype='complex64')

        # These are guaranteed to be misaligned (due to dtype size == 8)
        input_array = _input_array[:-1]
        output_array = _output_array[:-1]
        u_input_array = _input_array[1:]
        u_output_array = _output_array[1:]

        fft = FFTW(input_array, u_output_array)
        self.assertEqual(fft.flags, ('FFTW_MEASURE', 'FFTW_UNALIGNED'))

        fft = FFTW(u_input_array, output_array)
        self.assertEqual(fft.flags, ('FFTW_MEASURE', 'FFTW_UNALIGNED'))

        fft = FFTW(u_input_array, u_output_array)
        self.assertEqual(fft.flags, ('FFTW_MEASURE', 'FFTW_UNALIGNED'))

    def pre_process(self):
        in_pData = self.get_plugin_in_datasets()[0]
        sino_shape = list(in_pData.get_shape())

        width1 = sino_shape[1] + 2*self.pad
        height1 = sino_shape[0] + 2*self.pad

        v0 = np.abs(self.parameters['vvalue'])
        u0 = np.abs(self.parameters['uvalue'])
        n = np.abs(self.parameters['nvalue'])
        # Create filter
        centerx = np.ceil(width1/2.0)-1.0
        centery = np.int16(np.ceil(height1/2.0)-1)
        self.row1 = centery - v0
        self.row2 = centery + v0+1
        listx = np.arange(width1)-centerx
        filtershape = 1.0/(1.0 + np.power(listx/u0, 2*n))
        filtershapepad2d = np.zeros((self.row2 - self.row1, filtershape.size))
        filtershapepad2d[:] = np.float64(filtershape)
        self.filtercomplex = filtershapepad2d + filtershapepad2d*1j

        a = pyfftw.n_byte_align_empty((height1, width1), 16, 'complex128')
        b = pyfftw.n_byte_align_empty((height1, width1), 16, 'complex128')
        c = pyfftw.n_byte_align_empty((height1, width1), 16, 'complex128')
        d = pyfftw.n_byte_align_empty((height1, width1), 16, 'complex128')
        self.fft_object = pyfftw.FFTW(a, b, axes=(0, 1))
        self.ifft_object = pyfftw.FFTW(c, d, axes=(0, 1),
                                       direction='FFTW_BACKWARD')

   def __init__(self, nbeads, natoms):
      """Initializes nm_trans.

      Args:
         nbeads: The number of beads.
         natoms: The number of atoms.
      """

      self.nbeads = nbeads
      self.natoms = natoms
      try:
         import pyfftw
         info("Import of PyFFTW successful", verbosity.medium)
         self.qdummy = pyfftw.n_byte_align_empty((nbeads, 3*natoms), 16, 'float32')
         self.qnmdummy = pyfftw.n_byte_align_empty((nbeads//2+1, 3*natoms), 16, 'complex64')
         self.fft = pyfftw.FFTW(self.qdummy, self.qnmdummy, axes=(0,), direction='FFTW_FORWARD')
         self.ifft = pyfftw.FFTW(self.qnmdummy, self.qdummy, axes=(0,), direction='FFTW_BACKWARD')
      except ImportError: #Uses standard numpy fft library if nothing better
                          #is available
         info("Import of PyFFTW unsuccessful, using NumPy library instead", verbosity.medium)
         self.qdummy = np.zeros((nbeads,3*natoms), dtype='float32')
         self.qnmdummy = np.zeros((nbeads//2+1,3*natoms), dtype='complex64')
         def dummy_fft(self):
            self.qnmdummy = np.fft.rfft(self.qdummy, axis=0)
         def dummy_ifft(self):
            self.qdummy = np.fft.irfft(self.qnmdummy, n=self.nbeads, axis=0)
         self.fft = lambda: dummy_fft(self)
         self.ifft = lambda: dummy_ifft(self)

    def ShiftConvNumbaFFT(h, N, M, xdtype=np.complex_, powerof2=True):
        # implements Doppler filter:
        # y[n, p] = SUM_k (exp(2*pi*j*n*(k - (L-1))/N) * h[k]) * x[p - k]
        #         = SUM_k (exp(-2*pi*j*n*k/N) * s*[k]) * x[p - (L-1) + k]
        L = len(h)
        outlen = M + L - 1
        nfft = outlen
        if powerof2:
            nfft = pow2(nfft)

        dopplermat = np.exp(2*np.pi*1j*np.arange(N)[:, np.newaxis]*(np.arange(L) - (L - 1))/N)
        dopplermat.astype(np.result_type(h.dtype, np.complex64)) # cast to complex type with precision of h
        hbank = h*dopplermat
        # speed not critical here, just use numpy fft
        hbankpad = zero_pad(hbank, nfft)
        H = np.fft.fft(hbankpad) / nfft # divide by nfft b/c FFTW's ifft does not do this

        xcdtype = np.result_type(xdtype, np.complex64) # cast to complex type with precision of x
        xpad = pyfftw.n_byte_align(np.zeros(nfft, xcdtype), 16)
        X = pyfftw.n_byte_align(np.zeros(nfft, xcdtype), 16)
        xfft = pyfftw.FFTW(xpad, X, threads=_THREADS)

        ydtype = np.result_type(H.dtype, xcdtype)
        Y = pyfftw.n_byte_align_empty(H.shape, 16, ydtype)
        y = pyfftw.n_byte_align_empty(H.shape, 16, ydtype)
        ifft = pyfftw.FFTW(Y, y, direction='FFTW_BACKWARD', threads=_THREADS)

        xtype = numba.__getattribute__(str(np.dtype(xdtype)))

        #htype = numba.__getattribute__(str(H.dtype))
        #xctype = numba.__getattribute__(str(X.dtype))
        #ytype = numba.__getattribute__(str(Y.dtype))
        #@jit(argtypes=[htype[:, ::1], xctype[::1], ytype[:, ::1], xtype[::1]])
        #def fun(H, X, Y, x):
            #xpad[:M] = x
            #xfft.execute() # input is xpad, output is X
            #Y[:, :] = H*X # need expression optimized by numba but that writes into Y
            #ifft.execute() # input is Y, output is y

            #yc = np.array(y)[:, :outlen] # need a copy, which np.array provides
            #return yc

        #@dopplerbank_dec(h, N, M, nfft=nfft, H=H)
        #def shiftconv_numba_fft(x):
            #return fun(H, X, Y, x)

        #@jit(argtypes=[xtype[::1]])
        @jit
        def shiftconv_numba_fft(x):
            xpad[:M] = x
            xfft.execute() # input is xpad, output is X
            Y[:, :] = X*H # need expression optimized by numba but that writes into Y
            ifft.execute() # input is Y, output is y

            yc = np.array(y[:, :outlen]) # need a copy, which np.array provides
            return yc

        shiftconv_numba_fft = dopplerbank_dec(h, N, M, nfft=nfft, H=H)(shiftconv_numba_fft)

        return shiftconv_numba_fft

    def generate_wisdom(self):
        for each_dtype in (numpy.complex128, numpy.complex64, 
                numpy.clongdouble):

            a = n_byte_align_empty((1,1024), 16, each_dtype)
            b = n_byte_align_empty(a.shape, 16, dtype=a.dtype)
            fft = FFTW(a,b)

	def __init__(self, ctrlNs, printQueue, rawDataRingBuf, fftDataRingBuf):

		self.log = logging.getLogger("Main.FFTWorker")
		logSetup.initLogging(printQ = printQueue)

		self.log.info("FFT Worker Starting up")

		self.ctrlNs         = ctrlNs
		self.printQueue     = printQueue
		self.rawDataRingBuf = rawDataRingBuf
		self.fftDataRingBuf = fftDataRingBuf

		self.fftChunkSize = s.FFT_CHUNK_SIZE
		self.outputSize = self.fftChunkSize//2 + 1
		self.chunksPerAcq = int(SignalHound.rawSweepArrSize/self.fftChunkSize)
		self.overlap = s.FFT_OVERLAP
		self.window = np.hamming(self.fftChunkSize)
		inArr = pyfftw.n_byte_align_empty(self.fftChunkSize, 16, dtype=np.float32)
		outArr = pyfftw.n_byte_align_empty(self.outputSize, 16, dtype=np.complex64)

		self.log.info("Choosing maximally optimized transform")
		self.fftFunc = pyfftw.FFTW(inArr, outArr, flags=('FFTW_PATIENT', "FFTW_DESTROY_INPUT"))
		self.log.info("Optimized transform selected. Run starting")

		self.run()

def radial_average(tiltseries,kr_cutoffs):
    (Nx,Ny,Nproj) = tiltseries.shape

    f = pyfftw.n_byte_align_empty((Nx,Ny/2+1),16,dtype='complex128')
    r = pyfftw.n_byte_align_empty((Nx,Ny),16,dtype='float64')
    p_fftw_object = pyfftw.FFTW(r,f,axes=(0,1))
    Ir = np.zeros(kr_cutoffs.size); I = np.zeros(kr_cutoffs.size)

    kx = np.fft.fftfreq(Nx)
    ky = np.fft.fftfreq(Ny)
    ky = ky[0:int(np.ceil(Ny/2)+1)]

    kX,kY = np.meshgrid(ky,kx)
    kR = np.sqrt(kY**2+kX**2)

    for a in range(0,Nproj):
        r = tiltseries[:,:,a].copy().astype('float64')
        p_fftw_object.update_arrays(r,f); p_fftw_object.execute()
        shell = kR<=kr_cutoffs[0]
        I[0] = np.sum(np.absolute(f[shell]))
        I[0] = I[0]/np.sum(shell)
        for j in range(1,kr_cutoffs.size):
            shell = np.logical_and(kR>kr_cutoffs[j-1],kR<=kr_cutoffs[j])
            I[j] = np.sum(np.absolute(f[shell]))
            I[j] = I[j]/np.sum(shell)
        Ir = Ir + I
    Ir = Ir/Nproj
    return Ir

    def test_update_data_with_unaligned_original(self):
        in_shape = self.input_shapes['2d']
        out_shape = self.output_shapes['2d']

        input_dtype_alignment = self.get_input_dtype_alignment()
        
        axes=(-1,)
        a, b = self.create_test_arrays(in_shape, out_shape)

        # Offset from 16 byte aligned to guarantee it's not
        # 16 byte aligned
        a__ = n_byte_align_empty(
                numpy.prod(in_shape)*a.itemsize + input_dtype_alignment,
                16, dtype='int8')
        
        a_ = a__[input_dtype_alignment:].view(dtype=self.input_dtype).reshape(*in_shape)
        a_[:] = a
        
        b__ = n_byte_align_empty(
                numpy.prod(out_shape)*b.itemsize + input_dtype_alignment, 
                16, dtype='int8')
        
        b_ = b__[input_dtype_alignment:].view(dtype=self.output_dtype).reshape(*out_shape)
        b_[:] = b
        
        fft, ifft = self.run_validate_fft(a_, b_, axes, 
                force_unaligned_data=True)
        
        self.run_validate_fft(a, b_, axes, fft=fft, ifft=ifft)
        self.run_validate_fft(a_, b, axes, fft=fft, ifft=ifft)
        self.run_validate_fft(a_, b_, axes, fft=fft, ifft=ifft)

    def test_call_with_normalisation_precision(self):
        '''The normalisation should use a double precision scaling.
        '''
        # Should be the case for double inputs...
        _input_array = n_byte_align_empty((256, 512), 16,
                dtype='complex128')

        ifft = FFTW(self.output_array, _input_array, 
                direction='FFTW_BACKWARD')

        ref_output = ifft(normalise_idft=False).copy()/numpy.float64(ifft.N)
        test_output = ifft(normalise_idft=True).copy()

        self.assertTrue(numpy.alltrue(ref_output == test_output))

        # ... and single inputs.
        _input_array = n_byte_align_empty((256, 512), 16,
                dtype='complex64')

        ifft = FFTW(numpy.array(self.output_array, _input_array.dtype), 
                    _input_array, 
                    direction='FFTW_BACKWARD')

        ref_output = ifft(normalise_idft=False).copy()/numpy.float64(ifft.N)
        test_output = ifft(normalise_idft=True).copy()

        self.assertTrue(numpy.alltrue(ref_output == test_output))

 def __init__(self,
              n = 1364,
              N = 2048,
              iter0 = 138000,
              nfiles = 64,
              fft_threads = 32,
              out_threads = 4,
              src_dir = './',
              dst_dir = './',
              src_format = 'K{0:0>6}QNP{1:0>3}',
              dst_format = 'RMHD_{0}_t{1:0>4x}_z{2:0>7x}'):
     self.src_format = src_format
     self.dst_format = dst_format
     self.src_dir = src_dir
     self.dst_dir = dst_dir
     self.n = n
     self.N = N
     self.iter0 = iter0
     self.nfiles = nfiles
     self.kdata = pyfftw.n_byte_align_empty(
             (self.N//2+1, 2, self.N, self.N),
             pyfftw.simd_alignment,
             dtype = np.complex64)
     self.rrdata = pyfftw.n_byte_align_empty(
             (self.N, self.N, self.N, 2),
             pyfftw.simd_alignment,
             dtype = np.float32)
     self.rzdata = pyfftw.n_byte_align_empty(
             ((self.N//8) * (self.N//8) * (self.N//8),
              8*8*8*2),
             pyfftw.simd_alignment,
             dtype = np.float32)
     if type(self.dst_dir) == type([]):
         self.zdir = np.array(
             range(0,
                   self.rzdata.shape[0],
                   self.rzdata.shape[0] // len(self.dst_dir)))
     self.cubbies_per_file = self.rzdata.shape[0] // self.nfiles
     if (os.path.isfile('fftw_wisdom.pickle.gz')):
         pyfftw.import_wisdom(
             pickle.load(gzip.open('fftw_wisdom.pickle.gz', 'rb')))
     print('about to initialize the fftw plan, which can take a while')
     self.plan = pyfftw.FFTW(
             self.kdata.transpose(3, 2, 0, 1), self.rrdata,
             axes = (0, 1, 2),
             direction = 'FFTW_BACKWARD',
             flags = ('FFTW_MEASURE',
                      'FFTW_DESTROY_INPUT'),
             threads = fft_threads)
     print('finalized fftw initialization')
     bla = pyfftw.export_wisdom()
     pickle.dump(bla, gzip.open('fftw_wisdom.pickle.gz', 'wb'))
     self.fft_threads = fft_threads
     self.out_threads = out_threads
     self.shuffle_lib = np.ctypeslib.load_library(
         'libzshuffle.so',
         os.path.abspath(os.path.join(
             os.path.expanduser('~'), 'repos/RMHD_converter/C-shuffle')))
     return None

def do_it(params):
  DIROUT = params[0]
  FILEINs = params[1]
  t_now = params[2]
  NW = params[3]
  K = params[4]
  PVAL = params[5]
  FS = params[6]
  NFFT = params[7]
  NWINDOWS_PER_WORKER = params[8]
  tapers = params[9]
  worker = params[10]
  t_offset_tic = params[11]
  FILE_TYPE = params[12]
  FILE_LEN_TIC = params[13]

  NCHANNELS = len(FILEINs)
  sig=stats.f.ppf((1-PVAL/NFFT),2,2*K-2)

  fft_in = pyfftw.n_byte_align_empty((NFFT,K,NCHANNELS), 16, 'float32')
  fft_out = pyfftw.n_byte_align_empty((NFFT//2+1,K,NCHANNELS), 16, 'complex64')
  fft_plan = pyfftw.FFTW(fft_in, fft_out, axes=(0,), flags=('FFTW_PATIENT',))
  fft_in2 = pyfftw.n_byte_align_empty(NFFT, 16, 'float32')
  fft_out2 = pyfftw.n_byte_align_empty(NFFT//2+1, 16, 'complex64')
  fft_plan2 = pyfftw.FFTW(fft_in2, fft_out2, flags=('FFTW_PATIENT',))

  NSAMPLES = NFFT//2*(NWINDOWS_PER_WORKER+1)
  dd = np.empty([NSAMPLES, NCHANNELS], dtype=np.float32)
  for i in range(0,NCHANNELS):
    tmp=np.empty(0)
    tmp2=(t_now+worker*NWINDOWS_PER_WORKER)*(NFFT//2)+t_offset_tic
    filei=os.path.join(DIROUT,FILEINs[i])
    if tmp2<FILE_LEN_TIC:
      if FILE_TYPE==1:
        fid=open(filei,'rb')
        fid.seek(tmp2*4,0);
        tmp = fid.read(NSAMPLES*4)
        fid.close()
        tmp = np.array(struct.unpack(str(int(len(tmp)/4))+'f', tmp), dtype=np.float32)
      if FILE_TYPE==2:
        rate, fid = wavfile.read(filei,mmap=True)
        tmp=fid[(tmp2+1):min(tmp2+NSAMPLES, FILE_LEN_TIC)]
    if len(tmp) < NSAMPLES:
      tmp = np.concatenate((tmp, np.zeros(NSAMPLES-len(tmp), dtype=np.float32)))
    dd[:,i] = tmp

  idx=list()
  for j in range(0,NWINDOWS_PER_WORKER):
    ddd=dd[np.array(range(0,NFFT))+j*NFFT//2,:]
    #F, A, f, sig, sd = ftest(ddd, tapers, FS, PVAL, fft_in, fft_out, fft_plan)
    F = ftest(ddd, tapers, PVAL, fft_in, fft_out, fft_plan)
    for l in range(0,NCHANNELS):
      tmp=[i+1 for (i,v) in enumerate(F[1:-1,l]) if v>sig]
      for m in range(0,len(tmp)):
        freq,amp = brown_puckette(ddd[:,l],tmp[m],FS, fft_in2, fft_out2, fft_plan2)
        idx.append((j+worker*NWINDOWS_PER_WORKER, freq, amp, l))
  return idx

 def __init__(self,gridsize):
     self.gridsize = gridsize
     self.a = pyfftw.n_byte_align_empty(gridsize,16,'complex128')
     self.b = pyfftw.n_byte_align_empty(gridsize,16,'complex128')
     self.fftw = pyfftw.FFTW(self.a,self.b,direction='FFTW_FORWARD')
     
     self.c = pyfftw.n_byte_align_empty(gridsize,16,'complex128')
     self.d = pyfftw.n_byte_align_empty(gridsize,16,'complex128')
     self.ifftw = pyfftw.FFTW(self.c,self.d,direction='FFTW_BACKWARD')

def seqDatab_to_four(seqList,fftPlan,inverseVariable,method=1,n=0):		
	four_seq_datab = [];		
	for i in range(len(seqList)):
		if n==0:
			four_seq_datab.append(pyfftw.n_byte_align_empty(len(seqList[i]),1,'complex'));
		else:
			four_seq_datab.append(pyfftw.n_byte_align_empty(n,1,'complex'));					
		four_seq_datab[i][:] = seq_to_four(seqList[i],fftPlan,inverseVariable,method,n)
	return four_seq_datab

    def test_is_n_byte_aligned(self):
        a = n_byte_align_empty(100, 16)
        self.assertTrue(is_n_byte_aligned(a, 16))

        a = n_byte_align_empty(100, 5)
        self.assertTrue(is_n_byte_aligned(a, 5))

        a = n_byte_align_empty(100, 16, dtype="float32")[1:]
        self.assertFalse(is_n_byte_aligned(a, 16))
        self.assertTrue(is_n_byte_aligned(a, 4))