Main loop of OUTRES. Run for each object

@param s start dimension
@param subspace Current subspace
@param id Current object ID
@param kernel Kernel
@param cands neighbor candidates
@return Score
552215
public double outresScore(final int s, long[] subspace, DBIDRef id, KernelDensityEstimator kernel, DBIDs cands) {
    double score = 1.0; // Initial score is 1.0
    final SubspaceEuclideanDistanceFunction df = new SubspaceEuclideanDistanceFunction(subspace);
    MeanVariance meanv = new MeanVariance();
    ModifiableDoubleDBIDList neighcand = DBIDUtil.newDistanceDBIDList(cands.size());
    ModifiableDoubleDBIDList nn = DBIDUtil.newDistanceDBIDList(cands.size());

    for(int i = s; i < kernel.dim; i++) {
      assert !BitsUtil.get(subspace, i);
      BitsUtil.setI(subspace, i);
      df.setSelectedDimensions(subspace);
      final double adjustedEps = kernel.adjustedEps(kernel.dim);
      DoubleDBIDList neigh = initialRange(id, cands, df, adjustedEps * 2, kernel, neighcand);
      // Relevance test
      if(neigh.size() > 2) {
        if(relevantSubspace(subspace, neigh, kernel)) {
          final double density = kernel.subspaceDensity(subspace, neigh);
          // Compute mean and standard deviation for densities of neighbors.
          meanv.reset();
          for(DoubleDBIDListIter neighbor = neigh.iter(); neighbor.valid(); neighbor.advance()) {
            subsetNeighborhoodQuery(neighcand, neighbor, df, adjustedEps, kernel, nn);
            meanv.put(kernel.subspaceDensity(subspace, nn));
          }
          final double deviation = (meanv.getMean() - density) / (2. * meanv.getSampleStddev());
          // High deviation:
          if(deviation >= 1) {
            score *= density / deviation;
          }
          // Recursion
          score *= outresScore(i + 1, subspace, id, kernel, neighcand);
        }
      }
      BitsUtil.clearI(subspace, i);
    }
    return score;
  }
847024
private ScalesResult run(Relation<? extends NumberVector> rel) {
    final int dim = RelationUtil.dimensionality(rel);
    LinearScale[] scales = new LinearScale[dim];
    if(minmax == null) {
      DoubleMinMax mm = new DoubleMinMax();
      for(DBIDIter iditer = rel.iterDBIDs(); iditer.valid(); iditer.advance()) {
        NumberVector vec = rel.get(iditer);
        for(int d = 0; d < dim; d++) {
          final double val = vec.doubleValue(d);
          if(val != val) {
            continue; // NaN
          }
          mm.put(val);
        }
      }
      LinearScale scale = new LinearScale(mm.getMin(), mm.getMax());
      for(int i = 0; i < dim; i++) {
        scales[i] = scale;
      }
    }
    else {
      // Use predefined.
      LinearScale scale = new LinearScale(minmax[0], minmax[1]);
      for(int i = 0; i < dim; i++) {
        scales[i] = scale;
      }
    }
    ScalesResult res = new ScalesResult(scales);
    return res;
  }
1218012
public final Solution<T> cluster(final JavaRDD<T> input_data) {
        if (k <= 0) {
            throw new IllegalStateException("k must be > 0!");
        }

        if (similarity == null) {
            throw new IllegalStateException("Similarity is not defined!");
        }

        if (budget == null) {
            throw new IllegalStateException("No budget is defined!");
        }

        JavaRDD<T> data = input_data.cache();

        // Keep track of best solution
        Solution<T> solution = new Solution<>(data.count());
        logger.debug("Dataset contains {} objects", solution.getDatasetSize());

        neighbor_generator.init(k);
        points_supplier = new RandomPointsSupplier<>(
                data, solution.getDatasetSize());

        // Select random initial solution and evaluate
        ArrayList<T> random_medoids = points_supplier.pick(k);
        solution.setNewMedoids(random_medoids, evaluate(data, random_medoids));
        solution.incComputedSimilarities(k * solution.getDatasetSize());

        while (true) {
            logger.trace("Trial {}", solution.getTrials());

            // Select neighbor solution
            CountingSimilarity<T> counting_sim =
                    new CountingSimilarity<>(similarity);

            ArrayList<T> candidate;
            try {
                candidate = neighbor_generator.getNeighbor(
                        new NeighborGeneratorHelper<>(this),
                        solution,
                        counting_sim);
            } catch (NoNeighborFoundException ex) {
                solution.incComputedSimilarities(counting_sim.getCount());
                break;
            }
            solution.incComputedSimilarities(counting_sim.getCount());

            // Evaluate
            double candidate_similarity = evaluate(data, candidate);
            solution.incComputedSimilarities(k * solution.getDatasetSize());
            neighbor_generator.notifyCandidateSolutionCost(
                    candidate, candidate_similarity);

            if (candidate_similarity > solution.getTotalSimilarity()) {
                solution.setNewMedoids(candidate, candidate_similarity);
                solution.incIterations();
                neighbor_generator.notifyNewSolution(
                        candidate, candidate_similarity);
            }

            solution.incTrials();

            if (budget.isExhausted(solution)) {
                break;
            }
        }

        solution.end();
        logger.debug("Found solution {}", solution);
        return solution;

    }
260357
public OutlierResult run(Relation<V> relation) {
    final DBIDs ids = relation.getDBIDs();

    ArrayList<ArrayDBIDs> subspaceIndex = buildOneDimIndexes(relation);
    Set<HiCSSubspace> subspaces = calculateSubspaces(relation, subspaceIndex, rnd.getSingleThreadedRandom());

    if(LOG.isVerbose()) {
      LOG.verbose("Number of high-contrast subspaces: " + subspaces.size());
    }
    List<DoubleRelation> results = new ArrayList<>();
    FiniteProgress prog = LOG.isVerbose() ? new FiniteProgress("Calculating Outlier scores for high Contrast subspaces", subspaces.size(), LOG) : null;

    // run outlier detection and collect the result
    // TODO extend so that any outlierAlgorithm can be used (use materialized
    // relation instead of SubspaceEuclideanDistanceFunction?)
    for(HiCSSubspace dimset : subspaces) {
      if(LOG.isVerbose()) {
        LOG.verbose("Performing outlier detection in subspace " + dimset);
      }

      ProxyDatabase pdb = new ProxyDatabase(ids);
      pdb.addRelation(new ProjectedView<>(relation, new NumericalFeatureSelection<V>(dimset)));

      // run LOF and collect the result
      OutlierResult result = outlierAlgorithm.run(pdb);
      results.add(result.getScores());
      LOG.incrementProcessed(prog);
    }
    LOG.ensureCompleted(prog);

    WritableDoubleDataStore scores = DataStoreUtil.makeDoubleStorage(relation.getDBIDs(), DataStoreFactory.HINT_STATIC);
    DoubleMinMax minmax = new DoubleMinMax();

    for(DBIDIter iditer = relation.iterDBIDs(); iditer.valid(); iditer.advance()) {
      double sum = 0.0;
      for(DoubleRelation r : results) {
        final double s = r.doubleValue(iditer);
        if(!Double.isNaN(s)) {
          sum += s;
        }
      }
      scores.putDouble(iditer, sum);
      minmax.put(sum);
    }
    OutlierScoreMeta meta = new BasicOutlierScoreMeta(minmax.getMin(), minmax.getMax());
    DoubleRelation scoreres = new MaterializedDoubleRelation("HiCS", "HiCS-outlier", scores, relation.getDBIDs());

    return new OutlierResult(meta, scoreres);
  }
204793
public PointerHierarchyRepresentationResult run(Database db, Relation<O> relation) {
    if(SingleLinkage.class.isInstance(linkage)) {
      LOG.verbose("Notice: SLINK is a much faster algorithm for single-linkage clustering!");
    }
    DistanceQuery<O> dq = db.getDistanceQuery(relation, getDistanceFunction());
    final DBIDs ids = relation.getDBIDs();
    MatrixParadigm mat = new MatrixParadigm(ids);

    // Compute the initial (lower triangular) distance matrix.
    initializeDistanceMatrix(mat, dq, linkage);

    // Initialize space for result:
    PointerHierarchyRepresentationBuilder builder = new PointerHierarchyRepresentationBuilder(ids, dq.getDistanceFunction().isSquared());

    nnChainCore(mat, builder);
    return builder.complete();
  }
1279514
public OutlierResult run(Relation<? extends NumberVector> relation) {
    final DBIDs ids = relation.getDBIDs();
    WritableDoubleDataStore ranks = DataStoreUtil.makeDoubleStorage(ids, DataStoreFactory.HINT_STATIC);
    DoubleMinMax minmax = new DoubleMinMax();

    KernelDensityEstimator kernel = new KernelDensityEstimator(relation, eps);
    long[] subspace = BitsUtil.zero(kernel.dim);

    FiniteProgress progress = LOG.isVerbose() ? new FiniteProgress("OUTRES scores", ids.size(), LOG) : null;

    for(DBIDIter iditer = ids.iter(); iditer.valid(); iditer.advance()) {
      BitsUtil.zeroI(subspace);
      double score = outresScore(0, subspace, iditer, kernel, ids);
      ranks.putDouble(iditer, score);
      minmax.put(score);
      LOG.incrementProcessed(progress);
    }
    LOG.ensureCompleted(progress);

    OutlierScoreMeta meta = new InvertedOutlierScoreMeta(minmax.getMin(), minmax.getMax(), 0., 1., 1.);
    return new OutlierResult(meta, new MaterializedDoubleRelation("OUTRES", "outres-score", ranks, ids));
  }
990265
protected void computeOutlierScores(KNNQuery<O> knnq, final DBIDs ids, WritableDataStore<double[]> densities, WritableDoubleDataStore kdeos, DoubleMinMax minmax) {
    final int knum = kmax + 1 - kmin;
    FiniteProgress prog = LOG.isVerbose() ? new FiniteProgress("Computing KDEOS scores", ids.size(), LOG) : null;

    double[][] scratch = new double[knum][kmax + 5];
    MeanVariance mv = new MeanVariance();

    for(DBIDIter iter = ids.iter(); iter.valid(); iter.advance()) {
      double[] dens = densities.get(iter);
      KNNList neighbors = knnq.getKNNForDBID(iter, kmax + 1);
      if(scratch[0].length < neighbors.size()) {
        // Resize scratch. Add some extra margin again.
        scratch = new double[knum][neighbors.size() + 5];
      }
      { // Store density matrix of neighbors
        int i = 0;
        for(DoubleDBIDListIter neighbor = neighbors.iter(); neighbor.valid(); neighbor.advance(), i++) {
          double[] ndens = densities.get(neighbor);
          for(int k = 0; k < knum; k++) {
            scratch[k][i] = ndens[k];
          }
        }
        assert (i == neighbors.size());
      }
      // Compute means and stddevs for each k
      double score = 0.;
      for(int i = 0; i < knum; i++) {
        mv.reset();
        for(int j = 0; j < neighbors.size(); j++) {
          mv.put(scratch[i][j]);
        }
        final double mean = mv.getMean(), stddev = mv.getSampleStddev();
        if(stddev > 0.) {
          score += (mean - dens[i]) / stddev;
        }
      }
      score /= knum; // average
      score = NormalDistribution.standardNormalCDF(score);
      minmax.put(score);
      kdeos.put(iter, score);
      LOG.incrementProcessed(prog);
    }
    LOG.ensureCompleted(prog);
  }
1230782
def edge_cluster_angle(edge_dirs, subspaces1, subspaces2):
  ''''''
  QG = edge_dirs / np.linalg.norm(edge_dirs, ord=2, axis=1)[:,None]
  X1 = np.einsum('...ij,...i->...j', subspaces1, QG)
  X2 = np.einsum('...ij,...i->...j', subspaces2, QG)
  # TODO: check the math on this for more cases
  # angles = np.maximum(1-np.sum(X1**2, axis=1), 1-np.sum(X2**2, axis=1))
  C1 = np.linalg.svd(X1[:,:,None], compute_uv=False)
  C2 = np.linalg.svd(X2[:,:,None], compute_uv=False)
  angles = np.maximum(1-C1**2, 1-C2**2)[:,0]
  angles[np.isnan(angles)] = 0.0  # nan when edge length == 0
  return angles
159196
@Override
  public OutlierResult run(Database db, Relation<V> relation) {
    DBIDs ids = relation.getDBIDs();
    WritableDoubleDataStore abodvalues = DataStoreUtil.makeDoubleStorage(ids, DataStoreFactory.HINT_STATIC);
    DoubleMinMax minmaxabod = new DoubleMinMax();
    if(kernelFunction.getClass() == LinearKernelFunction.class) {
      if(!kNNABOD(db, relation, ids, abodvalues, minmaxabod)) {
        // Fallback, if we do not have an index.
        fastABOD(db, relation, ids, abodvalues, minmaxabod);
      }
    }
    else {
      fastABOD(db, relation, ids, abodvalues, minmaxabod);
    }

    // Build result representation.
    DoubleRelation scoreResult = new MaterializedDoubleRelation("Angle-Based Outlier Degree", "abod-outlier", abodvalues, relation.getDBIDs());
    OutlierScoreMeta scoreMeta = new InvertedOutlierScoreMeta(minmaxabod.getMin(), minmaxabod.getMax(), 0.0, Double.POSITIVE_INFINITY);
    return new OutlierResult(scoreMeta, scoreResult);
  }
1138707
public static int rank(DMatrixRMaj A , double threshold ) {
        SingularValueDecomposition_F64<DMatrixRMaj> svd = DecompositionFactory_DDRM.svd(A.numRows,A.numCols,false,false,true);

        if( svd.inputModified() )
            A = A.copy();

        if( !svd.decompose(A) )
            throw new RuntimeException("Decomposition failed");

        return SingularOps_DDRM.rank(svd, threshold);
    }
1201057
@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);
  }
1092413
private Clustering<SubspaceModel> computeClusters(Relation<V> database, DiSHClusterOrder clusterOrder) {
    final int dimensionality = RelationUtil.dimensionality(database);

    // extract clusters
    Object2ObjectOpenCustomHashMap<long[], List<ArrayModifiableDBIDs>> clustersMap = extractClusters(database, clusterOrder);
    logClusterSizes("Step 1: extract clusters", dimensionality, clustersMap);

    // check if there are clusters < minpts
    checkClusters(database, clustersMap);
    logClusterSizes("Step 2: check clusters", dimensionality, clustersMap);

    // sort the clusters
    List<Cluster<SubspaceModel>> clusters = sortClusters(database, clustersMap);
    if(LOG.isVerbose()) {
      StringBuilder msg = new StringBuilder("Step 3: sort clusters");
      for(Cluster<SubspaceModel> c : clusters) {
        msg.append('\n').append(BitsUtil.toStringLow(c.getModel().getSubspace().getDimensions(), dimensionality)).append(" ids ").append(c.size());
      }
      LOG.verbose(msg.toString());
    }

    // build the hierarchy
    Clustering<SubspaceModel> clustering = new Clustering<>("DiSH clustering", "dish-clustering");
    buildHierarchy(database, clustering, clusters, dimensionality);
    if(LOG.isVerbose()) {
      StringBuilder msg = new StringBuilder("Step 4: build hierarchy");
      for(Cluster<SubspaceModel> c : clusters) {
        msg.append('\n').append(BitsUtil.toStringLow(c.getModel().getSubspace().getDimensions(), dimensionality)).append(" ids ").append(c.size());
        for(It<Cluster<SubspaceModel>> iter = clustering.getClusterHierarchy().iterParents(c); iter.valid(); iter.advance()) {
          msg.append("\n   parent ").append(iter.get());
        }
        for(It<Cluster<SubspaceModel>> iter = clustering.getClusterHierarchy().iterChildren(c); iter.valid(); iter.advance()) {
          msg.append("\n   child ").append(iter.get());
        }
      }
      LOG.verbose(msg.toString());
    }

    // build result
    for(Cluster<SubspaceModel> c : clusters) {
      if(clustering.getClusterHierarchy().numParents(c) == 0) {
        clustering.addToplevelCluster(c);
      }
    }
    return clustering;
  }
984780
@Override
    public void sgesvd(byte jobu, byte jobvt, int M, int N, INDArray A, INDArray S, INDArray U, INDArray VT,
                    INDArray INFO) {
        INDArray superb = Nd4j.create(DataType.FLOAT, M < N ? M : N ) ;
        int status = LAPACKE_sgesvd(getColumnOrder(A), jobu, jobvt, M, N, 
                        (FloatPointer)A.data().addressPointer(), getLda(A),
                        (FloatPointer)S.data().addressPointer(), 
                        U == null ? null : (FloatPointer)U.data().addressPointer(), U == null ? 1 : getLda(U),
                        VT == null ? null : (FloatPointer)VT.data().addressPointer(), VT == null ? 1 : getLda(VT), 
                        (FloatPointer)superb.data().addressPointer() 
                        );
        if( status != 0 ) {
            throw new BlasException( "Failed to execute sgesvd", status ) ;
        }
    }
337466
public Clustering<Model> run(Database database, Relation<V> relation) {
    // current dimensionality associated with each seed
    int dim_c = RelationUtil.dimensionality(relation);

    if(dim_c < l) {
      throw new IllegalStateException("Dimensionality of data < parameter l! " + "(" + dim_c + " < " + l + ")");
    }

    // current number of seeds
    int k_c = Math.min(relation.size(), k_i * k);

    // pick k0 > k points from the database
    List<ORCLUSCluster> clusters = initialSeeds(relation, k_c);

    double beta = FastMath.exp(-FastMath.log(dim_c / (double) l) * FastMath.log(1 / alpha) / FastMath.log(k_c / (double) k));

    IndefiniteProgress cprogress = LOG.isVerbose() ? new IndefiniteProgress("Current number of clusters:", LOG) : null;

    while(k_c > k) {
      // find partitioning induced by the seeds of the clusters
      assign(relation, clusters);

      // determine current subspace associated with each cluster
      for(ORCLUSCluster cluster : clusters) {
        if(cluster.objectIDs.size() > 0) {
          cluster.basis = findBasis(relation, cluster, dim_c);
        }
      }

      // reduce number of seeds and dimensionality associated with
      // each seed
      k_c = (int) Math.max(k, k_c * alpha);
      dim_c = (int) Math.max(l, dim_c * beta);
      merge(relation, clusters, k_c, dim_c, cprogress);
      if(cprogress != null) {
        cprogress.setProcessed(clusters.size(), LOG);
      }
    }
    assign(relation, clusters);

    LOG.setCompleted(cprogress);

    // get the result
    Clustering<Model> r = new Clustering<>("ORCLUS clustering", "orclus-clustering");
    for(ORCLUSCluster c : clusters) {
      r.addToplevelCluster(new Cluster<Model>(c.objectIDs, ClusterModel.CLUSTER));
    }
    return r;
  }
1198453
protected double computeLOFScore(KNNQuery<O> knnq, DBIDRef cur, DoubleDataStore lrds) {
    final double lrdp = lrds.doubleValue(cur);
    if(Double.isInfinite(lrdp)) {
      return 1.0;
    }
    double sum = 0.;
    int count = 0;
    final KNNList neighbors = knnq.getKNNForDBID(cur, k);
    for(DBIDIter neighbor = neighbors.iter(); neighbor.valid(); neighbor.advance()) {
      // skip the point itself
      if(DBIDUtil.equal(cur, neighbor)) {
        continue;
      }
      sum += lrds.doubleValue(neighbor);
      ++count;
    }
    return sum / (lrdp * count);
  }
577386
private void calculateContrast(Relation<? extends NumberVector> relation, HiCSSubspace subspace, ArrayList<ArrayDBIDs> subspaceIndex, Random random) {
    final int card = subspace.cardinality();
    final double alpha1 = FastMath.pow(alpha, (1.0 / card));
    final int windowsize = (int) (relation.size() * alpha1);
    final FiniteProgress prog = LOG.isDebugging() ? new FiniteProgress("Monte-Carlo iterations", m, LOG) : null;

    int retries = 0;
    double deviationSum = 0.0;
    for(int i = 0; i < m; i++) {
      // Choose a random set bit.
      int chosen = -1;
      for(int tmp = random.nextInt(card); tmp >= 0; tmp--) {
        chosen = subspace.nextSetBit(chosen + 1);
      }
      // initialize sample
      DBIDs conditionalSample = relation.getDBIDs();

      for(int j = subspace.nextSetBit(0); j >= 0; j = subspace.nextSetBit(j + 1)) {
        if(j == chosen) {
          continue;
        }
        ArrayDBIDs sortedIndices = subspaceIndex.get(j);
        ArrayModifiableDBIDs indexBlock = DBIDUtil.newArray(windowsize);
        // initialize index block
        DBIDArrayIter iter = sortedIndices.iter();
        iter.seek(random.nextInt(relation.size() - windowsize));
        for(int k = 0; k < windowsize; k++, iter.advance()) {
          indexBlock.add(iter); // select index block
        }

        conditionalSample = DBIDUtil.intersection(conditionalSample, indexBlock);
      }
      if(conditionalSample.size() < 10) {
        retries++;
        if(LOG.isDebugging()) {
          LOG.debug("Sample size very small. Retry no. " + retries);
        }
        if(retries >= MAX_RETRIES) {
          LOG.warning("Too many retries, for small samples: " + retries);
        }
        else {
          i--;
          continue;
        }
      }
      // Project conditional set
      double[] sampleValues = new double[conditionalSample.size()];
      {
        int l = 0;
        for(DBIDIter iter = conditionalSample.iter(); iter.valid(); iter.advance()) {
          sampleValues[l++] = relation.get(iter).doubleValue(chosen);
        }
      }
      // Project full set
      double[] fullValues = new double[relation.size()];
      {
        int l = 0;
        for(DBIDIter iter = subspaceIndex.get(chosen).iter(); iter.valid(); iter.advance()) {
          fullValues[l++] = relation.get(iter).doubleValue(chosen);
        }
      }
      double contrast = statTest.deviation(fullValues, sampleValues);
      if(Double.isNaN(contrast)) {
        i--;
        LOG.warning("Contrast was NaN");
        continue;
      }
      deviationSum += contrast;
      LOG.incrementProcessed(prog);
    }
    LOG.ensureCompleted(prog);
    subspace.contrast = deviationSum / m;
  }
1150031
protected Cluster<SubspaceModel> makeCluster(Relation<V> relation, DBIDs C, long[] D) {
    DBIDs ids = DBIDUtil.newHashSet(C); // copy, also to lose distance values!
    Cluster<SubspaceModel> cluster = new Cluster<>(ids);
    cluster.setModel(new SubspaceModel(new Subspace(D), Centroid.make(relation, ids).getArrayRef()));
    return cluster;
  }