docs: def barycenter_graph(distance_matrix, X, reg=1e-3):
      """
      Computes the barycenter weighted graph for points in X
  
      Parameters
      ----------
      distance_matrix: sparse Ndarray, (N_obs, N_obs) pairwise distance matrix.
      X : Ndarray (N_obs, N_dim) observed data matrix.
      reg : float, optional
          Amount of regularization when solving the least-squares
          problem. Only relevant if mode='barycenter'. If None, use the
          default.
  
      Returns
      -------
      W : sparse matrix in CSR format, shape = [n_samples, n_samples]
          W[i, j] is assigned the weight of edge that connects i to j.
      """
      (N, d_in) = X.shape
      (rows, cols) = distance_matrix.nonzero()
      W = sparse.lil_matrix((N, N)) # best for W[i, nbrs_i] = w/np.sum(w)
      for i in range(N):
          nbrs_i = cols[rows == i]
          n_neighbors_i = len(nbrs_i)
          v = np.ones(n_neighbors_i, dtype=X.dtype)
          C = X[nbrs_i] - X[i]
          G = np.dot(C, C.T)
          trace = np.trace(G)
          if trace > 0:
              R = reg * trace
          else:
              R = reg
          G.flat[::n_neighbors_i + 1] += R
          w = solve(G, v, sym_pos = True)
          W[i, nbrs_i] = w / np.sum(w)
      return W