cluster

  auto dist = [n, &dists](size_t i, size_t j) { return dists[i * n + j]; } // dists is pre-calculated n x n distances matrix
  auto index = seq(n); //[0,...,n-1]

  std::vector<size_t> rho(n);
  for (size_t i=0; i<n; ++i) {
      rho[i] = std::count_if(index.begin(), index.end(), [&](size_t j) { return dist(i,j) < dist_cutoff; });
  }

  std::vector<D> delta(n);
  for (size_t i=0; i<n; ++i) {
      auto it = std::min_element_if(index.begin(), index.end(), [&](size_t j, size_t k) { return dist(i,j) < dist(i,k); }, [&](size_t j) { return rho[j] > rho[i]; });
      if (it == index.end())
          it = std::max_element(index.begin(), index.end(), [&](size_t j, size_t k) { return dist(i,j) < dist(i,k); });
      delta[i] = dist(i, *it);
  }