def test_init(self):
        import numpy as np
        import math
        import sys

        assert np.intp() == np.intp(0)
        assert np.intp("123") == np.intp(123)
        raises(TypeError, np.intp, None)
        assert np.float64() == np.float64(0)
        assert math.isnan(np.float64(None))
        assert np.bool_() == np.bool_(False)
        assert np.bool_("abc") == np.bool_(True)
        assert np.bool_(None) == np.bool_(False)
        assert np.complex_() == np.complex_(0)
        # raises(TypeError, np.complex_, '1+2j')
        assert math.isnan(np.complex_(None))
        for c in ["i", "I", "l", "L", "q", "Q"]:
            assert np.dtype(c).type().dtype.char == c
        for c in ["l", "q"]:
            assert np.dtype(c).type(sys.maxint) == sys.maxint
        for c in ["L", "Q"]:
            assert np.dtype(c).type(sys.maxint + 42) == sys.maxint + 42
        assert np.float32(np.array([True, False])).dtype == np.float32
        assert type(np.float32(np.array([True]))) is np.ndarray
        assert type(np.float32(1.0)) is np.float32
        a = np.array([True, False])
        assert np.bool_(a) is a

	def fun(self,t, args):
		"""Maximum likelihood function to be minimized.

		:param numpy_array t: t values.
		:param numpy_array args: first entry contains correlation counts, second the corresponding basis as string.
		:param numpy_array args: PSIs, all possible basis as Jones vectors.

		:return: Function value. See for further information D. F. V. James et al. Phys. Rev. A, 64, 052312 (2001).
		:rtype: numpy float
		"""
		corr_counts=args[0]
		basis=args[1]
		PSI=args[2]
		nbrOfElements=len(basis)

		rho_phys = self.rho_phys(t)

		#Estimate NormFactor
		estNormFactor=[]

		for i in range(len(corr_counts)):
			estNormFactor.append(np.dot(np.dot(np.conj(PSI[i]), rho_phys),PSI[i]))

		NormFactor=np.sum(corr_counts)/np.sum(estNormFactor)

		#Optimize density matrix
		BraRoh_physKet = np.complex_(np.zeros(nbrOfElements))

		for i in range(nbrOfElements):
			rhoket = np.array(np.dot(rho_phys,PSI[i]).flat) #Convert 2d to 1d array with flat
			BraRoh_physKet[i] = np.complex_(np.dot(np.conj(PSI[i]), rhoket))

		return np.real(np.sum((NormFactor*BraRoh_physKet-corr_counts)**2.0/(2.0*NormFactor*BraRoh_physKet)))

 def test_init(self):
     import numpy as np
     import math
     import sys
     assert np.intp() == np.intp(0)
     assert np.intp('123') == np.intp(123)
     raises(TypeError, np.intp, None)
     assert np.float64() == np.float64(0)
     assert math.isnan(np.float64(None))
     assert np.bool_() == np.bool_(False)
     assert np.bool_('abc') == np.bool_(True)
     assert np.bool_(None) == np.bool_(False)
     assert np.complex_() == np.complex_(0)
     #raises(TypeError, np.complex_, '1+2j')
     assert math.isnan(np.complex_(None))
     for c in ['i', 'I', 'l', 'L', 'q', 'Q']:
         assert np.dtype(c).type().dtype.char == c
     for c in ['l', 'q']:
         assert np.dtype(c).type(sys.maxint) == sys.maxint
     for c in ['L', 'Q']:
         assert np.dtype(c).type(sys.maxint + 42) == sys.maxint + 42
     assert np.float32(np.array([True, False])).dtype == np.float32
     assert type(np.float32(np.array([True]))) is np.ndarray
     assert type(np.float32(1.0)) is np.float32
     a = np.array([True, False])
     assert np.bool_(a) is a

def fft_deriv ( var ):
    #on_error, 2
                 
    n = numpy.size(var)

    F = old_div(numpy.fft.fft(var),n)  #different definition between IDL - python
                 
    imag = numpy.complex(0.0, 1.0)
    imag = numpy.complex_(imag)

    F[0] = 0.0
          
                
    if (n % 2) == 0 :
      # even number
        for i in range (1, old_div(n,2)) :
          a = imag*2.0*numpy.pi*numpy.float(i)/numpy.float(n)
          F[i] = F[i] * a         # positive frequencies
          F[n-i] = - F[n-i] * a   # negative frequencies
           
        F[old_div(n,2)] = F[old_div(n,2)] * (imag*numpy.pi)
    else:
      # odd number
        for i in range (1, old_div((n-1),2)+1) :
          a = imag*2.0*numpy.pi*numpy.float(i)/numpy.float(n)
          F[i] = F[i] * a
          F[n-i] = - F[n-i] * a 
        
    

    result = numpy.fft.ifft(F)*n  #different definition between IDL - python
    
      
    return result 

	def test_gamma(self,gamma):
		"""Test if :math:`\Gamma_i, {i=1,...,16}` matrices are properly defined.

		:param array GAMMA: Gamma matrices.

		Test for:

			:math:`\mathrm{Tr}(\Gamma_i\Gamma_j)= \delta_{i,j}`

		:return: True if equation fullfilled for all gamma matrices, False otherwise.
		:rtype: bool
		"""

		#initialize empty test matrix
		test_matrix = np.complex_(np.zeros((16,16)))

		for i in range(len(gamma)):
			for j in range(len(gamma)):
				test_matrix[i,j] = np.trace(np.dot(gamma[i], gamma[j]))

		diag_matrix = np.diag(np.ones(16))


		test_result = np.einsum('ij,ij',test_matrix - diag_matrix, test_matrix - diag_matrix)-16

		if np.abs(test_result) < 10**-6:
			return True
		else:
			False

 def __new__(cls, x=0):
     if isinstance(x, afnumpy.ndarray):
         return x.astype(cls)
     elif isinstance(x, numbers.Number):
         return numpy.complex_(x)
     else:
         return afnumpy.array(x).astype(cls)

    def test_against_cmath(self):
        import cmath, sys

        # cmath.asinh is broken in some versions of Python, see
        # http://bugs.python.org/issue1381
        broken_cmath_asinh = False
        if sys.version_info < (2, 6):
            broken_cmath_asinh = True

        points = [-1 - 1j, -1 + 1j, +1 - 1j, +1 + 1j]
        name_map = {
            "arcsin": "asin",
            "arccos": "acos",
            "arctan": "atan",
            "arcsinh": "asinh",
            "arccosh": "acosh",
            "arctanh": "atanh",
        }
        atol = 4 * np.finfo(np.complex).eps
        for func in self.funcs:
            fname = func.__name__.split(".")[-1]
            cname = name_map.get(fname, fname)
            try:
                cfunc = getattr(cmath, cname)
            except AttributeError:
                continue
            for p in points:
                a = complex(func(np.complex_(p)))
                b = cfunc(p)

                if cname == "asinh" and broken_cmath_asinh:
                    continue

                assert_(abs(a - b) < atol, "%s %s: %s; cmath: %s" % (fname, p, a, b))

def halo_transform(image):

    na=image.shape[0]
    nb=image.shape[1]
    ko= 5.336
    delta = (2.*np.sqrt(-2.*np.log(.75)))/ko

    x=np.arange(nb)
    y=np.arange(na)
    x,y=np.meshgrid(x,y)
    if (nb % 2) == 0:
        x=(1.*x - (nb)/2.)/nb
        shiftx = (nb)/2.
    else:
       x= (1.*x - (nb-1)/2.)/nb 
       shiftx=(nb-1.)/2.+1

    if (na % 2) == 0:
        y=(1.*y-(na/2.))/na
        shifty=(na)/2.
    else:
        y= (1.*y - (na-1)/2.)/ na
        shifty=(na-1.)/2.+1

    #--------------Spectral Logarithm--------------------
    M=int(np.log(nb)/delta)
    a2= np.zeros(M)
    a2[0]=np.log(nb)

    for i in range(M-1):
        a2[i+1]=a2[i]-delta

    a2=np.exp(a2)
    tab_k = 1. / (a2)
    wt= np.complex_(np.zeros(((na,nb,M))))


    a= ko*a2

    imageFT= np.fft.fft2(image)
    #imageFT=np.fft.fftshift(image)
    imageFT= np.roll(imageFT,int( shiftx), axis=1)
    imageFT= np.roll(imageFT,int(shifty), axis=0)

    for j in range(M):
        uv = 0

        uv=np.exp(-0.5*((abs(a[j]*np.sqrt(x**2.+y**2.))**2.  - abs(ko))**2.))

        uv = uv * a[j]

        W1FT=imageFT*(uv)
        W1F2=np.roll(W1FT,int(shiftx), axis =1)
        W1F2=np.roll(W1F2,int(shifty),axis=0)
       #W1F2=np.fft.ifftshift(W1FT)
        W1=np.fft.ifft2(W1F2)
        wt[:,:,j]=wt[:,:,j]+ W1

    return wt, tab_k

 def test_params_to_array_inconsistent_types(self):
     """
     Tests if an assertion error is raised when parameters of different
     types are passed in
     """
     guess_params_adj = self.guess_params
     guess_params_adj[-1] = np.complex_(guess_params_adj[-1])
     self.assertRaises(AssertionError, self.affine_obj.params_to_array,
                       guess_params_adj)

 def test_opt_gen_pred_coef_complex(self):
     """
     Tests if complex values are return properly
     """
     # DOC: Need to make sure that the arrays are also np.complex
     guess_params = [np.complex_(el) for el in self.guess_params]
     params = self.affine_obj.params_to_array(guess_params)
     arrays_gpc = self.affine_obj.opt_gen_pred_coef(*params)
     for array in arrays_gpc:
         self.assertEqual(array.dtype, np.complex_)

    def imdct(self,y,i):
        """imdct: inverse mdct
            y(np.array) -> mdct signal
            i -> index frame size
        """
        # frame size
        L = self.sizes[i]
        # signal size
        N = y.size
        # Number of frequency channels
        K = L/2
        # Number of frames
        P = N/K

        # Reshape y
        tmp = y.reshape(P,K)
        y = np.zeros( (P,L) )
        y[:,:K] = tmp

        # Pre-twidle
        y = np.complex_(y)
        y *= np.tile(np.exp(1j*2.*np.pi*np.arange(L)*(L/2+1)/2./L), (P,1))

        # IFFT
        x = np.fft.ifft(y);

        # Post-twidle
        x *= np.tile(np.exp((1./2.)*1j*2*np.pi*(np.arange(L)+((L/2+1)/2.))/L),(P,1));

        # Windowing
        win = self.window(np.arange(L))
        winL = np.copy(win)
        winL[:L/4] = 0
        winL[L/4:L/2] = 1
        winR = np.copy(win)
        winR[L/2:3*L/4] = 1
        winR[3*L/4:] = 0
        x[0,:] *= winL
        x[1:-1,:] *= np.tile(win,(P-2,1))
        x[-1,:] *= winR

        # Real part & scaling
        x = np.sqrt(2./K)*L*x.real

        # Overlap and add
        b = np.repeat(np.arange(P),L)
        a = np.tile(np.arange(L),P) + b*K
        x = sparse.coo_matrix((x.ravel(),(b,a)),shape=(P,N+K)).sum(axis=0).A1

        # Cut edges
        return x[K/2:-K/2]

	def construct_b_matrix(self, PSI, GAMMA):
		"""Construct B matrix as in D. F. V. James et al. Phys. Rev. A, 64, 052312 (2001).

		:param array PSI: :math:`\psi_\\nu` vector with :math:`\\nu=1,...,16`, computed in __init__
		:param array GAMMA: :math:`\Gamma` matrices, computed in __init__

		:return: :math:`B_{\\nu,\mu} = \\langle \psi_\\nu \\rvert  \Gamma_\mu  \\lvert \psi_\\nu \\rangle`
		:rtype: numpy array
		"""

		B = np.complex_(np.zeros((16,16)))

		for i in range(16):
			for j in range(16):
				B[i,j] = np.dot(np.conj(PSI[i]) , np.dot(GAMMA[j], PSI[i]))
		return B

	def from_csv(cls, csv, delimiter="\t"):
		'''Constructs class instance from the given csv table. Is used in testing mode. Method is able to construct only 2d space. Is supposed to read tables produced by 'to_csv' method
		
			csv string: path to csv table
			delimiter string: csv table delimiter
			
			Return Grid: class instance
		'''		
		def _converter(l):
			a = l.strip().split(delimiter)
			a = [x.replace("|", "+") + "j" for x in a]
			return delimiter.join(a);
			
		array = np.genfromtxt((_converter(x) for x in open(csv)), dtype=str, delimiter=delimiter)
		array = np.complex_(array)
		encoding_table = [dict([(x,x) for x in range(array.shape[0])]), dict([(x,x) for x in range(array.shape[1])])]
		return cls(array, encoding_table);

def besselj(n, z):
    """Bessel function of first kind of order n at kr.
    Wraps scipy.special.jn(n, z).

    Parameters
    ----------
    n : array_like
       Order
    z: array_like
       Argument

    Returns
    -------
    J : array_like
       Values of Bessel function of order n at position z
    """
    return scy.jv(n, _np.complex_(z))

    def test_against_cmath(self):
        import cmath, sys

        points = [-1-1j, -1+1j, +1-1j, +1+1j]
        name_map = {'arcsin': 'asin', 'arccos': 'acos', 'arctan': 'atan',
                    'arcsinh': 'asinh', 'arccosh': 'acosh', 'arctanh': 'atanh'}
        atol = 4*np.finfo(np.complex).eps
        for func in self.funcs:
            fname = func.__name__.split('.')[-1]
            cname = name_map.get(fname, fname)
            try:
                cfunc = getattr(cmath, cname)
            except AttributeError:
                continue
            for p in points:
                a = complex(func(np.complex_(p)))
                b = cfunc(p)
                assert_(abs(a - b) < atol, "%s %s: %s; cmath: %s"%(fname, p, a, b))

 def test_generic_roundtrip(self):
     values = [
         np.int_(1),
         np.int32(-2),
         np.float_(2.5),
         np.nan,
         -np.inf,
         np.inf,
         np.datetime64('2014-01-01'),
         np.str_('foo'),
         np.unicode_('bar'),
         np.object_({'a': 'b'}),
         np.complex_(1 - 2j)
     ]
     for value in values:
         decoded = self.roundtrip(value)
         assert_equal(decoded, value)
         self.assertTrue(isinstance(decoded, type(value)))

 def test_dataframe_roundtrip(self):
     if self.should_skip:
         return self.skip('pandas is not importable')
     df = pd.DataFrame({
         'an_int': np.int_([1, 2, 3]),
         'a_float': np.float_([2.5, 3.5, 4.5]),
         'a_nan': np.array([np.nan] * 3),
         'a_minus_inf': np.array([-np.inf] * 3),
         'an_inf': np.array([np.inf] * 3),
         'a_str': np.str_('foo'),
         'a_unicode': np.unicode_('bar'),
         'date': np.array([np.datetime64('2014-01-01')] * 3),
         'complex': np.complex_([1 - 2j, 2 - 1.2j, 3 - 1.3j]),
         # TODO: the following dtypes are not currently supported.
         # 'object': np.object_([{'a': 'b'}]*3),
     })
     decoded_df = self.roundtrip(df)
     assert_frame_equal(decoded_df, df)

 def test_series_roundtrip(self):
     if self.should_skip:
         return self.skip('pandas is not importable')
     ser = pd.Series({
         'an_int': np.int_(1),
         'a_float': np.float_(2.5),
         'a_nan': np.nan,
         'a_minus_inf': -np.inf,
         'an_inf': np.inf,
         'a_str': np.str_('foo'),
         'a_unicode': np.unicode_('bar'),
         'date': np.datetime64('2014-01-01'),
         'complex': np.complex_(1 - 2j),
         # TODO: the following dtypes are not currently supported.
         # 'object': np.object_({'a': 'b'}),
     })
     decoded_ser = self.roundtrip(ser)
     assert_series_equal(decoded_ser, ser)

def frequencies(mol, H):
  
  mol = mol.copy()
  mol.to_angstrom()

  # 3. build mass-weighted Hessian matrix
  natom = mol.natom
  m     = np.repeat(mol.masses, 3) ** -0.5
  M     = np.broadcast_to(m, (3*natom, 3*natom)).T
  tH    = M.T * H * M

  # 4. compute eigenvalues and eigenvectors of the mass-weighted Hessian matrix
  k, tQ  = np.linalg.eigh(tH)

  # 5. un-mass-weight the normal coordinates
  Q = M * tQ

  # 6. get wavenumbers from a.u. force constants
  hartree2J = 4.3597443e-18
  amu2kg = 1.6605389e-27
  bohr2m = 5.2917721e-11
  c = 29979245800.0 # speed of light in cm/s
  v = np.sqrt(np.complex_(k)) * np.sqrt(hartree2J/(amu2kg*bohr2m*bohr2m))/(c*2*np.pi)

  ret = ''
  fmt = '{:2s}' + '{: >15.10f}'*6 + '\n'
  for A in range(3*natom):
    if np.isreal(v[A]):
      ret += '{:d}\n{: >7.2f}  cm^-1\n'.format(natom, np.absolute(v[A]))
    else:
      ret += '{:d}\n{: >7.2f}i cm^-1\n'.format(natom, np.absolute(v[A]))
    for B in range(natom):
      label      = mol.labels[B]
      x ,  y,  z = mol.geom[B]
      dx, dy, dz = Q[3*B:3*B+3, A]
      ret += fmt.format(label, x, y, z, dx, dy, dz)
    ret += '\n'

  # 7. save the normal modes to a file
  outfile = open('modes.xyz', 'w+')
  outfile.write(ret)
  outfile.close()

  return v, Q

	def rho_phys(self, t):
		"""Positive semidefinite matrix based on t values.

		:param numpy_array t: tvalues

		:return: A positive semidefinite matrix which is an estimation of the actual density matrix.
		:rtype: numpy matrix
		"""
		T = np.complex_(np.matrix([
						[t[0], 				0, 				0, 				0],
						[t[4]+1j*t[5], 		t[1], 			0, 				0],
						[t[10]+1j*t[11], 	t[6]+1j*t[7], 	t[2], 			0],
						[t[14]+1j*t[15], 	t[12]+1j*t[13], t[8]+1j*t[9], 	t[3]]
						]))

		TdagT = np.dot(T.conj().T , T)
		norm = np.trace(TdagT)

		return TdagT/norm