[2] We use k nearest neighbor interpolation to get the bias for a given raw estimate
The choice of k = 6 is rather arbitrary. The best value of k could be determined experimentally,
but we found that the choice has only a minuscule influence.
762940
private long estimateBias(double e) {
        int i = 0;
        double[] rawEstimates = DenseHyperLogLogConstants.RAW_ESTIMATE_DATA[p - 4];
        double closestToZero = Math.abs(e - rawEstimates[0]);
        NavigableMap<Double, Integer> distances = new TreeMap<Double, Integer>();
        for (double est : rawEstimates) {
            double distance = e - est;
            distances.put(distance, i++);
            if (Math.abs(distance) < closestToZero) {
                closestToZero = distance;
            }
        }

        // abomination to compute kNN elements (we could ideally use a tree structure)
        int kNN = 6;
        double sum = 0;
        Iterator<Map.Entry<Double, Integer>> firstX = distances.descendingMap().tailMap(closestToZero).entrySet().iterator();
        Iterator<Map.Entry<Double, Integer>> lastX = distances.tailMap(closestToZero).entrySet().iterator();

        int kNNLeft = kNN;
        while (kNNLeft-- > kNN / 2 && firstX.hasNext()) {
            sum += DenseHyperLogLogConstants.BIAS_DATA[p - 4][firstX.next().getValue()];
        }

        while (kNNLeft-- >= 0 && lastX.hasNext()) {
            sum += DenseHyperLogLogConstants.BIAS_DATA[p - 4][lastX.next().getValue()];
        }

        return (long) (sum / kNN);
    }
1066294
def _get_interpolator_nd(method, **kwargs):
    '''
    '''
    valid_methods = ['linear', 'nearest']

    try:
        from scipy import interpolate
    except ImportError:
        raise ImportError(
            'Interpolation with method `%s` requires scipy' % method)

    if method in valid_methods:
        kwargs.update(method=method)
        interp_class = interpolate.interpn
    else:
        raise ValueError('%s is not a valid interpolator for interpolating '
                         'over multiple dimensions.' % method)

    return interp_class, kwargs
225817
@Override
    public double inverseCumulativeProbability(final double p) throws OutOfRangeException {
        /*
         * IMPLEMENTATION NOTES
         * --------------------
         * Where applicable, use is made of the one-sided Chebyshev inequality
         * to bracket the root. This inequality states that
         * P(X - mu >= k * sig) <= 1 / (1 + k^2),
         * mu: mean, sig: standard deviation. Equivalently
         * 1 - P(X < mu + k * sig) <= 1 / (1 + k^2),
         * F(mu + k * sig) >= k^2 / (1 + k^2).
         *
         * For k = sqrt(p / (1 - p)), we find
         * F(mu + k * sig) >= p,
         * and (mu + k * sig) is an upper-bound for the root.
         *
         * Then, introducing Y = -X, mean(Y) = -mu, sd(Y) = sig, and
         * P(Y >= -mu + k * sig) <= 1 / (1 + k^2),
         * P(-X >= -mu + k * sig) <= 1 / (1 + k^2),
         * P(X <= mu - k * sig) <= 1 / (1 + k^2),
         * F(mu - k * sig) <= 1 / (1 + k^2).
         *
         * For k = sqrt((1 - p) / p), we find
         * F(mu - k * sig) <= p,
         * and (mu - k * sig) is a lower-bound for the root.
         *
         * In cases where the Chebyshev inequality does not apply, geometric
         * progressions 1, 2, 4, ... and -1, -2, -4, ... are used to bracket
         * the root.
         */
        if (p < 0.0 || p > 1.0) {
            throw new OutOfRangeException(p, 0, 1);
        }

        double lowerBound = getSupportLowerBound();
        if (p == 0.0) {
            return lowerBound;
        }

        double upperBound = getSupportUpperBound();
        if (p == 1.0) {
            return upperBound;
        }

        final double mu = getNumericalMean();
        final double sig = FastMath.sqrt(getNumericalVariance());
        final boolean chebyshevApplies;
        chebyshevApplies = !(Double.isInfinite(mu) || Double.isNaN(mu) || Double.isInfinite(sig) || Double.isNaN(sig));

        if (lowerBound == Double.NEGATIVE_INFINITY) {
            if (chebyshevApplies) {
                lowerBound = mu - sig * FastMath.sqrt((1. - p) / p);
            } else {
                lowerBound = -1.0;
                while (cumulativeProbability(lowerBound) >= p) {
                    lowerBound *= 2.0;
                }
            }
        }

        if (upperBound == Double.POSITIVE_INFINITY) {
            if (chebyshevApplies) {
                upperBound = mu + sig * FastMath.sqrt(p / (1. - p));
            } else {
                upperBound = 1.0;
                while (cumulativeProbability(upperBound) < p) {
                    upperBound *= 2.0;
                }
            }
        }

        final UnivariateFunction toSolve = new UnivariateFunction() {

            public double value(final double x) {
                return cumulativeProbability(x) - p;
            }
        };

        double x = UnivariateSolverUtils.solve(toSolve, lowerBound, upperBound, getSolverAbsoluteAccuracy());

        if (!isSupportConnected()) {
            /* Test for plateau. */
            final double dx = getSolverAbsoluteAccuracy();
            if (x - dx >= getSupportLowerBound()) {
                double px = cumulativeProbability(x);
                if (cumulativeProbability(x - dx) == px) {
                    upperBound = x;
                    while (upperBound - lowerBound > dx) {
                        final double midPoint = 0.5 * (lowerBound + upperBound);
                        if (cumulativeProbability(midPoint) < px) {
                            lowerBound = midPoint;
                        } else {
                            upperBound = midPoint;
                        }
                    }
                    return upperBound;
                }
            }
        }
        return x;
    }
968974
def _K(self, R):
        """"""
        return self._ns * self._N / R / self._sin_alpha
253038
private boolean adjustParam(DenseMatrix64F X, DenseMatrix64F Y,
                                double prevCost) {
        // lambda adjusts how big of a step it takes
        double lambda = initialLambda;
        // the difference between the current and previous cost
        double difference = 1000;

        for( int iter = 0; iter < iter1 && difference > maxDifference ; iter++ ) {
            // compute some variables based on the gradient
            computeDandH(param,X,Y);

            // try various step sizes and see if any of them improve the
            // results over what has already been done
            boolean foundBetter = false;
            for( int i = 0; i < iter2; i++ ) {
                computeA(A,H,lambda);

                if( !solve(A,d,negDelta) ) {
                    return false;
                }
                // compute the candidate parameters
                subtract(param, negDelta, tempParam);

                double cost = cost(tempParam,X,Y);
                if( cost < prevCost ) {
                    // the candidate parameters produced better results so use it
                    foundBetter = true;
                    param.set(tempParam);
                    difference = prevCost - cost;
                    prevCost = cost;
                    lambda /= 10.0;
                } else {
                    lambda *= 10.0;
                }
            }

            // it reached a point where it can't improve so exit
            if( !foundBetter )
                break;
        }
        finalCost = prevCost;
        return true;
    }
1106506
def nonconformity(self, knns_labels):
    
    nb_data = knns_labels[self.layers[0]].shape[0]
    knns_not_in_class = np.zeros((nb_data, self.nb_classes), dtype=np.int32)
    for i in range(nb_data):
      # Compute number of nearest neighbors per class
      knns_in_class = np.zeros(
          (len(self.layers), self.nb_classes), dtype=np.int32)
      for layer_id, layer in enumerate(self.layers):
        knns_in_class[layer_id, :] = np.bincount(
            knns_labels[layer][i], minlength=self.nb_classes)

      # Compute number of knns in other class than class_id
      for class_id in range(self.nb_classes):
        knns_not_in_class[i, class_id] = np.sum(
            knns_in_class) - np.sum(knns_in_class[:, class_id])
    return knns_not_in_class
908668
def K(self):
        
        K = np.eye(3, dtype=np.float64)
        K[0, 0] = self.focal[0]
        K[1, 1] = self.focal[1]
        K[0, 2] = self.resolution[0] / 2.0
        K[1, 2] = self.resolution[1] / 2.0
        return K
425764
def _determine_nTrackIterations(self,nTrackIterations):
        """"""
        if not nTrackIterations is None:
            self.nTrackIterations= nTrackIterations
            return None
        if numpy.fabs(self.misalignment(quantity=False)) < 1./180.*numpy.pi:
            self.nTrackIterations= 0
        elif numpy.fabs(self.misalignment(quantity=False)) >= 1./180.*numpy.pi \
                and numpy.fabs(self.misalignment(quantity=False)) < 3./180.*numpy.pi:
            self.nTrackIterations= 1
        elif numpy.fabs(self.misalignment(quantity=False)) >= 3./180.*numpy.pi:
            self.nTrackIterations= 2
        return None
1337167
def median_kneighbour_distance(X, k=5):
    
    N_all = X.shape[0]
    k = min(k, N_all)
    N_subset = min(N_all, 2000)
    sample_idx_train = np.random.permutation(N_all)[:N_subset]
    nn = neighbors.NearestNeighbors(k)
    nn.fit(X[sample_idx_train, :])
    d, idx = nn.kneighbors(X[sample_idx_train, :])
    return np.median(d[:, -1])
286887
def infer_x(self, y, sigma=None, k=None, **kwargs):
        
        assert len(y) == self.fmodel.dim_y, "Wrong dimension for y. Expected %i, got %i" % (self.fmodel.dim_y, len(y))
        k = k or self.k
        sigma = sigma or self.sigma
        x_guess = self._guess_x(y, k = k)[0]
        dists, index = self.fmodel.dataset.nn_x(x_guess, k = k)
        w = self._weights(index, dists, sigma*sigma,  y)
        return [np.sum([wi*self.fmodel.dataset.get_x(idx)
                for wi, idx in zip(w, index)], axis = 0)]
486570
def Kdiag(self,X,target):
        """"""
        foo = np.zeros((X.shape[0],X.shape[0]))
        self.K(X,X,foo)
        target += np.diag(foo)
986725
private boolean kNNABOD(Database db, Relation<V> relation, DBIDs ids, WritableDoubleDataStore abodvalues, DoubleMinMax minmaxabod) {
    DistanceQuery<V> dq = db.getDistanceQuery(relation, SquaredEuclideanDistanceFunction.STATIC);
    KNNQuery<V> knnq = db.getKNNQuery(dq, DatabaseQuery.HINT_OPTIMIZED_ONLY);
    boolean squared = true;
    if(knnq == null) {
      dq = db.getDistanceQuery(relation, EuclideanDistanceFunction.STATIC);
      knnq = db.getKNNQuery(dq, DatabaseQuery.HINT_OPTIMIZED_ONLY);
      if(knnq == null) {
        return false;
      }
      squared = false;
    }
    SimilarityQuery<V> lk = db.getSimilarityQuery(relation, LinearKernelFunction.STATIC);
    int k1 = k + 1; // We will get the query point back by the knnq.

    MeanVariance s = new MeanVariance();
    for(DBIDIter pA = ids.iter(); pA.valid(); pA.advance()) {
      KNNList nl = knnq.getKNNForDBID(pA, k1);
      double simAA = lk.similarity(pA, pA);

      s.reset();
      DoubleDBIDListIter iB = nl.iter(), iC = nl.iter();
      for(; iB.valid(); iB.advance()) {
        double dAB = iB.doubleValue();
        double simAB = lk.similarity(pA, iB);
        if(!(dAB > 0.)) {
          continue;
        }
        for(iC.seek(iB.getOffset() + 1); iC.valid(); iC.advance()) {
          double dAC = iC.doubleValue();
          double simAC = lk.similarity(pA, iC);
          if(!(dAC > 0.)) {
            continue;
          }
          // Exploit bilinearity of scalar product:
          // <B-A, C-A> = <B, C-A> - <A,C-A>
          // = <B,C> - <B,A> - <A,C> + <A,A>
          double simBC = lk.similarity(iB, iC);
          double numerator = simBC - simAB - simAC + simAA;
          if(squared) {
            double div = 1. / (dAB * dAC);
            s.put(numerator * div, FastMath.sqrt(div));
          }
          else {
            double sqrtdiv = 1. / (dAB * dAC);
            s.put(numerator * sqrtdiv * sqrtdiv, sqrtdiv);
          }
        }
      }
      final double abof = s.getNaiveVariance();
      minmaxabod.put(abof);
      abodvalues.putDouble(pA, abof);
    }
    return true;
  }
761638
def vq_nearest_neighbor(x, hparams):
  """"""
  bottleneck_size = 2**hparams.bottleneck_bits
  means = hparams.means
  x_norm_sq = tf.reduce_sum(tf.square(x), axis=-1, keepdims=True)
  means_norm_sq = tf.reduce_sum(tf.square(means), axis=-1, keepdims=True)
  scalar_prod = tf.matmul(x, means, transpose_b=True)
  dist = x_norm_sq + tf.transpose(means_norm_sq) - 2 * scalar_prod
  if hparams.bottleneck_kind == "em":
    x_means_idx = tf.multinomial(-dist, num_samples=hparams.num_samples)
    x_means_hot = tf.one_hot(
        x_means_idx, depth=bottleneck_size)
    x_means_hot = tf.reduce_mean(x_means_hot, axis=1)
  else:
    x_means_idx = tf.argmax(-dist, axis=-1)
    x_means_hot = tf.one_hot(x_means_idx, depth=bottleneck_size)
  x_means = tf.matmul(x_means_hot, means)
  e_loss = tf.reduce_mean(tf.squared_difference(x, tf.stop_gradient(x_means)))
  return x_means_hot, e_loss
699394
private static double snorm(double[] x, int itol) {
        int n = x.length;

        if (itol <= 3) {
            double ans = 0.0;
            for (int i = 0; i < n; i++) {
                ans += x[i] * x[i];
            }
            return Math.sqrt(ans);
        } else {
            int isamax = 0;
            for (int i = 0; i < n; i++) {
                if (Math.abs(x[i]) > Math.abs(x[isamax])) {
                    isamax = i;
                }
            }

            return Math.abs(x[isamax]);
        }
    }
100401
private void calculateNearestNeighbors(int nn) {
        int n = this.dimension;
        if ((nn > 0) && (nn < n)) {
            nearestNeighbors = new int[n][nn];

            Neighbor[] neighbors = new Neighbor[n];
            for (int i = 0; i < n; i++) neighbors[i] = new Neighbor();

            for (int i = 0; i < n; i++) {
                for (int j = 0; j < n; j++) {
                    neighbors[j].id = j;
                    neighbors[j].distance = this.distances[i][j];
                }
                neighbors[i].distance = Integer.MAX_VALUE;
                Arrays.sort(neighbors);
                for (int j = 0; j < nn; j++) nearestNeighbors[i][j] = neighbors[j].id;
            }
        }
     }
1293319
def K(self, r, R, kwargs):
        
        u = r / R
        if np.min(u) < 1:
            raise ValueError("3d radius is smaller than projected radius! Does not make sense.")

        if self._type == 'const_wrong':
            beta = kwargs['beta']
            k = 1./2. * u**(2*beta - 1.) * ((3./2 - beta) * np.sqrt(np.pi) * special.gamma(beta - 1./2)/special.gamma(beta)
                        + beta * self._B(x=1./u**2, a=beta+1./2, b=1./2) - self._B(x=1./u**2, a=beta-1./2, b=1./2))
        elif self._type == 'const':
            beta = kwargs['beta']
            k = np.sqrt(1 - 1./u**2) / (1. - 2*beta) + np.sqrt(np.pi)/2 * special.gamma(beta - 1./2)/special.gamma(beta)\
                * (3./2 - beta) * u**(2*beta - 1.) * (1 - special.betainc(beta+1./2, 1./2, 1./u**2))
        elif self._type == 'isotropic':
            k = np.sqrt(1 - 1./u**2)
        elif self._type == 'radial':
            k = np.pi/4 * u - 1./2 * np.sqrt(1 - 1./u**2) - u/2. * np.arcsin(1./u)
        elif self._type == 'OsipkovMerritt':
            r_ani = kwargs['r_ani']
            ua = r_ani / R
            k = (ua**2 + 1./2) / (ua**2 + 1)**(3./2) * (u**2 + ua**2) / u * np.arctan(np.sqrt((u**2 - 1) / (ua**2 + 1))) \
                - 1./2/(ua**2 + 1) * np.sqrt(1 -1./u**2)
        elif self._type == 'Colin':
            r_ani = kwargs['r_ani']
            ua = r_ani / R
            if ua == 1:
                k = (1 + 1./u) * np.arccosh(u) - 1./6 * (8./u + 7) * np.sqrt((u-1.)/(u+1.))
            else:
                k = 0.5 / (ua**2 - 1) * np.sqrt(1 - 1./u**2) + (1. + ua/u) * np.cosh(u) - np.sign(ua - 1) * ua * \
                            (ua**2 - 0.5) / np.abs(ua**2-1)**(3./2) * (1. + ua/u) * np.arccosh((ua*u + 1)/(u + ua))
        else:
            raise ValueError('anisotropy type %s not supported!' % self._type)
        return k
935018
def _nbinom_ztrunc_p(mu, k_agg):
        

        p_eq = lambda p, mu, k_agg: (k_agg * p) / (1 - (1 + p)**-k_agg) - mu

        # The upper bound needs to be large. p will increase with increasing mu
        # and decreasing k_agg
        p = optim.brentq(p_eq, 1e-10, 1e10, args=(mu, k_agg))
        return p