docs: protected static double[][] computePij(double[][] dist, double perplexity, double initialScale) {
      final int size = dist.length;
      final double logPerp = FastMath.log(perplexity);
      double[][] pij = new double[size][size];
      FiniteProgress prog = LOG.isVerbose() ? new FiniteProgress("Optimizing perplexities", size, LOG) : null;
      Duration timer = LOG.isStatistics() ? LOG.newDuration(PerplexityAffinityMatrixBuilder.class.getName() + ".runtime.pijmatrix").begin() : null;
      MeanVariance mv = LOG.isStatistics() ? new MeanVariance() : null;
      for(int i = 0; i < size; i++) {
        double beta = computePi(i, dist[i], pij[i], perplexity, logPerp);
        if(mv != null) {
          mv.put(beta > 0 ? FastMath.sqrt(.5 / beta) : 0.); // Sigma
        }
        LOG.incrementProcessed(prog);
      }
      LOG.ensureCompleted(prog);
      if(LOG.isStatistics()) { // timer != null, mv != null
        LOG.statistics(timer.end());
        LOG.statistics(new DoubleStatistic(PerplexityAffinityMatrixBuilder.class.getName() + ".sigma.average", mv.getMean()));
        LOG.statistics(new DoubleStatistic(PerplexityAffinityMatrixBuilder.class.getName() + ".sigma.stddev", mv.getSampleStddev()));
      }
      // Scale pij to have the desired sum EARLY_EXAGGERATION
      double sum = 0.;
      for(int i = 1; i < size; i++) {
        final double[] pij_i = pij[i];
        for(int j = 0; j < i; j++) { // Nur über halbe Matrix!
          sum += (pij_i[j] += pij[j][i]); // Symmetrie herstellen
        }
      }
      // Scaling taken from original tSNE code:
      final double scale = initialScale / (2. * sum);
      for(int i = 1; i < size; i++) {
        final double[] pij_i = pij[i];
        for(int j = 0; j < i; j++) {
          pij_i[j] = pij[j][i] = MathUtil.max(pij_i[j] * scale, MIN_PIJ);
        }
      }
      return pij;
    }