docs: protected void estimateDensities(Relation<O> rel, KNNQuery<O> knnq, final DBIDs ids, WritableDataStore<double[]> densities) {
      final int dim = dimensionality(rel);
      final int knum = kmax + 1 - kmin;
      // Initialize storage:
      for(DBIDIter iter = ids.iter(); iter.valid(); iter.advance()) {
        densities.put(iter, new double[knum]);
      }
      // Distribute densities:
      FiniteProgress prog = LOG.isVerbose() ? new FiniteProgress("Computing densities", ids.size(), LOG) : null;
      double iminbw = (minBandwidth > 0.) ? 1. / (minBandwidth * scale) : Double.POSITIVE_INFINITY;
      for(DBIDIter iter = ids.iter(); iter.valid(); iter.advance()) {
        KNNList neighbors = knnq.getKNNForDBID(iter, kmax + 1);
        int k = 1, idx = 0;
        double sum = 0.;
        for(DoubleDBIDListIter kneighbor = neighbors.iter(); k <= kmax && kneighbor.valid(); kneighbor.advance(), k++) {
          sum += kneighbor.doubleValue();
          if(k < kmin) {
            continue;
          }
          final double ibw = Math.min(k / (sum * scale), iminbw);
          final double sca = MathUtil.powi(ibw, dim);
          for(DoubleDBIDListIter neighbor = neighbors.iter(); neighbor.valid(); neighbor.advance()) {
            final double dens;
            if(sca < Double.POSITIVE_INFINITY) { // NaNs with duplicate points!
              dens = sca * kernel.density(neighbor.doubleValue() * ibw);
            }
            else {
              dens = neighbor.doubleValue() == 0. ? 1. : 0.;
            }
            densities.get(neighbor)[idx] += dens;
            if(dens < CUTOFF) {
              break;
            }
          }
          ++idx; // Only if k >= kmin
        }
        LOG.incrementProcessed(prog);
      }
      LOG.ensureCompleted(prog);
    }