docs: private double[] findAndEvaluateThreshold(DoubleDynamicHistogram histogram) {
      int n = histogram.getNumBins();
      double[] p1 = new double[n];
      double[] p2 = new double[n];
      double[] mu1 = new double[n];
      double[] mu2 = new double[n];
      double[] sigma1 = new double[n];
      double[] sigma2 = new double[n];
      double[] jt = new double[n];
      // Forward pass
      {
        MeanVariance mv = new MeanVariance();
        DoubleHistogram.Iter forward = histogram.iter();
        for(int i = 0; forward.valid(); i++, forward.advance()) {
          p1[i] = forward.getValue() + ((i > 0) ? p1[i - 1] : 0);
          mv.put(i, forward.getValue());
          mu1[i] = mv.getMean();
          sigma1[i] = mv.getNaiveStddev();
        }
      }
      // Backwards pass
      {
        MeanVariance mv = new MeanVariance();
        DoubleHistogram.Iter backwards = histogram.iter();
        backwards.seek(histogram.getNumBins() - 1); // Seek to last
  
        for(int j = n - 1; backwards.valid(); j--, backwards.retract()) {
          p2[j] = backwards.getValue() + ((j + 1 < n) ? p2[j + 1] : 0);
          mv.put(j, backwards.getValue());
          mu2[j] = mv.getMean();
          sigma2[j] = mv.getNaiveStddev();
        }
      }
  
      for(int i = 0; i < n; i++) {
        jt[i] = 1.0 + 2 * (p1[i] * (FastMath.log(sigma1[i]) - FastMath.log(p1[i])) + p2[i] * (FastMath.log(sigma2[i]) - FastMath.log(p2[i])));
      }
  
      int bestpos = -1;
      double bestgoodness = Double.NEGATIVE_INFINITY;
  
      double devPrev = jt[1] - jt[0];
      for(int i = 1; i < jt.length - 1; i++) {
        double devCur = jt[i + 1] - jt[i];
        // Local minimum found - calculate depth
        if(devCur >= 0 && devPrev <= 0) {
          double lowestMaxima = Double.POSITIVE_INFINITY;
          for(int j = i - 1; j > 0; j--) {
            if(jt[j - 1] < jt[j]) {
              lowestMaxima = Math.min(lowestMaxima, jt[j]);
              break;
            }
          }
          for(int j = i + 1; j < n - 2; j++) {
            if(jt[j + 1] < jt[j]) {
              lowestMaxima = Math.min(lowestMaxima, jt[j]);
              break;
            }
          }
          double localDepth = lowestMaxima - jt[i];
  
          final double mud = mu1[i] - mu2[i];
          double discriminability = mud * mud / (sigma1[i] * sigma1[i] + sigma2[i] * sigma2[i]);
          if(Double.isNaN(discriminability)) {
            discriminability = -1;
          }
          double goodness = localDepth * discriminability;
          if(goodness > bestgoodness) {
            bestgoodness = goodness;
            bestpos = i;
          }
        }
        devPrev = devCur;
      }
      Iter iter = histogram.iter();
      iter.seek(bestpos);
      return new double[] { iter.getRight(), bestgoodness };
    }