if(progressPreproc != null) {
progressPreproc.ensureCompleted(logger);
}
// LOCI main step
FiniteProgress progressLOCI = logger.isVerbose() ? new FiniteProgress("LOCI scores", relation.size(), logger) : null;
WritableDoubleDataStore mdef_norm = DataStoreUtil.makeDoubleStorage(relation.getDBIDs(), DataStoreFactory.HINT_STATIC);
WritableDoubleDataStore mdef_radius = DataStoreUtil.makeDoubleStorage(relation.getDBIDs(), DataStoreFactory.HINT_STATIC);
DoubleMinMax minmax = new DoubleMinMax();
for(DBID id : relation.iterDBIDs()) {
final List<DoubleIntPair> cdist = interestingDistances.get(id);
final double maxdist = cdist.get(cdist.size() - 1).first;
final int maxneig = cdist.get(cdist.size() - 1).second;
double maxmdefnorm = 0.0;
double maxnormr = 0;
if(maxneig >= nmin) {
D range = distFunc.getDistanceFactory().fromDouble(maxdist);
// Compute the largest neighborhood we will need.
List<DistanceResultPair<D>> maxneighbors = rangeQuery.getRangeForDBID(id, range);
// Ensure the set is sorted. Should be a no-op with most indexes.
Collections.sort(maxneighbors);
// For any critical distance, compute the normalized MDEF score.
for(DoubleIntPair c : cdist) {
// Only start when minimum size is fulfilled
if (c.second < nmin) {
continue;
}
final double r = c.first;
final double alpha_r = alpha * r;
// compute n(p_i, \alpha * r) from list (note: alpha_r is different from c!)
final int n_alphar = elementsAtRadius(cdist, alpha_r);
// compute \hat{n}(p_i, r, \alpha) and the corresponding \simga_{MDEF}
MeanVariance mv_n_r_alpha = new MeanVariance();
for(DistanceResultPair<D> ne : maxneighbors) {
// Stop at radius r
if(ne.getDistance().doubleValue() > r) {
break;
}
int rn_alphar = elementsAtRadius(interestingDistances.get(ne.getDBID()), alpha_r);
mv_n_r_alpha.put(rn_alphar);
}
// We only use the average and standard deviation
final double nhat_r_alpha = mv_n_r_alpha.getMean();
final double sigma_nhat_r_alpha = mv_n_r_alpha.getNaiveStddev();
// Redundant divisions removed.
final double mdef = (nhat_r_alpha - n_alphar); // / nhat_r_alpha;
final double sigmamdef = sigma_nhat_r_alpha; // / nhat_r_alpha;
final double mdefnorm = mdef / sigmamdef;
if(mdefnorm > maxmdefnorm) {
maxmdefnorm = mdefnorm;
maxnormr = r;
}
}
}
else {
// FIXME: when nmin was not fulfilled - what is the proper value then?
maxmdefnorm = 1.0;
maxnormr = maxdist;
}
mdef_norm.putDouble(id, maxmdefnorm);
mdef_radius.putDouble(id, maxnormr);
minmax.put(maxmdefnorm);
if(progressLOCI != null) {
progressLOCI.incrementProcessed(logger);
}
}