def order1(self, ps, vm_in):
        """
        Optical propagation of O(1) periodic input
        
        Here we take the input as an phasor in the optical domain
        and the output is returned as a phasor in the optical domain
        
        vm_out = [a₀, a₁, a₂ ... aᵥ, a₋ᵥ ... a₋₂, a₋₁, a₋₀]

        where the optical signal at optical carrier frequency wc is:
        E(t) = sum_{k=-v}^v a_k exp(j k w0 t) exp(j wc t)

        using the positive frequency convention (as in Rubiola)
        """

        from scipy.integrate import ode

  
        fiber = self.fibers[0]
        alpha = np.complex128(fiber.alpha(self.wl))
        gamma = np.complex128(fiber.gamma(self.wl))
        betaw = fiber.beta(self.wl, ps.omega).astype(np.complex128)

        r = ode(df_nls_order1)
        r.set_integrator('zvode', nsteps=1000, atol=1e-10, rtol=1e-10)
        r.set_initial_value(vm_in, 0)
        r.set_f_params(alpha,gamma,betaw)

        #Solve
        r.integrate(fiber.length)

        return ps, r.y

    def order1_alt(self, ps, v0_in):
        fiber = self.fibers[0]
        alpha = np.complex128(fiber.alpha(self.wl))
        gamma = np.complex128(fiber.gamma(self.wl))
        betaw = fiber.beta(self.wl, ps.omega).astype(np.complex128)

        # r = odesolve.ode_rk45(df_nls_order1)
        r = integrate.ode(df_nls_order1)
        r.set_integrator("zvode", method="adams", nsteps=5000, atol=1e-12, rtol=1e-12)
        r.set_initial_value(v0_in, 0)
        r.set_f_params(alpha, gamma, betaw)

        # Solve
        Nz = 100
        dz = fiber.length / (Nz - 1.0)
        self.O1_z = np.arange(0, Nz) * dz
        self.O1_y = np.zeros(v0_in.shape + (Nz,), np.complex_)

        # Solve to each step in the given range
        self.O1_y[..., 0] = v0_in
        for ii in xrange(1, Nz):
            y_out = r.integrate(self.O1_z[ii])
            self.O1_y[..., ii] = y_out

        # Flag O(1) run as done
        self.O1_run_complete = True
        return ps, y_out

    def eta_fil(self, x, V_app, apprx=(0, 0, 0, 0)):
        m_eff = self.m_r * const.electron_mass

        mpmath.mp.dps = 20
        x0 = Symbol('x0')  # eta_fil
        x1 = Symbol('x1')  # eta_ac
        x2 = Symbol('x2')  # eta_hop
        x3 = Symbol('x3')  # V_tunnel

        f0 = const.Boltzmann * self.T / (1 - self.alpha) / const.elementary_charge / self.z * \
             ln(self.A_fil/self.A_ac*(exp(- self.alpha * const.elementary_charge * self.z / const.Boltzmann / self.T * x0) - 1) + 1) - x1# eta_ac = f(eta_fil) x1 = f(x0)
        f1 = x*2*const.Boltzmann*self.T/self.a/self.z/const.elementary_charge*\
             asinh(self.j_0et/self.j_0hop*(exp(- self.alpha * const.elementary_charge * self.z / const.Boltzmann / self.T * x0) - 1)) - x2# eta_hop = f(eta_fil)
        f2 = x1 - x0 + x2 - x3

        f3 = -V_app + ((self.C * 3 * sqrt(2 * m_eff * ((4+x3/2)*const.elementary_charge)) / 2 / x * (const.elementary_charge / const.Planck)**2 * \
             exp(- 4 * const.pi * x / const.Planck * sqrt(2 * m_eff * ((4+x3/2)*const.elementary_charge))) * self.A_fil*x3)
                       + (self.j_0et*self.A_fil*(exp(-self.alpha*const.elementary_charge*self.z*x0/const.Boltzmann/self.T) - 1))) * (self.R_el + self.R_S + self.rho_fil*(self.L - x) / self.A_fil) \
             + x3

        eta_fil, eta_ac, eta_hop, V_tunnel = nsolve((f0, f1, f2, f3), [x0, x1, x2, x3], apprx)
        eta_fil = np.real(np.complex128(eta_fil))
        eta_ac = np.real(np.complex128(eta_ac))
        eta_hop = np.real(np.complex128(eta_hop))
        V_tunnel = np.real(np.complex128(V_tunnel))
        current = ((self.C * 3 * sqrt(2 * m_eff * ((4+V_tunnel)*const.elementary_charge)) / 2 / x * (const.elementary_charge / const.Planck)**2 * \
            exp(- 4 * const.pi * x / const.Planck * sqrt(2 * m_eff * ((4+V_tunnel)*const.elementary_charge))) * self.A_fil*V_tunnel)
                       + (self.j_0et*self.A_fil*(exp(-self.alpha*const.elementary_charge*self.z*eta_fil/const.Boltzmann/self.T) - 1)))
        print(eta_fil, eta_ac, eta_hop, V_tunnel)
        # print(eta_ac - eta_fil + eta_hop - V_tunnel)
        return eta_fil, eta_ac, eta_hop, V_tunnel, current

    def test_cublasZgemmBatched(self):
        l, m, k, n = 11, 7, 5, 3
        A = (np.random.rand(l, m, k)+1j*np.random.rand(l, m, k)).astype(np.complex128)
        B = (np.random.rand(l, k, n)+1j*np.random.rand(l, k, n)).astype(np.complex128)

        C_res = np.einsum('nij,njk->nik', A, B)

        a_gpu = gpuarray.to_gpu(A)
        b_gpu = gpuarray.to_gpu(B)
        c_gpu = gpuarray.empty((l, m, n), np.complex128)

        alpha = np.complex128(1.0)
        beta = np.complex128(0.0)

        a_arr = bptrs(a_gpu)
        b_arr = bptrs(b_gpu)
        c_arr = bptrs(c_gpu)

        cublas.cublasZgemmBatched(self.cublas_handle, 'n','n',
                                  n, m, k, alpha,
                                  b_arr.gpudata, n,
                                  a_arr.gpudata, k,
                                  beta, c_arr.gpudata, n, l)

        assert np.allclose(C_res, c_gpu.get())

    def test_numpy(self):
        assert chash(np.bool_(True)) == chash(np.bool_(True))

        assert chash(np.int8(1)) == chash(np.int8(1))
        assert chash(np.int16(1))
        assert chash(np.int32(1))
        assert chash(np.int64(1))

        assert chash(np.uint8(1))
        assert chash(np.uint16(1))
        assert chash(np.uint32(1))
        assert chash(np.uint64(1))

        assert chash(np.float32(1)) == chash(np.float32(1))
        assert chash(np.float64(1)) == chash(np.float64(1))
        assert chash(np.float128(1)) == chash(np.float128(1))

        assert chash(np.complex64(1+1j)) == chash(np.complex64(1+1j))
        assert chash(np.complex128(1+1j)) == chash(np.complex128(1+1j))
        assert chash(np.complex256(1+1j)) == chash(np.complex256(1+1j))

        assert chash(np.datetime64('2000-01-01')) == chash(np.datetime64('2000-01-01'))
        assert chash(np.timedelta64(1,'W')) == chash(np.timedelta64(1,'W'))

        self.assertRaises(ValueError, chash, np.object())

        assert chash(np.array([[1, 2], [3, 4]])) == \
            chash(np.array([[1, 2], [3, 4]]))
        assert chash(np.array([[1, 2], [3, 4]])) != \
            chash(np.array([[1, 2], [3, 4]]).T)
        assert chash(np.array([1, 2, 3])) == chash(np.array([1, 2, 3]))
        assert chash(np.array([1, 2, 3], dtype=np.int32)) != \
            chash(np.array([1, 2, 3], dtype=np.int64))

def distance_matrix(K, N, dim, obs, sampl, method):
	
	dist = np.zeros((K, N), dtype=np.complex128)
	
	if method == "eucl":
		for i in range(0, N): 
		
			for j in range(0, K): 
			
			  zero = np.complex128(0)
			  
			  for h in range(0, dim): 
				
				zero += (obs[j, h] - sampl[i, h])**2  
			  
			  dist[j,i] = ma.sqrt(zero)
			  
	if method == "city":
		for i in range(0, N): 
		
			for j in range(0, K): 
			
			  zero = np.complex128(0)
			  
			  for h in range(0, dim): 
				
				zero += abs(obs[j, h] - sampl[i, h])  
			  
			  dist[j,i] = ma.sqrt(zero)
	
	
	
	
	
	return dist 

    def ordere(self, ms, v1_in):
        """
        Modulator response to O(ε) slow-time input
        linearized about an O(1) behaviour
        """
        # Check for a prior O(1) solution
        assert self.O1_complete, "O(1) solve must be performed before calling O(ε) solve"

        # Fiber parameters
        fiber = self.fibers[0]
        alpha = np.complex128(fiber.alpha(self.wl))
        gamma = np.complex128(fiber.gamma(self.wl))
        betaw = fiber.beta(self.wl, ms.omega).astype(np.complex128)

        # Solver only deals with flat arrays .. how sad
        y_in = ms.fft(v1_in).flatten()

        # Setup solver
        # r = odesolve.ode_rk45(df_nls_ordere, atol=1e-10, rtol=1e-10)
        r = ode(df_nls_ordere)
        r.set_integrator("zvode", method="adams", nsteps=5000, atol=1e-8, rtol=1e-8)
        r.set_initial_value(y_in, 0)
        r.set_f_params(self.O1_y, self.O1_z, alpha, gamma, betaw, v1_in.shape[0], v1_in.shape[1])

        # Solve to the end - no need for intermediate solutions
        v1_out = r.integrate(fiber.length)
        v1_out = ms.ifft(v1_out.reshape(v1_in.shape))

        return ms, v1_out

 def test_dict_numpy_complex(self):
     x = {'foo': np.complex128(1.0 + 1.0j),
          'bar': np.complex128(2.0 + 2.0j)}
     x_rec = self.encode_decode(x)
     self.assertEqual(x, x_rec)
     for key in x:
         self.assertEqual(type(x[key]), type(x_rec[key]))

 def test_dict_numpy_complex(self):
     x = {"foo": np.complex128(1.0 + 1.0j), "bar": np.complex128(2.0 + 2.0j)}
     x_rec = self.encode_decode(x)
     self.assert_(
         all(map(lambda x, y: x == y, x.values(), x_rec.values()))
         and all(map(lambda x, y: type(x) == type(y), x.values(), x_rec.values()))
     )

    def test_dict_numpy_complex(self):
        x = {'foo': np.complex128(1.0 + 1.0j),
             'bar': np.complex128(2.0 + 2.0j)}
        x_rec = self.encode_decode(x)
        tm.assert_dict_equal(x, x_rec)

        for key in x:
            tm.assert_class_equal(x[key], x_rec[key], obj="numpy complex128")

def PV(fns,n,p,cols,rows,bands):
    '''Return p-values for change indices R^ell_j'''
    j = np.float64(len(fns))
    eps = sys.float_info.min
    k = 0.0; a = 0.0; rho = 0.0; xsi = 0.0; b = 0.0; zeta = 0.0
    for fn in fns:
        result = getmat(fn,cols,rows,bands)
        if p==3:
            k1,a1,rho1,xsi1,b1,zeta1 = result
            k1 = n*np.float64(k1)
            a1 = n*np.complex128(a1)
            rho1 = n*np.complex128(rho1)
            xsi1 = n*np.float64(xsi1)
            b1 = n*np.complex128(b1)
            zeta1 = n*np.float64(zeta1)
            k += k1; a += a1; rho += rho1; xsi += xsi1; b += b1; zeta += zeta1  
        elif p==2:
            k1,a1,xsi1 = result
            k1 = n*np.float64(k1)
            a1 = n*np.complex128(a1)
            xsi1 = n*np.float64(xsi1)
            k += k1; a += a1; xsi += xsi1
        elif p==1:
            k1 = n*np.float64(result[0])
            k += k1              
    if p==3: 
        detsumj = k*xsi*zeta + 2*np.real(a*b*np.conj(rho)) - xsi*(abs(rho)**2) - k*(abs(b)**2) - zeta*(abs(a)**2) 
        k -= k1; a -= a1; rho -= rho1; xsi -= xsi1; b -= b1; zeta -= zeta1 
        detsumj1 = k*xsi*zeta + 2*np.real(a*b*np.conj(rho)) - xsi*(abs(rho)**2) - k*(abs(b)**2) - zeta*(abs(a)**2)
        detj = k1*xsi1*zeta1 + 2*np.real(a1*b1*np.conj(rho1)) - xsi1*(abs(rho1)**2) - k1*(abs(b1)**2) - zeta1*(abs(a1)**2)
    elif p==2:
        detsumj = k*xsi - abs(a)**2 
        k -= k1; a -= a1; xsi -= xsi1
        detsumj1 = k*xsi - abs(a)**2
        detj = k1*xsi1 - abs(a1)**2
    elif p==1:
        detsumj = k  
        k -= k1
        detsumj1 = k
        detj = k1    
    detsumj = np.nan_to_num(detsumj)    
    detsumj = np.where(detsumj <= eps,eps,detsumj)
    logdetsumj = np.log(detsumj)
    detsumj1 = np.nan_to_num(detsumj1)
    detsumj1 = np.where(detsumj1 <= eps,eps,detsumj1)
    logdetsumj1 = np.log(detsumj1)
    detj = np.nan_to_num(detj)
    detj = np.where(detj <= eps,eps,detj)
    logdetj = np.log(detj)
#  test statistic    
    lnRj = n*( p*( j*np.log(j)-(j-1)*np.log(j-1.) ) + (j-1)*logdetsumj1 + logdetj - j*logdetsumj )   
    f =p**2
    rhoj = 1 - (2.*p**2 - 1)*(1. + 1./(j*(j-1)))/(6.*p*n)
    omega2j = -(p*p/4.)*(1.-1./rhoj)**2 + (1./(24.*n*n))*p*p*(p*p-1)*(1+(2.*j-1)/(j*(j-1))**2)/rhoj**2   
#  return p-values  
    Z = -2*rhoj*lnRj
    return 1.0 - ((1.-omega2j)*stats.chi2.cdf(Z,[f])+omega2j*stats.chi2.cdf(Z,[f+4]))

def add_vdot(M, v, out, beta=0.0, transM='N', handle=None):
    if handle is None:
        handle = scm._global_cublas_handle

    assert M.strides[1] <= M.strides[0], 'only C-order arrays supported'

    transM = transM.lower()
    if transM == 'n':
        trans = 't'
        m = M.shape[1]
        n = M.shape[0]
        alpha = 1.0
        lda = M.strides[0] // M.dtype.itemsize
        if v.shape[0] != M.shape[1] or out.shape[0] != M.shape[0]:
            raise ValueError('dimension mismatch: %s %s %s' %
                             (M.shape, v.shape, out.shape))
    elif transM == 't':
        trans = 'n'
        m = M.shape[1]
        n = M.shape[0]
        alpha = 1.0
        lda = M.strides[0] // M.dtype.itemsize
        if v.shape[0] != M.shape[0] or out.shape[0] != M.shape[1]:
            raise ValueError('dimension mismatch: %s %s %s' %
                             (M.shape, v.shape, out.shape))
    else:
        raise ValueError('transM must be n or t')

    if (M.dtype == np.complex64 and v.dtype == np.complex64):
        cublas_func = scikits.cuda.cublas.cublasCgemv
        alpha = np.complex64(alpha)
        beta = np.complex64(beta)
    elif (M.dtype == np.float32 and v.dtype == np.float32):
        cublas_func = scikits.cuda.cublas.cublasSgemv
        alpha = np.float32(alpha)
        beta = np.float32(beta)
    elif (M.dtype == np.complex128 and v.dtype == np.complex128):
        cublas_func = scikits.cuda.cublas.cublasZgemv
        alpha = np.complex128(alpha)
        beta = np.complex128(beta)
    elif (M.dtype == np.float64 and v.dtype == np.float64):
        cublas_func = scikits.cuda.cublas.cublasDgemv
        alpha = np.float64(alpha)
        beta = np.float64(beta)
    else:
        raise ValueError('unsupported combination of input types')

    incx = 1
    incy = 1
    cublas_func(handle,
                trans, m, n,
                alpha,
                M.gpudata, lda,
                v.gpudata, incx,
                beta,
                out.gpudata, incy)

 def test_dict_numpy_complex(self):
     x = {b'foo': np.complex128(1.0+1.0j), b'bar': np.complex128(2.0+2.0j)}
     x_rec = self.encode_decode(x)
     assert_array_equal(sorted(x.values(), key=np.linalg.norm),
                        sorted(x_rec.values(), key=np.linalg.norm))
     assert_array_equal([type(e) for e in sorted(x.values(), key=np.linalg.norm)],
                        [type(e) for e in sorted(x_rec.values(), key=np.linalg.norm)])
     assert_array_equal(sorted(x.keys()), sorted(x_rec.keys()))
     assert_array_equal([type(e) for e in sorted(x.keys())],
                        [type(e) for e in sorted(x_rec.keys())])

 def test_prod_sum_fn (self):
     compiled_prod_sum_fn = self.jit(argtypes = [complex128, complex128, complex128],
                                     restype = complex128)(prod_sum_fn)
     rng = numpy.arange(-1., 1.1, 0.5)
     for ar, ai, xr, xi, br, bi in itertools.product(rng, rng, rng, rng, rng,
                                                     rng):
         a = numpy.complex128(ar + ai * 1j)
         x = numpy.complex128(xr + xi * 1j)
         b = numpy.complex128(br + bi * 1j)
         self.assertEqual(prod_sum_fn(a, x, b),
                          compiled_prod_sum_fn(a, x, b))

    def test_legacy_mode_scalars(self):
        # in legacy mode, str of floats get truncated, and complex scalars
        # use * for non-finite imaginary part
        np.set_printoptions(legacy='1.13')
        assert_equal(str(np.float64(1.123456789123456789)), '1.12345678912')
        assert_equal(str(np.complex128(complex(1, np.nan))), '(1+nan*j)')

        np.set_printoptions(legacy=False)
        assert_equal(str(np.float64(1.123456789123456789)),
                     '1.1234567891234568')
        assert_equal(str(np.complex128(complex(1, np.nan))), '(1+nanj)')

 def test_cublasZgemv(self):
     a = (np.random.rand(2, 3)+1j*np.random.rand(2, 3)).astype(np.complex128)
     x = (np.random.rand(3, 1)+1j*np.random.rand(3, 1)).astype(np.complex128)
     a_gpu = gpuarray.to_gpu(a.T.copy())
     x_gpu = gpuarray.to_gpu(x)
     y_gpu = gpuarray.empty((2, 1), np.complex128)
     alpha = np.complex128(1.0)
     beta = np.complex128(0.0)
     cublas.cublasZgemv('n', 2, 3, alpha, a_gpu.gpudata, 2, x_gpu.gpudata,
                        1, beta, y_gpu.gpudata, 1)
     assert np.allclose(y_gpu.get(), np.dot(a, x))

 def test_cublasXtCgemm(self):
     a = (np.random.rand(4, 4)+1j*np.random.rand(4, 4)).astype(np.complex128)
     b = (np.random.rand(4, 4)+1j*np.random.rand(4, 4)).astype(np.complex128)
     c = np.zeros((4, 4), np.complex128)
     
     cublasxt.cublasXtCgemm(self.handle,
                            'N', 'N', 4, 4, 4,           
                            np.complex128(1.0),
                            a.ctypes.data, 4, b.ctypes.data, 4, np.complex128(0.0),
                            c.ctypes.data, 4)
     np.allclose(np.dot(b.T, a.T).T, c)

 def test_prod_sum_fn (self):
     compiled_prod_sum_fn = numba_compile(arg_types = ['D', 'D', 'D'],
                                          ret_type = 'D')(prod_sum_fn)
     rng = numpy.arange(-1., 1.1, 0.5)
     for ar, ai, xr, xi, br, bi in itertools.product(rng, rng, rng, rng, rng,
                                                     rng):
         a = numpy.complex128(ar + ai * 1j)
         x = numpy.complex128(xr + xi * 1j)
         b = numpy.complex128(br + bi * 1j)
         self.assertEqual(prod_sum_fn(a, x, b),
                          compiled_prod_sum_fn(a, x, b))

 def test_cublasXtCgemm(self):
     a = (np.random.rand(4, 4)+1j*np.random.rand(4, 4)).astype(np.complex128)
     b = (np.random.rand(4, 4)+1j*np.random.rand(4, 4)).astype(np.complex128)
     c = np.zeros((4, 4), np.complex128)
     
     cublasxt.cublasXtDeviceSelect(handle, 2, np.array([0, 1], np.int32).ctypes.data)
     cublasxt.cublasXtCgemm(self.handle, cublasxt._CUBLAS_OP['N'],
                            cublasxt._CUBLAS_OP['N'], 4, 4, 4, np.complex128(1.0),
                            a.ctypes.data, 4, b.ctypes.data, 4, np.complex128(0.0),
                            c.ctypes.data, 4)
     np.allclose(np.dot(b.T, a.T).T, c)

    def get_wavefunction(self, ibzk, n, check_focc=True, spin=0):
        if (self.calc.wfs.world.size == 1 or self.calc.wfs.gd.comm.size != 1 
            or self.calc.input_parameters['mode'] == 'lcao'):
            if not check_focc:
                return
            else:
                psit_G = self.calc.wfs.get_wave_function_array(n, ibzk, spin)
            
                if self.calc.wfs.world.size == 1:
                    return np.complex128(psit_G)
                
                if self.calc.wfs.world.rank != 0:
                    psit_G = self.calc.wfs.gd.empty(dtype=self.calc.wfs.dtype,
                                                    global_array=True)
                self.calc.wfs.world.broadcast(psit_G, 0)
        
                return np.complex128(psit_G)
        else:
            # support ground state calculation with kpoint and band parallelization
            # but domain decomposition must = 1
            kpt_rank, u = self.calc.wfs.kd.get_rank_and_index(0, ibzk)
            bzkpt_rank = self.kcomm.rank
            band_rank, myn = self.calc.wfs.bd.who_has(n)
            assert self.calc.wfs.gd.comm.size == 1
            world_rank = (kpt_rank * self.calc.wfs.band_comm.size + band_rank)

            # in the following, kpt_rank is assigned to world_rank
            klist = np.array([world_rank, u, bzkpt_rank, myn])
            klist_kcomm = np.zeros((self.kcomm.size, 4), dtype=int)            
            self.kcomm.all_gather(klist, klist_kcomm)

            check_focc_global = np.zeros(self.kcomm.size, dtype=bool)
            self.kcomm.all_gather(np.array([check_focc]), check_focc_global)

            psit_G = self.calc.wfs.gd.empty(dtype=self.calc.wfs.dtype)

            for i in range(self.kcomm.size):
                if check_focc_global[i]:
                    kpt_rank, u, bzkpt_rank, nlocal = klist_kcomm[i]
                    if kpt_rank == bzkpt_rank:
                        if rank == kpt_rank:
                            psit_G = self.calc.wfs.kpt_u[u].psit_nG[nlocal]
                    else:
                        if rank == kpt_rank:
                            world.send(self.calc.wfs.kpt_u[u].psit_nG[nlocal],
                                       bzkpt_rank, 1300+bzkpt_rank)
                        if rank == bzkpt_rank:
                            psit_G = self.calc.wfs.gd.empty(dtype=self.calc.wfs.dtype)
                            world.receive(psit_G, kpt_rank, 1300+bzkpt_rank)
                    
            self.wScomm.broadcast(psit_G, 0)

            return psit_G