docs: public Clustering<MeanModel> run(Database database, Relation<V> relation) {
      final DistanceQuery<V> distq = database.getDistanceQuery(relation, getDistanceFunction());
      final RangeQuery<V> rangeq = database.getRangeQuery(distq);
      final NumberVector.Factory<V> factory = RelationUtil.getNumberVectorFactory(relation);
      final int dim = RelationUtil.dimensionality(relation);
  
      // Stopping threshold
      final double threshold = bandwidth * 1E-10;
  
      // Result store:
      ArrayList<Pair<V, ModifiableDBIDs>> clusters = new ArrayList<>();
  
      ModifiableDBIDs noise = DBIDUtil.newArray();
  
      FiniteProgress prog = LOG.isVerbose() ? new FiniteProgress("Mean-shift clustering", relation.size(), LOG) : null;
  
      for(DBIDIter iter = relation.iterDBIDs(); iter.valid(); iter.advance()) {
        // Initial position:
        V position = relation.get(iter);
        iterations: for(int j = 1; ; j++) {
          // Compute new position:
          V newvec = null;
          {
            DoubleDBIDList neigh = rangeq.getRangeForObject(position, bandwidth);
            boolean okay = (neigh.size() > 1) || (neigh.size() >= 1 && j > 1);
            if(okay) {
              Centroid newpos = new Centroid(dim);
              for(DoubleDBIDListIter niter = neigh.iter(); niter.valid(); niter.advance()) {
                final double weight = kernel.density(niter.doubleValue() / bandwidth);
                newpos.put(relation.get(niter), weight);
              }
              newvec = factory.newNumberVector(newpos.getArrayRef());
              // TODO: detect 0 weight!
            }
            if(!okay) {
              noise.add(iter);
              break iterations;
            }
          }
          // Test if we are close to one of the known clusters:
          double bestd = Double.POSITIVE_INFINITY;
          Pair<V, ModifiableDBIDs> bestp = null;
          for(Pair<V, ModifiableDBIDs> pair : clusters) {
            final double merged = distq.distance(newvec, pair.first);
            if(merged < bestd) {
              bestd = merged;
              bestp = pair;
            }
          }
          // Check for convergence:
          double delta = distq.distance(position, newvec);
          if(bestd < 10 * threshold || bestd * 2 < delta) {
            assert(bestp != null);
            bestp.second.add(iter);
            break iterations;
          }
          if(Double.isNaN(delta)) {
            LOG.warning("Encountered NaN distance. Invalid center vector? " + newvec.toString());
            break iterations;
          }
          if(j == MAXITER || delta < threshold) {
            if(j == MAXITER) {
              LOG.warning("No convergence after " + MAXITER + " iterations. Distance: " + delta);
            }
            if(LOG.isDebuggingFine()) {
              LOG.debugFine("New cluster:" + newvec + " delta: " + delta + " threshold: " + threshold + " bestd: " + bestd);
            }
            ArrayModifiableDBIDs cids = DBIDUtil.newArray(1);
            cids.add(iter);
            clusters.add(new Pair<V, ModifiableDBIDs>(newvec, cids));
            break iterations;
          }
          position = newvec;
        }
        LOG.incrementProcessed(prog);
      }
      LOG.ensureCompleted(prog);
  
      ArrayList<Cluster<MeanModel>> cs = new ArrayList<>(clusters.size());
      for(Pair<V, ModifiableDBIDs> pair : clusters) {
        cs.add(new Cluster<>(pair.second, new MeanModel(pair.first.toArray())));
      }
      if(noise.size() > 0) {
        cs.add(new Cluster<MeanModel>(noise, true));
      }
      return new Clustering<>("Mean-shift Clustering", "mean-shift-clustering", cs);
    }