docs: protected void computePij(DBIDRange ids, KNNQuery<?> knnq, boolean square, int numberOfNeighbours, double[][] pij, int[][] indices, double initialScale) {
      Duration timer = LOG.isStatistics() ? LOG.newDuration(this.getClass().getName() + ".runtime.neighborspijmatrix").begin() : null;
      final double logPerp = FastMath.log(perplexity);
      // Scratch arrays, resizable
      DoubleArray dists = new DoubleArray(numberOfNeighbours + 10);
      IntegerArray inds = new IntegerArray(numberOfNeighbours + 10);
      // Compute nearest-neighbor sparse affinity matrix
      FiniteProgress prog = LOG.isVerbose() ? new FiniteProgress("Finding neighbors and optimizing perplexity", ids.size(), LOG) : null;
      MeanVariance mv = LOG.isStatistics() ? new MeanVariance() : null;
      Mean mid = LOG.isStatistics() ? new Mean() : null;
      for(DBIDArrayIter ix = ids.iter(); ix.valid(); ix.advance()) {
        dists.clear();
        inds.clear();
        KNNList neighbours = knnq.getKNNForDBID(ix, numberOfNeighbours + 1);
        convertNeighbors(ids, ix, square, neighbours, dists, inds, mid);
        double beta = computeSigma(ix.getOffset(), dists, perplexity, logPerp, //
            pij[ix.getOffset()] = new double[dists.size()]);
        if(mv != null) {
          mv.put(beta > 0 ? FastMath.sqrt(.5 / beta) : 0.); // Sigma
        }
        indices[ix.getOffset()] = inds.toArray();
        LOG.incrementProcessed(prog);
      }
      LOG.ensureCompleted(prog);
      if(mid != null) {
        LOG.statistics(new DoubleStatistic(getClass() + ".average-original-id", mid.getMean()));
      }
      // Sum of the sparse affinity matrix:
      double sum = 0.;
      for(int i = 0; i < pij.length; i++) {
        final double[] pij_i = pij[i];
        for(int offi = 0; offi < pij_i.length; offi++) {
          int j = indices[i][offi];
          if(j > i) {
            continue; // Exploit symmetry.
          }
          assert (i != j);
          int offj = containsIndex(indices[j], i);
          if(offj >= 0) { // Found
            sum += FastMath.sqrt(pij_i[offi] * pij[j][offj]);
          }
        }
      }
      final double scale = initialScale / (2 * sum);
      for(int i = 0; i < pij.length; i++) {
        final double[] pij_i = pij[i];
        for(int offi = 0; offi < pij_i.length; offi++) {
          int j = indices[i][offi];
          assert (i != j);
          int offj = containsIndex(indices[j], i);
          if(offj >= 0) { // Found
            assert (indices[j][offj] == i);
            // Exploit symmetry:
            if(i < j) {
              final double val = FastMath.sqrt(pij_i[offi] * pij[j][offj]); // Symmetrize
              pij_i[offi] = pij[j][offj] = MathUtil.max(val * scale, MIN_PIJ);
            }
          }
          else { // Not found, so zero.
            pij_i[offi] = 0;
          }
        }
      }
      if(LOG.isStatistics()) { // timer != null, mv != null
        LOG.statistics(timer.end());
        LOG.statistics(new DoubleStatistic(NearestNeighborAffinityMatrixBuilder.class.getName() + ".sigma.average", mv.getMean()));
        LOG.statistics(new DoubleStatistic(NearestNeighborAffinityMatrixBuilder.class.getName() + ".sigma.stddev", mv.getSampleStddev()));
      }
    }