def test_betti0_periodic(horiz_complex, vert_complex):
    """
    Verify that the 0-form Hodge Laplacian has kernel of dimension
    equal to the 0th Betti number of the periodic extruded interval,
    i.e. 1.  Also verify that the 0-form Hodge Laplacian with
    Dirichlet boundary conditions has kernel of dimension equal to the
    2nd Betti number of the extruded mesh, i.e. 0.
    """
    U0, U1 = horiz_complex
    V0, V1 = vert_complex

    m = PeriodicUnitIntervalMesh(5)
    mesh = ExtrudedMesh(m, layers=4, layer_height=0.25)
    U0 = FiniteElement(U0[0], "interval", U0[1])
    V0 = FiniteElement(V0[0], "interval", V0[1])

    W0_elt = TensorProductElement(U0, V0)
    W0 = FunctionSpace(mesh, W0_elt)

    u = TrialFunction(W0)
    v = TestFunction(W0)

    L = assemble(inner(grad(u), grad(v))*dx)
    uvecs, s, vvecs = linalg.svd(L.M.values)
    nharmonic = sum(s < 1.0e-5)
    assert(nharmonic == 1)

    bcs = [DirichletBC(W0, 0., x) for x in ["top", "bottom"]]
    L = assemble(inner(grad(u), grad(v))*dx, bcs=bcs)
    uvecs, s, vvecs = linalg.svd(L.M.values)
    nharmonic = sum(s < 1.0e-5)
    assert(nharmonic == 0)

def test_betti0(space, mesh):
    """
    Verify that the 0-form Hodge Laplacian with strong Dirichlet
    boundary conditions has kernel of dimension equal to the 2nd Betti
    number of the annulus mesh, i.e. 0.
    """
    V0tag, V1tag, V2tag = space

    if(len(V0tag) == 2):
        V0 = FunctionSpace(mesh, V0tag[0], V0tag[1])
    else:
        V0a = FiniteElement(V0tag[0], "triangle", V0tag[1])
        V0b = FiniteElement(V0tag[2], "triangle", V0tag[3])
        V0 = FunctionSpace(mesh, V0a + V0b)
    # V0 Hodge Laplacian
    u = TrialFunction(V0)
    v = TestFunction(V0)

    L = assemble(inner(nabla_grad(u), nabla_grad(v))*dx)

    bc0 = DirichletBC(V0, Constant(0.0), 9)
    L0 = assemble(inner(nabla_grad(u), nabla_grad(v))*dx, bcs=[bc0])

    u, s, v = linalg.svd(L.M.values)
    nharmonic = sum(s < 1.0e-5)
    assert(nharmonic == 1)

    u, s, v = linalg.svd(L0.M.values)
    nharmonic = sum(s < 1.0e-5)
    assert(nharmonic == 0)

def admira(r, b, m, n, iter, A, A_star):
	if 2*r > min(m,n):
		r_prime = min(m,n)
	else:
		r_prime = 2*r

	# initialization
	X_hat = np.random.randn(m,n) # step 1
	Psi_hatU = np.matrix([])
	Psi_hatV = np.matrix([])
	for i in range(iter):
		Y = A_star(b - A(X_hat))
		(U, s, Vt) = svd(Y)
		Psi_primeU = U[:, 0:r_prime]
		Psi_primeV = Vt.T[:, 0:r_prime]
		if i > 0:
			Psi_tildeU = np.bmat([Psi_primeU, Psi_hatU])
			Psi_tildeV = np.bmat([Psi_primeV, Psi_hatV])
		else:
			Psi_tildeU = Psi_primeU
			Psi_tildeV = Psi_primeV
		AP = lambda b: APsiUV(b, A, Psi_tildeU, Psi_tildeV)
		APt = lambda s: APsitUV(s, A_star, Psi_tildeU, Psi_tildeV)
		ALS = lambda b: APt(AP(b))
		(s, res, iter) = cgsolve(ALS, APt(b), 1e-6, 100, False)
		X_tilde = Psi_tildeU*np.matrix(np.diag(np.array(s).reshape(-1)))*Psi_tildeV.T
		(U, s, Vt) = svd(X_tilde)
		Psi_hatU = U[:, 0:r]
		Psi_hatV = Vt.T[:, 0:r]
		X_hat = Psi_hatU*np.diag(s[0:r])*Psi_hatV.T

	return X_hat

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,data):
     self.data = data
     self.N,self.p = self.data.shape
     self.covMatrix = covmatrix(self.data)
     self.corMatrix = corrmatrix(self.data)
     _,self.covlambda,_ = svd(self.covMatrix, full_matrices=False)
     _,self.corlambda,_ = svd(self.corMatrix, full_matrices=False)

def frequent_directions(A, ell):
  """A matrix "A" should be 256x7291
  """

  m = 256
  n = 7291

  if A.shape[0] != m or A.shape[1] != n: raise ValueError('Error: incorrect matrix size')

  start = time.clock()

  B = np.hstack((A[:, :(ell-1)], np.zeros((m, 1))))

  for i in range(ell-1, n):

    # new matrix is just a single vector (i-th column of A)
    B[:, ell-1] = A[:, i]
    U, s, V = ln.svd(B, full_matrices=False)

    delta = s[-1] ** 2 # squared smallest singular value

    B = np.dot(U, np.diag(np.sqrt(abs(s ** 2 - delta))))

  U, s, V = ln.svd(B, full_matrices=False)

  elapsed_time = time.clock() - start
  print 'time:', elapsed_time

  return U, s, V

def AND(C, B):
	
	dim, col = C.shape
	tolerance = 1e-14

	UC, SC, UtC = svd(C)
	UB, SB, UtB = svd(B)

	diag_SC = diag(SC)
	diag_SB = diag(SB)

	# sum up how many elements on diagonal 
	# are bigger than tolerance
	numRankC =  (1.0 * (diag_SC > tolerance)).sum()
	numRankB =  (1.0 * (diag_SB > tolerance)).sum()

	UC0 = matrix(UC[:, numRankC:])
	UB0 = matrix(UB[:, numRankB:])
	W, Sigma, Wt = svd(UC0 * UC0.transpose() + UB0 * UB0.transpose())
	numRankSigma =  (1.0 * (diag(Sigma) > tolerance)).sum()
	Wgk = matrix(W[:, numRankSigma:])
	I = matrix(identity(dim))
	CandB = \
	  Wgk * inv(Wgk.transpose() *  \
	  ( pinv(C, tolerance) + pinv(B, tolerance) - \
	    I) * Wgk) *Wgk.transpose()
	return CandB

def get_eigenvectors(spikes, nfeats, nspikes):
    """Calculate eigenvectors of spike waveforms

    spikes: resampled and aligned spike waveforms, dimensions (nspikes, nsamples)
    nfeats: the number of the most significant eigenvectors to return
    nspikes: the number of spikes to use

    Returns eigenvectors, dimension (nsamples, nfeats). Does not need to be
    transposed to calculate projections.

    The call to svd may "fail to converge", which just means dgesdd (a faster
    algorithm) didn't work. In this case, the algorithm tries to decompose the
    transpose. (see
    http://r.789695.n4.nabble.com/Observations-on-SVD-linpack-errors-and-a-workaround-td837282.html)

    """
    from numpy.linalg import svd, LinAlgError
    # center data
    data = spikes[:nspikes] - spikes[:nspikes].mean(0)
    try:
        u, s, v = svd(data, full_matrices=0)
        return v[:nfeats].T.copy()
    except LinAlgError:
        u, s, v = svd(data.T, full_matrices=0)
        return u[:, :nfeats].copy()

def PCA(Y, components):
	"""
	run PCA, retrieving the first (components) principle components
	return [s0, eig, w0]
	s0: factors
	w0: weights
	"""

	N,D = Y.shape
	sv = linalg.svd(Y, full_matrices=0);
	[s0, w0] = [sv[0][:, 0:components], np.dot(np.diag(sv[1]), sv[2]).T[:, 0:components]]
	v = s0.std(axis=0)
	s0 /= v;
	w0 *= v;
	return [s0, w0]

	if N>D:
		sv = linalg.svd(Y, full_matrices=0);
		[s0, w0] = [sv[0][:, 0:components], np.dot(np.diag(sv[1]), sv[2]).T[:, 0:components]]
		v = s0.std(axis=0)
		s0 /= v;
		w0 *= v;
		return [s0, w0]
	else:
		K=np.cov(Y)
		sv = linalg.eigh(K)
		std_var = np.sqrt(sv[0])
		pc = sv[1]*std_var[np.newaxis(),0]
		#
		#ipdb.set_trace()
		return [pc,std_var]

    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 stdNorm(self, U1, U2):
        print "U1"
        print U1
        print "U2"
        print U2

        mat1 = np.matrix(U1).T
        print mat1
        print mat1.mean(axis=1)
        mat1 = mat1 - mat1.mean(axis=1)
        print mat1
        mat1cov = np.cov(mat1)
        print mat1cov
        p1,l1,p1t = NLA.svd(mat1cov)
        print p1
        print l1
        print p1t
        l1sq = SLA.sqrtm(SLA.inv(np.diag(l1))) 
        snU1 =  np.dot(np.dot(l1sq, p1.T), mat1)

        mat2 = np.matrix(U2).T
        mat2 = mat2 - mat2.mean(axis=1)
        mat2cov = np.cov(mat2)
        p2,l2,p2t = NLA.svd(mat2cov)
        l2sq = SLA.sqrtm(SLA.inv(np.diag(l2))) 
        snU2 =  np.dot(np.dot(l2sq, p2.T), mat2)

        print "cov:"
        print np.cov(snU1)
        print np.cov(snU2)

        return snU1, snU2

	def planar_ransac(self, pc):
		#fit to model Ax  +By + Cz + D = 0 (a plane)
		sample_iter = 20
		pc = pc[np.nonzero(np.nansum(pc, axis=1)>0)[0], :]
		pc = pc[::sample_iter,:]
		pc_len = pc.shape[0]
		n=int(0.1*pc_len)		#size of random sample
#		print "Points: ", n
		k=15; 		#number of iteration
		err_thresh=0.01				#deviation - meters
		min_points=.5*pc_len	 	#minime amount of points within deviation		
		
		iter_ = 0
		best_model = None
		best_consensus_set = None
		best_error = np.inf
		best_offset = None
		
		while best_model == None:
			while iter_ < k:

				maybe_inliers = pc[np.random.randint(0, pc_len, (n))] #get n random points from pc
				offset = np.mean(maybe_inliers, axis=0)
				maybe_inliers -= offset
				# Find model
				_,_,Vs = svd(maybe_inliers)
				V = Vs.T.conj() 
				Normal = V[:,2]
				maybe_model = Normal
				consensus_set = maybe_inliers

				err = np.sum((pc-offset) * maybe_model, axis=1)
				consensus_set = pc[np.nonzero(err < err_thresh)[0]]
#				print consensus_set.shape, pc_len
			
				if consensus_set.shape[0] > min_points:
					_,_,Vs = svd(maybe_inliers)
					V = Vs.T.conj() 
					Normal = V[:,2]
					better_model = Normal
				
					offset = np.mean(consensus_set, axis=0)
					new_error = np.sum(np.sum((consensus_set-offset) * maybe_model, axis=1))
#					print new_error
				
					if abs(new_error) < best_error:
	#					print "iter: ", iter_
	#					print "e: ", best_error					
						best_model = better_model
						best_consensus_set = consensus_set
						best_error = abs(new_error)
						best_offset = offset
					
				iter_ += 1
			#if there is no model, decrease the number of points required for 'best' model
			min_points = min_points*3/4 
			iter_ = 0
			
		return best_model, best_offset

 def check(dtype):
     x = np.array([[1, 0.5], [0.5, 1]], dtype=dtype)
     u, s, vh = linalg.svd(x)
     assert_equal(u.dtype, dtype)
     assert_equal(s.dtype, get_real_dtype(dtype))
     assert_equal(vh.dtype, dtype)
     s = linalg.svd(x, compute_uv=False)
     assert_equal(s.dtype, get_real_dtype(dtype))

def nk_singular(a,k, opt):
  (size,_) = a.shape
  nseq = rdft.all_leading_sequence(size)[-1]
  seq  = rdft.all_leading_sequence(k)[-1]
  if opt == 0:
    _, sings, _ = linalg.svd(a[np.ix_(nseq,seq)])
  else:
    _, sings, _ = linalg.svd(a[np.ix_(seq,nseq)])
  return sings

def conceptor_similarity (a, b): 
      U_a, S_a, V_a = svd(a);
      U_b, S_b, V_b = svd(b);
      S_a = diag(S_a)
      S_b = diag(S_b)
      # similarity with previous conceptor
      return pow(norm((sqrt(S_a) * U_a.transpose() * \
              U_b * sqrt(S_b))),2) / \
              (norm(a) * norm(b));

def test_betti1(space):
    """
    Verify that the 1-form Hodge Laplacian with strong Dirichlet
    boundary conditions has kernel of dimension equal to the 1st Betti
    number of the annulus mesh, i.e. 1.
    """
    mesh = Mesh(join(cwd, "annulus.msh"))
    V0tag, V1tag, V2tag = space

    if(len(V0tag) == 2):
        V0 = FunctionSpace(mesh, V0tag[0], V0tag[1])
    else:
        V0a = FiniteElement(V0tag[0], "triangle", V0tag[1])
        V0b = FiniteElement(V0tag[2], "triangle", V0tag[3])
        V0 = FunctionSpace(mesh, V0a + V0b)

    V1 = FunctionSpace(mesh, V1tag[0], V1tag[1])

    W = V0*V1
    sigma, u = TrialFunctions(W)
    tau, v = TestFunctions(W)
    L = assemble((sigma*tau - inner(rot(tau), u) + inner(rot(sigma), v) +
                  div(u)*div(v))*dx)

    bc0 = DirichletBC(W.sub(0), 0., 9)
    bc1 = DirichletBC(W.sub(1), Expression(("0.0", "0.0")), 9)
    L0 = assemble((sigma*tau - inner(rot(tau), u) + inner(rot(sigma), v) +
                   div(u)*div(v))*dx, bcs=[bc0, bc1])

    dV0 = V0.dof_count
    dV1 = V1.dof_count

    A = numpy.zeros((dV0+dV1, dV0+dV1))
    A[:dV0, :dV0] = L.M[0, 0].values
    A[:dV0, dV0:dV0+dV1] = L.M[0, 1].values
    A[dV0:dV0+dV1, :dV0] = L.M[1, 0].values
    A[dV0:dV0+dV1, dV0:dV0+dV1] = L.M[1, 1].values

    u, s, v = linalg.svd(A)

    nharmonic = sum(s < 1.0e-5)
    assert(nharmonic == 1)

    dV0 = V0.dof_count
    dV1 = V1.dof_count

    A0 = numpy.zeros((dV0+dV1, dV0+dV1))
    A0[:dV0, :dV0] = L0.M[0, 0].values
    A0[:dV0, dV0:dV0+dV1] = L0.M[0, 1].values
    A0[dV0:dV0+dV1, :dV0] = L0.M[1, 0].values
    A0[dV0:dV0+dV1, dV0:dV0+dV1] = L0.M[1, 1].values

    u, s, v = linalg.svd(A0)

    nharmonic = sum(s < 1.0e-5)
    assert(nharmonic == 1)

def get_singular(a,fra):
  (size,_) = a.shape
  a_fra_subsings = []
  _, a_sings, _ = linalg.svd(a)
  i = 0
  for seq in rdft.all_leading_sequence(size):
    _, fra_subsing, _ = linalg.svd(fra[np.ix_(seq,seq)])
    a_fra_subsings.append( (a_sings[i], fra_subsing[-1]) )
    i = i + 1
  return a_fra_subsings

def incremental_svd(A, qr_flg=False):
  """A matrix "A" should be 256x7291
  """

  m = 256
  n = 7291

  n0 = 256

  if A.shape[0] != m or A.shape[1] != n: raise ValueError('Error: incorrect matrix size')

  start = time.clock()

  A0 = A[:, :n0]
  U, s, V = ln.svd(A0, full_matrices=False)

  # NOTE: s is a vector; np.diag(s) will produce a diagonal matrix
  for i in range(n0, n):

    # new matrix is just a single vector (i-th column of A)
    A1 = np.matrix(A[:, i]).T

    if qr_flg:
      J, K = ln.qr(A1 - np.dot(np.dot(U, U.T), A1))
      U_, s_, V_ = ln.svd(
          np.vstack((
            np.hstack((np.diag(s), np.dot(U.T, A1))),
            np.hstack((np.zeros((K.shape[0], s.shape[0])), K))
          )),
          full_matrices=False)

      # update the result of SVD
      U = np.dot(np.hstack((U, J)), U_)

    else:
      U_, s_, V_ = ln.svd(np.hstack((np.diag(s), np.dot(U.T, A1))), full_matrices=False)
      U = np.dot(U, U_)

    s = s_

    # NOTE: V from svd on NumPy is already transposed
    V = np.dot(V_,
          np.vstack((
            np.hstack((V, np.zeros((V.shape[0], i+1-V.shape[1])))),
            np.hstack((np.zeros((V_.shape[1]-V.shape[0], V.shape[1])), np.eye(V_.shape[1]-V.shape[0], i+1-V.shape[1])))
          ))
        )

    # for next computation, update A0
    A0 = np.hstack((A0, A1))

  elapsed_time = time.clock() - start
  print 'time:', elapsed_time

  return U, s, V

def PCA(X):
    # tansform X to be centered data
    # cauz the feature dimention is larger than data number, so we should calculate the SVD of X's transpose
    V,D,UT = svd(X.T)
    U1,D1,VT1 = svd(X)
    U = array(matrix(UT).T)
    results = []
    results.append(U[:,0])
    results.append(U[:,1])
    results = dot(results,X) 
    return results

 def G(self):
     if self.pos in (0,len(self)):
         B=self[self.pos] if self.pos==0 else self[self.pos-1]
         U,S,V=svd(B.V.dot(B.U).transpose().conjugate()+np.diag(B.S))
         U,V=B.U.dot(U).transpose().conjugate(),V.dot(B.V).transpose().conjugate()
         return np.einsum('ij,j,jk->ik',V,1.0/S,U)
     else:
         BR,BL=self[self.pos-1],self[self.pos]
         U,S,V=svd(BL.V.dot(BR.U).transpose().conjugate()+np.einsum('i,ij,jk,k->ik',BR.S,BR.V,BL.U,BL.S))
         U,V=BR.U.dot(U).transpose().conjugate(),V.dot(BL.V).transpose().conjugate()
         return np.einsum('ij,j,jk->ik',V,1.0/S,U)