def fromEBmodes(cls,fourier_E,fourier_B,angle=3.14*deg):

		"""
		This class method allows to build a shear map specifying its E and B mode components

		:param fourier_E: E mode of the shear map in fourier space
		:type fourier_E: numpy 2D array, must be of type np.complex128 and must have a shape that is appropriate for a real fourier transform, i.e. (N,N/2 + 1); N should be a power of 2

		:param fourier_B: B mode of the shear map in fourier space
		:type fourier_B: numpy 2D array, must be of type np.complex128 and must have a shape that is appropriate for a real fourier transform, i.e. (N,N/2 + 1); N should be a power of 2

		:param angle: Side angle of the real space map in degrees
		:type angle: float.

		:returns: the corresponding ShearMap instance

		:raises: AssertionErrors for inappropriate inputs

		"""

		assert fourier_E.dtype == np.complex128 and fourier_B.dtype == np.complex128
		assert fourier_E.shape[1] == fourier_E.shape[0]/2 + 1
		assert fourier_B.shape[1] == fourier_B.shape[0]/2 + 1
		assert fourier_E.shape == fourier_B.shape

		#Compute frequencies
		lx = rfftfreq(fourier_E.shape[0])
		ly = fftfreq(fourier_E.shape[0])

		#Safety check
		assert len(lx)==fourier_E.shape[1]
		assert len(ly)==fourier_E.shape[0]

		#Compute sines and cosines of rotation angles
		l_squared = lx[np.newaxis,:]**2 + ly[:,np.newaxis]**2
		l_squared[0,0] = 1.0

		sin_2_phi = 2.0 * lx[np.newaxis,:] * ly[:,np.newaxis] / l_squared
		cos_2_phi = (lx[np.newaxis,:]**2 - ly[:,np.newaxis]**2) / l_squared

		sin_2_phi[0,0] = 0.0
		cos_2_phi[0,0] = 0.0

		#Invert E/B modes and find the components of the shear
		ft_data1 = cos_2_phi * fourier_E - sin_2_phi * fourier_B
		ft_data2 = sin_2_phi * fourier_E + cos_2_phi * fourier_B

		#Invert Fourier transforms
		data1 = irfft2(ft_data1)
		data2 = irfft2(ft_data2)

		#Instantiate new shear map class
		new = cls(np.array([data1,data2]),angle)
		setattr(new,"fourier_E",fourier_E)
		setattr(new,"fourier_B",fourier_B)

		return new

    def output(self):
        """ """

        # One dimension
        if len(self._source.shape) == 1:
            source = self._actual_source
            # Use FFT convolution
            if self._fft:
                if not self._toric:
                    P = rfft(source,self._fft_shape[0])*self._fft_weights
                    R = irfft(P, self._fft_shape[0]).real
                    R = R[self._fft_indices]
                else:
                    P = rfft(source)*self._fft_weights
                    R = irfft(P,source.shape[0]).real

                # if self._toric:
                #     R  = ifft(fft(source)*self._fft_weights).real
                # else:
                #     n = source.shape[0]
                #     self._src_holder[n//2:n//2+n] = source
                #     R = ifft(fft(self._src_holder)*self._fft_weights)
                #     R = R.real[n//2:n//2+n]
            # Use regular convolution
            else:
                R = convolve1d(source, self._weights[::-1], self._toric)
            if self._src_rows is not None:
                R = R[self._src_rows]
            return R.reshape(self._target.shape)
        # Two dimensions
        else:
            source = self._actual_source
            # Use FFT convolution
            if self._fft:
                if not self._toric:
                    P = rfft2(source,self._fft_shape)*self._fft_weights
                    R = irfft2(P, self._fft_shape).real
                    R = R[self._fft_indices]
                else:
                    P = rfft2(source)*self._fft_weights
                    R = irfft2(P,source.shape).real

            # Use SVD convolution
            else:
                R = convolve2d(source, self._weights, self._USV, self._toric)
            if self._src_rows is not None and self._src_cols is not None:
                R = R[self._src_rows, self._src_cols]
        return R.reshape(self._target.shape)

def _correlate2d_fft(data0, kernel0, output=None, mode="nearest", cval=0.0):
    """_correlate2d_fft does 2d correlation of 'data' with 'kernel', storing
    the result in 'output' using the FFT to perform the correlation.

    supported 'mode's include:
        'nearest'   elements beyond boundary come from nearest edge pixel.
        'wrap'      elements beyond boundary come from the opposite array edge.
        'reflect'   elements beyond boundary come from reflection on same array
                    edge.
        'constant'  elements beyond boundary are set to 'cval'
    """

    shape = data0.shape
    kshape = kernel0.shape
    oversized = (np.array(shape) + np.array(kshape))

    dy = kshape[0] // 2
    dx = kshape[1] // 2

    kernel = np.zeros(oversized, dtype=np.float64)
    # convolution <-> correlation
    kernel[:kshape[0], :kshape[1]] = kernel0[::-1,::-1]
    data = iraf_frame.frame(data0, oversized, mode=mode, cval=cval)

    complex_result = (isinstance(data, np.complexfloating) or
                      isinstance(kernel, np.complexfloating))

    Fdata = dft.fft2(data)
    del data

    Fkernel = dft.fft2(kernel)
    del kernel

    np.multiply(Fdata, Fkernel, Fdata)
    del Fkernel

    if complex_result:
        convolved = dft.irfft2(Fdata, s=oversized)
    else:
        convolved = dft.irfft2(Fdata, s=oversized)

    ks0, ks1 = kshape[0] - 1, kshape[1] - 1
    result = convolved[ks0:shape[0] + ks0, ks1:shape[1] + ks1]

    if output is not None:
        output._copyFrom(result)
    else:
        return result

def _filter(filter, array):
    X = fft.rfft2(array)
    X = fft.fftshift(X)
    numpy.multiply(X, filter, out=X)
    X = fft.ifftshift(X)
    x = fft.irfft2(X, s=array.shape)
    return x

def _BandPassFilter(image, len_noise, len_object):
    """
    bandpass filter implementation.
    Source: http://physics-server.uoregon.edu/~raghu/particle_tracking.html
    """
    b = len_noise
    w = round(len_object)
    N = 2 * w + 1

    # Gaussian Convolution Kernel
    sm = numpy.arange(0, N, dtype=numpy.float)
    r = (sm - w) / (2 * b)
    gx = numpy.power(math.e, -r ** 2) / (2 * b * math.sqrt(math.pi))
    gx = numpy.reshape(gx, (gx.shape[0], 1))
    gy = gx.conj().transpose()

    # Boxcar average kernel, background
    bx = numpy.zeros((1, N), numpy.float) + 1 / N
    by = bx.conj().transpose()

    # Convolution with the matrix and kernels
    gxy = gx * gy
    bxy = bx * by
    kernel = fft.rfft2(gxy - bxy, image.shape)

    res = fft.irfft2(fft.rfft2(image) * kernel)
    arr_out = numpy.zeros((image.shape))
    arr_out[w:-w, w:-w] = res[2 * w:, 2 * w:]
    res = numpy.maximum(arr_out, 0)
    return res

def matches_exist(template, image, tolerance=1):
    # just taken from convolution def
    expected = numpy.count_nonzero(template)
    ih, iw = image.shape
    th, tw = template.shape

    # Pad image to even dimensions
    if ih % 2 or iw % 2:
        if ih % 2:
            ih += 1
        if iw % 2:
            iw += 1
        bin_image = pad_bin_image_to_shape(image, (ih, iw))
    if expected == 0:
        return []

    # Calculate the convolution of the FFT's of the image & template
    convolution_freqs = rfft2(image) * rfft2(template[::-1, ::-1],
                                                 image.shape)
    convolution_image = irfft2(convolution_freqs)
    found_bitmap = convolution_image > (expected - tolerance)
    if True in found_bitmap:
        return True
    else:
        return False

def ifftn_mpi(fu, u):
    """ifft in three directions using mpi.
    Need to do ifft in reversed order of fft
    """
    if num_processes == 1:
        #u[:] = irfft(ifft(ifft(fu, axis=0), axis=2), axis=1)
        u[:] = irfftn(fu, axes=(0,2,1))
        return
    
    # Do first owned direction
    Uc_hat[:] = ifft(fu, axis=0)
    
    # Communicate all values
    comm.Alltoall([Uc_hat, MPI.DOUBLE_COMPLEX], [U_mpi, MPI.DOUBLE_COMPLEX])
    for i in range(num_processes): 
        Uc_hatT[:, :, i*Np:(i+1)*Np] = U_mpi[i]

    #for i in range(num_processes):
    #    if not i == rank:
    #        comm.Sendrecv_replace([Uc_send[i], MPI.DOUBLE_COMPLEX], i, 0, i, 0)   
    #    Uc_hatT[:, :, i*Np:(i+1)*Np] = Uc_send[i]
           
    # Do last two directions
    #u[:] = irfft(ifft(Uc_hatT, axis=2), axis=1)
    u[:] = irfft2(Uc_hatT, axes=(2,1))

def ifftn_mpi(fu, u):
    Uc_hat[:] = ifft(fu, axis=0)
    comm.Alltoall([Uc_hat, MPI.DOUBLE_COMPLEX], [U_mpi, MPI.DOUBLE_COMPLEX])
    for i in range(num_processes):
        Uc_hatT[:, i*Np:(i+1)*Np] = U_mpi[i]
    u[:] = irfft2(Uc_hatT, axes=(1,2))
    return u

def convolution(bin_template, bin_image, tollerance=0.5):
    expected = numpy.count_nonzero(bin_template)
    ih, iw = bin_image.shape
    th, tw = bin_template.shape

    # Padd image to even dimensions
    if ih % 2 or iw % 2:
        if ih % 2:
            ih += 1
        if iw % 2:
            iw += 1
        bin_image = pad_bin_image_to_shape(bin_image, (ih, iw))
    if expected == 0:
        return []

    # Calculate the convolution of the FFT's of the image & template
    convolution_freqs = rfft2(bin_image) * rfft2(bin_template[::-1, ::-1],
                                                 bin_image.shape)
    # Reverse the FFT to find the result image
    convolution_image = irfft2(convolution_freqs)
    # At this point, the maximum point in convolution_image should be the
    # bottom right (why?) of the area of greatest match

    # The areas in the result image within expected +- tollerance are where we
    # saw matches
    found_bitmap = ((convolution_image > (expected - tollerance)) &
                    (convolution_image < (expected + tollerance)))

    match_points = numpy.transpose(numpy.nonzero(found_bitmap))  # bottom right

    # Find the top left point from the template (remember match_point is
    # inside the template (hence -1)
    return [((fx - (tw - 1)), (fy - (th - 1))) for (fy, fx) in match_points]

 def getDisplacements2D(self, Z=None, window=False):
     """
     Use phase correlation to find the relative displacement between
     each time step
     """
     if Z is None:
         Z = self.getNbPixelsPerFrame()/self.getNbPixelsPerSlice()/2
     shape = np.asarray(self.get2DShape())
     if window:
         ham = np.hamming(shape[1])*np.atleast_2d(np.hamming(shape[0])).T
     else:
         ham = 1.0
     displs = np.zeros((self.getNbFrames(),2))
     a = rfft2(self.get2DSlice(T=0, Z=Z)*ham)
     for t in range(1,self.getNbFrames()):
         b = rfft2(self.get2DSlice(T=t, Z=Z)*ham)
         #calculate the normalized cross-power spectrum
         #R = numexpr.evaluate(
         #    'a*complex(real(b), -imag(b)/abs(a*complex(real(b), -imag(b))))'
         #    )
         R = a*b.conj()
         Ra = np.abs(a*b.conj())
         R[Ra>0] /= Ra[Ra>0]
         r = irfft2(R)
         #Get the periodic position of the peak
         l = r.argmax()
         displs[t] = np.unravel_index(l, r.shape)
         #prepare next step
         a = b
     return np.where(displs<shape/2, displs, displs-shape)

	def fromConvPower(self,power_func,seed=0,**kwargs):

		"""
		This method uses a supplied power spectrum to generate correlated noise maps in real space via FFTs

		:param power_func: function that given a numpy array of l's returns a numpy array with the according Pl's (this is the input power spectrum); alternatively you can pass an array (l,Pl) and the power spectrum will be calculated with scipy's interpolation routines
		:type power_func: function with the above specifications, or numpy array (l,Pl) of shape (2,n) 

		:param seed: seed of the random generator 
		:type seed: int.

		:param kwargs: keyword arguments to be passed to power_func, or to the interpolate.interp1d routine

		:returns: ConvergenceMap instance of the same exact shape as the one used as blueprint

		"""
		assert self.label == "convergence"

		#Initialize random number generator
		np.random.seed(seed)

		#Generate a random Fourier realization and invert it
		ft_map = self._fourierMap(power_func,**kwargs)
		noise_map = irfft2(ft_map)

		return ConvergenceMap(noise_map,self.side_angle)

def convolve_dcr_image(flux_arr, x_loc, y_loc, bandpass=None, x_size=None, y_size=None, seed=None,
                       psf=None, pad_image=1.5, pixel_scale=None, kernel_radius=None,
                       oversample_image=1, photon_noise=False, sky_noise=0.0, verbose=True, **kwargs):
    """Wrapper to call fast_dft with multiple DCR planes."""
    x_size_use = int(x_size * pad_image)
    y_size_use = int(y_size * pad_image)
    oversample_image = int(oversample_image)
    pixel_scale_use = pixel_scale / oversample_image
    x0 = oversample_image * ((x_size_use - x_size) // 2)
    x1 = x0 + x_size * oversample_image
    y0 = oversample_image * ((y_size_use - y_size) // 2)
    y1 = y0 + y_size * oversample_image
    x_loc_use = x_loc * oversample_image + x0
    y_loc_use = y_loc * oversample_image + y0
    x_size_use *= oversample_image
    y_size_use *= oversample_image
    timing_model = -time.time()
    source_image = fast_dft(flux_arr, x_loc_use, y_loc_use, x_size=x_size_use, y_size=y_size_use,
                            kernel_radius=kernel_radius, **kwargs)
    timing_model += time.time()
    n_star = len(x_loc)
    if oversample_image > 1:
        bright_star = "bright "
    else:
        bright_star = ""
    if verbose:
        if n_star == 1:
            print("Time to model %i %sstar: [%0.3fs]"
                  % (n_star, bright_star, timing_model))
        else:
            print("Time to model %i %sstars: [%0.3fs | %0.5fs per star]"
                  % (n_star, bright_star, timing_model, timing_model / n_star))
    rand_gen = np.random
    if seed is not None:
        rand_gen.seed(seed - 1)
    # The images are purely real, so we can save time by using the real FFT,
    # which uses only half of the complex plane
    convol = np.zeros((y_size_use, x_size_use // 2 + 1), dtype='complex64')
    dcr_gen = dcr_generator(bandpass, pixel_scale=pixel_scale_use, **kwargs)
    timing_fft = -time.time()

    for _i, offset in enumerate(dcr_gen):
        source_image_use = source_image[_i]

        psf_image = psf.drawImage(scale=pixel_scale_use, method='fft', offset=offset,
                                  nx=x_size_use, ny=y_size_use, use_true_center=False)
        if photon_noise:
            base_noise = np.random.normal(scale=1.0, size=(y_size_use, x_size_use))
            base_noise *= np.sqrt(np.abs(source_image_use) / photons_per_adu)
            source_image_use += base_noise
        if sky_noise > 0:
            source_image_use += (rand_gen.normal(scale=sky_noise, size=(y_size_use, x_size_use))
                                 / np.sqrt(bandpass_nstep(bandpass)))
        convol += rfft2(source_image_use) * rfft2(psf_image.array)
    return_image = np.real(fftshift(irfft2(convol)))
    timing_fft += time.time()
    if verbose:
        print("FFT timing for %i DCR planes: [%0.3fs | %0.3fs per plane]"
              % (_i, timing_fft, timing_fft / _i))
    return(return_image[y0:y1:oversample_image, x0:x1:oversample_image] * oversample_image**2)

def fft_convolve(in1, in2, times):
    def _centered(arr, newsize):
        # Return the center newsize portion of the array.
        currsize = np.array(arr.shape)
        startind = (currsize - newsize) // 2
        endind = startind + newsize
        myslice = [slice(startind[k], endind[k]) for k in range(len(endind))]
        return arr[tuple(myslice)]

    if times == 0:
        return in1.copy()


    s1 = np.array(in1.shape)
    s2 = np.array(in2.shape)
    shape = s1 + (s2 - 1) * times

    # Speed up FFT by padding to optimal size for FFTPACK
    fshape = [next_fast_len(int(d)) for d in shape]
    fslice = tuple([slice(0, int(sz)) for sz in shape])

    resfft = fast_power(rfft2(in2, fshape), times)
    resfft = resfft * rfft2(in1, fshape)
    ret = irfft2(resfft, fshape)[fslice].copy()
    ret = ret.real

    return _centered(ret, s1)

def linespectra (arr, freqs, sigma=4, channelWidth=20, kms=False, source_speed=0): #nb sigma is given in px (can be fractional)
    """arr should be an array of shape (x,:,>pix,>pix)
freqs an array or list of nums of length x"""
    shifts=[int(round((freqs[-1]-freqs[i])*299792458/(channelWidth*freqs[-1]))) for i in xrange(len(freqs))]
#    print shifts
    x=[[] for _ in xrange(arr.shape[0])]
    mid=arr.shape[2]/2.0-0.5

    gauss_mask=garray(arr.shape[-2:],sigma)
    s=[y*2 for y in gauss_mask.shape]
    ftg=rfft2(gauss_mask, s)

    for i in xrange(len(x)):
        for j in xrange(arr.shape[1]):
            convolved=irfft2(rfft2(arr[i,j,:,:],s)*ftg)
            x[i].append(convolved[s[0]/2,s[1]/2])

    padding=abs(max(shifts))
    padded=[0 for _ in xrange(arr.shape[1]+padding*2+2)]
    for i in xrange(len(x[0])):
        for j in xrange(len(x)):
            try:
                padded[i+shifts[j]+padding]+=x[j][i]
            except IndexError : 
                print j,i,len(x),len(x[j])
                None
    if kms: return [((i-150)*20/1000.0,x) for i,x in enumerate(padded)]
    else:   return [((i-150)*20,x)        for i,x in enumerate(padded)]

 def getDispl2DImage(self, t0=0, t1=1, Z=0):
     ham = np.hamming(self.get2DShape()[1])*np.atleast_2d(np.hamming(self.get2DShape()[0])).T
     a = rfft2(self.get2DSlice(T=t0, Z=Z)*ham)
     b = rfft2(self.get2DSlice(T=t1, Z=Z)*ham)
     R = numexpr.evaluate(
         'a*complex(real(b), -imag(b)/abs(a*complex(real(b), -imag(b))'
         )
     return irfft2(R)

def register_imgs(imgs,template):
    "save some time by only taking fft of template once"
    rfft2_template_conj = rfft2(template).conj()
    shifts = []
    for img in imgs:
        corr = irfft2(rfft2(img)*rfft2_template_conj)
        shifts.append(balanced_mod(np.unravel_index(corr.argmax(),corr.shape),corr.shape))
    return shifts

def ifct(fu, u):
    """Inverse Cheb transform of x-direction, Fourier in y and z"""
    Uc_hat3[:] = ST.ifct(fu, Uc_hat3)
    comm.Alltoall([Uc_hat3, MPI.DOUBLE_COMPLEX], [U_mpi, MPI.DOUBLE_COMPLEX])
    n0 = U_mpi.shape[2]
    for i in range(num_processes):
        Uc_hatT[:, i*n0:(i+1)*n0] = U_mpi[i]
    u[:] = irfft2(Uc_hatT, axes=(1,2))
    return u

def beam_convolve(arr, sigma):
    "convoles a 2D image with a gaussian profile with sigma in px"
    if len(arr.shape)!=2 or 3*sigma > max(arr.shape): raise ValueError ("arr is not 2d or beam is too wide")
    else: 
        shape=arr.shape
        gauss_mask=garray(shape,sigma)
        s=[y*2 for y in gauss_mask.shape]
        ftg=rfft2(gauss_mask, s)
        return irfft2(rfft2(arr,s)*ftg)

def dnf_response( n, Rn, stimulus, w, we, wi, time, dt ):
		alpha, tau = 0.1, 1.0
		U  = np.random.random((n,n)) * .01
		V  = np.random.random((n,n)) * .01

		V_shape = np.array(V.shape)

		# Computes field input accordingly
		D = (( np.abs( w - stimulus )).sum(axis=-1))/float(Rn*Rn)
		I = ( 1.0 - D.reshape(n,n) ) * alpha

		for j in range( int(time/dt) ):
				Z = rfft2( V * alpha )
				Le = irfft2( Z * we, V_shape).real
				Li = irfft2( Z * wi, V_shape).real
				U += ( -U + ( Le - Li ) + I )* dt * tau
				V = np.maximum( U, 0.0 )
		return V

def conv(im, ker):
    ''' Convolves image im with kernel ker 
        Both image and kernel's dimensions should be even: ker.shape % 2 == 0
    '''
    sy,sx = array(ker.shape)/2
    y0,x0 = array(im.shape)/2
    big_ker = zeros(im.shape)
    big_ker[y0-sy:y0+sy,x0-sx:x0+sx] = ker
    return irfft2(rfft2(im)*rfft2(fftshift(big_ker)))