build_snn_graph.hpp

Go to the documentation of this file.

#ifndef SCRAN_GRAPH_CLUSTER_BUILD_SNN_GRAPH_HPP
#define SCRAN_GRAPH_CLUSTER_BUILD_SNN_GRAPH_HPP
 
#include <vector>
#include <algorithm>
#include <memory>
#include <cstddef>
 
#include "knncolle/knncolle.hpp"
#include "sanisizer/sanisizer.hpp"
 
#if __has_include("igraph.h")
#include "raiigraph/raiigraph.hpp"
#include "edges_to_graph.hpp"
#endif
 
namespace scran_graph_cluster {
 
enum class SnnWeightScheme : char { RANKED, NUMBER, JACCARD };
 
struct BuildSnnGraphOptions {
    int num_neighbors = 10;
 
    SnnWeightScheme weighting_scheme = SnnWeightScheme::RANKED;
 
    int num_threads = 1;
};
 
template<typename Node_, typename Weight_>
struct BuildSnnGraphResults {
    Node_ num_cells;
 
    std::vector<Node_> edges;
 
    std::vector<Weight_> weights;
};
 
#if __has_include("igraph.h")
typedef igraph_integer_t DefaultNode;
typedef igraph_real_t DefaultWeight;
#else
typedef int DefaultNode;
 
typedef double DefaultWeight;
#endif
 
template<typename Input_>
std::remove_cv_t<std::remove_reference_t<Input_> > I(Input_ x) {
    return x;
}
template<typename Node_ = DefaultNode, typename Weight_ = DefaultWeight, typename Index_, class GetNeighbors_, class GetIndex_>
void build_snn_graph(const Index_ num_cells, const GetNeighbors_ get_neighbors, const GetIndex_ get_index, const BuildSnnGraphOptions& options, BuildSnnGraphResults<Node_, Weight_>& output) {
    // Reverse mapping is not parallel-frendly, so we don't construct this with the neighbor search.
    std::vector<std::vector<Node_> > simple_hosts;
    std::vector<std::vector<std::pair<Node_, Weight_> > > ranked_hosts;
 
    // Check that all implicit casts from Index_ to Node_ will be safe.
    sanisizer::cast<Node_>(num_cells);
 
    if (options.weighting_scheme == SnnWeightScheme::RANKED) {
        sanisizer::resize(ranked_hosts, num_cells);
        for (Index_ i = 0; i < num_cells; ++i) {
            ranked_hosts[i].emplace_back(i, 0); // each point is its own 0-th nearest neighbor
            const auto& current = get_neighbors(i);
            Weight_ rank = 1;
            for (const auto& x : current) {
                ranked_hosts[get_index(x)].emplace_back(i, rank);
                ++rank;
            }
        }
    } else {
        sanisizer::resize(simple_hosts, num_cells);
        for (Index_ i = 0; i < num_cells; ++i) {
            simple_hosts[i].emplace_back(i); // each point is its own 0-th nearest neighbor
            const auto& current = get_neighbors(i);
            for (const auto& x : current) {
                simple_hosts[get_index(x)].emplace_back(i);
            }
        }
    }
 
    // Constructing the shared neighbor graph.
    auto edge_stores = sanisizer::create<std::vector<std::vector<Node_> > >(num_cells);
    auto weight_stores = sanisizer::create<std::vector<std::vector<Weight_> > >(num_cells);
 
    knncolle::parallelize(options.num_threads, num_cells, [&](const int, const Index_ start, const Index_ length) -> void {
        auto current_score = sanisizer::create<std::vector<Weight_> >(num_cells);
        std::vector<Node_> current_added;
        current_added.reserve(num_cells);
 
        for (Index_ j = start, end = start + length; j < end; ++j) {
            const auto& current_neighbors = get_neighbors(j);
            const auto nneighbors = current_neighbors.size();
 
            for (decltype(I(nneighbors)) i = 0; i <= nneighbors; ++i) {
                // First iteration treats 'j' as the zero-th neighbor.
                // Remaining iterations go through the neighbors of 'j'.
                const Node_ cur_neighbor = (i == 0 ? j : get_index(current_neighbors[i-1]));
 
                // Going through all observations 'h' for which 'cur_neighbor'
                // is a nearest neighbor, a.k.a., 'cur_neighbor' is a shared
                // neighbor of both 'h' and 'j'.
                if (options.weighting_scheme == SnnWeightScheme::RANKED) {
                    for (const auto& h : ranked_hosts[cur_neighbor]) {
                        const auto othernode = h.first;
                        if (othernode < static_cast<Node_>(j)) { // avoid duplicates from symmetry in the SNN calculations.
                            auto& existing_other = current_score[othernode];
 
                            // Recording the lowest combined rank per neighbor. 
                            const Weight_ currank = h.second + static_cast<Weight_>(i);
                            if (existing_other == 0) { 
                                existing_other = currank;
                                current_added.push_back(othernode);
                            } else if (existing_other > currank) {
                                existing_other = currank;
                            }
                        }
                    }
 
                } else {
                    for (const auto& othernode : simple_hosts[cur_neighbor]) {
                        if (othernode < static_cast<Node_>(j)) { // avoid duplicates from symmetry in the SNN calculations.
                            auto& existing_other = current_score[othernode];
 
                            // Recording the number of shared neighbors.
                            if (existing_other == 0) { 
                                current_added.push_back(othernode);
                            } 
                            ++existing_other;
                        }
                    }
                }
            }
           
            // Converting to edges.
            auto& current_edges = edge_stores[j];
            current_edges.reserve(current_added.size() * 2);
            auto& current_weights = weight_stores[j];
            current_weights.reserve(current_added.size());
 
            for (const auto othernode : current_added) {
                current_edges.push_back(j);
                current_edges.push_back(othernode);
 
                Weight_& otherscore = current_score[othernode];
                current_weights.push_back([&]{
                    if (options.weighting_scheme == SnnWeightScheme::RANKED) {
                        const Weight_ preliminary = static_cast<Weight_>(nneighbors) - otherscore / 2;
                        return std::max(preliminary, static_cast<Weight_>(1e-6)); // Ensuring that an edge with a positive weight is always reported.
                    } else if (options.weighting_scheme == SnnWeightScheme::JACCARD) {
                        return otherscore / (2 * (static_cast<Weight_>(nneighbors) + 1) - otherscore);
                    } else {
                        return otherscore;
                    }
                }());
 
                // Resetting all those added to zero.
                otherscore = 0;
            }
            current_added.clear();
        }
    });
 
    // Collating the total number of edges.
    std::size_t nedges = 0;
    for (const auto& w : weight_stores) {
        nedges = sanisizer::sum<std::size_t>(nedges, w.size());
    }
 
    output.num_cells = num_cells;
    output.weights.clear();
    output.weights.reserve(nedges);
    for (const auto& w : weight_stores) {
        output.weights.insert(output.weights.end(), w.begin(), w.end());
    }
    weight_stores.clear();
    weight_stores.shrink_to_fit(); // forcibly release memory so that we have some more space for edges.
 
    output.edges.clear();
    output.edges.reserve(nedges * 2);
    for (const auto& e : edge_stores) {
        output.edges.insert(output.edges.end(), e.begin(), e.end());
    }
 
    return;
}
 
template<typename Node_ = DefaultNode, typename Weight_ = DefaultWeight, typename Index_, typename Distance_>
BuildSnnGraphResults<Node_, Weight_> build_snn_graph(const knncolle::NeighborList<Index_, Distance_>& neighbors, const BuildSnnGraphOptions& options) {
    BuildSnnGraphResults<Node_, Weight_> output;
    build_snn_graph(
        sanisizer::cast<Index_>(neighbors.size()), 
        [&](const Index_ i) -> const std::vector<std::pair<Index_, Distance_> >& { return neighbors[i]; }, 
        [](const std::pair<Index_, Distance_>& x) -> Index_ { return x.first; }, 
        options,
        output
    );
    return output;
}
 
template<typename Node_ = int, typename Weight_ = double, typename Index_>
BuildSnnGraphResults<Node_, Weight_> build_snn_graph(const std::vector<std::vector<Index_> >& neighbors, const BuildSnnGraphOptions& options) {
    BuildSnnGraphResults<Node_, Weight_> output;
    build_snn_graph(
        sanisizer::cast<Index_>(neighbors.size()),
        [&](const Index_ i) -> const std::vector<Index_>& { return neighbors[i]; }, 
        [](const Index_ x) -> Index_ { return x; }, 
        options,
        output
    );
    return output;
}
 
template<typename Node_ = DefaultNode, typename Weight_ = DefaultWeight, typename Index_, typename Input_, typename Distance_>
BuildSnnGraphResults<Node_, Weight_> build_snn_graph(const knncolle::Prebuilt<Index_, Input_, Distance_>& prebuilt, const BuildSnnGraphOptions& options) {
    const auto neighbors = knncolle::find_nearest_neighbors_index_only(prebuilt, options.num_neighbors, options.num_threads);
    return build_snn_graph<Node_, Weight_>(neighbors, options);
}
 
template<typename Node_ = DefaultNode, typename Weight_ = DefaultWeight, typename Index_, typename Input_, typename Distance_, class Matrix_ = knncolle::Matrix<Index_, Input_> >
BuildSnnGraphResults<Node_, Weight_> build_snn_graph(
    const std::size_t num_dims, 
    const Index_ num_cells, 
    const Input_* data, 
    const knncolle::Builder<Index_, Input_, Distance_, Matrix_>& knn_method,
    const BuildSnnGraphOptions& options) 
{
    const auto prebuilt = knn_method.build_unique(knncolle::SimpleMatrix(num_dims, num_cells, data));
    return build_snn_graph<Node_, Weight_>(*prebuilt, options);
}
 
#if __has_include("igraph.h")
template<typename Node_ = DefaultNode, typename Weight_ = DefaultWeight>
raiigraph::Graph convert_to_graph(const BuildSnnGraphResults<Node_, Weight_>& result) {
    return edges_to_graph(result.edges.size(), result.edges.data(), result.num_cells, IGRAPH_UNDIRECTED);
}
#endif
 
}
 
#endif