def test_not_last_axis_success(self):
        ar, ai = np.random.random((2, 16, 8, 32))
        a = ar + 1j*ai

        axes = (-2,)

        # Should not raise error
        fft.irfftn(a, axes=axes)

def getCovMatrix(IQdata, lags=100, hp=False):
	# 0: <I1I2> # 1: <Q1Q2> # 2: <I1Q2> # 3: <Q1I2> # 4: <Squeezing> Magnitude # 5: <Squeezing> Phase
	lags = int(lags)
	I1 = np.asarray(IQdata[0])
	Q1 = np.asarray(IQdata[1])
	I2 = np.asarray(IQdata[2])
	Q2 = np.asarray(IQdata[3])
	CovMat = np.zeros([10, lags * 2 + 1])
	start = len(I1) - lags - 1
	stop = len(I1) + lags
	sI1 = np.array(I1.shape)
	shape0 = sI1 * 2 - 1
	fshape = [_next_regular(int(d)) for d in shape0]  # padding to optimal size for FFTPACK
	fslice = tuple([slice(0, int(sz)) for sz in shape0])
	# Do FFTs and get Cov Matrix
	fftI1 = rfftn(I1, fshape)
	fftQ1 = rfftn(Q1, fshape)
	fftI2 = rfftn(I2, fshape)
	fftQ2 = rfftn(Q2, fshape)
	rfftI1 = rfftn(I1[::-1], fshape)  # there should be a simple relationship to fftI1
	rfftQ1 = rfftn(Q1[::-1], fshape)
	rfftI2 = rfftn(I2[::-1], fshape)
	rfftQ2 = rfftn(Q2[::-1], fshape)
	CovMat[0, :] = irfftn((fftI1 * rfftI2))[fslice].copy()[start:stop] / fshape
	CovMat[1, :] = irfftn((fftQ1 * rfftQ2))[fslice].copy()[start:stop] / fshape
	CovMat[2, :] = irfftn((fftI1 * rfftQ2))[fslice].copy()[start:stop] / fshape
	CovMat[3, :] = irfftn((fftQ1 * rfftI2))[fslice].copy()[start:stop] / fshape
	psi = (1j * (CovMat[2, :] + CovMat[3, :]) + (CovMat[0, :] - CovMat[1, :]))
	CovMat[4, :] = abs(psi)
	CovMat[5, :] = np.angle(psi)
	CovMat[6, :] = irfftn((fftI1 * rfftI1))[fslice].copy()[start:stop] / fshape
	CovMat[7, :] = irfftn((fftQ1 * rfftQ1))[fslice].copy()[start:stop] / fshape
	CovMat[8, :] = irfftn((fftI2 * rfftI2))[fslice].copy()[start:stop] / fshape
	CovMat[9, :] = irfftn((fftQ2 * rfftQ2))[fslice].copy()[start:stop] / fshape
	return CovMat

 def c2f(self, x_h):
     'see class doc'
     assert x_h.shape == (self.n*2, self.n*2, self.n+1)
     x_h[:,:,[0,-1]] *= sqrt(2)
     x = irfftn(self.extend(x_h)) * self.n32**3
     x_h[:,:,[0,-1]] /= sqrt(2)
     return x

def computeDisplacement(arry1, arry2):
    """
    Compute an optimal displacement between two ndarrays.

    Finds the displacement between two ndimensional arrays. Arrays must be
    of the same size. Algorithm uses a cross correlation, computed efficiently
    through an n-dimensional fft.

    Parameters
    ----------
    arry1 : ndarray
        The first array

    arry2 : ndarray
        The second array
    """

    from numpy.fft import rfftn, irfftn
    from numpy import unravel_index, argmax

    # compute real-valued cross-correlation in fourier domain
    s = arry1.shape
    f = rfftn(arry1)
    f *= rfftn(arry2).conjugate()
    c = abs(irfftn(f, s))

    # find location of maximum
    inds = unravel_index(argmax(c), s)

    # fix displacements that are greater than half the total size
    pairs = zip(inds, arry1.shape)
    # cast to basic python int for serialization
    adjusted = [int(d - n) if d > n // 2 else int(d) for (d, n) in pairs]

    return adjusted

def hg_conv(f, g):
	ub = tuple(array(f.shape) + array(g.shape) - 1);
	fh = rfftn(cpad(f,ub));
	gh = rfftn(cpad(g,ub));
	res = ifftshift(irfftn(fh * sp.conjugate(gh)));
	del fh, gh;
	return res;

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 customfftconvolve(in1, in2, mode="full", types=('','')):
  """ Pretty much the same as original fftconvolve, but supports
      having operands as fft already 
  """

  in1 = asarray(in1)
  in2 = asarray(in2)

  if in1.ndim == in2.ndim == 0:  # scalar inputs
    return in1 * in2
  elif not in1.ndim == in2.ndim:
    raise ValueError("in1 and in2 should have the same dimensionality")
  elif in1.size == 0 or in2.size == 0:  # empty arrays
    return array([])

  s1 = array(in1.shape)
  s2 = array(in2.shape)
  complex_result = False
  #complex_result = (np.issubdtype(in1.dtype, np.complex) or
  #                  np.issubdtype(in2.dtype, np.complex))
  shape = s1 + s2 - 1
  
  if mode == "valid":
    _check_valid_mode_shapes(s1, s2)

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

  if not complex_result:
    if types[0] == 'fft':
      fin1 = in1#_unfold_fft(in1, fshape)
    else:
      fin1 = rfftn(in1, fshape)

    if types[1] == 'fft':
      fin2 = in2#_unfold_fft(in2, fshape)
    else:
      fin2 = rfftn(in2, fshape)
    ret = irfftn(fin1 * fin2, fshape)[fslice].copy()
  else:
    if types[0] == 'fft':
      fin1 = _unfold_fft(in1, fshape)
    else:
      fin1 = fftn(in1, fshape)
    if types[1] == 'fft':
      fin2 = _unfold_fft(in2, fshape)
    else:
      fin2 = fftn(in2, fshape)
    ret = ifftn(fin1 * fin2)[fslice].copy()

  if mode == "full":
    return ret
  elif mode == "same":
    return _centered(ret, s1)
  elif mode == "valid":
    return _centered(ret, s1 - s2 + 1)
  else:
    raise ValueError("Acceptable mode flags are 'valid',"
                     " 'same', or 'full'.")

def GaussianRandomInitializer(gridShape, sigma=0.2, seed=None, slipSystem=None, slipPlanes=None, slipDirections=None, vacancy=None, smectic=None):

    oldgrid = copy.copy(gridShape)
   
    if len(gridShape) == 1:
	    gridShape = (128,)
    if len(gridShape) == 2:
	    gridShape = (128,128)
    if len(gridShape) == 3:
	    gridShape = (128,128,128)

    """ Returns a random initial set of fields of class type PlasticityState """
    if slipSystem=='gamma':
        state = SlipSystemState.SlipSystemState(gridShape,slipPlanes=slipPlanes,slipDirections=slipDirections)
    elif slipSystem=='betaP':
        state = SlipSystemBetaPState.SlipSystemState(gridShape,slipPlanes=slipPlanes,slipDirections=slipDirections)
    else:
        if vacancy is not None:
            state = VacancyState.VacancyState(gridShape,alpha=vacancy)
        elif smectic is not None:
            state = SmecticState.SmecticState(gridShape)
        else:
            state = PlasticityState.PlasticityState(gridShape)

    field = state.GetOrderParameterField()
    Ksq_prime = FourierSpaceTools.FourierSpaceTools(gridShape).kSq * (-sigma**2/4.)

    if seed is None:
        seed = 0
    n = 0
    random.seed(seed)

    Ksq = FourierSpaceTools.FourierSpaceTools(gridShape).kSq.numpy_array()

    for component in field.components:
        temp = random.normal(scale=gridShape[0],size=gridShape)
        ktemp = fft.rfftn(temp)*(sqrt(pi)*sigma)**len(gridShape)*exp(-Ksq*sigma**2/4.)
        field[component] = numpy.real(fft.irfftn(ktemp))
        #field[component] = GenerateGaussianRandomArray(gridShape, temp ,sigma)
        n += 1

    """
    t, s = LoadState("2dstate32.save", 0)
    for component in field.components:
        for j in range(0,32):
            field[component][:,:,j] = s.betaP[component].numpy_array()
    """

    ## To make seed consistent across grid sizes and convergence comparison
    gridShape = copy.copy(oldgrid)
    if gridShape[0] != 128:
        state = ResizeState(state,gridShape[0],Dim=len(gridShape))

    state = ReformatState(state)
    state.ktools = FourierSpaceTools.FourierSpaceTools(gridShape)
    
    return state 

def getCovMatrix(I1, Q1, I2, Q2, lags=20):
    # calc <I1I2>, <I1Q2>, Q1I2, Q1Q2
    lags = int(lags)
    I1 = np.asarray(I1)
    Q1 = np.asarray(Q1)
    I2 = np.asarray(I2)
    Q2 = np.asarray(Q2)
    CovMat = np.zeros([6, lags*2-1])
    start = len(I1*2-1)-lags
    stop = len(I1*2-1)-1+lags
    sI1 = np.array(I1.shape)
    sQ2 = np.array(Q2.shape)
    shape = sI1 + sQ2 - 1
    HPfilt = (int(sI1/(lags*4)))  # smallest features visible is lamda/4
    fshape = [_next_regular(int(d)) for d in shape]  # padding to optimal size for FFTPACK
    fslice = tuple([slice(0, int(sz)) for sz in shape])
    # Do FFTs and get Cov Matrix
    fftI1 = rfftn(I1, fshape)
    fftQ1 = rfftn(Q1, fshape)
    fftI2 = rfftn(I2, fshape)
    fftQ2 = rfftn(Q2, fshape)
    rfftI1 = rfftn(I1[::-1], fshape)
    rfftQ1 = rfftn(Q1[::-1], fshape)
    rfftI2 = rfftn(I2[::-1], fshape)
    rfftQ2 = rfftn(Q2[::-1], fshape)
    # filter frequencies outside the lags range
    fftI1 = np.concatenate((np.zeros(HPfilt), fftI1[HPfilt:]))
    fftQ1 = np.concatenate((np.zeros(HPfilt), fftQ1[HPfilt:]))
    fftI2 = np.concatenate((np.zeros(HPfilt), fftI2[HPfilt:]))
    fftQ2 = np.concatenate((np.zeros(HPfilt), fftQ2[HPfilt:]))
    # filter frequencies outside the lags range
    rfftI1 = np.concatenate((np.zeros(HPfilt), rfftI1[HPfilt:]))
    rfftQ1 = np.concatenate((np.zeros(HPfilt), rfftQ1[HPfilt:]))
    rfftI2 = np.concatenate((np.zeros(HPfilt), rfftI2[HPfilt:]))
    rfftQ2 = np.concatenate((np.zeros(HPfilt), rfftQ2[HPfilt:]))
    CovMat[0, :] = (irfftn((fftI1*rfftI2))[fslice].copy()[start:stop] / len(fftI1))  # 0: <I1I2>
    CovMat[1, :] = (irfftn((fftQ1*rfftQ2))[fslice].copy()[start:stop] / len(fftI1))  # 1: <Q1Q2>
    CovMat[2, :] = (irfftn((fftI1*rfftQ2))[fslice].copy()[start:stop] / len(fftI1))  # 2: <I1Q2>
    CovMat[3, :] = (irfftn((fftQ1*rfftI2))[fslice].copy()[start:stop] / len(fftI1))  # 3: <Q1I2>
    CovMat[4, :] = (abs(1j*(CovMat[2, :]+CovMat[3, :]) + (CovMat[0, :] - CovMat[1, :])))  # 4: <Squeezing> Magnitude
    CovMat[5, :] = np.angle(1j*(CovMat[2, :]+CovMat[3, :]) + (CovMat[0, :] - CovMat[1, :]))  # 5: <Squeezing> Angle
    return CovMat

def GenerateGaussianRandomArray(gridShape, temp, sigma):
    dimension = len(gridShape)
    if dimension == 1:
        kfactor = fromfunction(lambda kz: exp(-0.5*(sigma*kz)**2),[gridShape[0]/2+1,])
        ktemp = fft.rfft(temp)
        ktemp *= kfactor
        data = fft.irfft(ktemp)
    elif dimension == 2:
        X,Y = gridShape
        kfactor = fromfunction(lambda kx,ky: exp(-0.5*sigma**2*((kx*(kx<X/2)+(X-kx)*(kx>=X/2))**2+ky**2)),[X,Y/2+1])
        ktemp = fft.rfftn(temp)
        ktemp *= kfactor
        data = fft.irfftn(ktemp)
    elif dimension == 3:
        X,Y,Z = gridShape
        kfactor = fromfunction(lambda kx,ky,kz: exp(-0.5*sigma**2*( (kx*(kx<X/2)+(X-kx)*(kx>=X/2))**2 + \
                                                (ky*(ky<Y/2)+(Y-ky)*(ky>=Y/2))**2 + kz**2)),[X,Y,Z/2+1])
        ktemp = fft.rfftn(temp)
        ktemp *= kfactor
        data = fft.irfftn(ktemp)
    return data 

def zeroextract2d(N,filename):
    a = fromfile(filename)
    a = a.reshape(9,N,N)
    b = numpy.zeros((9,N/2,N/2),float)
    for i in range(9):
        ka = fft.rfftn(a[i])
        kb = numpy.zeros((N/2,N/4+1),complex)
        kb[:N/4,:]=ka[:N/4,:N/4+1]
        kb[-N/4:,:]=ka[-N/4:,:N/4+1]
        b[i] = fft.irfftn(kb)
    b /= 4.
    b.tofile(filename.replace(str(N),str(N/2)))

def SpatialCorrelationFunctionA(Field1, Field2):
    """
    Designed for Periodic Boundary Condition.

    Corr_12(r) = <Phi_1(r)Phi_2(0)>
    Corr_12(k) = Phi_1(k)* Phi_2(k)/V 
    """
    V = float(numpy.array(Field1.shape).prod())
    KField1 = fft.rfftn(Field1).conj()
    KField2 = fft.rfftn(Field2)
    KCorr = KField1 * KField2 / V
    Corr = fft.irfftn(KCorr)
    return Corr

def get_covariance_submatrix(IQdata, lags, q):
    logging.debug('Calculating Submatrix')
    I1 = np.asarray(IQdata[0])
    Q1 = np.asarray(IQdata[1])
    I2 = np.asarray(IQdata[2])
    Q2 = np.asarray(IQdata[3])
    lags = int(lags)
    start = len(I1) - lags - 1
    stop = len(I1) + lags
    sub_matrix = np.zeros([4, lags * 2 + 1])
    sI1 = np.array(I1.shape)
    shape0 = sI1*2 - 1
    fshape = [_next_regular(int(d)) for d in shape0]  # padding to optimal size for FFTPACK
    fslice = tuple([slice(0, int(sz)) for sz in shape0])  # remove padding later
    fftI1 = rfftn(I1, fshape)/fshape
    fftQ1 = rfftn(Q1, fshape)/fshape
    rfftI2 = rfftn(I2[::-1], fshape)/fshape
    rfftQ2 = rfftn(Q2[::-1], fshape)/fshape
    sub_matrix[0] = irfftn((fftI1 * rfftI2))[fslice].copy()[start:stop]  # <II>
    sub_matrix[1] = irfftn((fftI1 * rfftQ2))[fslice].copy()[start:stop]  # <IQ>
    sub_matrix[2] = irfftn((fftQ1 * rfftI2))[fslice].copy()[start:stop]  # <QI>
    sub_matrix[3] = irfftn((fftQ1 * rfftQ2))[fslice].copy()[start:stop]  # <QQ>
    q.put(sub_matrix)

def fftconvolve(in1, in2):
    """Convolve two N-dimensional arrays using FFT.

    This is a modified version of the scipy.signal.fftconvolve.
    The new feature is derived from the fftconvolve algorithm used in the IDL package.

    Parameters
    ----------
    in1 : array_like
        First input.
    in2 : array_like
        Second input. Should have the same number of dimensions as `in1`;
        if sizes of `in1` and `in2` are not equal then `in1` has to be the
        larger array.

    Returns
    -------
    out : array
        An N-dimensional array containing a subset of the discrete linear
        convolution of `in1` with `in2`.

    """
    in1 = asarray(in1)
    in2 = asarray(in2)

    if matrix_rank(in1) == matrix_rank(in2) == 0:  # scalar inputs
        return in1 * in2
    elif not in1.ndim == in2.ndim:
        raise ValueError("in1 and in2 should have the same rank")
    elif in1.size == 0 or in2.size == 0:  # empty arrays
        return array([])

    s1 = np.array(in1.shape)
    s2 = np.array(in2.shape)
    complex_result = (np.issubdtype(in1.dtype, np.complex) or
                      np.issubdtype(in2.dtype, np.complex))

    fsize = s1

    fslice = tuple([slice(0, int(sz)) for sz in fsize])
    if not complex_result:
        ret = irfftn(rfftn(in1, fsize) *
                     rfftn(in2, fsize), fsize)[fslice].copy()
        ret = ret.real
    else:
        ret = ifftn(fftn(in1, fsize) * fftn(in2, fsize))[fslice].copy()

    shift = array([int(floor(fsize[0]/2.0)), int(floor(fsize[1]/2.0))])
    ret   = roll(roll(ret, -shift[0], axis=0), -shift[1], axis=1)
    return ret

def KspaceInitializer(gridShape, fileDictionary, state = None):
    """
    Initialize plasticity state by reading from files given by the file
    dictionary.

    File dictionary must be a dictionary(hash) with component of field as
    its keys and filename pair(for real and imaginary) as its values.
   
    State must first be initialized and passed in for non default plasticity
    states.
 
    example:
    fileDictionary = {('x','z') : \
        ("InitialConditions/InitFourierCoeffXZ_real256.dat", \
         "InitialConditions/InitFourierCoeffXZ_im256.dat"),\
        ('y','z') : ("InitialConditions/InitFourierCoeffYZ_real256.dat", \
         "InitialConditions/InitFourierCoeffYZ_im256.dat"),\
        ('z','x') : ("InitialConditions/InitFourierCoeffZX_real256.dat", \
         "InitialConditions/InitFourierCoeffZX_im256.dat")}
    """
    if state is None:
        state = PlasticityState.PlasticityState(gridShape)
    field = state.GetOrderParameterField() 
 
    kGridShape = list(gridShape)
    kGridShape[-1] = int(kGridShape[-1]/2)+1 
    kGridShape = tuple(kGridShape)
    totalSize = 1
    for sz in kGridShape:
        totalSize *= sz  

    for component in fileDictionary: 
        rePart = numpy.fromfile(fileDictionary[component][0])
        imPart = numpy.fromfile(fileDictionary[component][1])
        kSpaceArray = rePart + 1.0j*imPart
       
        """
        Strip down only first half rows as this is the only data that gets
        used. Notice that we are actually taking real fourier transform on
        last axis, nevertheless we only take half data from the top. i.e.
        this may not very intuitive, but is coded this way for compatibility
        with Yor's old C++ version. i.e., strip down only half rows, and re-
        arrange so that column is halved. (That is, second half of first row
        becomes the second row and first half of second row becomes the third.
        """
        kSpaceArray = kSpaceArray[:totalSize].reshape(kGridShape)
       
        field[component] = fft.irfftn(kSpaceArray)
    return state

    def _cpu_search(self):

        d = self.data
        c = self.cpu_data

        time0 = _time()
        for n in xrange(c['rotmat'].shape[0]):
            # rotate ligand image
            rotate_image3d(c['im_lsurf'], c['vlength'], 
                    np.linalg.inv(c['rotmat'][n]), d['im_center'], c['lsurf'])

            c['ft_lsurf'] = rfftn(c['lsurf']).conj()
            c['clashvol'] = irfftn(c['ft_lsurf'] * c['ft_rcore'], s=c['shape'])
            c['intervol'] = irfftn(c['ft_lsurf'] * c['ft_rsurf'], s=c['shape'])

            np.logical_and(c['clashvol'] < c['max_clash'],
                           c['intervol'] > c['min_interaction'],
                           c['interspace'])


	    print('Number of complexes to analyze: ', c['interspace'].sum())
            c['chi2'].fill(0)
            calc_chi2(c['interspace'], c['q'], c['base_Iq'], 
                    c['rind'], c['rxyz'], c['lind'], (np.mat(c['rotmat'][n])*np.mat(c['lxyz']).T).T, 
                    c['origin'], self.voxelspacing, 
                    c['fifj'], c['targetIq'], c['sq'], c['chi2'])

            ind = c['chi2'] > c['best_chi2']
            c['best_chi2'][ind] = c['chi2'][ind]
            c['rot_ind'][ind] = n

            if _stdout.isatty():
                self._print_progress(n, c['nrot'], time0)

        d['best_chi2'] = c['best_chi2']
        d['rot_ind'] = c['rot_ind']

def SpatialCorrelationFunction(Field1,Field2):
    """
	Designed for Periodic Boundary Condition.
	
	.. math::
		C_{12}(r) = <F_1(r) F_2(0)>
	
		C_{12}(k) = F_1(k)* F_2(k)/V 
    """ 
    V = float(numpy.array(Field1.shape).prod())
    KField1 = fft.rfftn(Field1).conj()
    KField2 = fft.rfftn(Field2) 
    KCorr = KField1*KField2/V 
    Corr  = fft.irfftn(KCorr)
    return Corr

def multichannel_overlap_add_fftconvolve(x, h, mode='valid'):
    """Given a signal x compute the convolution with h using the overlap-add algorithm.

    This is an fft based convolution algorithm optimized to work on a signal x 
    and filter, h, where x is much longer than h. 

    Input:
    x: float array with dimensions (N, num_channel)
    h: float array with dimensions (filter_len, num_channel)
    mode: only accepts valid - same profile scipy.convolve 

    Returns:
    float array of length N-filter_len+1 (for mode = valid)
    """
    assert mode == 'valid'
    
    # pad x so that the boundaries are dealt with correctly
    x_len = x.shape[0]
    num_channels = h.shape[1]
    h_len = h.shape[0]
    assert x.shape[1] == num_channels
    assert x_len >= h_len, \
        "The signal needs to be at least as long as the filter"
    
    #x = np.vstack((np.zeros((h_len, num_channels)), 
    #               x, 
    #               np.zeros((h_len, num_channels))))    
    # make sure that the desired block size is long enough to capture the motif
    block_size = max(2**OVERLAP_ADD_BLOCK_POWER, h_len)
    N = int(2**math.ceil(np.log2(block_size+h_len-1)))
    step_size = N-h_len+1
    
    H = rfftn(h,(N,num_channels))
    n_blocks = int(math.ceil(float(len(x))/step_size))
    y = np.zeros((n_blocks+1)*step_size)
    for block_index in xrange(n_blocks):
        start = block_index*step_size
        yt = irfftn( rfftn(x[start:start+step_size,:],(N, num_channels))*H, 
                     (N, num_channels) )
        y[start:start+N] += yt[:,num_channels-1]

    #y = y[h_len:2*h_len+x_len-1]
    if mode == 'full':
        return y
    elif mode == 'valid':
        return y[h_len-1:x_len]
    elif mode == 'same':
        raise NotImplementedError, "'same' mode is not implemented"

def cross_correlation(seqs):
    # deal with the shape, and upcast to the next reasonable shape
    shape = np.array(seqs.shape[1:]) + np.array(seqs.shape[1:]) - 1
    fshape = [next_good_fshape(x) for x in shape]
    fslice = tuple([slice(0, int(sz)) for sz in shape])
    flipped_seqs_fft = np.zeros([seqs.shape[0],] + fshape[:-1] + [fshape[-1]//2+1,], dtype='complex')
    for i in xrange(seqs.shape[0]):
        rev_slice = tuple([i,] + [slice(None, None, -1) for sz in shape])
        flipped_seqs_fft[i] = rfftn(seqs[rev_slice], fshape)
    rv = np.zeros((seqs.shape[0], seqs.shape[0]), dtype='float32')
    for i in xrange(seqs.shape[0]):
        fft_seq = rfftn(seqs[i], fshape)
        for j in xrange(i+1, seqs.shape[0]):
            rv[i,j] = irfftn(fft_seq*flipped_seqs_fft[j], fshape)[fslice].max()
            #print rv[i,j], correlate(seqs[i], seqs[j]).max()
    return rv

def get_g2(P1, P2, lags=20):
    ''' Returns the Top part of the G2 equation (<P1P2> - <P1><P2>)'''
    lags = int(lags)
    P1 = np.asarray(P1)
    P2 = np.asarray(P2)
    # G2 = np.zeros([lags*2-1])

    start = len(P1*2-1)-lags
    stop = len(P1*2-1)-1+lags

    # assume I1 Q1 have the same shape
    sP1 = np.array(P1.shape)
    complex_result = np.issubdtype(P1.dtype, np.complex)
    shape = sP1 - 1
    HPfilt = (int(sP1/(lags*4)))  # smallest features visible is lamda/4

    # Speed up FFT by padding to optimal size for FFTPACK
    fshape = [_next_regular(int(d)) for d in shape]
    fslice = tuple([slice(0, int(sz)) for sz in shape])
    # Pre-1.9 NumPy FFT routines are not threadsafe.  For older NumPys, make
    # sure we only call rfftn/irfftn from one thread at a time.
    if not complex_result and _rfft_lock.acquire(False):
        try:
            fftP1 = rfftn(P1, fshape)
            rfftP2 = rfftn(P2[::-1], fshape)
            fftP1 = np.concatenate((np.zeros(HPfilt), fftP1[HPfilt:]))
            rfftP2 = np.concatenate((np.zeros(HPfilt), rfftP2[HPfilt:]))
            G2 = irfftn((fftP1*rfftP2))[fslice].copy()[start:stop]/len(fftP1)
            return 

        finally:
            _rfft_lock.release()

    else:
        # If we're here, it's either because we need a complex result, or we
        # failed to acquire _rfft_lock (meaning rfftn isn't threadsafe and
        # is already in use by another thread).  In either case, use the
        # (threadsafe but slower) SciPy complex-FFT routines instead.
        # ret = ifftn(fftn(in1, fshape) * fftn(in2, fshape))[fslice].copy()
        print 'Abort, reason:complex input or Multithreaded FFT not available'

        if not complex_result:
            pass  # ret = ret.real

    P12var = np.var(P1)*np.var(P2)
    return G2-P12var