def test_plan_call():
    for shape in tested_shapes:
        plan = Plan(
            input_array=ranf_unit_complex(shape),
            output_array=numpy.empty(shape, dtype=numpy.complex),
            direction=Direction.forward,
        )
        testing.assert_allclose(
            plan(),
            fft.fftn(plan.input_array)
        )
        testing.assert_allclose(
            plan(normalize=True),
            fft.fftn(plan.input_array) / plan.input_array.size
        )
        plan = Plan(
            input_array=ranf_unit_complex(shape),
            output_array=numpy.empty(shape, dtype=numpy.complex),
            direction=Direction.backward
        )
        testing.assert_allclose(
            plan(),
            fft.ifftn(plan.input_array) * plan.input_array.size
        )
        testing.assert_allclose(
            plan(normalize=True),
            fft.ifftn(plan.input_array)
        )

	def compute_init_displacement(self, Dens):
		"""compute Zeldovich displacement from initial density"""
		fDens = fft.fftn(Dens)
		fPot  = fDens / self.km2 * self.mass_res**2
		vx = fft.ifftn(fPot * -1j * np.sin(self.Km[0])).real / self.mass_res
		vy = fft.ifftn(fPot * -1j * np.sin(self.Km[1])).real / self.mass_res
		return np.array([vx,vy])

    def fft_g2r(self, fg, fg_ishifted=False):
        """
        FFT of array ``fg`` given in G-space.
        """
        ndim, shape = fg.ndim, fg.shape

        if ndim == 1:
            fg = np.reshape(fg, self.shape)
            return self.fft_g2r(fg, fg_ishifted=fg_ishifted).flatten()

        if ndim == 3:
            assert self.size == np.prod(shape[-3:])
            if fg_ishifted: fg = ifftshift(fg)
            fr = ifftn(fg)

        elif ndim > 3:
            assert self.size == np.prod(shape[-3:])
            axes = np.arange(ndim)[-3:]
            if fg_ishifted: fg = ifftshift(fg, axes=axes)
            fr = ifftn(fg, axes=axes)

        else:
            raise NotImplementedError("ndim < 3 are not supported")

        return fr * self.size

def garfield(B, P, T = Identity(), seed = None):
	if seed != None:
		random.seed(seed)
	wn = random.normal(0, 1, B.shape)
	f  = fft.ifftn(fft.fftn(wn) * np.sqrt(P(B.K))).real
	#f /= f.std()
	return fft.ifftn(fft.fftn(f) * T(B.K)).real

def Au(U,GF,EpsArr,NX,NY,NZ):
    """Returns the result of matrix-vector multiplication
       by the system matrix A=I-GX
    """
    # reshaping input vector into 4-D array
    Uarr=sci.reshape(U,(NX,NY,NZ,3))
    # extended zero-padded arrays
    Uext=sci.zeros((2*NX,2*NY,2*NZ,3),complex)
    Vext=sci.zeros((2*NX,2*NY,2*NZ,3),complex)
    Jext=sci.zeros((2*NX,2*NY,2*NZ,3),complex)
    JFext=sci.zeros((2*NX,2*NY,2*NZ,3),complex)
    Uext[0:NX,0:NY,0:NZ,:]=Uarr
    # contrast current array
    s=0
    while s<=2:
        Jext[0:NX,0:NY,0:NZ,s]=Uext[0:NX,0:NY,0:NZ,s]*(EpsArr[0:NX,0:NY,0:NZ]-1.0)
        JFext[:,:,:,s]=fft.fftn(sci.squeeze(Jext[:,:,:,s]))
        s=s+1
    Vext[:,:,:,0]=Uext[:,:,:,0]-\
    fft.ifftn(sci.squeeze(sci.multiply(GF[:,:,:,0,0],JFext[:,:,:,0])+\
                          sci.multiply(GF[:,:,:,0,1],JFext[:,:,:,1])+\
                          sci.multiply(GF[:,:,:,0,2],JFext[:,:,:,2])))
    Vext[:,:,:,1]=Uext[:,:,:,1]-\
    fft.ifftn(sci.squeeze(sci.multiply(GF[:,:,:,1,0],JFext[:,:,:,0])+\
                          sci.multiply(GF[:,:,:,1,1],JFext[:,:,:,1])+\
                          sci.multiply(GF[:,:,:,1,2],JFext[:,:,:,2])))
    Vext[:,:,:,2]=Uext[:,:,:,2]-\
    fft.ifftn(sci.squeeze(sci.multiply(GF[:,:,:,2,0],JFext[:,:,:,0])+\
                          sci.multiply(GF[:,:,:,2,1],JFext[:,:,:,1])+\
                          sci.multiply(GF[:,:,:,2,2],JFext[:,:,:,2])))
    # reshaping output into column vector
    V=sci.reshape(Vext[0:NX,0:NY,0:NZ,:],(NX*NY*NZ*3,1))

    return V

    def solve_neutral(self, phi_g, rho_g, eps=None):
        """Solve Poissons equation for a neutral and periodic charge density.

        Parameters
        ----------
        phi_g: ndarray
            Potential (output array).
        rho_g: ndarray
            Charge distribution (in units of -e).

        """

        assert phi_g.dtype == self.dtype
        assert rho_g.dtype == self.dtype
        
        if self.gd.comm.size == 1:
            # Note, implicit downcast from complex to float when the dtype of
            # phi_g is float
            phi_g[:] = ifftn(fftn(rho_g) * 4.0 * pi / self.k2_Q)
        else:
            rho_g = self.gd.collect(rho_g)
            if self.gd.comm.rank == 0:
                globalphi_g = ifftn(fftn(rho_g) * 4.0 * pi / self.k2_Q)
            else:
                globalphi_g = None
            # What happens here if globalphi is complex and phi is real ??????
            self.gd.distribute(globalphi_g, phi_g)
            
        return 1

    def source_terms(self, mara, retphi=False):
        from numpy.fft import fftfreq, fftn, ifftn
        ng = mara.number_guard_zones()
        G = self.G
        L = 1.0
        Nx, Ny, Nz = mara.fluid.shape
        Nx -= 2*ng
        Ny -= 2*ng
        Nz -= 2*ng
        P = mara.fluid.primitive[ng:-ng,ng:-ng,ng:-ng]
        rho = P[...,0]
        vx = P[...,2]
        vy = P[...,3]
        vz = P[...,4]

        K = [fftfreq(Nx)[:,np.newaxis,np.newaxis] * (2*np.pi*Nx/L),
             fftfreq(Ny)[np.newaxis,:,np.newaxis] * (2*np.pi*Ny/L),
             fftfreq(Nz)[np.newaxis,np.newaxis,:] * (2*np.pi*Nz/L)]
        delsq = -(K[0]**2 + K[1]**2 + K[2]**2)
        delsq[0,0,0] = 1.0 # prevent division by 0

        rhohat = fftn(rho)
        phihat = (4*np.pi*G) * rhohat / delsq
        fx = -ifftn(1.j * K[0] * phihat).real
        fy = -ifftn(1.j * K[1] * phihat).real
        fz = -ifftn(1.j * K[2] * phihat).real

        S = np.zeros(mara.fluid.shape + (5,))
        S[ng:-ng,ng:-ng,ng:-ng,0] = 0.0
        S[ng:-ng,ng:-ng,ng:-ng,1] = rho * (fx*vx + fy*vy + fz*vz)
        S[ng:-ng,ng:-ng,ng:-ng,2] = rho * fx
        S[ng:-ng,ng:-ng,ng:-ng,3] = rho * fy
        S[ng:-ng,ng:-ng,ng:-ng,4] = rho * fz
        return (S, ifftn(phihat).real) if retphi else S

def nccfft(s, p, room, fact=1):
    """Used for all Patterns that do not fall other categories.

    Cross correlates normalized Source and Pattern images while
    taking advantage of FFTs for the convolution step. 
    
    ----------
    s, p : Image
       Pattern and Source images for comparison.
       
    fact : int, optional
       Factor by which both Source and Pattern are
       scaled down.

    ----------
    out1 : ndarray[float]
       Confidence matrix for matches.

    out2 : float
       Threshold for deciding if a match has been found.

    out3 : float
       Mean of the confidence matrix.
    """
    
    # subtract mean from Pattern
    pmm = p - p.mean()
    pstd = p.std()
    n = p.size
    
    # make matrix of ones the same size as pattern
    u = np.ones(p.shape)

    # pad matrices (necessary for convolution)
    s = pad_by(s, room)

    upad = pad_to_size_of(u, s)
    pmmpad = pad_to_size_of(pmm, s)
    
    # compute neccessary ffts
    fftppad = fftn(pmmpad)
    ffts = fftn(s)
    fftss = fftn(s**2)
    fftu = fftn(upad)

    # compute conjugates
    cfppad = np.conj(fftppad)
    cfu = np.conj(fftu)

    # do multiplications and ifft's
    top = ifftn(cfppad * ffts)
    bot1 = n * ifftn(cfu * fftss)
    bot2 = ifftn(cfu * ffts) ** 2

    # finish it off!
    bottom = pstd * np.sqrt(bot1 - bot2)
    full = top / bottom

    return np.where(full.real.max() == full.real)

def FieldIFFT(a_hat, a):
	"""Calculate the component-wise 2D or 3D inverse FFT of a vector field a_hat, and stores it in a."""

	if DimOfVectorFieldDomain(a) == 2:
		a[:,:,0], a[:,:,1] = real(fft.ifftn(a_hat[:,:,0])), real(fft.ifftn(a_hat[:,:,1]))
	else:
		a[...,0], a[...,1], a[...,2] = real(fft.ifftn(a_hat[...,0])), real(fft.ifftn(a_hat[...,1])), real(fft.ifftn(a_hat[...,2]))

	return a

 def run(self, image):
     nslice, n_pe, n_fe = image.shape[-3:]
     mask = checkercube(nslice, n_pe, n_fe)
     from recon.tools import Recon
     if Recon._FAST_ARRAY:
         image[:] = mask*ifftn(mask*image[:], axes=[-3,-2,-1])
     else:
         for vol in image:
             vol[:] = mask*ifftn(mask*vol[:])

 def solve_neutral(self, phi_g, rho_g, eps=None):
     if self.gd.comm.size == 1:
         phi_g[:] = ifftn(fftn(rho_g) * 4.0 * pi / self.k2_Q).real
     else:
         rho_g = self.gd.collect(rho_g)
         if self.gd.comm.rank == 0:
             globalphi_g = ifftn(fftn(rho_g) * 4.0 * pi / self.k2_Q).real
         else:
             globalphi_g = None
         self.gd.distribute(globalphi_g, phi_g)
     return 1

def ifftd(I, dims=None):

    # Compute fft
    if dims is None:
        X = ifftn(I)
    elif dims == 2:
        X = ifft2(I, axes=(0, 1))
    else:
        X = ifftn(I, axes=tuple(range(dims)))

    return X

    def propagate(self):
        r"""Given the wavefunction values :math:`\Psi(\Gamma)` at time :math:`t`, calculate
        new values :math:`\Psi^\prime(\Gamma)` at time :math:`t + \tau`. We perform exactly
        one single timestep of size :math:`\tau` within this function.
        """
        # How many components does Psi have
        N = self._psi.get_number_components()

        # Unpack the values from the current WaveFunction
        values = self._psi.get_values()

        # First step with the potential
        tmp = [zeros(value.shape, dtype=complexfloating) for value in values]
        for row in range(0, N):
            for col in range(0, N):
                tmp[row] = tmp[row] + self._VE[row * N + col] * values[col]

        # Go to Fourier space
        tmp = [fftn(component) for component in tmp]

        # First step with the kinetic operator
        tmp = [self._TE * component for component in tmp]

        # Go back to real space
        tmp = [ifftn(component) for component in tmp]

        # Central step with V-tilde
        tmp2 = [zeros(value.shape, dtype=complexfloating) for value in values]
        for row in range(0, N):
            for col in range(0, N):
                tmp2[row] = tmp2[row] + self._VEtilde[row * N + col] * tmp[col]

        # Go to Fourier space
        tmp = [fftn(component) for component in tmp2]

        # Second step with the kinetic operator
        tmp = [self._TE * component for component in tmp]

        # Go back to real space
        tmp = [ifftn(component) for component in tmp]

        # Second step with the potential
        values = [zeros(component.shape, dtype=complexfloating) for component in tmp]
        for row in range(0, N):
            for col in range(0, N):
                values[row] = values[row] + self._VE[row * N + col] * tmp[col]

        # Pack values back to WaveFunction object
        # TODO: Consider squeeze(.) of data before repacking
        self._psi.set_values(values)

def spec2d(u,v,w,dims=(1.0,1.0)):
	ny = u.shape[0]
	nx = u.shape[1]

	lens = np.array(dims)
	res = np.array(u.shape)
	minres = res.min()
	limiter = (res/lens).argmin()
	factors = lens / lens[limiter]

	fu = np.abs(fft.ifftn(u))**2 + np.abs(fft.ifftn(v))**2 + np.abs(fft.ifftn(w))**2
	fu = circleAvg(fu[0:ny/2,0:nx/2],limiter,factors)

	return fu

  def f2_solve(self, rhs, y, t, dt, f2, **kwargs):
    """Solve and evaluate implicit piece."""

    # solve (rebuild operator every time, as dt may change)
    invop = 1.0 / (1.0 - self.nu*dt*self.laplacian)

    z = fft.fftn(rhs)
    z = invop * z

    y[...] = np.real(fft.ifftn(z))

    # evaluate
    z = self.nu * self.laplacian * z

    f2[...] = np.real(fft.ifftn(z))

    def icwt2d(self, da=0.25):
        '''
        Inverse bi-dimensional continuous wavelet transform as in Wang and
        Lu (2010), equation [5].

        Parameters
        ----------
        da : float, optional
            Spacing in the frequency axis.
        '''
        if self.Wf is None:
            raise TypeError("Run cwt2D before icwt2D")
        m0, l0, k0 = self.Wf.shape

        if m0 != self.scales.size:
            raise Warning('Scale parameter array shape does not match\
                           wavelet transform array shape.')
        # Calculates the zonal and meridional wave numters.
        L, K = 2 ** int(np.ceil(np.log2(l0))), 2 ** int(np.ceil(np.log2(k0)))
        # Calculates the zonal and meridional wave numbers.
        l, k = fftfreq(L, self.dy), fftfreq(K, self.dx)
        # Creates empty inverse wavelet transform array and fills it for every
        # discrete scale using the convolution theorem.
        self.iWf = np.zeros((m0, L, K), 'complex')
        for i, an in enumerate(self.scales):
            psi_ft_bar = an * self.wavelet.psi_ft(an * k, an * l)
            W_ft = fftn(self.Wf[i, :, :], s=(L, K))
            self.iWf[i, :, :] = ifftn(W_ft * psi_ft_bar, s=(L, K)) *\
                da / an ** 2.

        self.iWf = self.iWf[:, :l0, :k0].real.sum(axis=0) / self.wavelet.cpsi

        return self

def interpolate(yF, yG, fevalF=None, fevalG=None,
                dim=1, xrat=2, interpolation_order=-1, **kwargs):
  """Interpolate yG to yF."""

  if interpolation_order == -1:

    zG = fft.fftn(yG)
    zF = np.zeros(fevalF.shape, zG.dtype)

    zF[fevalF.half] = zG[fevalG.full]

    yF[...] = np.real(2**dim*fft.ifftn(zF))

  elif interpolation_order == 2:

    if dim != 1:
      raise NotImplementedError

    yF[0::xrat] = yG
    yF[1::xrat] = (yG + np.roll(yG, -1)) / 2.0

  elif interpolation_order == 4:

    if dim != 1:
      raise NotImplementedError

    yF[0::xrat] = yG
    yF[1::xrat] = ( - np.roll(yG,1)
                    + 9.0*yG
                    + 9.0*np.roll(yG,-1)
                    - np.roll(yG,-2) ) / 16.0

  else:
    raise ValueError, 'interpolation order must be -1, 2 or 4'

def fftconvolve(in1, in2, mode='same'):
    """Convolve two N-dimensional arrays using FFT. See convolve.

    """
    s1 = array(in1.shape)
    s2 = array(in2.shape)
    complex_result = (np.issubdtype(in1.dtype, np.complex) or
                      np.issubdtype(in2.dtype, np.complex))
    size = s1 + s2 - 1

    # Always use 2**n-sized FFT
    fsize = (2 ** np.ceil(np.log2(size))).astype('int')
    IN1 = fftn(in1, fsize)
    IN1 *= fftn(in2, fsize)
    fslice = tuple([slice(0, int(sz)) for sz in size])
    ret = ifftn(IN1)[fslice].copy()
    del IN1 
    if not complex_result:
        ret = ret.real
    if mode == "full":
        return ret
    elif mode == "same":
        if np.product(s1, axis=0) > np.product(s2, axis=0):
            osize = s1
        else:
            osize = s2
        return _centered(ret, osize)
    elif mode == "valid":
        return _centered(ret, abs(s2 - s1) + 1)

    return conv[:s[0], :s[1], :s[2]]

 def cost_closure(x, k):
     if k is None:
         return lambda: x.ravel().T.dot(x.ravel())
     else:
         kx = ifftn(k[..., None] * fftn(x, axes=(-2, -1)),
                    axes=(-2, -1))
         return lambda: x.ravel().T.dot(kx.ravel())

def ktoi(data,axis=-1):
    if (axis == -1):
        ax = fth.arange(0,data.ndim)
    else:
        ax = axis

    return fth.fftshift(ft.ifftn(fth.ifftshift(data,axes=ax),axes=ax),axes=ax)