docs: public OutlierResult run(Relation<V> relation) {
      final DBIDs ids = relation.getDBIDs();
      KNNQuery<V> knnQuery = QueryUtil.getKNNQuery(relation, getDistanceFunction(), k + 1);
  
      final int dim = RelationUtil.dimensionality(relation);
      if(k <= dim + 1) {
        LOG.warning("PCA is underspecified with a too low k! k should be at much larger than " + dim);
      }
  
      WritableDoubleDataStore cop_score = DataStoreUtil.makeDoubleStorage(ids, DataStoreFactory.HINT_HOT | DataStoreFactory.HINT_STATIC);
      WritableDataStore<double[]> cop_err_v = models ? DataStoreUtil.makeStorage(ids, DataStoreFactory.HINT_HOT | DataStoreFactory.HINT_STATIC, double[].class) : null;
      WritableIntegerDataStore cop_dim = models ? DataStoreUtil.makeIntegerStorage(ids, DataStoreFactory.HINT_HOT | DataStoreFactory.HINT_STATIC, -1) : null;
      // compute neighbors of each db object
      FiniteProgress prog = LOG.isVerbose() ? new FiniteProgress("Correlation Outlier Probabilities", relation.size(), LOG) : null;
  
      double[] centroid = new double[dim];
      double[] scores = new double[dim];
      ModifiableDBIDs nids = DBIDUtil.newHashSet(k + 10);
      for(DBIDIter id = ids.iter(); id.valid(); id.advance()) {
        nids.clear();
        nids.addDBIDs(knnQuery.getKNNForDBID(id, k + 1));
        nids.remove(id); // Do not use query object
  
        computeCentroid(centroid, relation, nids);
        PCAResult pcares = pca.processIds(nids, relation);
        double[][] tevecs = pcares.getEigenvectors();
        double[] evs = pcares.getEigenvalues();
        double[] projected = times(tevecs, minusEquals(relation.get(id).toArray(), centroid));
  
        if(dist == DistanceDist.CHISQUARED) {
          double sqdevs = 0;
          for(int d = 0; d < dim; d++) {
            double dev = projected[d];
            // Scale with variance and accumulate
            sqdevs += dev * dev / evs[d];
            scores[d] = 1 - ChiSquaredDistribution.cdf(sqdevs, d + 1);
          }
        }
        else {
          assert (dist == DistanceDist.GAMMA);
          double[][] dists = new double[dim][nids.size()];
          int j = 0;
          double[] srel = new double[dim];
          for(DBIDIter s = nids.iter(); s.valid() && j < nids.size(); s.advance(), j++) {
            V vec = relation.get(s);
            for(int d = 0; d < dim; d++) {
              srel[d] = vec.doubleValue(d) - centroid[d];
            }
            double sqdist = 0.0;
            for(int d = 0; d < dim; d++) {
              double serrd = transposeTimes(tevecs[d], srel);
              dists[d][j] = (sqdist += serrd * serrd / evs[d]);
            }
          }
          double sqdevs = 0;
          for(int d = 0; d < dim; d++) {
            // Scale with Stddev
            final double dev = projected[d];
            // Accumulate
            sqdevs += dev * dev / evs[d];
            // Sort, so we can trim the top 15% below.
            Arrays.sort(dists[d]);
            // Evaluate
            scores[d] = 1 - GammaChoiWetteEstimator.STATIC.estimate(dists[d], SHORTENED_ARRAY).cdf(sqdevs);
          }
        }
        // Find best score
        double min = Double.POSITIVE_INFINITY;
        int vdim = dim - 1;
        for(int d = 0; d < dim; d++) {
          double v = scores[d];
          if(v < min) {
            min = v;
            vdim = d;
          }
        }
        // Normalize the value
        final double prob = expect * (1 - min) / (expect + min);
        cop_score.putDouble(id, prob);
        if(models) {
          // Construct the error vector:
          Arrays.fill(projected, vdim + 1, dim, 0.);
          cop_err_v.put(id, timesEquals(transposeTimes(tevecs, projected), -prob));
          cop_dim.putInt(id, dim - vdim);
        }
        LOG.incrementProcessed(prog);
      }
      LOG.ensureCompleted(prog);
  
      // combine results.
      DoubleRelation scoreResult = new MaterializedDoubleRelation("Correlation Outlier Probabilities", COP_SCORES, cop_score, ids);
      OutlierScoreMeta scoreMeta = new ProbabilisticOutlierScore();
      OutlierResult result = new OutlierResult(scoreMeta, scoreResult);
      if(models) {
        result.addChildResult(new MaterializedRelation<>("Local Dimensionality", COP_DIM, TypeUtil.INTEGER, cop_dim, ids));
        result.addChildResult(new MaterializedRelation<>("Error vectors", COP_ERRORVEC, TypeUtil.DOUBLE_ARRAY, cop_err_v, ids));
      }
      return result;
    }