def idft(dft_vec, dt=1.0):
        """
        computes the inverse DFT of vec
        takes in the one-sided spectrum
        """
        N = len(dft_vec) ### if N is even, then n is even
                         ### if N is odd, then n is odd

        if N%2: ### if N is odd, n is odd
                n = 2*N-1
        else: ### if N is even, n is even
                n = 2*N

        seglen = n*dt ### length of time series

        vec = np.empty((n,), complex)
        vec[:N] = dft_vec
        if n%2: ### odd number of points
                vec[N:] = np.conjugate(dft_vec[1:])[::-1]
        else: ### even number of points
                vec[N:] = np.conjugate(dft_vec)[::-1]

        vec = np.fft.ifft( vec ) / seglen
        time = np.arange(0, seglen, dt)

        return vec, time

 def _npBatchMatmul(self, x, y, adjoint_a, adjoint_b):
   # output's shape depends on adj[0] and adj[1]
   if adjoint_a:
     x = np.conjugate(np.swapaxes(x, -1, -2))
   if adjoint_b:
     y = np.conjugate(np.swapaxes(y, -1, -2))
   return np.matmul(x, y)

def FFT_Correlation(x,y):
    """
    FFT-based correlation, much faster than numpy autocorr.
    x and y are row-based vectors of arbitrary lengths.
    This is a vectorized implementation of O(N*log(N)) flops.
    """

    lengthx = x.shape[0]
    lengthy = y.shape[0]

    x = np.reshape(x,(1,lengthx))
    y = np.reshape(y,(1,lengthy))

    length = np.array([lengthx, lengthy]).min()
    
    x = x[:length]
    y = y[:length]
    
    fftx = fft(x, 2 * length - 1, axis=1) #pad with zeros
    ffty = fft(y, 2 * length - 1, axis=1)

    corr_xy = fft.ifft(fftx * np.conjugate(ffty), axis=1)
    corr_xy = np.real(fft.fftshift(corr_xy, axes=1)) #should be no imaginary part

    corr_yx = fft.ifft(ffty * np.conjugate(fftx), axis=1)
    corr_yx = np.real(fft.fftshift(corr_yx, axes=1))

    corr = 0.5 * (corr_xy[:,length:] + corr_yx[:,length:]) / range(1,length)[::-1]
    return np.reshape(corr,corr.shape[1])

def make_topo(a, A, B, C, D, M, E, t, sigma_0,mixing,iterations,N, n_layers):
	"""
	Creates a topological insulator from given parameters.
	"""	
	E_s = C + M - 4*((B+D)/(a**2))
	E_p = C - M - 4*((D-B)/(a**2))	
	V_ss = (B+D)/(a**2)
	V_pp = (D-B)/(a**2)
	V_sp = (-1)*(1.j)*(A)/(2*a)
	h_topo = np.diag([E_s,E_p,E_s,E_p])
	t_y = np.zeros((4,4), dtype=np.complex_)
	t_y[0,0] = V_ss
	t_y[0,1] = (1.j)*V_sp
	t_y[1,0] = (1.j)*np.conjugate(V_sp)
	t_y[1,1] = V_pp
	t_y[2,2] = V_ss
	t_y[2,3] = (-1)*t_y[1,0]
	t_y[3,2] = (-1)*t_y[0,1]
	t_y[3,3] = V_pp
	t_x = np.zeros((4,4), dtype=np.complex_)
	t_x[0,0] = V_ss
	t_x[0,1] = 1*V_sp
	t_x[1,0] = (-1)*np.conjugate(V_sp)
	t_x[1,1] = V_pp
	t_x[2,2] = V_ss
	t_x[2,3] = 1*np.conjugate(V_sp)
	t_x[3,2] = (-1)*V_sp
	t_x[3,3] = V_pp
	for arr in [h_topo,t_x,t_y]:
		arr = swap_cols(arr,1,2)
		arr = swap_rows(arr,1,2)
	topo = constructor.Constructor(E,h_topo,t_y,sigma_0,mixing,iterations,N, n_layers, t_x)
	return (topo, t_x)

    def mix_parameters(self, Pibra, Piket):
        r"""Mix the two parameter sets :math:`\Pi_i` and :math:`\Pi_j`
        from the 'bra' and the 'ket' wavepackets :math:`\Phi\left[\Pi_i\right]`
        and :math:`\Phi^\prime\left[\Pi_j\right]`.

        :param Pibra: The parameter set :math:`\Pi_i` from the bra part wavepacket.
        :param Piket: The parameter set :math:`\Pi_j` from the ket part wavepacket.
        :return: The mixed parameters :math:`q_0` and :math:`Q_S`. (See the theory for details.)
        """
        # <Pibra | ... | Piket>
        qr, pr, Qr, Pr = Pibra
        qc, pc, Qc, Pc = Piket

        # Mix the parameters
        Gr = dot(Pr, inv(Qr))
        Gc = dot(Pc, inv(Qc))

        r = imag(Gc - conjugate(Gr.T))
        s = imag(dot(Gc, qc) - dot(conjugate(Gr.T), qr))

        q0 = dot(inv(r), s)
        Q0 = 0.5 * r

        # Here we can not avoid the matrix root by using svd
        Qs = inv(sqrtm(Q0))

        return (q0, Qs)

    def exact_result_ground(self, Pibra, Piket, eps):
        r"""Compute the overlap integral :math:`\langle \phi_0 | \phi_0 \rangle` of
        the groundstate :math:`\phi_0` by using the symbolic formula:

        .. math::
            \langle \phi_0 | \phi_0 \rangle =
            \sqrt{\frac{-2 i}{Q_2 \overline{P_1} - P_2 \overline{Q_1}}} \cdot
              \exp \Biggl(
                \frac{i}{2 \varepsilon^2}
                \frac{Q_2 \overline{Q_1} \left(p_2-p_1\right)^2 + P_2 \overline{P_1} \left(q_2-q_1\right)^2}
                      {\left(Q_2 \overline{P_1} - P_2 \overline{Q_1}\right)}
              \\
              -\frac{i}{\varepsilon^2}
              \frac{\left(q_2-q_1\right) \left( Q_2 \overline{P_1} p_2 - P_2 \overline{Q_1} p_1\right)}
                   {\left(Q_2 \overline{P_1} - P_2 \overline{Q_1}\right)}
              \Biggr)

        Note that this is an internal method and usually there is no
        reason to call it from outside.

        :param Pibra: The parameter set :math:`\Pi = \{q_1,p_1,Q_1,P_1\}` of the bra :math:`\langle \phi_0 |`.
        :param Piket: The parameter set :math:`\Pi^\prime = \{q_2,p_2,Q_2,P_2\}` of the ket :math:`| \phi_0 \rangle`.
        :param eps: The semi-classical scaling parameter :math:`\varepsilon`.
        :return: The value of the integral :math:`\langle \phi_0 | \phi_0 \rangle`.
        """
        q1, p1, Q1, P1 = Pibra
        q2, p2, Q2, P2 = Piket
        hbar = eps**2
        X = Q2*conjugate(P1) - P2*conjugate(Q1)
        I = sqrt(-2.0j/X) * exp( 1.0j/(2*hbar) * (Q2*conjugate(Q1)*(p2 - p1)**2 + P2*conjugate(P1)*(q2 - q1)**2) / X
                                -1.0j/hbar *     ((q2 - q1)*(Q2*conjugate(P1)*p2 - P2*conjugate(Q1)*p1)) / X
                               )
        return I

def get_moments(v,m,n=100,use_fortran=use_fortran):
  """ Get the first n moments of a certain vector
  using the Chebychev recursion relations"""
  if use_fortran:
    from kpmf90 import get_momentsf90 # fortran routine
    mo = coo_matrix(m) # convert to coo matrix
    vo = v.todense() # convert to conventional vector
    vo = np.array([vo[i,0] for i in range(len(vo))])
# call the fortran routine
    mus = get_momentsf90(mo.row+1,mo.col+1,mo.data,vo,n) 
    return mus # return fortran result
  else:
    mus = np.array([0.0j for i in range(2*n)]) # empty arrray for the moments
    a = v.copy() # first vector
    am = v.copy() # zero vector
    a = m*v  # vector number 1
    bk = (np.transpose(np.conjugate(v))*v)[0,0] # scalar product
    bk1 = (np.transpose(np.conjugate(v))*a)[0,0] # scalar product
    mus[0] = bk  # mu0
    mus[1] = bk1 # mu1
    for i in range(1,n): 
      ap = 2*m*a - am # recursion relation
      bk = (np.transpose(np.conjugate(a))*a)[0,0] # scalar product
      bk1 = (np.transpose(np.conjugate(ap))*a)[0,0] # scalar product
      mus[2*i] = 2.*bk
      mus[2*i+1] = 2.*bk1
      am = a +0. # new variables
      a = ap+0. # new variables
    mu0 = mus[0] # first
    mu1 = mus[1] # second
    for i in range(1,n): 
      mus[2*i] +=  - mu0
      mus[2*i+1] += -mu1 
    return mus

    def mfunc(uv, p, d):
        crd,t,(i,j) = p
        p1,p2 = a.miriad.pol2str[uv['pol']]
        #if i == j and (p1,p2) == ('y','x'): return p, None, None

#        if is_run1(t):
        ni = rewire_run1[i][p1]
        nj = rewire_run1[j][p2]
        if t!= curtime:
            aa.set_jultime(t)
            uvo['lst'] = aa.sidereal_time()
            uvo['ra'] = aa.sidereal_time()
            uvo['obsra'] = aa.sidereal_time()
#        else: return p, None
        if ni > nj:
            ni,nj = nj,ni
            
            if (p1,p2) != ('y','x'): d = n.conjugate(d)
        elif ni < nj and (p1,p2) == ('y','x'):
            d = n.conjugate(d)

        p = crd,t,(ni,nj)
#        print t,i,j,a.miriad.pol2str[uv['pol']],'->',ni,nj,'xx'
        d[orbchan].mask = 1
        return p,d

 def testDelta_simmetric(self):
     m=1;
     c1 = 2; n1 = 4
     c2 = 3; n2 = 4;
     func = lambda x: pro_ang1(m, n1, c1, x)[0]  * numpy.conjugate(pro_ang1(m, n2, c2, x)[0])
     func2 = lambda x: pro_ang1(m, n2, c2, x)[0]  * numpy.conjugate(pro_ang1(m, n1, c1, x)[0])
     self.assertAlmostEqual(quad(func, -1, 1), quad(func2, -1, 1), places = 7)

    def perform_quadrature(self, row, col):
        r"""Evaluates the integral :math:`\langle \Phi_i | \Phi^\prime_j \rangle`
        by an exact symbolic formula.

        .. warning:: This method does only take into account the ground state
                     basis components :math:`\phi_{\underline{0}}` from both,
                     the 'bra' and the 'ket'. If the wavepacket :math:`\Phi`
                     contains higher order basis functions :math:`\phi_{\underline{k}}`
                     with non-zero coefficients :math:`c_{\underline{k}}`, the inner products
                     computed are wrong! There is also no warning about that.

        :param row: The index :math:`i` of the component :math:`\Phi_i` of :math:`\Psi`.
        :param row: The index :math:`j` of the component :math:`\Phi^\prime_j` of :math:`\Psi^\prime`.
        :return: A single complex floating point number.
        """
        eps = self._packet.get_eps()
        D = self._packet.get_dimension()

        Pibra = self._pacbra.get_parameters(component=row)
        Piket = self._packet.get_parameters(component=col)
        cbra = self._pacbra.get_coefficient_vector(component=row)
        cket = self._packet.get_coefficient_vector(component=col)
        Kbra = self._pacbra.get_basis_shapes(component=row)
        Kket = self._packet.get_basis_shapes(component=col)

        phase = exp(1.0j/eps**2 * (Piket[4]-conjugate(Pibra[4])))

        z = tuple(D*[0])
        cr = cbra[Kbra[z],0]
        cc = cket[Kket[z],0]
        i = self.exact_result_gauss(Pibra[:4], Piket[:4], D, eps)
        result = phase * conjugate(cr) * cc * i

        return result

    def quadrature(self, lcket, operator=None, component=None):
        r"""Delegates the evaluation of :math:`\langle\Upsilon|f|\Upsilon\rangle` for a general
        function :math:`f(x)` with :math:`x \in \mathbb{R}^D`.

        :param lcket: The linear combination :math:`\Upsilon` with :math:`J` summands :math:`\Psi_j`.
        :param operator: A matrix-valued function :math:`f(x): \mathbb{R}^D \rightarrow \mathbb{R}^{N \times N}`.
        :return: The value of :math:`\langle\Upsilon|f|\Upsilon\rangle`.
        :type: An :py:class:`ndarray`.
        """
        J = lcket.get_number_packets()
        packets = lcket.get_wavepackets()

        M = zeros((J, J), dtype=complexfloating)

        # Elements below the diagonal
        for row, pacbra in enumerate(packets):
            for col, packet in enumerate(packets[:row]):
                if self._obey_oracle:
                    if self._oracle.is_not_zero(pacbra, packet):
                        # TODO: Handle multi-component packets
                        M[row, col] = self._quad.quadrature(pacbra, packet, operator=operator, component=0)
                else:
                    # TODO: Handle multi-component packets
                    M[row, col] = self._quad.quadrature(pacbra, packet, operator=operator, component=0)

        M = M + conjugate(transpose(M))

        # Diagonal Elements
        for d, packet in enumerate(packets):
            # TODO: Handle multi-component packets
            M[d, d] = self._quad.quadrature(packet, packet, operator=operator, component=0)

        c = lcket.get_coefficients()

        return dot(conjugate(transpose(c)), dot(M, c))

def verify_gaus_sum_complex_conj(N):

    ALMOST_ZERO = 0.001

    # get dirich_char_mat, which is a 2D matrix
    dirich_char_mat, phi_n, coprime_list, prim_char_stat_all = dirichi_char_for_n(N)

    # number of row. Note: 1st row is principle
    for i in range(0, phi_n):

        chi = dirich_char_mat[i]

        gaus_sum = gaus_sum_for_dirich_char(chi)
        chi_conj = np.conjugate(chi)
        gaus_sum_for_chi_conj = gaus_sum_for_dirich_char(chi_conj)

        gaus_sum_conj = np.conjugate(gaus_sum)

        # chi(-1) = chi(N-1)
        # It seems this value is always be 1 or -1, why ?
        chi_minus_1 = chi[N-1]

        tmp = chi_minus_1*gaus_sum_conj

        assert(np.absolute(gaus_sum_for_chi_conj - tmp) < ALMOST_ZERO)
        # print 'chi      = ', chi
        # print 'chi_conj = ', chi_conj
        # print 'N = %d, i = %d'%(N, i), ', gaus_sum_for_chi_conj = ', gaus_sum_for_chi_conj, ', tmp = ', tmp, ', chi_minus_1 = ', chi_minus_1

    print 'verify_gaus_sum_complex_conj() pass !'

  def update_expectation_values(self):
    """Calculate the expectation values of the different operators"""
    # this conjugate comes from being inconsistent
    # in the routines to calculate exectation values
    voccs = np.conjugate(self.wavefunctions) # get wavefunctions
    ks = self.kvectors # kpoints
    mode = self.correlator_mode # 
#    mode = "1by1"
    if mode=="plain": # conventional mode
      for v in self.interactions:
        v.vav = (voccs*v.a*voccs.H).trace()[0,0]/self.kfac # <vAv>
        v.vbv = (voccs*v.b*voccs.H).trace()[0,0]/self.kfac # <vBv>
    elif mode=="1by1": # conventional mode
      for v in self.interactions:
        phis = [self.hamiltonian.geometry.bloch_phase(v.dir,k*0.) for k in ks]
        v.vav = meanfield.expectation_value(voccs,v.a,np.conjugate(phis))/self.kfac # <vAv>
        v.vbv = meanfield.expectation_value(voccs,v.b,phis)/self.kfac # <vBv>
      self.v2cij() # update the v vector
    elif mode=="multicorrelator": # multicorrelator mode
      numc = len(self.interactions)*2 # number of correlators
      if self.bloch_multicorrelator:
        cs = multicorrelator_bloch(voccs,ks,self.lamb,self.ijk,self.dir,numc)
      else: cs = multicorrelator(voccs,self.lamb,self.ijk,numc)
      self.cij = cs/self.kfac # store in the object, already normalized
      self.cij2v() # update the expectation values
    else: raise

    def calclogI(self):
        """
        The logarithm intensity function.

        Returns:
        Numpy.complex data type representing the value of the logarithm of the intensity function.
        """
        ret=numpy.complex(0.,0.)
        for n in range(0,len(self.alphaList)-1,1):    
            argret=numpy.complex(0.,0.)            
            for wave1 in self.waves:
                for wave2 in self.waves:
                    if len(self.productionAmplitudes)!=0:
                                #logarithmic domain error
                        arg = self.productionAmplitudes[self.waves.index(wave1)]*numpy.conjugate(self.productionAmplitudes[self.waves.index(wave2)])*wave1.complexamplitudes[n]*numpy.conjugate(wave2.complexamplitudes[n])*spinDensity(self.beamPolarization,self.alphaList[n])[wave1.epsilon,wave2.epsilon]
                        argret+=arg
            argret=argret.real
            if self.debugPrinting==1:                        
                print"loop#",n,"="*10
                print"argval:",arg
                print"argtype:",type(arg)
                print"productionAmps1:",self.productionAmplitudes[self.waves.index(wave1)]
                print"productionAmps2*:",numpy.conjugate(self.productionAmplitudes[self.waves.index(wave2)])
                print"spinDensityValue:",spinDensity(self.beamPolarization,self.alphaList[n])[wave1.epsilon,wave2.epsilon]
                print"A1:",wave1.complexamplitudes[n]                        
                print"A2*:",numpy.conjugate(wave2.complexamplitudes[n])
            if argret > 0.:                        
                ret+=log(argret)
            
            self.iList.append(argret)                           
        return ret

def dft2d(img, flags):
    if flags == 1:
        return Fdft2d(img)
    elif flags == -1:
        res = np.conjugate(img)
        #return np.conjugate(Fdft2d(img))
        return np.conjugate(Fdft2d(img))

 def myzpk2tf(self, z, p, k):
         z = np.atleast_1d(z)
         k = np.atleast_1d(k)
         if len(z.shape) > 1:
                 temp = np.poly(z[0])
                 b = np.zeros((z.shape[0], z.shape[1] + 1), temp.dtype.char)
                 if len(k) == 1:
                         k = [k[0]] * z.shape[0]
                 for i in range(z.shape[0]):
                         b[i] = k[i] * poly(z[i])
         else:
                 b = k * np.poly(z)
         a = np.atleast_1d(np.poly(p))
         # Use real output if possible. Copied from numpy.poly, since
         # we can't depend on a specific version of numpy.
         if issubclass(b.dtype.type, np.complexfloating):
                 # if complex roots are all complex conjugates, the roots are real.
                 roots = np.asarray(z, complex)
                 pos_roots = np.compress(roots.imag > 0, roots)
                 neg_roots = np.conjugate(np.compress(roots.imag < 0, roots))
                 if len(pos_roots) == len(neg_roots):
                         if np.all(np.sort_complex(neg_roots) == np.sort_complex(pos_roots)):
                                 b = b.real.copy()
         if issubclass(a.dtype.type, np.complexfloating):
                 # if complex roots are all complex conjugates, the roots are real.
                 roots = np.asarray(p, complex)
                 pos_roots = np.compress(roots.imag > 0, roots)
                 neg_roots = np.conjugate(np.compress(roots.imag < 0, roots))
                 if len(pos_roots) == len(neg_roots):
                         if np.all(np.sort_complex(neg_roots) == np.sort_complex(pos_roots)):
                                 a = a.real.copy()
         return b, a

def MakeEpsilonScreenFromList(EpsilonList, N):
    epsilon = np.zeros((N,N),dtype=np.complex)
    #There are (N^2-1)/2 real elements followed by (N^2-1)/2 complex elements

    #The first (N-1)/2 are the top row
    N_re = (N*N-1)/2
    i = 0
    for x in range(1,(N+1)/2):
        epsilon[0][x] = EpsilonList[i] + 1j * EpsilonList[i+N_re]
        epsilon[0][N-x] = np.conjugate(epsilon[0][x])
        i=i+1

    #The next N(N-1)/2 are filling the next N rows
    for y in range(1,(N+1)/2):
        for x in range(N):
            epsilon[y][x] = EpsilonList[i] + 1j * EpsilonList[i+N_re]

            x2 = N - x
            y2 = N - y
            if x2 == N:
                x2 = 0
            if y2 == N:
                y2 = 0

            epsilon[y2][x2] = np.conjugate(epsilon[y][x])
            i=i+1    

    if no_linear_shift == True:
        epsilon[0,0] = 0
        epsilon[1,0] = 0
        epsilon[0,1] = 0
        epsilon[-1,0] = 0
        epsilon[0,-1] = 0

    return epsilon

def xcorr(x,y,**kwargs):
	"""cross correlation by rfft"""
	x = np.asarray(x)
	y = np.asarray(y)
	if np.ndim(x) == np.ndim(y):
		shape=kwargs.get('shape',np.max((x.shape, y.shape), axis = 0))
		return np.fft.irfftn(np.conjugate(np.fft.rfftn(x,s=shape))*np.fft.rfftn(y,s=shape))
	elif np.ndim(y) == 1:
		axis = kwargs.get('axis', 0)
		shape=kwargs.get('shape', max(x.shape[axis], len(y)))
		shape+=shape%2
		outshape = np.array(x.shape[:])
		outshape[axis] = shape
		out = np.zeros(outshape)
		y = np.fft.ifftshift(np.pad(y, pad_width = (int((shape-len(y)+1)/2), int((shape-len(y))/2)), mode = 'constant'))
		y_fft = np.fft.rfft(y, n=shape)
		x_fft = np.fft.rfft(x, n=shape, axis=axis)
		if axis == 0:
			for ii in range(len(x_fft[0])):
				out[:,ii] = np.fft.irfft(x_fft[:,ii]*np.conjugate(y_fft))
		else:
			for ii in range(len(x_fft)):
				out[ii] = np.fft.irfft(x_fft[ii]*np.conjugate(y_fft))
		return out
	else:
		raise ValueError('Only inputs with dimensions of 1 or 2 can be processed.')

def MakeEpsilonScreen(Nx, Ny, rngseed = 0):
    if rngseed != 0:
        np.random.seed( rngseed )

    epsilon = np.random.normal(loc=0.0, scale=1.0/math.sqrt(2), size=(Nx,Ny)) + 1j * np.random.normal(loc=0.0, scale=1.0/math.sqrt(2), size=(Nx,Ny))
    epsilon[0][0] = 0.0

    #Now let's ensure that it has the necessary conjugation symmetry
    for x in range(Nx):
        if x > (Nx-1)/2:
            epsilon[0][x] = np.conjugate(epsilon[0][Nx-x])
        for y in range((Ny-1)/2, Ny):
            x2 = Nx - x
            y2 = Ny - y
            if x2 == Nx:
                x2 = 0
            if y2 == Ny:
                y2 = 0
            epsilon[y][x] = np.conjugate(epsilon[y2][x2])

    if no_linear_shift == True:
        epsilon[0,0] = 0
        epsilon[1,0] = 0
        epsilon[0,1] = 0
        epsilon[-1,0] = 0
        epsilon[0,-1] = 0

    return epsilon

    def perform_quadrature(self, row, col):
        r"""Evaluates the integral :math:`\langle \Phi_i | \Phi^\prime_j \rangle`
        by an exact symbolic formula.

        :param row: The index :math:`i` of the component :math:`\Phi_i` of :math:`\Psi`.
        :param row: The index :math:`j` of the component :math:`\Phi^\prime_j` of :math:`\Psi^\prime`.
        :return: A single complex floating point number.
        """
        eps = self._packet.get_eps()

        Pibra = self._pacbra.get_parameters(component=row)
        Piket = self._packet.get_parameters(component=col)
        cbra = self._pacbra.get_coefficient_vector(component=row)
        cket = self._packet.get_coefficient_vector(component=col)
        Kbra = self._pacbra.get_basis_shapes(component=row)
        Kket = self._packet.get_basis_shapes(component=col)

        self._cache_factors(Pibra[:4], Piket[:4], Kbra, Kket, eps)

        result = array([[0.0j]], dtype=complexfloating)

        for r in Kbra:
            for c in Kket:
                cr = cbra[Kbra[r],0]
                cc = cket[Kket[c],0]
                i = self.exact_result_higher(Pibra[:4], Piket[:4], eps, r[0], c[0])
                result = result + conjugate(cr) * cc * i

        phase = exp(1.0j/eps**2 * (Piket[4]-conjugate(Pibra[4])))
        return phase * result