def lowrank_to_qr(A, B):
	"""
	Convert a low-rank factorization C appr A*B into an equivalent QR 
	C approx QR 


	Parameters:
	-----------
	A:	(m,k) ndarray
		First of the two low-rank matrices
	B:	(k,n) ndarray
		Second of the two low-rank matrices
		
	Returns:
	--------
	Q:	(m,k) ndarray
		
		
		
	
	Notes:
	------
	
	"""

	qa, ra = qr(A, mode = 'economic')
	d = dot(ra, B.conj().T)
	qd, r = qr(d, mode = 'economic')

	q = dot(qa,qd)

	return q,r 

def _overlap_projector(data_int, data_res, corr):
    """Calculate projector for removal of subspace intersection in tSSS"""
    # corr necessary to deal with noise when finding identical signal
    # directions in the subspace. See the end of the Results section in [2]_

    # Note that the procedure here is an updated version of [2]_ (and used in
    # Elekta's tSSS) that uses residuals instead of internal/external spaces
    # directly. This provides more degrees of freedom when analyzing for
    # intersections between internal and external spaces.

    # Normalize data, then compute orth to get temporal bases. Matrices
    # must have shape (n_samps x effective_rank) when passed into svd
    # computation
    Q_int = linalg.qr(_orth_overwrite((data_int / np.linalg.norm(data_int)).T),
                      overwrite_a=True, mode='economic', **check_disable)[0].T
    Q_res = linalg.qr(_orth_overwrite((data_res / np.linalg.norm(data_res)).T),
                      overwrite_a=True, mode='economic', **check_disable)[0]
    assert data_int.shape[1] > 0
    C_mat = np.dot(Q_int, Q_res)
    del Q_int

    # Compute angles between subspace and which bases to keep
    S_intersect, Vh_intersect = linalg.svd(C_mat, overwrite_a=True,
                                           full_matrices=False,
                                           **check_disable)[1:]
    del C_mat
    intersect_mask = (S_intersect >= corr)
    del S_intersect

    # Compute projection operator as (I-LL_T) Eq. 12 in [2]_
    # V_principal should be shape (n_time_pts x n_retained_inds)
    Vh_intersect = Vh_intersect[intersect_mask].T
    V_principal = np.dot(Q_res, Vh_intersect)
    return V_principal

    def _compute_xDAWN_filters(self, X, D):
        # Compute xDAWN spatial filters

        # QR decompositions of X and D
        if map(int, __import__("scipy").__version__.split('.')) >= [0,9,0]:
            # NOTE: mode='economy'required since otherwise the memory
            # consumption is excessive
            Qx, Rx = qr(X, overwrite_a=True, mode='economic')
            Qd, Rd = qr(D, overwrite_a=True, mode='economic')
        else:
            # NOTE: econ=True required since otherwise the memory consumption
            # is excessive
            Qx, Rx = qr(X, overwrite_a=True, econ=True)
            Qd, Rd = qr(D, overwrite_a=True, econ=True)

        # Singular value decomposition of Qd.T Qx
        # NOTE: full_matrices=True required since otherwise we do not get
        #       num_channels filters.
        Phi, Lambda, Psi = numpy.linalg.svd(numpy.dot(Qd.T, Qx),
                                           full_matrices=True)
        Psi = Psi.T

        # Construct the spatial filters
        for i in range(Psi.shape[1]):
            # Construct spatial filter with index i as Rx^-1*Psi_i
            ui = numpy.dot(numpy.linalg.inv(Rx), Psi[:,i])
            if i == 0:
                filters = numpy.atleast_2d(ui).T
            else:
                filters = numpy.hstack((filters, numpy.atleast_2d(ui).T))

        return filters

def random_non_singular(p, sing_min=1., sing_max=2., random_state=0):
    """Generate a random nonsingular matrix.

    Parameters
    ----------
    p : int
        The first dimension of the array.

    sing_min : float, optional (default to 1.)
        Minimal singular value.

    sing_max : float, optional (default to 2.)
        Maximal singular value.

    random_state : int or numpy.random.RandomState instance, optional
        random number generator, or seed.

    Returns
    -------
    output : numpy.ndarray, shape (p, p)
        A nonsingular matrix with the given minimal and maximal singular
        values.
    """
    random_state = check_random_state(random_state)
    diag = random_diagonal(p, v_min=sing_min, v_max=sing_max,
                           random_state=random_state)
    mat1 = random_state.randn(p, p)
    mat2 = random_state.randn(p, p)
    unitary1, _ = linalg.qr(mat1)
    unitary2, _ = linalg.qr(mat2)
    return unitary1.dot(diag).dot(unitary2.T)

def  addblock_svd_update( Uarg, Sarg, Varg, Aarg, force_orth = False):
  U = Varg
  V = Uarg
  S = np.eye(len(Sarg),len(Sarg))*Sarg
  A = Aarg.T
  
  current_rank = U.shape[1]
  m = np.dot(U.T,A)
  p = A - np.dot(U,m)
  P = lin.orth(p)
  Ra = np.dot(P.T,p)
  z = np.zeros(m.shape)
  K = np.vstack(( np.hstack((S,m)), np.hstack((z.T,Ra)) ))
  tUp,tSp,tVp = lin.svd(K);
  tUp = tUp[:,:current_rank]
  tSp = np.diag(tSp[:current_rank])
  tVp = tVp[:,:current_rank]
  Sp = tSp
  Up = np.dot(np.hstack((U,P)),tUp)
  Vp = np.dot(V,tVp[:current_rank,:])
  Vp = np.vstack((Vp, tVp[current_rank:tVp.shape[0], :]))
  
  if force_orth:
    UQ,UR = lin.qr(Up,mode='economic')
    VQ,VR = lin.qr(Vp,mode='economic')
    tUp,tSp,tVp = lin.svd( np.dot(np.dot(UR,Sp),VR.T));
    tSp = np.diag(tSp)
    Up = np.dot(UQ,tUp)
    Vp = np.dot(VQ,tVp)
    Sp = tSp;

  Up1 = Vp;
  Vp1 = Up;
    
  return Up1,Sp,Vp1

def eig(A, normal = False, iter = 100):
	'''Finds eigenvalues of an nxn array A. If A is normal, QRalg.eig 
	may also return eigenvectors.
	
	Parameters
	----------
	A :  nxn array
	     May be real or complex
	normal : bool, optional
		     Set to True if A is normal and you want to calculate
		     the eigenvectors.
	iter : positive integer, optional
			
	Returns
	-------
	v : 1xn array of eigenvectors, may be real or complex
	Q : (only returned if normal = True) 
		nxn array whose columns are eigenvectors, s.t. A*Q = Q*diag(v)
		real if A is real, complex if A is complex
	
	For more on the QR algorithm, see Eigenvalue Solvers lab.
	'''
	def getSchurEig(A):
		#Find the eigenvalues of a Schur form matrix. These are the 
		#elements on the main diagonal, except where there's a 2x2 
		#block on the main diagonal. Then we have to find the 
		#eigenvalues of that block.
		D = sp.diag(A).astype(complex)
		#Find all the 2x2 blocks:
		LD = sp.diag(A,-1)
		index = sp.nonzero(abs(LD)>.01)[0] #is this a good tolerance?
		#Find the eigenvalues of those blocks:
		a = 1
		b = -D[index]-D[index+1]
		c = D[index]*D[index+1] - A[index,index+1]*LD[index]
		discr = sp.sqrt(b**2-4*a*c)
		#Fill in vector D with those eigenvalues
		D[index] = (-b + discr)/(2*a)
		D[index+1] = (-b - discr)/(2*a)
		return D

	n,n = A.shape
	I = sp.eye(n)
	A,Q = hessenberg(A,True)
	if normal == False:
		for i in sp.arange(iter):
			s = A[n-1,n-1].copy()
			Qi,R = la.qr(A-s*I)
			A = sp.dot(R,Qi) + s*I
		v = getSchurEig(A)
		return v
	
	elif normal == True:
		for i in sp.arange(iter):
			s = A[n-1,n-1].copy()
			Qi,R = la.qr(A-s*I)
			A = sp.dot(R,Qi) + s*I
			Q = sp.dot(Q,Qi)
		v = sp.diag(A)
		return v,Q

    def normalise_site_qr(self, site, left=True):
        '''
        Normalise the MP at a given site. Use the QR method,
        arXiv:1008.3477, sec. 4.1.3
        '''
        if (left and site == L) or (not left and site == 1):
            self.normalise_extremal_site(site, left=left)
            return
        
        M_asb = self.mps[site] # the matrix we have to work on, M_site
        r, s, c = M_asb.shape  # row, spin, column

        if left:
            M_aab = M_asb.reshape((r * s, c)) # merge first two indices
            

            Q, R = la.qr(M_aab, mode='economic') # economic = thin QR in
                                                     # 1008.3477
            self.mps[site] = Q.reshape((r, s, Q.size // (r * s))) # new A_site

            # Contract R with matrix at site + 1
            self.mps[site+1] = np.tensordot(R, self.mps[site+1], ((1),(0)))
        else:
            M_aab = M_asb.reshape((r, s * c)) # merge last two indices

            # QR = M^+ ==> M = R^+ Q^+
            Q, R = la.qr(M_aab.transpose())
            Q = Q.transpose()
            R = R.transpose()
            self.mps[site] = Q.reshape((Q.size // (s * c), s, c)) # new B_site

            # Contract R with matrix at site - 1
            self.mps[site-1] = np.tensordot(self.mps[site-1], R, ((2),(0)))
        return

def eigenm(B):
    c,d= householder(B)
    a = diag(c,0)+diag(d,1)+diag(d,-1)
    q,r = ln.qr(a)
    for i in range(100):
        q,r = ln.qr(a);
        a = dot(r,q);
    return sort(diag(a,0))

def random_non_singular(shape):
    """Generates random non singular matrix"""
    d = random_diagonal_spd(shape)
    ran1 = np.random.rand(shape, shape)
    ran2 = np.random.rand(shape, shape)
    u, _ = linalg.qr(ran1)
    v, _ = linalg.qr(ran2)
    return u.dot(d).dot(v.T)

def lowrank_to_svd(A, B):
	"""M = AB^T"""

	qa, ra = qr(A, mode = 'economic')
	qb, rb = qr(B, mode = 'economic')

	mat = dot(ra,rb.conj().T)
	u, s, vh = svd(mat)

	return dot(qa, u), s, dot(qb,vh.conj().T)

def randomized_pca(A, n_components, n_oversamples=10, n_iter="auto",
                   flip_sign=True, random_state=0):
    """Compute the randomized PCA decomposition of a given matrix.

    This method differs from the scikit-learn implementation in that it supports
    and handles sparse matrices well.

    """
    if n_iter == "auto":
        # Checks if the number of iterations is explicitly specified
        # Adjust n_iter. 7 was found a good compromise for PCA. See sklearn #5299
        n_iter = 7 if n_components < .1 * min(A.shape) else 4

    n_samples, n_features = A.shape

    c = np.atleast_2d(A.mean(axis=0))

    if n_samples >= n_features:
        Q = random_state.normal(size=(n_features, n_components + n_oversamples))
        Q = safe_sparse_dot(A, Q) - safe_sparse_dot(c, Q)

        # Normalized power iterations
        for _ in range(n_iter):
            Q = safe_sparse_dot(A.T, Q) - safe_sparse_dot(c.T, Q.sum(axis=0)[None, :])
            Q, _ = lu(Q, permute_l=True)
            Q = safe_sparse_dot(A, Q) - safe_sparse_dot(c, Q)
            Q, _ = lu(Q, permute_l=True)

        Q, _ = qr(Q, mode="economic")

        QA = safe_sparse_dot(A.T, Q) - safe_sparse_dot(c.T, Q.sum(axis=0)[None, :])
        R, s, V = svd(QA.T, full_matrices=False)
        U = Q.dot(R)

    else:  # n_features > n_samples
        Q = random_state.normal(size=(n_samples, n_components + n_oversamples))
        Q = safe_sparse_dot(A.T, Q) - safe_sparse_dot(c.T, Q.sum(axis=0)[None, :])

        # Normalized power iterations
        for _ in range(n_iter):
            Q = safe_sparse_dot(A, Q) - safe_sparse_dot(c, Q)
            Q, _ = lu(Q, permute_l=True)
            Q = safe_sparse_dot(A.T, Q) - safe_sparse_dot(c.T, Q.sum(axis=0)[None, :])
            Q, _ = lu(Q, permute_l=True)

        Q, _ = qr(Q, mode="economic")

        QA = safe_sparse_dot(A, Q) - safe_sparse_dot(c, Q)
        U, s, R = svd(QA, full_matrices=False)
        V = R.dot(Q.T)

    if flip_sign:
        U, V = svd_flip(U, V)

    return U[:, :n_components], s[:n_components], V[:n_components, :]

    def _stop_training(self, debug=False):
        # The following if statement is needed only to account for
        # different versions of scipy
        if map(int, __import__("scipy").__version__.split('.')) >= [0, 9, 0]:
            # NOTE: mode='economy'required since otherwise
            # the memory consumption is excessive;
            # QR decompositions of X
            Qx, Rx = qr(self.X, overwrite_a=True, mode='economic')
            # QR decompositions of D
            Qd, Rd = qr(self.D, overwrite_a=True, mode='economic')
        else:
            # NOTE: econ=True required since otherwise
            #       the memory consumption is excessive
            # QR decompositions of X
            Qx, Rx = qr(self.X, overwrite_a=True, econ=True)
            # QR decompositions of D
            Qd, Rd = qr(self.D, overwrite_a=True, econ=True)

        # Singular value decomposition of Qd.T Qx
        # NOTE: full_matrices=True required since otherwise we do not get 
        #       num_channels filters. 
        self.Phi, self.Lambda, self.Psi = \
                    numpy.linalg.svd(numpy.dot(Qd.T, Qx), full_matrices=True)
        self.Psi = self.Psi.T
       
        SNR = numpy.zeros(self.X.shape[1])
        # Construct the spatial filters
        for i in range(self.Psi.shape[1]):
            # Construct spatial filter with index i as Rx^-1*Psi_i
            ui = numpy.dot(numpy.linalg.inv(Rx), self.Psi[:,i])
            wi = numpy.dot(Rx.T, self.Psi[:,i]) 
            if i < self.Phi.shape[1]:
                ai = numpy.dot(numpy.dot(numpy.linalg.inv(Rd), self.Phi[:,i]),
                               self.Lambda[i])
            if i == 0:
                self.filters = numpy.atleast_2d(ui).T
                self.wi = numpy.atleast_2d(wi)
                self.ai = numpy.atleast_2d(ai)
            else:
                self.filters = numpy.hstack((self.filters,
                                             numpy.atleast_2d(ui).T))
                self.wi = numpy.vstack((self.wi, numpy.atleast_2d(wi)))
                if i < self.Phi.shape[1]:
                    self.ai = numpy.vstack((self.ai, numpy.atleast_2d(ai)))
            a = numpy.dot(self.D, ai.T)
            b = numpy.dot(self.X, ui)
#            b.view(numpy.ndarray)
#            bb = numpy.dot(b.T, b)
#            aa = numpy.dot(a.T, a)
            SNR[i] = numpy.dot(a.T, a)/numpy.dot(b.T, b)

        self.SNR = SNR
        self.D = None
        self.X = None

def QRalg(A):
	#initial step
	Q, R = cacca.qr(A)
	eigenVec = Q
	Ap = R.dot(Q)

	for i in range(10):
		Q, R = cacca.qr(Ap)
		Ap = R.dot(Q)
		eigenVec = eigenvec.dot(Q)
	
	eigenVal = np.diag(Ap)
	return eigenVal, eigenVec

def test_outcome_dependent_data():
    np.random.seed(10)
    m = 1000
    max_terms = 100
    y = np.random.normal(size=m)
    w = np.random.normal(size=m) ** 2
    weight = SingleWeightDependentData.alloc(w, m, max_terms, 1e-16)
    data = SingleOutcomeDependentData.alloc(y, weight, m, max_terms)

    # Test updating
    B = np.empty(shape=(m, max_terms))
    for k in range(max_terms):
        b = np.random.normal(size=m)
        B[:, k] = b
        code = weight.update_from_array(b)
        if k >= 99:
            1 + 1
        data.update()
        assert_equal(code, 0)
        assert_almost_equal(
            np.dot(weight.Q_t[:k + 1, :], np.transpose(weight.Q_t[:k + 1, :])),
            np.eye(k + 1))
    assert_equal(weight.update_from_array(b), -1)
#     data.update(1e-16)

    # Test downdating
    q = np.array(weight.Q_t).copy()
    theta = np.array(data.theta[:max_terms]).copy()
    weight.downdate()
    data.downdate()
    weight.update_from_array(b)
    data.update()
    assert_almost_equal(q, np.array(weight.Q_t))
    assert_almost_equal(theta, np.array(data.theta[:max_terms]))
    assert_almost_equal(
        np.array(data.theta[:max_terms]), np.dot(weight.Q_t, w * y))
    wB = B * w[:, None]
    Q, _ = qr(wB, pivoting=False, mode='economic')
    assert_almost_equal(np.abs(np.dot(weight.Q_t, Q)), np.eye(max_terms))

    # Test that reweighting works
    assert_equal(data.k, max_terms)
    w2 = np.random.normal(size=m) ** 2
    weight.reweight(w2, B, max_terms)
    data.synchronize()
    assert_equal(data.k, max_terms)
    w2B = B * w2[:, None]
    Q2, _ = qr(w2B, pivoting=False, mode='economic')
    assert_almost_equal(np.abs(np.dot(weight.Q_t, Q2)), np.eye(max_terms))
    assert_almost_equal(
        np.array(data.theta[:max_terms]), np.dot(weight.Q_t, w2 * y))

def xDAWN(X,tau,N_e,remain = 5):
    '''
    xDAWN spatial filter for enhancing event-related potentials.
    
    xDAWN tries to construct spatial filters such that the 
    signal-to-signal plus noise ratio is maximized. This spatial filter is 
    particularly suited for paradigms where classification is based on 
    event-related potentials.
    
    For more details on xDAWN, please refer to 
    http://www.icp.inpg.fr/~rivetber/Publications/references/Rivet2009a.pdf

    this code is inspired by 'xdawn.py' in pySPACE:
    https://github.com/pyspace/pyspace/blob/master/pySPACE/missions/nodes/spatial_filtering/xdawn.py
    ##################
    use the same notations as in the paper linked above (../River2009a.pdf)
    N_t:the number of temporal samples. (over 100,000 per session)
    N_s:the number of sensors. (56 EEG sensors)
    @param:
        X: input EEG signal (N_t,N_s)
        tau: index list where stimulus onset
        N_e: number of temporal samples of the ERP (<1.3sec)
    return:

    '''
    N_t,N_s= X.shape
    D = build_toe_mat(N_t,N_e,tau)#construct Toeplitz matrix
    #print X.shape
    Qx, Rx = qr(X, overwrite_a = True, mode='economic')
    Qd, Rd = qr(D, overwrite_a = True, mode='economic')  
    Phi,Lambda,Psi = svd(np.dot(Qd.T,Qx),full_matrices = True)
    Psi = Psi.T
    SNR = []
    U = None
    A = None
    for i in range(remain):
        ui = np.dot(np.linalg.inv(Rx),Psi[:,i])
        ai = np.dot(np.dot(np.linalg.inv(Rd),Phi[:,i]),Lambda[i])
        if U == None:
            U = np.atleast_2d(ui).T
        else:
            U = np.hstack((U,np.atleast_2d(ui).T))
        if A == None:
            A = np.atleast_2d(ai)
        else:
            A = np.vstack((A,np.atleast_2d(ai)))
        tmp_a = np.dot(D, ai.T)
        tmp_b = np.dot(X, ui)
        #print np.dot(tmp_a.T,tmp_a)/np.dot(tmp_b.T,tmp_b)
    return U,A

def random_spd(p, eig_min, cond, random_state=0):
    """Generate a random symmetric positive definite matrix.

    Parameters
    ----------
    p : int
        The first dimension of the array.

    eig_min : float
        Minimal eigenvalue.

    cond : float
        Condition number, defined as the ratio of the maximum eigenvalue to the
        minimum one.

    random_state : int or numpy.random.RandomState instance, optional
        random number generator, or seed.

    Returns
    -------
    ouput : numpy.ndarray, shape (p, p)
        A symmetric positive definite matrix with the given minimal eigenvalue
        and condition number.
    """
    random_state = check_random_state(random_state)
    mat = random_state.randn(p, p)
    unitary, _ = linalg.qr(mat)
    diag = random_diagonal(p, v_min=eig_min, v_max=cond * eig_min,
                           random_state=random_state)
    return unitary.dot(diag).dot(unitary.T)

def _qr_economic_old(A, **kwargs):
    """
    Compat function for the QR-decomposition in economic mode
    Scipy 0.9 changed the keyword econ=True to mode='economic'
    """
    with warnings.catch_warnings(record=True):
        return linalg.qr(A, econ=True, **kwargs)

def _ar_model_qr(data, p=1):
    """QR factorization for a (multivariate) zero-mean AR model

    :Parameters:
        data : ndarray
            data with observations on the rows and variables on the columns
        p : int or list
            the model order, how many samples to regress over
    """

    # inits
    n, m = data.shape            # observations, channels
    ne = n - p                   # number of block equations of size m
    np = m * p                   # number of parameter vectors of size m
    K = N.zeros((ne, np + m))  # the lag shifted data matrix

    # compute predictors
    for i in xrange(p):
        K[:, m * i:m * (i + 1)] = data[p - i - 1:n - i - 1, :]
    K[:, np:np + m] = data[p:n, :]

    # contition the matrix and factorize
    scale = N.sqrt(((np + m) ** 2 + np + m + 1) * EPS)
    R = NL.qr(
        N.concatenate((
            K,
            scale * N.diag([NL.norm(K[:, i]) for i in xrange(K.shape[1])])
        )),
        mode='r'
    )

    # return
    del K
    return R

def groupInverse(M):
    """
    Computes the group inverse of stochastic matrix using the algorithm
    given by Golub and Meyer in:
    G. H. Golub and C. D. Meyer, Jr, SIAM J. Alg. Disc. Meth. 7, 273-
    281 (1986)

    Parameters
    ----------
        M : ndarray
            A square matrix with index 1.
    
    Returns
    -------
        grpInvM : ndarray
            The group inverse of M.
    """
    L=np.shape(M)[1]
    q,r=qr(M)
    piDist=q[:,L-1]
    piDist=(1/np.sum(piDist))*piDist
    specProjector=np.identity(L)-np.outer(np.ones(L),piDist)
    u=r[0:(L-1),0:(L-1)]#remember 0:(L-1) actually means 0 to L-2!
    uInv= inv(u)#REPLACE W. lapack, invert triangular matrix ROUTINE
    uInv = np.real(uInv)
    grpInvM=np.zeros((L,L))
    grpInvM[0:(L-1),0:(L-1)]=uInv
    grpInvM=np.dot(specProjector,np.dot(grpInvM,np.dot(q.transpose(),specProjector)))
    return grpInvM

 def test_random_complex(self):
     n = 20
     for k in range(2):
         a = random([n,n])+1j*random([n,n])
         q,r = qr(a)
         assert_array_almost_equal(dot(conj(transpose(q)),q),identity(n))
         assert_array_almost_equal(dot(q,r),a)