def principal_angles(A, B):
    '''Compute the principal angles between subspaces A and B.

    The algorithm for computing the principal angles is described in :
    A. V. Knyazev and M. E. Argentati,
    Principal Angles between Subspaces in an A-Based Scalar Product: 
    Algorithms and Perturbation Estimates. SIAM Journal on Scientific Computing, 
    23 (2002), no. 6, 2009-2041.
    http://epubs.siam.org/sam-bin/dbq/article/37733
    '''    
    # eps = np.finfo(np.float64).eps**.981
    # for i in range(A.shape[1]):
    #     normi = la.norm(A[:,i],np.inf)
    #     if normi > eps: A[:,i] = A[:,i]/normi
    # for i in range(B.shape[1]):
    #     normi = la.norm(B[:,i],np.inf)
    #     if normi > eps: B[:,i] = B[:,i]/normi
    QA = sl.orth(A)
    QB = sl.orth(B)
    _, s, Zs = svd(QA.T.dot(QB), full_matrices=False)
    s = np.minimum(s, ones_like(s))
    theta = np.maximum(np.arccos(s), np.zeros_like(s))
    V = QB.dot(Zs)
    idxSmall = s > np.sqrt(2.)/2.
    if np.any(idxSmall):
        RB = V[:,idxSmall]
        _, x, _ = svd(RB-QA.dot(QA.T.dot(RB)),full_matrices=False)
        thetaSmall = np.flipud(np.maximum(arcsin(np.minimum(x, ones_like(x))), zeros_like(x)))
        theta[idxSmall] = thetaSmall
    return theta

 def __init__(self, basef, 
              translate=True, 
              rotate=False, 
              conditioning=None, 
              asymmetry=None,
              oscillate=False, 
              penalize=None,
              ):
     FunctionEnvironment.__init__(self, basef.xdim, basef.xopt)
     self.desiredValue = basef.desiredValue            
     self.toBeMinimized = basef.toBeMinimized
     
     if translate:            
         self.xopt = (rand(self.xdim) - 0.5) * 9.8
         
     self._diags = eye(self.xdim)            
     self._R = eye(self.xdim)            
     self._Q = eye(self.xdim)            
     
     if conditioning is not None:
         self._diags = generateDiags(conditioning, self.xdim)
     if rotate:
         self._R = orth(rand(basef.xdim, basef.xdim))        
         if conditioning:
             self._Q = orth(rand(basef.xdim, basef.xdim))
                     
     tmp = lambda x: dot(self._Q, dot(self._diags, dot(self._R, x-self.xopt)))
     if asymmetry is not None:
         tmp2 = tmp
         tmp = lambda x: asymmetrify(tmp2(x), asymmetry)
     if oscillate:
         tmp3 = tmp
         tmp = lambda x: oscillatify(tmp3(x))
     
     self.f = lambda x: basef.f(tmp(x))

def subspace(a, b, deg=True):
    """
    Angle between two subspaces specified by the columns of a and b
    Ported from MATLAB 'subspace' function

    Parameters
    ----------
    a : matrix
    b : matrix
    deg : bool
        return degree or radian

    Returns
    -------
    double
        angle
    """
    warnings.warn(
        "Deprecated. Use scipy.linalg.subspace_angles instead.", FutureWarning
    )
    oa = linalg.orth(a)
    ob = linalg.orth(b)
    if oa.shape[1] < ob.shape[1]:
        oa, ob = ob.copy(), oa.copy()
    ob -= oa @ (oa.T @ ob)
    rad = np.arcsin(min(1, linalg.norm(ob, ord=2)))
    return np.degrees(rad) if deg else rad

 def __init__(self, *args, **kwargs):
     MultiModalFunction.__init__(self, *args, **kwargs)
     self._mu0 = 2.5
     self._s = 1 - 1 / (2 * sqrt(self.xdim + 20) - 8.2)
     self._mu1 = -sqrt((self._mu0 ** 2 - 1) / self._s)
     self._signs = sign(randn(self.xdim))
     self._R1 = orth(rand(self.xdim, self.xdim))
     self._R2 = orth(rand(self.xdim, self.xdim))
     self._diags = generateDiags(100, self.xdim)

def fubini_study(A, B):
    '''
    fubini_study(A, B) Compute the Fubini-Study distance
    Compute the Fubini-Study distance based on principal angles between A and B
    as d=\acos{ \prod_i \theta_i}
    '''
    if A.shape != B.shape:
        raise ValueError('Atoms have different dim (', A.shape, ' and ', B.shape,'). Error raised in fubini_study(A, B)')
    if np.allclose(A, B): return 0.
    return arccos(det(sl.orth(A).T.dot(sl.orth(B))))

def chordal(A, B):
    '''
    chordal(A, B) Compute the chordal distance
    Compute the chordal distance between A and B
    as d=\sqrt{K - ||\bar{A}^T\bar{B}||_F^2}
    where K is the rank of A and B, || . ||_F is the Frobenius norm,
    \bar{A} is the orthogonal basis associated with A and the same goes for B.
    '''
    if A.shape != B.shape:
        raise ValueError('Atoms have not the same dimension (', A.shape, ' and ', B.shape,'). Error raised in chordal(A, B)')
    
    if np.allclose(A, B): return 0.
    else: 
        d2 = A.shape[1] - norm(sl.orth(A).T.dot(sl.orth(B)), 'fro')**2
        if d2 < 0.: return sqrt(abs(d2))
        else: return sqrt(d2)

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 fastica_defl(X, nIC=None, guess=None,
             nonlinfn = pow3nonlin,
             termtol = 5e-7, maxiters = 2e3):
    nPC, siglen = X.shape
    nIC = nIC or nPC-1
    guess = guess or randn(nPC,nIC)

    if _orth_loaded:
        guess = orth(guess)

    B = zeros(guess.shape, np.float64)

    errvec = []
    icc = 0
    while icc < nIC:
        w = randn(nPC,1) - 0.5
        w -= dot(dot(B, transp(B)), w)
        w /= norm(w)

        wprev = zeros(w.shape)
        for i in xrange(long(maxiters) +1):
            w -= dot(dot(B, transp(B)), w)
            w /= norm(w)
            #wprev = w.copy()
            if (norm(w-wprev) < termtol) or (norm(w + wprev) < termtol):
                B[:,icc]  = transp(w)
                icc += 1
                break
            wprev = w.copy()
    return B.real, errvec

def bench_lobpcg_mikota():
    print()
    print('                 lobpcg benchmark using mikota pairs')
    print('==============================================================')
    print('      shape      | blocksize |    operation   |   time   ')
    print('                                              | (seconds)')
    print('--------------------------------------------------------------')
    fmt = ' %15s |   %3d     |     %6s     | %6.2f '

    m = 10
    for n in 128, 256, 512, 1024, 2048:
        shape = (n, n)
        A, B = _mikota_pair(n)
        desired_evs = np.square(np.arange(1, m+1))

        tt = time.clock()
        X = rand(n, m)
        X = orth(X)
        LorU, lower = cho_factor(A, lower=0, overwrite_a=0)
        M = LinearOperator(shape,
                matvec=partial(_precond, LorU, lower),
                matmat=partial(_precond, LorU, lower))
        eigs, vecs = lobpcg(A, X, B, M, tol=1e-4, maxiter=40)
        eigs = sorted(eigs)
        elapsed = time.clock() - tt
        assert_allclose(eigs, desired_evs)
        print(fmt % (shape, m, 'lobpcg', elapsed))

        tt = time.clock()
        w = eigh(A, B, eigvals_only=True, eigvals=(0, m-1))
        elapsed = time.clock() - tt
        assert_allclose(w, desired_evs)
        print(fmt % (shape, m, 'eigh', elapsed))

    def _create_SDP(self):
        """ Creates the SDP knockoff of X"""
 
        # Check for rank deficiency (will add later).
 
        # SVD and come up with perpendicular matrix
        U, d, V = nplin.svd(self.X,full_matrices=True) 
        d[d<0] = 0
        U_perp = U[:,self.p:(2*self.p)]
        if self.randomize:
            U_perp = np.dot(U_perp,splin.orth(npran.randn(self.p,self.p)))
 
        # Compute the Gram matrix and its (pseudo)inverse.
        G     = np.dot(V.T * d**2 ,V)
        G_inv = np.dot(V.T * d**-2,V)
 
        # Optimize the parameter s of Equation 1.3 using SDP.
        self.s = solve_sdp(G)
        self.s[s <= self.zerotol] = 0
 
        # Construct the knockoff according to Equation 1.4:
        C_U,C_d,C_V = nplin.svd(2*np.diag(s) - (self.s * G_inv.T).T * self.s)
        C_d[C_d < 0] = 0
        X_ko = self.X - np.dot(self.X,G_inv*s) + np.dot(U_perp*np.sqrt(C_d),C_V)
        self.X_lrg = np.concatenate((self.X,X_ko), axis=1)

def random_walk(G, initial_prob, subspace_dim=3, walk_steps=3):
    assert type(initial_prob) == np.ndarray, "Initial probability distribution is not a numpy array"

    # Transform the adjacent matrix to a laplacian matrix P
    P = adj_to_laplacian(G)

    prob_matrix = np.zeros((G.shape[0], subspace_dim))
    prob_matrix[:, 0] = initial_prob
    for i in range(1, subspace_dim):
        prob_matrix[:, i] = np.dot(prob_matrix[:, i - 1], P)

    orth_prob_matrix = splin.orth(prob_matrix)

    for i in range(walk_steps):
        temp = np.dot(orth_prob_matrix.T, P)
        orth_prob_matrix = splin.orth(temp.T)
    return orth_prob_matrix

def calc_subspace_proj_error(U, U_hat, ortho=False):
    """Calculate the normalized projection error between two orthogonal subspaces.
    Keyword arguments:
    U: ground truth subspace
    U_hat: estimated subspace
    """
    if not ortho:
        U = splinalg.orth(U)
        U_hat = splinalg.orth(U_hat)

    I = np.identity(U.shape[0])
    top = np.linalg.norm((I - U_hat @ U_hat.T) @ U, ord="fro")
    bottom = np.linalg.norm(U, ord="fro")

    error = float(top) / float(bottom)

    return error

def condex(n, k=4, theta=100):
    """
    CONDEX   `Counterexamples' to matrix condition number estimators.
         CONDEX(N, K, THETA) is a `counterexample' matrix to a condition
         estimator.  It has order N and scalar parameter THETA (default 100).
         If N is not equal to the `natural' size of the matrix then
         the matrix is padded out with an identity matrix to order N.
         The matrix, its natural size, and the estimator to which it applies
         are specified by K (default K = 4) as follows:
             K = 1:   4-by-4,     LINPACK (RCOND)
             K = 2:   3-by-3,     LINPACK (RCOND)
             K = 3:   arbitrary,  LINPACK (RCOND) (independent of THETA)
             K = 4:   N >= 4,     SONEST (Higham 1988)
         (Note that in practice the K = 4 matrix is not usually a
          counterexample because of the rounding errors in forming it.)

         References:
         A.K. Cline and R.K. Rew, A set of counter-examples to three
            condition number estimators, SIAM J. Sci. Stat. Comput.,
            4 (1983), pp. 602-611.
         N.J. Higham, FORTRAN codes for estimating the one-norm of a real or
            complex matrix, with applications to condition estimation
            (Algorithm 674), ACM Trans. Math. Soft., 14 (1988), pp. 381-396.
    """

    if k == 1:  # Cline and Rew (1983), Example B.

        a = np.array([[1, -1, -2 * theta, 0], [0, 1, theta, -theta], [0, 1, 1 + theta, -(theta + 1)], [0, 0, 0, theta]])

    elif k == 2:  # Cline and Rew (1983), Example C.

        a = np.array([[1, 1 - 2 / theta ** 2, -2], [0, 1 / theta, -1 / theta], [0, 0, 1]])

    elif k == 3:  # Cline and Rew (1983), Example D.

        a = rogues.triw(n, -1).T
        a[-1, -1] = -1

    elif k == 4:  # Higham (1988), p. 390.

        x = np.ones((n, 3))  # First col is e
        x[1:n, 1] = np.zeros(n - 1)  # Second col is e(1)

        # Third col is special vector b in SONEST
        x[:, 2] = ((-1) ** np.arange(n)) * (1 + np.arange(n) / (n - 1))

        # Q*Q' is now the orthogonal projector onto span(e(1),e,b)).
        q = sl.orth(x)
        p = np.eye(n) - np.asmatrix(q) * np.asmatrix(q.T)
        a = np.eye(n) + theta * p

    # Pad out with identity as necessary.
    m, m = a.shape
    if m < n:
        for i in range(n - 1, m, -1):
            a[i, i] = 1

    return a

 def __init__(self, *args, **kwargs):
     MultiModalFunction.__init__(self, *args, **kwargs)
     self._opts = (rand((self.numPeaks, self.xdim)) - 0.5) * 9.8
     self._opts[0] = (rand(self.xdim) - 0.5) * 8
     alphas = [power(self.maxCond, 2 * i / float(self.numPeaks - 2)) for i in range(self.numPeaks - 1)]
     shuffle(alphas)
     self._covs = [generateDiags(alpha, self.xdim, shuffled=True) / power(alpha, 0.25) for alpha in [self.optCond] + alphas]
     self._R = orth(rand(self.xdim, self.xdim))
     self._ws = [10] + [1.1 + 8 * i / float(self.numPeaks - 2) for i in range(self.numPeaks - 1)]

    def fit(self, data):
        """
        Fit independent components using an iterative fixed-point algorithm

        Parameters
        ----------
        data: RDD of (tuple, array) pairs, or RowMatrix
            Data to estimate independent components from

        Returns
        ----------
        self : returns an instance of self.
        """

        if type(data) is not RowMatrix:
            data = RowMatrix(data)

        # reduce dimensionality
        svd = SVD(k=self.k, method=self.svdmethod).calc(data)

        # whiten data
        whtmat = real(dot(inv(diag(svd.s/sqrt(data.nrows))), svd.v))
        unwhtmat = real(dot(transpose(svd.v), diag(svd.s/sqrt(data.nrows))))
        wht = data.times(whtmat.T)

        # do multiple independent component extraction
        if self.seed != 0:
            random.seed(self.seed)
        b = orth(random.randn(self.k, self.c))
        b_old = zeros((self.k, self.c))
        iter = 0
        minabscos = 0
        errvec = zeros(self.maxiter)

        while (iter < self.maxiter) & ((1 - minabscos) > self.tol):
            iter += 1
            # update rule for pow3 non-linearity (TODO: add others)
            b = wht.rows().map(lambda x: outer(x, dot(x, b) ** 3)).sum() / wht.nrows - 3 * b
            # make orthogonal
            b = dot(b, real(sqrtm(inv(dot(transpose(b), b)))))
            # evaluate error
            minabscos = min(abs(diag(dot(transpose(b), b_old))))
            # store results
            b_old = b
            errvec[iter-1] = (1 - minabscos)

        # get un-mixing matrix
        w = dot(transpose(b), whtmat)

        # get components
        sigs = data.times(w.T).rdd

        self.w = w
        self.sigs = sigs

        return self

 def rvs(self, n=1):
     """
     It ignores the dimension ``n``.
     """
     D = self.D
     d = self.d
     k = self.k
     W_sub = orth(randn(D - k, d - k))
     return np.vstack([np.hstack([W_sub, np.zeros((D - k, k))]),
                     np.hstack([np.zeros((k, d - k)), np.eye(k)])]).flatten()

def init(X, init, ncomp):
    N, K = X[0].shape[0], len(X)
    if init == 'random':
        A = orth(rand(N, ncomp))
    elif init == 'nvecs':
        S = zeros(N, N)
        for k in range(K):
            S = S + X[k] + X[k].T
        _, A = eigsh(S, ncomp)
    return A

 def __init__(self, basef, rotMat = None):
     """ by default the rotation matrix is random. """
     FunctionEnvironment.__init__(self, basef.xdim, basef.xopt)
     if rotMat == None:
         # make a random orthogonal rotation matrix
         self.M = orth(rand(basef.xdim, basef.xdim))
     else:
         self.M = rotMat
     if isinstance(basef, FunctionEnvironment):
         self.desiredValue = basef.desiredValue
     self.xopt = dot(inv(self.M), self.xopt)
     self.f = lambda x: basef.f(dot(x,self.M))

def svd3(data, k, meanSubtract=1):

    n = data.count()
    d = len(data.first())

    if meanSubtract == 1:
        data = data.map(lambda x: x - mean(x))

    def outerProd(x):
        return outer(x, x)

    def outerSum(iterator):
        yield sum(outer(x, x) for x in iterator)

    def outerSum2(iterator, other1, other2):
        yield sum(outer(x, dot(dot(x, other1), other2)) for x in iterator)

    C = random.rand(k, d)
    iterNum = 0
    iterMax = 10
    error = 100
    tol = 0.000001

    while (iterNum < iterMax) & (error > tol):
        Cold = C
        Cinv = dot(transpose(C), inv(dot(C, transpose(C))))
        preMult1 = data.context.broadcast(Cinv)
        # X = data.times(preMult1.value)
        # XX' = X.cov()
        XX = data.map(lambda x: outerProd(dot(x, preMult1.value))).reduce(lambda x, y: x + y)
        XXinv = inv(XX)
        preMult2 = data.context.broadcast(dot(Cinv, XXinv))
        # data1 = data.times(dot(Cinv, inv(XX'))
        # C = data.times(data1)
        C = data.map(lambda x: outer(x, dot(x, preMult2.value))).reduce(lambda x, y: x + y)
        C = transpose(C)

        error = sum(sum((C-Cold) ** 2))
        iterNum += 1

    C = transpose(orth(transpose(C)))
    # cov = data.times(transpose(C)).cov()
    cov = data.map(lambda x: dot(x, transpose(C))).mapPartitions(outerSum).reduce(
        lambda x, y: x + y) / n
    w, v = eig(cov)
    w = real(w)
    v = real(v)
    inds = argsort(w)[::-1]
    latent = w[inds[0:k]]
    comps = dot(transpose(v[:, inds[0:k]]), C)
    scores = data.map(lambda x: inner(x, comps))

    return comps, latent, scores

 def time_mikota(self, n, solver):
     m = 10
     if solver == 'lobpcg':
         X = rand(n, m)
         X = orth(X)
         LorU, lower = cho_factor(self.A, lower=0, overwrite_a=0)
         M = LinearOperator(self.shape,
                 matvec=partial(_precond, LorU, lower),
                 matmat=partial(_precond, LorU, lower))
         eigs, vecs = lobpcg(self.A, X, self.B, M, tol=1e-4, maxiter=40)
     else:
         w = eigh(self.A, self.B, eigvals_only=True, eigvals=(0, m-1))