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

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)