def saffron(X, q=32, k=4, tangent_dim=1, curv_thresh=0.95, decay_rate=0.9,
            max_iter=15, verbose=False):
  '''
  SAFFRON graph construction method.

    X : (n,d)-array of coordinates
    q : int, median number of candidate friends per vertex
    k : int, number of friends to select per vertex, k < q
    tangent_dim : int, dimensionality of manifold tangent space
    curv_thresh : float, tolerance to curvature, lambda in the paper
    decay_rate : float, controls step size per iteration, between 0 and 1
    max_iter : int, cap on number of iterations
    verbose : bool, print goodness measure per iteration when True

  From "Tangent Space Guided Intelligent Neighbor Finding",
    by Gashler & Martinez, 2011.
  See http://axon.cs.byu.edu/papers/gashler2011ijcnn1.pdf
  '''
  n = len(X)
  dist = pairwise_distances(X)
  idx = np.argpartition(dist, q)[:, q]
  # radius for finding candidate friends: median distance to qth neighbor
  r = np.median(dist[np.arange(n), idx])

  # make candidate graph + weights
  W = neighbor_graph(dist, precomputed=True, epsilon=r).matrix('csr')
  # NOTE: this differs from the paper, where W.data[:] = 1 initially
  W.data[:] = 1 / W.data
  # row normalize
  normalize(W, norm='l1', axis=1, copy=False)
  # XXX: hacky densify
  W = W.toarray()

  # iterate to learn optimal weights
  prev_goodness = 1e-12
  for it in range(max_iter):
    goodness = 0
    S = _estimate_tangent_spaces(X, W, tangent_dim)
    # find aligned candidates
    for i, row in enumerate(W):
      nbrs = row.nonzero()[-1]

      # compute alignment scores
      edges = X[nbrs] - X[i]
      edge_norms = (edges**2).sum(axis=1)
      a1 = (edges.dot(S[i])**2).sum(axis=1) / edge_norms
      a2 = (np.einsum('ij,ijk->ik', edges, S[nbrs])**2).sum(axis=1) / edge_norms
      a3 = _principal_angle(S[i], S[nbrs]) ** 2
      x = (np.minimum(curv_thresh, a1) *
           np.minimum(curv_thresh, a2) *
           np.minimum(curv_thresh, a3))

      # decay weight of least-aligned candidates
      excess = x.shape[0] - k
      if excess > 0:
        bad_idx = np.argpartition(x, excess-1)[:excess]
        W[i, nbrs[bad_idx]] *= decay_rate
        W[i] /= W[i].sum()

      # update goodness measure (weighted alignment)
      goodness += x.dot(W[i,nbrs])

    if verbose:  # pragma: no cover
      goodness /= n
      print(it, goodness, goodness / prev_goodness)
    if goodness / prev_goodness <= 1.0001:
      break
    prev_goodness = goodness
  else:
    warnings.warn('Failed to converge after %d iterations.' % max_iter)

  # use the largest k weights for each row of W, weighted by original distance
  indptr, indices, data = [0], [], []
  for i, row in enumerate(W):
    nbrs = row.nonzero()[-1]
    if len(nbrs) > k:
      nbrs = nbrs[np.argpartition(row[nbrs], len(nbrs)-k)[-k:]]
    indices.extend(nbrs)
    indptr.append(len(nbrs))
    data.extend(dist[i, nbrs])
  indptr = np.cumsum(indptr)
  data = np.array(data)
  indices = np.array(indices)
  W = ss.csr_matrix((data, indices, indptr), shape=W.shape)
  return Graph.from_adj_matrix(W)
def saffron ( X , q = 32 , k = 4 , tangent_dim = 1 , curv_thresh = 0.95 , decay_rate = 0.9 , max_iter = 15 , verbose = False ) : n = len ( X ) dist = pairwise_distances ( X ) idx = np . argpartition ( dist , q ) [ : , q ] r = np . median ( dist [ np . arange ( n ) , idx ] ) W = neighbor_graph ( dist , precomputed = True , epsilon = r ) . matrix ( 'csr' ) W . data [ : ] = 1 / W . data normalize ( W , norm = 'l1' , axis = 1 , copy = False ) W = W . toarray ( ) prev_goodness = 1e-12 for it in range ( max_iter ) : goodness = 0 S = _estimate_tangent_spaces ( X , W , tangent_dim ) for i , row in enumerate ( W ) : nbrs = row . nonzero ( ) [ - 1 ] edges = X [ nbrs ] - X [ i ] edge_norms = ( edges ** 2 ) . sum ( axis = 1 ) a1 = ( edges . dot ( S [ i ] ) ** 2 ) . sum ( axis = 1 ) / edge_norms a2 = ( np . einsum ( 'ij,ijk->ik' , edges , S [ nbrs ] ) ** 2 ) . sum ( axis = 1 ) / edge_norms a3 = _principal_angle ( S [ i ] , S [ nbrs ] ) ** 2 x = ( np . minimum ( curv_thresh , a1 ) * np . minimum ( curv_thresh , a2 ) * np . minimum ( curv_thresh , a3 ) ) excess = x . shape [ 0 ] - k if excess > 0 : bad_idx = np . argpartition ( x , excess - 1 ) [ : excess ] W [ i , nbrs [ bad_idx ] ] *= decay_rate W [ i ] /= W [ i ] . sum ( ) goodness += x . dot ( W [ i , nbrs ] ) if verbose : goodness /= n print ( it , goodness , goodness / prev_goodness ) if goodness / prev_goodness <= 1.0001 : break prev_goodness = goodness else : warnings . warn ( 'Failed to converge after %d iterations.' % max_iter ) indptr , indices , data = [ 0 ] , [ ] , [ ] for i , row in enumerate ( W ) : nbrs = row . nonzero ( ) [ - 1 ] if len ( nbrs ) > k : nbrs = nbrs [ np . argpartition ( row [ nbrs ] , len ( nbrs ) - k ) [ - k : ] ] indices . extend ( nbrs ) indptr . append ( len ( nbrs ) ) data . extend ( dist [ i , nbrs ] ) indptr = np . cumsum ( indptr ) data = np . array ( data ) indices = np . array ( indices ) W = ss . csr_matrix ( ( data , indices , indptr ) , shape = W . shape ) return Graph . from_adj_matrix ( W )
SAFFRON graph construction method.

    X : (n,d)-array of coordinates
    q : int, median number of candidate friends per vertex
    k : int, number of friends to select per vertex, k < q
    tangent_dim : int, dimensionality of manifold tangent space
    curv_thresh : float, tolerance to curvature, lambda in the paper
    decay_rate : float, controls step size per iteration, between 0 and 1
    max_iter : int, cap on number of iterations
    verbose : bool, print goodness measure per iteration when True

  From "Tangent Space Guided Intelligent Neighbor Finding",
    by Gashler & Martinez, 2011.
  See http://axon.cs.byu.edu/papers/gashler2011ijcnn1.pdf