def lqmn(m,n,z):
    """Associated Legendre functions of the second kind, Qmn(z) and its
    derivative, ``Qmn'(z)`` of order m and degree n.  Returns two
    arrays of size ``(m+1, n+1)`` containing ``Qmn(z)`` and ``Qmn'(z)`` for
    all orders from ``0..m`` and degrees from ``0..n``.

    z can be complex.
    """
    if not isscalar(m) or (m<0):
        raise ValueError("m must be a non-negative integer.")
    if not isscalar(n) or (n<0):
        raise ValueError("n must be a non-negative integer.")
    if not isscalar(z):
        raise ValueError("z must be scalar.")
    m = int(m)
    n = int(n)

    # Ensure neither m nor n == 0
    mm = max(1,m)
    nn = max(1,n)

    if iscomplex(z):
        q,qd = specfun.clqmn(mm,nn,z)
    else:
        q,qd = specfun.lqmn(mm,nn,z)
    return q[:(m+1),:(n+1)],qd[:(m+1),:(n+1)]

def is_spd(M, decimal=15):
    """Assert that input matrix is real symmetric positive definite.

    M must be symmetric down to specified decimal places and with no complex
    entry.
    The check is performed by checking that all eigenvalues are positive.

    Parameters
    ==========
    M: numpy.ndarray
        matrix.

    Returns
    =======
    answer: boolean
        True if matrix is symmetric real positive definite, False otherwise.
    """
    if not np.allclose(M, M.T, atol=0.1 ** decimal):
        print ("matrix not symmetric to {0} decimals".format(decimal))
        return False
    if np.all(np.iscomplex(M)):
        print ("matrix has a non real value {0}".format(M[np.iscomplex(M)][0]))
    eigvalsh = np.linalg.eigvalsh(M)
    ispd = eigvalsh.min() > 0
    if not ispd:
        print ("matrix has a negative eigenvalue: %.3f" % eigvalsh.min())
    return ispd

def average_structure(X):
    """
    Calculate an average structure from an ensemble of structures
    (i.e. X is a rank-3 tensor: X[i] is a (N,3) configuration matrix).

    @param X: m x n x 3 input vector
    @type X: numpy array
    
    @return: average structure
    @rtype: (n,3) numpy.array
    """
    from numpy.linalg import eigh

    B = csb.numeric.gower_matrix(X)
    v, U = eigh(B)
    if numpy.iscomplex(v).any():
        v = v.real
    if numpy.iscomplex(U).any():
        U = U.real

    indices = numpy.argsort(v)[-3:]
    v = numpy.take(v, indices, 0)
    U = numpy.take(U, indices, 1)
        
    x = U * numpy.sqrt(v)
    i = 0
    while is_mirror_image(x, X[0]) and i < 2:
        x[:, i] *= -1
        i += 1
    return x

    def __init__(self, qpoint, wpts, gsphere, wggmat, inord="C"):
        """"
        Args:
            qpoint: Q-point object
            wpts: Frequency points in Ha.
            wggmat: numpy array of shape [nw, ng, ng]
            inord: storage order of wggmat. If inord == "F", wggmat in 
                in Fortran column-major order. Default: "C" i.e. C row-major order
        """
        self.qpoint = qpoint
        self.wpts = wpts
        self.gsphere = gsphere
        self.wggmat = np.reshape(wggmat, (self.nw, self.ng, self.ng))

        if inord == "F": 
            # Fortran to C.
            for iw in range(len(wpts)):
                self.wggmat[iw] = self.wggmat[iw].T

        for i in (1, 2):
            assert len(gsphere) == wggmat.shape[-i]
        assert len(self.wpts) == len(self.wggmat)

        # Find number of real/imaginary frequencies
        self.nrew = self.nw; self.nimw = 0
        for i, w in enumerate(self.wpts):
            if np.iscomplex(w):
                self.nrew = i
                break

        self.nimw = self.nw - self.nrew
        if self.nimw and not np.all(np.iscomplex(self.wpts[self.nrew+1:])):
            raise ValueError("wpts should contained real points packed in the first positions\n"
                "followed by imaginary points but got: %s" % str(self.wpts))

    def test_random_like(self):
        """
        Test that the random_like function produces sensible data
        """

        # Try for floats and complex data
        for dtype in [np.float32, np.float64, np.complex64, np.complex128]:
            # Test random array creation with same
            # shape and type as existing array
            shape = (np.random.randint(1, 50), np.random.randint(1, 50))
            ary = np.empty(shape=shape, dtype=dtype)    
            random_ary = mbu.random_like(ary)

            # Test that that the shape and type is correct
            self.assertTrue(random_ary.shape == ary.shape)
            self.assertTrue(random_ary.dtype == dtype)

            # Test that we're getting complex data out
            if np.issubdtype(dtype, np.complexfloating):
                proportion_cplx = np.sum(np.iscomplex(random_ary)) / random_ary.size
                self.assertTrue(proportion_cplx > 0.9)

            # Test random array creation with supplied shape and type
            shape = (np.random.randint(1, 50), np.random.randint(1, 50))
            random_ary = mbu.random_like(shape=shape, dtype=dtype)

            # Test that that the shape and type is correct
            self.assertTrue(random_ary.shape == shape)
            self.assertTrue(random_ary.dtype == dtype)

            # Test that we're getting complex data out
            if np.issubdtype(dtype, np.complexfloating):
                proportion_cplx = np.sum(np.iscomplex(random_ary)) / random_ary.size
                self.assertTrue(proportion_cplx > 0.9)

def min(X,Y=[],axis=1):
    axis -= 1
    tX, tY = X, Y
    if _N.iscomplex(X.flat[0]): tX = abs(X)
    if len(tY) > 0:
        if _N.iscomplex(Y.flat[0]): tY = abs(Y)
        return _N.minimum(tX,tY)
    else:
        nargout = _get_nargout()
        print nargout
        if nargout == 1:
            return _N.min(tX,axis)
        elif nargout == 2:
            # slow
            i = _N.argmin(tX,axis)
            return _N.min(tX,axis), i
#             i = _N.argmin(tX,axis)
#             sh = X.shape
#             index = [ slice(0,x,1) for x in sh ]
#             if axis == 0:
#                 index[1] = range(sh[1])
#             else:
#                 index[0] = range(sh[0])
#             index[axis] = i
#             return _N.ndarray.__getslice__(index)
        else:
            raise Exception('too many output vals')

def trapz2d(x_gpu, dx=1.0, dy=1.0, handle=None):
    """
    2D trapezoidal integration.

    Parameters
    ----------
    x_gpu : pycuda.gpuarray.GPUArray
        Input matrix to integrate.
    dx : float
        X-axis spacing.
    dy : float
        Y-axis spacing
    handle : int
        CUBLAS context. If no context is specified, the default handle from
        `skcuda.misc._global_cublas_handle` is used.

    Returns
    -------
    result : float
        Definite double integral as approximated by the trapezoidal rule.

    Examples
    --------
    >>> import pycuda.autoinit
    >>> import pycuda.gpuarray
    >>> import numpy as np
    >>> import integrate
    >>> integrate.init()
    >>> x = np.asarray(np.random.rand(10, 10), np.float32)
    >>> x_gpu = gpuarray.to_gpu(x)
    >>> z = integrate.trapz2d(x_gpu)
    >>> np.allclose(np.trapz(np.trapz(x)), z)
    True
    """

    if handle is None:
        handle = misc._global_cublas_handle

    if len(x_gpu.shape) != 2:
        raise ValueError('input array must be 2D')
    if np.iscomplex(dx) or np.iscomplex(dy):
        raise ValueError('dx and dy must be real')

    float_type = x_gpu.dtype.type
    if float_type == np.complex64:
        cublas_func = cublas.cublasCdotu
    elif float_type == np.float32:
        cublas_func = cublas.cublasSdot
    elif float_type == np.complex128:
        cublas_func = cublas.cublasZdotu
    elif float_type == np.float64:
        cublas_func = cublas.cublasDdot
    else:
        raise ValueError('unsupported input type')

    trapz_mult_gpu = gen_trapz2d_mult(x_gpu.shape, float_type)
    result = cublas_func(handle, x_gpu.size, x_gpu.gpudata, 1,
                         trapz_mult_gpu.gpudata, 1)

    return float_type(dx)*float_type(dy)*result

    def __init__(self, wpoints, gsphere, wggmat, inord="C"):
        """"
        Args:
            gsphere: |GSphere| with G-vectors and k-point object.
            wpoints: Complex frequency points in Hartree.
            wggmat: [nw, ng, ng] complex array.
            inord: storage order of ``wggmat``. If inord == "F", ``wggmat`` is in
                in Fortran column-major order. Default: "C" i.e. C row-major order.
        """
        self.wpoints = np.array(wpoints, dtype=np.complex)
        self.gsphere = gsphere
        self.wggmat = np.reshape(wggmat, (self.nw, self.ng, self.ng))

        if inord.lower() == "f":
            # Fortran to C.
            for iw, _ in enumerate(wpoints):
                self.wggmat[iw] = self.wggmat[iw].T.copy()

        for i in (1, 2):
            assert len(gsphere) == wggmat.shape[-i]
        assert len(self.wpoints) == len(self.wggmat)

        # Find number of real/imaginary frequencies.
        self.nrew = self.nw
        self.nimw = 0
        for i, w in enumerate(self.wpoints):
            if np.iscomplex(w):
                self.nrew = i
                break

        self.nimw = self.nw - self.nrew
        if self.nimw and not np.all(np.iscomplex(self.wpoints[self.nrew+1:])):
            raise ValueError("wpoints should contained real points packed in the first positions\n"
                "followed by imaginary points but got: %s" % str(self.wpoints))

    def fitToData(self, data):
        '''
        param data: numpy array where [:,0] is x and [:,1] is y
        '''
        x = data[:, 0][:, np.newaxis]
        y = data[:, 1][:, np.newaxis]
        D = np.hstack((x*x, x*y, y*y, x, y, np.ones_like(x)))
        S = np.dot(D.T, D)
        C = np.zeros([6, 6])
        C[0, 2] = C[2, 0] = 2; C[1, 1] = -1
        E, V = eig(np.dot(inv(S), C))
        n = np.argmax(np.abs(E))
        self.parameters = V[:, n]

        axes = self.ellipse_axis_length()
        self.a = axes[0]
        self.b = axes[1]
        self.angle = self.ellipse_angle_of_rotation()

        if not self.a or not self.b or self.parameters == None or np.iscomplexobj(self.parameters) or \
           math.isnan(self.a) or math.isnan(self.b) or math.isnan(self.ellipse_center()[0]) or \
           np.iscomplex(self.ellipse_center()[0]) or np.iscomplex(self.a) or np.iscomplex(self.b) or \
           np.iscomplexobj(self.angle):
            self.a = 0
            self.b = 0
            self.parameters = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
            self.angle = 0
            self.error = True

def solution_not_acceptable(P = -5e8, T = 293.15):
    """
    This function raises a flag if the newly calculated values of P or T are
    problematic (either complex, positive or not calculable)
    """
    a = np.any(np.isnan(P)) or np.any(np.iscomplex(P)) or np.any(P>0)
    b = np.any(np.isnan(T)) or np.any(np.iscomplex(T))
    return a or b

 def _check(self, res, ref):
     if hasattr(res, "get_x"):
         x = res.get_x()
         for k in list(res.keys()):
             if np.all(res[k] == x):
                 continue
             elif np.any(np.iscomplex(res[k])) or np.any(np.iscomplex(ref[k])):
                 # Interpolate Re and Im of the results to compare.
                 x = x.reshape((-1,))
                 refx = ref[ref.x].reshape((-1,))
                 d1 = InterpolatedUnivariateSpline(x, np.real(res[k]).reshape((-1,)))
                 d2 = InterpolatedUnivariateSpline(refx, np.real(ref[k]).reshape((-1,)))
                 ok(d1(x), d2(x), rtol=self.er, atol=self.ea, msg=("Test %s FAILED (Re)" % self.test_id))
                 d1 = InterpolatedUnivariateSpline(x, np.imag(res[k]).reshape((-1,)))
                 d2 = InterpolatedUnivariateSpline(refx, np.imag(ref[k]).reshape((-1,)))
                 ok(d1(x), d2(x), rtol=self.er, atol=self.ea, msg=("Test %s FAILED (Im)" % self.test_id))
             else:
                 # Interpolate the results to compare.
                 x = x.reshape((-1,))
                 refx = ref[ref.x].reshape((-1,))
                 d1 = InterpolatedUnivariateSpline(x, np.real_if_close(res[k]).reshape((-1,)))
                 d2 = InterpolatedUnivariateSpline(refx, np.real_if_close(ref[k]).reshape((-1,)))
                 ok(d1(x), d2(x), rtol=self.er, atol=self.ea, msg=("Test %s FAILED" % self.test_id))
     elif isinstance(res, results.op_solution):
         for k in list(res.keys()):
             assert k in ref
             ok(res[k], ref[k], rtol=self.er, atol=self.ea, msg=("Test %s FAILED" % self.test_id))
     elif isinstance(res, results.pz_solution):
         # recover the reference signularities from Re/Im data
         ref_sing_keys = list(ref.keys())[:]
         ref_sing_keys.sort()
         assert len(ref_sing_keys) % 2 == 0
         ref_sing = [
             ref[ref_sing_keys[int(len(ref_sing_keys) / 2) + k]] + ref[ref_sing_keys[k]] * 1j
             for k in range(int(len(ref_sing_keys) / 2))
         ]
         ref_poles_num = len([k for k in ref.keys() if k[:4] == "Re(p"])
         poles_ref, zeros_ref = ref_sing[:ref_poles_num], ref_sing[ref_poles_num:]
         assert len(poles_ref) == len(res.poles)
         pz._check_singularities(res.poles, poles_ref)
         assert len(zeros_ref) == len(res.zeros)
         pz._check_singularities(res.zeros, zeros_ref)
     else:
         if isinstance(res, list) or isinstance(res, tuple):
             for i, j in zip(res, ref):
                 self._check(i, j)
         elif res is not None:
             for k in list(res.keys()):
                 assert k in ref
                 if isinstance(res[k], dict):  # hence ref[k] will be a dict too
                     self._check(res[k], ref[k])
                 elif isinstance(ref[k], sympy.Basic) and isinstance(res[k], sympy.Basic):
                     # get rid of assumptions. Evaluate only expression
                     rf = parse_expr(str(ref[k]))
                     rs = parse_expr(str(res[k]))
                     assert (rs == rf) or (sympy.simplify(rf / rs) == 1)
                 else:
                     assert res[k] == ref[k]

def sph_yn(n, z):
    idx = np.isreal(z)
    out =  _sph_yn_bessel(n, z)
    if np.any(idx):
        # Ascending recurrence is more accurate for real z
        out[idx] = _sph_yn_a_recur(n[idx], z[idx])
    if np.any(np.iscomplex(out)):
        out[np.logical_and(np.isnan(out), np.iscomplex(out))] = np.inf*(1+1j)
    return out

 def test_simple(self):
     a = [[8,12,3],[2,9,3],[10,3,6]]
     t,z = schur(a)
     assert_array_almost_equal(dot(dot(z,t),transp(conj(z))),a)
     tc,zc = schur(a,'complex')
     assert_(any(ravel(iscomplex(zc))) and any(ravel(iscomplex(tc))))
     assert_array_almost_equal(dot(dot(zc,tc),transp(conj(zc))),a)
     tc2,zc2 = rsf2csf(tc,zc)
     assert_array_almost_equal(dot(dot(zc2,tc2),transp(conj(zc2))),a)

def allsortedclose(a, b, atol=1e-3, rtol=1e-3):
    if np.iscomplex(a).any():
        a = np.sort_complex(a)
    else:
        a = np.sort(a)
    if np.iscomplex(b).any():
        b = np.sort_complex(b)
    else:
        b = np.sort(b)
    return np.allclose(a, b, rtol=rtol, atol=atol)

    def _iter_initialize(self):
        """
        Perform any necessary pre-processing operations.

        Returns
        -------
        float
            Initial relative error in the user-specified residuals.
        float
            Initial absolute error in the user-specified residuals.
        """
        system = self._system
        if self.options['debug_print']:
            self._err_cache['inputs'] = self._system._inputs._copy_views()
            self._err_cache['outputs'] = self._system._outputs._copy_views()

        # Convert local storage if we are under complex step.
        if system.under_complex_step:
            if np.iscomplex(self.xm[0]):
                self.Gm = self.Gm.astype(np.complex)
                self.xm = self.xm.astype(np.complex)
                self.fxm = self.fxm.astype(np.complex)
        elif np.iscomplex(self.xm[0]):
            self.Gm = self.Gm.real
            self.xm = self.xm.real
            self.fxm = self.fxm.real

        self._converge_failures = 0
        self._computed_jacobians = 0

        # Execute guess_nonlinear if specified.
        system._guess_nonlinear()

        # When under a complex step from higher in the hierarchy, sometimes the step is too small
        # to trigger reconvergence, so nudge the outputs slightly so that we always get at least
        # one iteration of Broyden.
        if system.under_complex_step and self.options['cs_reconverge']:
            system._outputs._data += np.linalg.norm(self._system._outputs._data) * 1e-10

        # Start with initial states.
        self.xm = self.get_states()

        with Recording('Broyden', 0, self):
            self._solver_info.append_solver()

            # should call the subsystems solve before computing the first residual
            self._gs_iter()

            self._solver_info.pop()

        self._run_apply()
        norm = self._iter_get_norm()

        norm0 = norm if norm != 0.0 else 1.0
        return norm0, norm

def time_correlations_direct(P, pi, obs1, obs2=None, times=[1]):
    r"""Compute time-correlations of obs1, or time-cross-correlation with obs2.
    
    The time-correlation at time=k is computed by the matrix-vector expression: 
    cor(k) = obs1' diag(pi) P^k obs2
    
    
    Parameters
    ----------
    P : ndarray, shape=(n, n) or scipy.sparse matrix
        Transition matrix
    obs1 : ndarray, shape=(n)
        Vector representing observable 1 on discrete states
    obs2 : ndarray, shape=(n)
        Vector representing observable 2 on discrete states. If not given,
        the autocorrelation of obs1 will be computed
    pi : ndarray, shape=(n)
        stationary distribution vector. Will be computed if not given
    times : array-like, shape(n_t)
        Vector of time points at which the (auto)correlation will be evaluated 
    
    Returns
    -------
    
    """
    n_t = len(times)
    times = np.sort(times)  # sort it to use caching of previously computed correlations
    f = np.zeros(n_t)

    # maximum time > number of rows?
    if times[-1] > P.shape[0]:
        use_diagonalization = True
        R, D, L = rdl_decomposition(P)
        # discard imaginary part, if all elements i=0
        if not np.any(np.iscomplex(R)):
            R = np.real(R)
        if not np.any(np.iscomplex(D)):
            D = np.real(D)
        if not np.any(np.iscomplex(L)):
            L = np.real(L)
        rdl = (R, D, L)

    if use_diagonalization:
        for i in xrange(n_t):
            f[i] = time_correlation_by_diagonalization(P, pi, obs1, obs2, times[i], rdl)
    else:
        start_values = None
        for i in xrange(n_t):
            f[i], start_values = \
                time_correlation_direct_by_mtx_vec_prod(P, pi, obs1, obs2,
                                                        times[i], start_values, True)
    return f

 def __dot__(self,other):
     r1 = self.r
     r2 = other.r
     d = self.d
     if ( np.iscomplex(self.core).any() or np.iscomplex(other.core).any()):
         dt = np.zeros(r1[0]*r2[0]*r1[d]*r2[d],dtype=np.complex)
         dt = tt_f90.tt_f90.ztt_dotprod(self.n,r1,r2,self.ps,other.ps,self.core+0j,other.core+0j,dt.size)
     else:
         dt = np.zeros(r1[0]*r2[0]*r1[d]*r2[d])
         dt = tt_f90.tt_f90.dtt_dotprod(self.n,r1,r2,self.ps,other.ps,np.real(self.core),np.real(other.core),dt.size)
     if dt.size is 1:
         dt = dt[0]
     return dt

def testConnectLapEig(g):
	begin = time.time()
	print ''
	print '-----'
	print 'Computing eigenvalues of L..'
	n = len(g.nodes())
	L = np.zeros((n,n))
	for x,i in g.edges():
		L[x,i] = -1
		L[i,x] = -1
	for x in g.nodes():
		if (x,x) in g.edges():
			L[x,x] = g.degree(x)-2
		else:
			L[x,x] = g.degree(x)
	w, v = LA.eig(L)
	w = sorted(list(w))
	print ''
	print 'elapsed time:', time.time() - begin,' s'
	print ''
	print 'the eigenvalues of L are:'
	c = 0
	for x in w:
		if np.iscomplex(x):
			print 'Complex eigenvalue:',x
		else:
			print float(np.where(x < 1e-10, 0, x))
		c = c + 1
		if c == 4:
			print 'and more..'
			break

	if np.iscomplex(w[1]):
		print ''
		print 'the second smallest eigenvalue is complex:', w[1]
		print ''
		seconSmallestEig = np.real(w[1])
	else:
		seconSmallestEig = float(np.where(w[1] < 1e-10, 0, w[1]))
		print ''
		print 'the second smallest eigenvalue is:', seconSmallestEig
		print ''
	
	if seconSmallestEig > 0 :
		print 'which is positive: the graph is connected'
		print '-----'
		return True
	else:
		print 'which is negative: the graph is disconnected'
		print '-----'
		return False

	def get_projected_coordinates(self):
		self.compute_polynomials()
		# returns projected x,y coordinates based on calculated parabolas

		print(self.z_x_poly.r)
		print(self.z_y_poly.r)

		x_coord = self.z_x_poly.r[0]
		y_coord = self.z_y_poly.r[0]

		if np.iscomplex(x_coord) or np.iscomplex(y_coord):
			return (None, None)

		return (x_coord, y_coord)

def time_relaxations_direct(P, p0, obs, times = [1]):
    r"""Compute time-relaxations of obs with respect of given initial distribution.
    
    relaxation(k) = p0 P^k obs
    
    Parameters
    ----------
    P : ndarray, shape=(n, n) or scipy.sparse matrix
        Transition matrix
    p0 : ndarray, shape=(n)
        initial distribution
    obs : ndarray, shape=(n)
        Vector representing observable on discrete states. 
    times : array-like, shape(n_t)
        Vector of time points at which the (auto)correlation will be evaluated 
    
    Returns
    -------
    relaxations : ndarray, shape(n_t)
    """
    n_t = len(times)
    times = np.sort(times)
    
    # maximum time > number of rows?
    if times[-1] > P.shape[0]:
        use_diagonalization = True
        R, D, L = rdl_decomposition(P)
        # discard imaginary part, if all elements i=0
        if not np.any(np.iscomplex(R)):
            R = np.real(R)
        if not np.any(np.iscomplex(D)):
            D = np.real(D)
        if not np.any(np.iscomplex(L)):
            L = np.real(L)
        rdl = (R, D, L)
    
    f = np.empty(n_t, dtype=D.dtype)
    
    if use_diagonalization:
        for i in xrange(n_t):
            f[i] = time_relaxation_direct_by_diagonalization(
                                        P, p0, obs, times[i], rdl)
    else:
        start_values = None
        for i in xrange(n_t):
            f[i], start_values = time_relaxation_direct_by_mtx_vec_prod(
                                        P, p0, obs, times[i], start_values, True)
    return f