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;
  }
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. ) ; } LOG . incrementProcessed ( prog ) ; } LOG . ensureCompleted ( prog ) ; if ( LOG . isStatistics ( ) ) { 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 ( ) ) ) ; } double sum = 0. ; for ( int i = 1 ; i < size ; i ++ ) { final double [ ] pij_i = pij [ i ] ; for ( int j = 0 ; j < i ; j ++ ) { sum += ( pij_i [ j ] += pij [ j ] [ i ] ) ; } } 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 ; }
Compute the pij from the distance matrix.

@param dist Distance matrix.
@param perplexity Desired perplexity
@param initialScale Initial scale
@return Affinity matrix pij