def find_rank_loop(sdpRelaxation, x_mat, base_level=0):
    """Helper function to detect rank loop in the solution matrix.

    :param sdpRelaxation: The SDP relaxation.
    :type sdpRelaxation: :class:`ncpol2sdpa.SdpRelaxation`.
    :param x_mat: The solution of the moment matrix.
    :type x_mat: :class:`numpy.array`.
    :param base_level: Optional parameter for specifying the lower level
                       relaxation for which the rank loop should be tested
                       against.
    :type base_level: int.
    :returns: list of int -- the ranks of the solution matrix with in the
                             order of increasing degree.
    """
    ranks = []
    from numpy.linalg import matrix_rank
    if sdpRelaxation.hierarchy != "npa":
        raise Exception("The detection of rank loop is only implemented for \
                         the NPA hierarchy")
    if base_level == 0:
        levels = range(1, sdpRelaxation.level + 1)
    else:
        levels = [base_level]
    for level in levels:
        base_monomials = \
          pick_monomials_up_to_degree(sdpRelaxation.monomial_sets[0], level)
        ranks.append(matrix_rank(x_mat[:len(base_monomials),
                                       :len(base_monomials)]))

    if x_mat.shape != (len(base_monomials), len(base_monomials)):
        ranks.append(matrix_rank(x_mat))
    return ranks

def binaryMatrixTest(n,e,M,Q): # e séquence binaire 
    epsilon= np.array([int(x) for x in list(e)]) # conversion string (entrée) en vecteur epsilon 
    print(epsilon)
    N= n//(Q*M)
    slice=Q*M # taille choisie pour le découpage de la séquence binaire 
    Ranks=[]
    Ranks += [matrix_rank(np.reshape(epsilon[0:slice],(M,Q)))] 
    # rang des matrices formées par les blocs du découpage rangés dans une liste
    for k in range(1,N):
        Ranks +=[matrix_rank(np.reshape(epsilon[k*slice +1 :(k+1)*slice+1],(M,Q)))]
#    print(Ranks) 
    
    FM=0
    FM_1=0
    #comptage des rangs valant M et M-1
    for i in range(len(Ranks)):
        if(Ranks[i]==M):
            FM +=1
        if(Ranks[i]==M-1):
            FM_1 +=1
            
    ki_carre= ((FM-0.2888*N)**2)/(0.2888*N)+(FM_1-0.5776*N)**2/(0.5776*N)+((N-FM-FM_1-0.1336*N)**2)/(0.1336*N)
    P_value= np.exp((-1)*ki_carre/2)
#    print(P_value)
    return P_value 

def test():
    from numpy.linalg import matrix_rank
    from mfpy.materials.linearelastic import LinearElastic

    # Test create
    nodes = [array((0,0)), array((1,0)), array((1,1)), array((0,1))]
    enm = [0,1,2,3]
    mat = LinearElastic(lmbda=0,mu=1.0,rho=1)
    elem = QuadRI(nodes, enm, mat, thickness=1)

    # Test internal force
    u = array([-1,0,
               +1,0,
               -1,0,
               +1,0])
    fint = elem.calc_internal_force([], u)
    K = elem.calc_linear_stiffness([], u)
    print("Rank =", matrix_rank(K))
    print(fint)
    print(K.dot(u))


    elem = Quad(nodes, enm, mat, thickness=1)
    fint = elem.calc_internal_force([], u)
    K = elem.calc_linear_stiffness([], u)
    print("Rank =", matrix_rank(K))
    print(fint)
    print(K.dot(u))

    # Test lumped mass matrix
    M = elem.calc_lumped_mass()

def test_reduced_rank():
    # Test matrices with reduced rank
    rng = np.random.RandomState(20120714)
    for i in range(100):
        # Make a rank deficient matrix
        X = rng.normal(size=(40, 10))
        X[:, 0] = X[:, 1] + X[:, 2]
        # Assert that matrix_rank detected deficiency
        assert_equal(matrix_rank(X), 9)
        X[:, 3] = X[:, 4] + X[:, 5]
        assert_equal(matrix_rank(X), 8)

def read_train_images(_data,_class):
    # Read each image
    for i in _data:
        img = cv2.imread(i,0)
        img_flat = img.flatten()
        # Build our train_array list
        train_array.append(img_flat)
    # Calculate the mean of the set
    mean = [sum(i) for i in zip(*train_array)]
    length = len(_data)
    mean = [float(i) / length for i in mean]
    # Substract mean to each image (this can be negative so python
    # transform each vector to int)
    # This is Omega
    list_norm_train = [i - mean for i in train_array]
    # Transform the list to an array
    norm_train_array = numpy.asarray(list_norm_train)
    # Compute the covariance matrix of the array
    cov = numpy.dot(norm_train_array, norm_train_array.T)
    # This line will help us to define how many eigenvectors we
    # can calculate. # of eigs = rank -1 
    print matrix_rank(cov)
    eigval, eigvec = lin.eigs(cov, 38)
    print "size of the eigen vector " + str(len(eigvec[:, 0]))
    print "size of the eigen vector matrix" + str(len(eigvec))
    print "number of eigenvalues" + str(len(eigval))
    # Each eigvec[:,i] is an eigenvector of size 40
    # We need to find u_i = A * eigvec[:,i]
    A = norm_train_array.T
    for i in range(1,(len(eigvec[0]+1))):
        u.append(numpy.dot(A,eigvec[:,i]))
    for i, val in enumerate(u):
        u[i] = u[i] / numpy.linalg.norm(u[i])
    # We're only keeping 75% of the number of eigenvector u[i]
    # This will correspond to the largest eigenvalues
    for i in range(1,(int(0.75*len(u))+4)):
        u_reduced.append(u[i])
    # u_reduced[i] are called Eigenfaces
    # Now lets represent each face in this basis
    sigma = []
    omega = []
    for i, val in enumerate(list_norm_train):
        for j, val in enumerate(u_reduced):
            w = numpy.dot(u_reduced[j].T,list_norm_train[i])
            sigma.append(w.real)
            sigma_array = numpy.asarray(sigma)
        omega.append(sigma_array)
    #print omega
    print "size of eigenvector" + str(len(omega[0]))
    print "size of omega" + str(len(omega))
    #return eigval, eigvec
    return omega

def reachable_form(xsys):
    """Convert a system into reachable canonical form

    Parameters
    ----------
    xsys : StateSpace object
        System to be transformed, with state `x`

    Returns
    -------
    zsys : StateSpace object
        System in reachable canonical form, with state `z`
    T : matrix
        Coordinate transformation: z = T * x
    """
    # Check to make sure we have a SISO system
    if not issiso(xsys):
        raise ControlNotImplemented(
            "Canonical forms for MIMO systems not yet supported")

    # Create a new system, starting with a copy of the old one
    zsys = StateSpace(xsys)

    # Generate the system matrices for the desired canonical form
    zsys.B = zeros(shape(xsys.B))
    zsys.B[0, 0] = 1.0
    zsys.A = zeros(shape(xsys.A))
    Apoly = poly(xsys.A)                # characteristic polynomial
    for i in range(0, xsys.states):
        zsys.A[0, i] = -Apoly[i+1] / Apoly[0]
        if (i+1 < xsys.states):
            zsys.A[i+1, i] = 1.0

    # Compute the reachability matrices for each set of states
    Wrx = ctrb(xsys.A, xsys.B)
    Wrz = ctrb(zsys.A, zsys.B)

    if matrix_rank(Wrx) != xsys.states:
        raise ValueError("System not controllable to working precision.")

    # Transformation from one form to another
    Tzx = solve(Wrx.T, Wrz.T).T  # matrix right division, Tzx = Wrz * inv(Wrx)

    # Check to make sure inversion was OK.  Note that since we are inverting
    # Wrx and we already checked its rank, this exception should never occur
    if matrix_rank(Tzx) != xsys.states:         # pragma: no cover
        raise ValueError("Transformation matrix singular to working precision.")

    # Finally, compute the output matrix
    zsys.C = solve(Tzx.T, xsys.C.T).T  # matrix right division, zsys.C = xsys.C * inv(Tzx)

    return zsys, Tzx

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 geneigh(A,B,tol=1e-12):
    """
    Solves the generalized eigenvalue problem also in the case where A and B share a common
    null-space. The eigenvalues corresponding to the null-space are given a Nan value.
    The null-space is defined with the tolereance tol.
    """
    # first check if there is a null-space issue
    if lg.matrix_rank(B,tol)==np.shape(B)[0]:
        return eigh(A,B)
    # first diagonalize the overlap matrix B
    Be,Bv=eigh(B)
    # rewrite the A matrix in the B-matrix eigenspace
    At=np.dot(np.conj(Bv.T),np.dot(A,Bv))
    Bt=np.diag(Be)
    # detect shared null-space. that is given by the first n null eigenvalues of B
    try:
        idx=next(i for i,v in enumerate(Be) if v>tol)
    except StopIteration:
        raise(RuntimeError('geneigh: Rank of B < B.shape[0] but null-space could not be found!'))
    # check that the B matrix null-space is shared by A.
    m=np.amax(abs(At[0:idx,:].flatten()))
    if m>tol:
        warnings.warn('Maximum non-diagonal element in A written in B null-space is bigger than the tolerance \''+str(tol)+'\'.',UserWarning)
    # diagonalize the non-null-space part of the problem
    Et,Vt=eigh(At[idx:,idx:],Bt[idx:,idx:])
    # define Ut, the change of basis in the non-truncated space
    Ut=np.zeros(np.shape(A),A.dtype)
    Ut[0:idx,0:idx]=np.eye(idx)
    Ut[idx:,idx:]=Vt
    U=np.dot(Bv,Ut)
    E=np.concatenate((float('NaN')*np.ones(idx),Et))
    return E,U

def f_test(V, R, beta, r, df_d):
    """Arbitrary F test.

    Args:
        V (array): K-by-K variance-covariance matrix.
        R (array): K-by-K Test matrix.
        beta (array): Length-K vector of coefficient estimates.
        r (array): Length-K vector of null hypotheses.
        df_d (int): Denominator degrees of freedom.

    Returns:
        tuple: A tuple containing:
            - **F** (float): F-stat.
            - **pF** (float): p-score for ``F``.
    """
    Rbr = (R.dot(beta) - r)
    if Rbr.ndim == 1:
        Rbr = Rbr.reshape(-1, 1)

    middle = la.inv(R.dot(V).dot(R.T))
    df_n = matrix_rank(R)
    # Can't just squeeze, or we get a 0-d array
    F = (Rbr.T.dot(middle).dot(Rbr)/df_n).flatten()[0]
    pF = 1 - stats.f.cdf(F, df_n, df_d)
    return F, pF

def glm_diagnostics(B_4d, design, data_4d):
    """
    Return a tuple of the MRSS in 3 dimensions, fitted values in 4 
    dimensions, and residuals in 4 dimensions. 
    
    Parameters
    ----------
    B_4d: numpy array of 4 dimensions
        The estimated coefficients
    design: numpy array
        The design matrix used to get the estimated coefficients
    data_4d: numpy array of 4 dimensions 
        The corresponding image data
    
    Returns
    -------
    diagnostics : tuple
        MRSS (3d), fitted values (4d), and residuals (4d).
    """
    B_2d = np.reshape(B_4d, (-1, B_4d.shape[-1])).T
    data_2d = np.reshape(data_4d, (-1, data_4d.shape[-1]))
    
    fitted = design.dot(B_2d)
    residuals = data_2d.T - fitted
    df = design.shape[0] - npl.matrix_rank(design)
    MRSS = (residuals**2).sum(0)/df
    
    MRSS_3d = np.reshape(MRSS.T, data_4d.shape[:-1])
    fitted_4d = np.reshape(fitted.T, data_4d.shape)
    residuals_4d = np.reshape(residuals.T, data_4d.shape)
    return MRSS_3d, fitted_4d, residuals_4d

 def t_stat(self):
     """ betas, t statistic and significance test given data,
     design matix, contrast
     This is OLS estimation; we assume the errors to have independent
     and identical normal distributions around zero for each $i$ in
     $\e_i$ (i.i.d).
     """
     if self.design is None:
         self.get_design_matrix()
     if self.t_values is None:
         y = self.data.T
         X = self.design
         c = [0, 0, 1]
         c = np.atleast_2d(c).T
         beta = npl.pinv(X).dot(y)
         fitted = X.dot(beta)
         errors = y - fitted
         RSS = (errors**2).sum(axis=0)
         df = X.shape[0] - npl.matrix_rank(X)
         MRSS = RSS / df
         SE = np.sqrt(MRSS * c.T.dot(npl.pinv(X.T.dot(X)).dot(c)))
         try:
             SE[SE == 0] = np.amin(SE[SE != 0])
         except ValueError:
             pass
         t = c.T.dot(beta) / SE
         self.t_values = abs(t[0])
     self.t_indices = np.array(self.t_values).argsort(
     )[::-1][:self.t_values.size]
     return self.t_indices

def significant(X,Y,beta):
    """
    Calculates t statistic for the first two entries of given beta estimates. 
    Particularly, this is a function to calculate t values for beta gain and 
    beta loss for a single voxel.

    Parameters:
    -----------
    X: Design matrix
    Y: Data matrix for a single voxel
    beta: beta gain/loss estimates from OLS regression of a single voxel, 
    1-d array of length = 2 
    
    Returns:
    --------
    t1, t2: t value for beta gain, t value for beta loss, type: double

    Example use for ith voxel: significant(X, Y[:,i], beta[:,i]) 
    """

    y_hat = X.dot(beta)
    residuals = Y - y_hat
    RSS = np.sum(residuals ** 2)
    df = X.shape[0] - npl.matrix_rank(X)
    MRSS = RSS / df
    s2 = MRSS
    v_cov = s2 * npl.inv(X.T.dot(X))
    numerator1 = beta[0]
    denominator1 = np.sqrt(v_cov[0, 0])
    t1= numerator1 / denominator1
    numerator2 = beta[1]
    denominator2 = np.sqrt(v_cov[1, 1])
    t2= numerator2 / denominator2
    return t1,t2

def t_stat(data, X_matrix):
    """
    Return the estimated betas, t-values, degrees of freedom, and p-values for the glm_multi regression
    
    Parameters
    ----------
    data_4d: numpy array of 4 dimensions 
             The image data of one subject, one run
    X_matrix: numpy array 
       The design matrix for glm_multi
    Note that the fourth dimension of `data_4d` (time or the number 
    of volumes) must be the same as the number of rows that X has. 
    
    Returns
    -------
    beta: estimated beta values
    
    t: t-values of the betas
    
    df: degrees of freedom
    
    p: p-values corresponding to the t-values and degrees of freedom
    """

    beta = glm_beta(data, X_matrix)

    # Calculate the parameters - b hat
    beta = np.reshape(beta, (-1, beta.shape[-1])).T

    fitted = X_matrix.dot(beta)
    # Residual error
    y = np.reshape(data, (-1, data.shape[-1]))
    errors = y.T - fitted
    # Residual sum of squares
    RSS = (errors**2).sum(axis=0)
 
    df = X_matrix.shape[0] - npl.matrix_rank(X_matrix)
    # Mean residual sum of squares
    MRSS = RSS / df
    # calculate bottom half of t statistic
    Cov_beta=npl.pinv(X_matrix.T.dot(X_matrix))

    SE =np.zeros(beta.shape)
    for i in range(X_matrix.shape[-1]):
        c = np.zeros(X_matrix.shape[-1])
        c[i]=1
        c = np.atleast_2d(c).T
        SE[i,:]= np.sqrt(MRSS* c.T.dot(npl.pinv(X_matrix.T.dot(X_matrix)).dot(c)))


    zeros = np.where(SE==0)
    SE[zeros] = 1
    t = beta / SE

    t[:,zeros] =0
    # Get p value for t value using CDF of t didstribution
    ltp = t_dist.cdf(abs(t), df)
    p = 1 - ltp # upper tail
    
    return beta.T, t, df, p

def test_glm():
    # Read in the image data.
    img = nib.load(pathtoclassdata + "ds114_sub009_t2r1.nii")
    data = img.get_data()[..., 4:]
    # Read in the convolutions. 
    convolved = np.loadtxt(pathtoclassdata + "ds114_sub009_t2r1_conv.txt")[4:]
    # Create design matrix. 
    actual_design = np.ones((len(convolved), 2))
    actual_design[:, 1] = convolved
    
    # Calculate betas, copied from the exercise. 
    data_2d = np.reshape(data, (-1, data.shape[-1]))
    actual_B = npl.pinv(actual_design).dot(data_2d.T)
    actual_B_4d = np.reshape(actual_B.T, img.shape[:-1] + (-1,))
    
    # Run function.
    exp_B_4d, exp_design = glm(data, convolved)
    assert_almost_equal(actual_B_4d, exp_B_4d)
    assert_almost_equal(actual_design, exp_design)

    # Pick a single voxel to check diagnostics. 
    # Calculate actual fitted values, residuals, and MRSS of voxel.
    actual_fitted = actual_design.dot(actual_B_4d[42, 32, 19])
    actual_residuals = data[42, 32, 19] - actual_fitted
    actual_MRSS = np.sum(actual_residuals**2)/(actual_design.shape[0] - npl.matrix_rank(actual_design))
    
    # Calculate using glm_diagnostics function.
    exp_MRSS, exp_fitted, exp_residuals = glm_diagnostics(exp_B_4d, exp_design, data)
    assert_almost_equal(actual_fitted, exp_fitted[42, 32, 19])
    assert_almost_equal(actual_residuals, exp_residuals[42, 32, 19])
    assert_almost_equal(actual_MRSS, exp_MRSS[42, 32, 19])

def __LDL__(A, combined=False):
    import numpy.linalg as nplinalg
    assert(A.shape[0] == A.shape[1])

    L = np.zeros(A.shape)
    D = np.zeros(A.shape)
    n = A.shape[0]

    for i in xrange(n):
        for j in xrange(n):
            if i == j:
                D[i, i] = A[i, i]
                for k in xrange(i):
                    D[i, i] -= (L[i, k] ** 2) * D[k, k]
                L[i, i] = 1
            elif j <= i:
                L[i, j] = A[i, j]
                for k in xrange(j):
                    L[i, j] -= L[i, k] * D[k, k] * L[j, k]
                L[i, j] *= 1/D[j, j]

    if combined:
        return np.dot(L, np.sqrt(D[:,:nplinalg.matrix_rank(A)]))
    else:
        return L, D

def geneigh(A,B,tol=1e-12):
    """
    Solves the generalized eigenvalue problem also in the case where A and B share a common
    null-space. The eigenvalues corresponding to the null-space are given a Nan value.
    The null-space is defined with the tolereance tol.
    """
    # first check if there is a null-space issue
    if matrix_rank(B,tol)==shape(B)[0]:
        return eigh(A,B)
    # first diagonalize the overlap matrix B
    Be,Bv=eigh(B)
    # rewrite the A matrix in the B-matrix eigenspace
    At=dot(conj(Bv.T),dot(A,Bv))
    Bt=diag(Be)
    # detect shared null-space. that is given by the first n null eigenvalues of B
    idx=find(Be>tol)
    idx=idx[0]
    # check that the B matrix null-space is shared by A.
    m=amax(abs(At[0:idx,:].flatten()))
    if m>tol:
        warnings.warn('Maximum non-diagonal element in A written in B null-space is bigger than the tolerance \''+str(tol)+'\'.',UserWarning)
    # diagonalize the non-null-space part of the problem
    Et,Vt=eigh(At[idx:,idx:],Bt[idx:,idx:])
    # define Ut, the change of basis in the non-truncated space
    Ut=zeros(shape(A),A.dtype)
    Ut[0:idx,0:idx]=eye(idx)
    Ut[idx:,idx:]=Vt
    U=dot(Bv,Ut)
    E=append(float('NaN')*ones(idx),Et)
    return E,U

def Regression_Calculation():
    Y = matrix( Get_Y_Matrix(con), dtype = float )
    X = matrix( Get_X_Matrix(con), dtype = float )
    if matrix_rank(X) == 32:
        return (X.T * X).I * X.T * Y.T
    else:
        return numpy.linalg.pinv(X) * Y.T

def simplex(graph, wp_trajs, withODs=False):
    """Build simplex constraints from waypoint trajectories wp_trajs
    wp_trajs is given by WP.get_wp_trajs()[1]
    
    Parameters:
    -----------
    graph: Graph object
    wp_trajs: list of waypoint trajectories with paths along this trajectory [(wp_traj, path_list, flow)]
    """
    if wp_trajs is None:
        return None, None
    n = len(wp_trajs)
    I, J, r, i = [], [], matrix(0.0, (n,1)), 0
    for wp_traj, path_ids, flow in wp_trajs:
        r[i] = flow
        for id in path_ids:
            I.append(i)
            J.append(graph.indpaths[id])
        i += 1
    U = spmatrix(1.0, I, J, (n, graph.numpaths))
    if not withODs: return U, r
    else:
        U1, r1 = path.simplex(graph)
        U, r = matrix([U, U1]), matrix([r, r1])
        if la.matrix_rank(U) < U.size[0]:
            logging.info('Remove redundant constraint(s)'); ind = find_basis(U.trans())
            return U[ind,:], r[ind]
        return U, r

def comp(M):
    """Returns a basis for the space orthogonal
    to the range of M
    """
    I = eye(M.shape[0])
    Q,R = qr(concatenate((M,I),axis=1))
    return Q[:,matrix_rank(M):]

def train(x,y):
	"""
		Build the linear least weight vector W
		:param x: NxD matrix containing N attributes vectors for training
		:param y: NxK matrix containing N class vectors for training
	"""
	# D = Number of attributes
	D = x.shape[1] + 1
	# K = Number of classes
	K = y.shape[1]
	
	# Build the sums of xi*xi' and xi*yi'
	sum1 = np.zeros((D,D)) # init placeholder
	sum2 = np.zeros((D,K))
	i = 0
	for x_i in x:						# loop over all vectors
		x_i = np.append(1, x_i) 		# augment vector with a 1 
		y_i = y[i]						
		sum1 += np.outer(x_i, x_i)		# find xi*xi'
		sum2 += np.outer(x_i, y_i)		# find xi*yi'
		i += 1
	
	# Check that condition number is finite
	# and therefore sum1 is nonsingular (invertable)
	while matrix_rank(sum1) != D:
		# Naive choice of sigma.
		# Could cause inaccuracies when sum1 has small values
		# However, in most cases the matrix WILL be invertable
		sum1 = sum1 + 0.001 * np.eye(D) 
	
	# Return weight vector
	# Weight vector multiplies sums and inverse of sum1
	return np.dot(inv(sum1),sum2)