simple.hpp

Go to the documentation of this file.

#ifndef MUMOSA_SIMPLE_HPP
#define MUMOSA_SIMPLE_HPP
 
#include <vector>
#include <stdexcept>
#include <cmath>
#include <algorithm>
#include <limits>
#include <cstddef>
#include <type_traits>
 
#include "knncolle/knncolle.hpp"
#include "quickstats/quickstats.hpp"
#include "sanisizer/sanisizer.hpp"
 
#include "utils.hpp"
 
namespace mumosa {
 
struct Options {
    int num_neighbors = 20;
 
    int num_threads = 1;
};
 
template<typename Index_, typename Distance_>
std::pair<Distance_, Distance_> compute_distance(const Index_ num_cells, Distance_* const distances) {
    if (num_cells == 0) {
        return std::pair<Distance_, Distance_>(0, 0);
    }
 
    const Distance_ med = quickstats::median(num_cells, distances);
    Distance_ rmsd = 0;
    for (Index_ i = 0; i < num_cells; ++i) {
        const auto d = distances[i];
        rmsd += d * d;
    }
    rmsd = std::sqrt(rmsd / num_cells);
    return std::make_pair(med, rmsd);
}
 
template<typename Index_, typename Input_, typename Distance_>
std::pair<Distance_, Distance_> compute_distance(
    const knncolle::Prebuilt<Index_, Input_, Distance_>& prebuilt,
    Distance_* const buffer,
    const Options& options
) {
    const Index_ nobs = prebuilt.num_observations();
    const auto capped_k = knncolle::cap_k(options.num_neighbors, nobs);
 
    knncolle::parallelize(options.num_threads, nobs, [&](const int, const Index_ start, const Index_ length) -> void {
        const auto searcher = prebuilt.initialize();
        std::vector<Distance_> cur_distances;
        for (Index_ i = start, end = start + length; i < end; ++i) {
            searcher->search(i, capped_k, NULL, &cur_distances);
            if (cur_distances.size()) {
                buffer[i] = cur_distances.back();
            } else {
                buffer[i] = 0; // i.e., only distance is that to itself.
            }
        }
    });
 
    return compute_distance(nobs, buffer);
}
 
template<typename Index_, typename Input_, typename Distance_, class Matrix_ = knncolle::Matrix<Index_, Input_> >
std::pair<Distance_, Distance_> compute_distance(
    const std::size_t num_dim,
    const Index_ num_cells,
    const Input_* const data,
    const knncolle::Builder<Index_, Input_, Distance_, Matrix_>& builder,
    Distance_* const buffer,
    const Options& options
) {
    const auto prebuilt = builder.build_unique(knncolle::SimpleMatrix(num_dim, num_cells, data));
    return compute_distance(*prebuilt, buffer, options);
}
 
}
 
#endif