docs: @Override
    public OutlierResult run(Database db, Relation<V> relation) {
      ArrayDBIDs ids = DBIDUtil.ensureArray(relation.getDBIDs());
      DBIDArrayIter pB = ids.iter(), pC = ids.iter();
      SimilarityQuery<V> sq = db.getSimilarityQuery(relation, kernelFunction);
      KernelMatrix kernelMatrix = new KernelMatrix(sq, relation, ids);
  
      // Output storage.
      WritableDoubleDataStore abodvalues = DataStoreUtil.makeDoubleStorage(ids, DataStoreFactory.HINT_STATIC);
      DoubleMinMax minmaxabod = new DoubleMinMax();
      double max = 0.;
  
      // Storage for squared distances (will be reused!)
      WritableDoubleDataStore sqDists = DataStoreUtil.makeDoubleStorage(ids, DataStoreFactory.HINT_TEMP | DataStoreFactory.HINT_HOT);
      // Nearest neighbor heap (will be reused!)
      KNNHeap nn = DBIDUtil.newHeap(k);
  
      // Priority queue for candidates
      ModifiableDoubleDBIDList candidates = DBIDUtil.newDistanceDBIDList(relation.size());
      // get Candidate Ranking
      for(DBIDIter pA = relation.iterDBIDs(); pA.valid(); pA.advance()) {
        // Compute nearest neighbors and distances.
        nn.clear();
        double simAA = kernelMatrix.getSimilarity(pA, pA);
        // Sum of 1./(|AB|) and 1./(|AB|^2); for computing R2.
        double sumid = 0., sumisqd = 0.;
        for(pB.seek(0); pB.valid(); pB.advance()) {
          if(DBIDUtil.equal(pB, pA)) {
            continue;
          }
          double simBB = kernelMatrix.getSimilarity(pB, pB);
          double simAB = kernelMatrix.getSimilarity(pA, pB);
          double sqdAB = simAA + simBB - simAB - simAB;
          sqDists.putDouble(pB, sqdAB);
          final double isqdAB = 1. / sqdAB;
          sumid += FastMath.sqrt(isqdAB);
          sumisqd += isqdAB;
          // Update heap
          nn.insert(sqdAB, pB);
        }
  
        // Compute FastABOD approximation, adjust for lower bound.
        // LB-ABOF is defined via a numerically unstable formula.
        // Variance as E(X^2)-E(X)^2 suffers from catastrophic cancellation!
        // TODO: ensure numerical precision!
        double nnsum = 0., nnsumsq = 0., nnsumisqd = 0.;
        KNNList nl = nn.toKNNList();
        DoubleDBIDListIter iB = nl.iter(), iC = nl.iter();
        for(; iB.valid(); iB.advance()) {
          double sqdAB = iB.doubleValue();
          double simAB = kernelMatrix.getSimilarity(pA, iB);
          if(!(sqdAB > 0.)) {
            continue;
          }
          for(iC.seek(iB.getOffset() + 1); iC.valid(); iC.advance()) {
            double sqdAC = iC.doubleValue();
            double simAC = kernelMatrix.getSimilarity(pA, iC);
            if(!(sqdAC > 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 = kernelMatrix.getSimilarity(iB, iC);
            double numerator = simBC - simAB - simAC + simAA;
            double sqweight = 1. / (sqdAB * sqdAC);
            double weight = FastMath.sqrt(sqweight);
            double val = numerator * sqweight;
            nnsum += val * weight;
            nnsumsq += val * val * weight;
            nnsumisqd += sqweight;
          }
        }
        // Remaining weight, term R2:
        double r2 = sumisqd * sumisqd - 2. * nnsumisqd;
        double tmp = (2. * nnsum + r2) / (sumid * sumid);
        double lbabof = 2. * nnsumsq / (sumid * sumid) - tmp * tmp;
  
        // Track maximum?
        if(lbabof > max) {
          max = lbabof;
        }
        abodvalues.putDouble(pA, lbabof);
        candidates.add(lbabof, pA);
      }
      minmaxabod.put(max); // Put maximum from approximate values.
      candidates.sort();
  
      // refine Candidates
      int refinements = 0;
      DoubleMinHeap topscores = new DoubleMinHeap(l);
      MeanVariance s = new MeanVariance();
      for(DoubleDBIDListIter pA = candidates.iter(); pA.valid(); pA.advance()) {
        // Stop refining
        if(topscores.size() >= k && pA.doubleValue() > topscores.peek()) {
          break;
        }
        final double abof = computeABOF(kernelMatrix, pA, pB, pC, s);
        // Store refined score:
        abodvalues.putDouble(pA, abof);
        minmaxabod.put(abof);
        // Update the heap tracking the top scores.
        if(topscores.size() < k) {
          topscores.add(abof);
        }
        else {
          if(topscores.peek() > abof) {
            topscores.replaceTopElement(abof);
          }
        }
        refinements += 1;
      }
      if(LOG.isStatistics()) {
        LoggingConfiguration.setVerbose(Level.VERYVERBOSE);
        LOG.statistics(new LongStatistic("lb-abod.refinements", refinements));
      }
      // Build result representation.
      DoubleRelation scoreResult = new MaterializedDoubleRelation("Angle-based Outlier Detection", "abod-outlier", abodvalues, ids);
      OutlierScoreMeta scoreMeta = new InvertedOutlierScoreMeta(minmaxabod.getMin(), minmaxabod.getMax(), 0.0, Double.POSITIVE_INFINITY);
      return new OutlierResult(scoreMeta, scoreResult);
    }