docs: protected static double computeH(final int i, double[] dist_i, double[] pij_i, double mbeta) {
      double sumP = 0.;
      // Skip point "i", break loop in two:
      for(int j = 0; j < i; j++) {
        sumP += (pij_i[j] = FastMath.exp(dist_i[j] * mbeta));
      }
      for(int j = i + 1; j < dist_i.length; j++) {
        sumP += (pij_i[j] = FastMath.exp(dist_i[j] * mbeta));
      }
      if(!(sumP > 0)) {
        // All pij are zero. Bad news.
        return Double.NEGATIVE_INFINITY;
      }
      final double s = 1. / sumP; // Scaling factor
      double sum = 0.;
      // While we could skip pi[i], it should be 0 anyway.
      for(int j = 0; j < dist_i.length; j++) {
        sum += dist_i[j] * (pij_i[j] *= s);
      }
      return FastMath.log(sumP) - mbeta * sum;
    }