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 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 _cpu_init(self):

        self.cpu_data = {}
        c = self.cpu_data
        d = self.data

        c['rcore'] = d['rcore'].array
        c['rsurf'] = d['rsurf'].array
        c['im_lsurf'] = d['lsurf'].array

        c['lsurf'] = np.zeros_like(c['rcore'])
        c['clashvol'] = np.zeros_like(c['rcore'])
        c['intervol'] = np.zeros_like(c['rcore'])
        c['interspace'] = np.zeros_like(c['rcore'], dtype=np.int64)

        # complex arrays
        c['ft_shape'] = list(d['shape'])
        c['ft_shape'][-1] = d['shape'][-1]//2 + 1
        c['ft_lsurf'] = np.zeros(c['ft_shape'], dtype=np.complex128)
        c['ft_rcore'] = np.zeros(c['ft_shape'], dtype=np.complex128)
        c['ft_rsurf'] = np.zeros(c['ft_shape'], dtype=np.complex128)

        # initial calculations
        c['ft_rcore'] = rfftn(c['rcore'])
        c['ft_rsurf'] = rfftn(c['rsurf'])
        c['rotmat'] = np.asarray(self.rotations, dtype=np.float64)
        c['weights'] = np.asarray(self.weights, dtype=np.float64)

        c['nrot'] = d['nrot']
        c['shape'] = d['shape']
        c['max_clash'] = d['max_clash']
        c['min_interaction'] = d['min_interaction']
        c['vlength'] = int(np.linalg.norm(self.ligand.coor - \
                self.ligand.center, axis=1).max() + \
                self.interaction_radius + 1.5)/self.voxelspacing
        c['origin'] = d['origin']

        # SAXS arrays
	c['q'] = d['q']
        c['targetIq'] = d['targetIq']
	c['sq'] = d['sq']
        c['base_Iq'] = d['base_Iq']
        c['fifj'] = d['fifj']
        c['rind'] = d['rind']
        c['lind'] = d['lind']
        c['rxyz'] = d['rxyz']
        c['lxyz'] = d['lxyz']

        c['chi2'] = d['chi2']
        c['best_chi2'] = d['best_chi2']
        c['rot_ind'] = np.zeros(d['shape'], dtype=np.int32)

        c['Iq'] = np.zeros_like(c['targetIq'])
        c['tmplxyz'] = np.zeros_like(c['lxyz'])

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 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 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 _cpu_init(self):
        """Initialize all the arrays and data required for a CPU search"""

        self.cpu_data = {}
        c = self.cpu_data
        d = self.data

        # create views of data in cpu_data
        c['rcore'] = d['rcore'].array
        c['rsurf'] = d['rsurf'].array
        c['im_lsurf'] = d['lsurf'].array
        c['restraints'] = d['restraints']

        # allocate arrays used for search
        # real arrays
        c['lsurf'] = np.zeros_like(c['rcore'])
        c['clashvol'] = np.zeros_like(c['rcore'])
        c['intervol'] = np.zeros_like(c['rcore'])
        c['interspace'] = np.zeros_like(c['rcore'], dtype=np.int32)
        c['access_interspace'] = np.zeros_like(c['rcore'], dtype=np.int32)
        c['restspace'] = np.zeros_like(c['rcore'], dtype=np.int32)
        c['violations'] = np.zeros((d['nrestraints'], d['nrestraints']), dtype=np.float64)

        # complex arrays
        c['ft_shape'] = list(d['shape'])
        c['ft_shape'][-1] = d['shape'][-1]//2 + 1
        c['ft_lsurf'] = np.zeros(c['ft_shape'], dtype=np.complex128)
        c['ft_rcore'] = np.zeros(c['ft_shape'], dtype=np.complex128)
        c['ft_rsurf'] = np.zeros(c['ft_shape'], dtype=np.complex128)

        c['rotmat'] = np.asarray(self.rotations, dtype=np.float64)
        c['weights'] = np.asarray(self.weights, dtype=np.float64)

        c['nrot'] = d['nrot']
        c['shape'] = d['shape']
        c['max_clash'] = d['max_clash']
        c['min_interaction'] = d['min_interaction']

        # the vlenght is the longest distance from the ligand center to another
        # atom. This makes the rotation of the ligand object cheaper by only
        # rotation the inner part of the array where there is density
        c['vlength'] = int(np.linalg.norm(self.ligand.coor - \
                self.ligand.center, axis=1).max() + \
                self.interaction_radius + 1.5)/self.voxelspacing

        # initial calculations. Calculate the FFT of the fixed chain objects.
        # This only needs to be done once before the search.
        c['ft_rcore'] = rfftn(c['rcore'])
        c['ft_rsurf'] = rfftn(c['rsurf'])

def multichannel_fftconvolve(x, h, mode='valid'):
    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"

    assert mode == 'valid'
    fshape = (int(2**math.ceil(np.log2((x_len + h_len - 1)))), num_channels)
    x_fft = rfftn(x, fshape)
    h_fft = rfftn(h, fshape)
    ret = _transformed_fft_convolve(x_fft, h_fft)
    
    return ret[h_len-1:x_len, num_channels-1]

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 _setup_kernel(self):
        if not isinstance(self.coordmap, AffineTransform):
            raise ValueError, 'for FFT smoothing, need a regular (affine) coordmap'

        voxels = np.indices(self.bshape).astype(np.float64)

        center = np.asarray(self.bshape)/2
        center = self.coordmap([center[i] for i in range(len(self.bshape))])

        voxels.shape = (voxels.shape[0], np.product(voxels.shape[1:]))
        X = (self.coordmap(voxels.T) - center).T
        X.shape = (self.coordmap.ndim[0],) + tuple(self.bshape)
        kernel = self(X)
        
        kernel = _crop(kernel)
        self.norms = {'l2':np.sqrt((kernel**2).sum()),
                      'l1':np.fabs(kernel).sum(),
                      'l1sum':kernel.sum()}

        self._kernel = kernel

        self.shape = (np.ceil((np.asarray(self.bshape) +
                              np.asarray(kernel.shape))/2)*2+2)
        self.fkernel = np.zeros(self.shape)
        slices = [slice(0, kernel.shape[i]) for i in range(len(kernel.shape))]
        self.fkernel[slices] = kernel
        self.fkernel = fft.rfftn(self.fkernel)

        return kernel

    def _setup_kernel(self):
        # voxel indices of array implied by shape
        voxels = np.indices(self._bshape, dtype=np.float64)

        # coordinates of physical center.  XXX - why the 'floor' here?
        vox_center = np.floor((np.array(self._bshape) - 1) / 2.0)
        phys_center = get_physical_coords(self._affine, vox_center)

        # reshape to (N coordinates, -1).  We appear to need to assign
        # to shape instead of doing a reshape, in order to avoid memory
        # copies
        voxels.shape = (voxels.shape[0], np.product(voxels.shape[1:]))

        # physical coordinates relative to center
        X = get_physical_coords(self._affine,
                                                        voxels) - phys_center

        X.shape = (self._ndims[1],) + tuple(self._bshape)

        # compute kernel from these positions
        kernel = self(X, axis=0)
        kernel = _crop(kernel)

        # compute kernel norm
        self._norm = _get_kernel_norm(kernel, self._normalization)

        self._kernel = kernel
        self._shape = (np.ceil((np.asarray(self._bshape) +
                              np.asarray(kernel.shape)) / 2) * 2 + 2)
        self.fkernel = np.zeros(self._shape)
        slices = [slice(0, kernel.shape[i]) for i in range(kernel.ndim)]
        self.fkernel[slices] = kernel
        self.fkernel = npfft.rfftn(self.fkernel)

        return kernel

def fftn_mpi(u, fu):
    """fft in three directions using mpi
    """
    if num_processes == 1:
        #fu[:] = fft(fft(rfft(u, axis=1), axis=2), axis=0)  
        fu[:] = rfftn(u, axes=(0,2,1))
        return
    
    # Do 2 ffts in y-z directions on owned data
    #ft = fu.transpose(2,1,0)
    #ft[:] = fft(rfft(u, axis=1), axis=2)
    Uc_hatT[:] = rfft2(u, axes=(2,1))
    
    ## Communicating intermediate result 
    ##rstack(ft, Uc_hatT, Np, num_processes)       
    #fu_send = fu.reshape((num_processes, Np, Nf, Np))
    #for i in range(num_processes):
        #if not i == rank:
           #comm.Sendrecv_replace([fu_send[i], MPI.DOUBLE_COMPLEX], i, 0, i, 0)   
    #fu_send[:] = fu_send.transpose(0,3,2,1)
      
    # Transform data to align with x-direction  
    for i in range(num_processes): 
        #U_mpi[i] = ft[:, :, i*Np:(i+1)*Np]
        U_mpi[i] = Uc_hatT[:, :, i*Np:(i+1)*Np]
        
    # Communicate all values
    comm.Alltoall([U_mpi, MPI.DOUBLE_COMPLEX], [fu, MPI.DOUBLE_COMPLEX])  
                
    # Do fft for last direction 
    fu[:] = fft(fu, axis=0)

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 _setup_kernel(self):
     if not isinstance(self.coordmap, AffineTransform):
         raise ValueError('for FFT smoothing, we need a '
                          'regular (affine) coordmap')
     # voxel indices of array implied by shape
     voxels = np.indices(self.bshape).astype(np.float64)
     # coordinates of physical center.  XXX - why the 'floor' here?
     vox_center = np.floor((np.array(self.bshape) - 1) / 2.0)
     phys_center = self.coordmap(vox_center)
     # reshape to (N coordinates, -1).  We appear to need to assign
     # to shape instead of doing a reshape, in order to avoid memory
     # copies
     voxels.shape = (voxels.shape[0], np.product(voxels.shape[1:]))
     # physical coordinates relative to center
     X = (self.coordmap(voxels.T) - phys_center).T
     X.shape = (self.coordmap.ndims[0],) + tuple(self.bshape)
     # compute kernel from these positions
     kernel = self(X, axis=0)
     kernel = _crop(kernel)
     self.norms = {'l2':np.sqrt((kernel**2).sum()),
                   'l1':np.fabs(kernel).sum(),
                   'l1sum':kernel.sum()}
     self._kernel = kernel
     self.shape = (np.ceil((np.asarray(self.bshape) +
                           np.asarray(kernel.shape))/2)*2+2)
     self.fkernel = np.zeros(self.shape)
     slices = [slice(0, kernel.shape[i]) for i in range(len(kernel.shape))]
     self.fkernel[slices] = kernel
     self.fkernel = fft.rfftn(self.fkernel)
     return kernel

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

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 SpatialCorrelationFunctionA(Field1,Field2):
    """
    Corr_12(r) = <Phi_1(r)Phi_2(0)>

    Corr_12(k) = Phi_1(k)* Phi_2(k)/V 
    """ 
    dim = len(Field1.shape)
    if dim == 1:
        V=float(Field1.shape[0])
    elif dim == 2:
        V=float(Field1.shape[0]*Field1.shape[1])
    elif dim == 3:
        V=float(Field1.shape[0]*Field1.shape[1]*Field1.shape[2])
    KField1 = fft.rfftn(Field1).conj()
    KField2 = fft.rfftn(Field2) 
    KCorr = KField1*KField2/V 
    Corr  = fft.irfftn(KCorr)
    return Corr

 def compute_fp(self):
     #n_seg = matrix.shape[1] / 2**16
     #idx = np.arange(2**16)
     n_bands = self.processed.shape[0]
     #result = np.zeros((n_seg,n_bands,61))
     #take 61 instead 60 frequency bins as next step will discard last bin
     fs_model_weights = self.fs_model_curve((11000./2**16) * np.arange(60))
     #use rfftn
     self.processed = np.c_[self.processed,np.zeros(n_bands)]                #calculate fluctuation patter for each frame
     self.processed = fs_model_weights * np.abs(rfftn(self.processed.reshape(n_bands,self.frame_size*2,-1).transpose(2,0,1), axes=(2,))[:,:,:60])