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"
 
#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 Node_ = DefaultNode, typename Weight_ = DefaultWeight, typename Index_, class GetNeighbors_, class GetIndex_>
void build_snn_graph(Index_ num_cells, GetNeighbors_ get_neighbors, 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;
 
    if (options.weighting_scheme == SnnWeightScheme::RANKED) {
        ranked_hosts.resize(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 (auto x : current) {
                ranked_hosts[get_index(x)].emplace_back(i, rank);
                ++rank;
            }
        }
    } else {
        simple_hosts.resize(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 (auto x : current) {
                simple_hosts[get_index(x)].emplace_back(i);
            }
        }
    }
 
    // Constructing the shared neighbor graph.
    std::vector<std::vector<Node_> > edge_stores(num_cells);
    std::vector<std::vector<Weight_> > weight_stores(num_cells);
 
    knncolle::parallelize(options.num_threads, num_cells, [&](int, Index_ start, Index_ length) -> void {
        std::vector<Weight_> current_score(num_cells);
        std::vector<Node_> current_added;
        current_added.reserve(num_cells);
 
        for (Index_ j = start, end = start + length; j < end; ++j) {
            const Node_ j_cast = j;
 
            const auto& current_neighbors = get_neighbors(j);
            auto nneighbors = current_neighbors.size();
 
            for (decltype(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) {
                    Weight_ i_cast = i;
                    for (const auto& h : ranked_hosts[cur_neighbor]) {
                        auto othernode = h.first;
                        if (othernode < j_cast) { // avoid duplicates from symmetry in the SNN calculations.
                            auto& existing_other = current_score[othernode];
 
                            // Recording the lowest combined rank per neighbor (casting to avoid overflow on Node_).
                            Weight_ currank = h.second + i_cast;
                            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 < j_cast) { // 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 (auto othernode : current_added) {
                Weight_& otherscore = current_score[othernode];
                Weight_ finalscore;
                if (options.weighting_scheme == SnnWeightScheme::RANKED) {
                    finalscore = static_cast<Weight_>(nneighbors) - 0.5 * static_cast<Weight_>(otherscore);
                } else {
                    finalscore = otherscore;
                    if (options.weighting_scheme == SnnWeightScheme::JACCARD) {
                        finalscore = finalscore / (2 * (nneighbors + 1) - finalscore);
                    }
                }
 
                current_edges.push_back(j);
                current_edges.push_back(othernode);
                current_weights.push_back(std::max(finalscore, 1e-6)); // Ensuring that an edge with a positive weight is always reported.
 
                // 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 += 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(
        static_cast<Index_>(neighbors.size()), 
        [&](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(
        static_cast<Index_>(neighbors.size()),
        [&](Index_ i) -> const std::vector<Index_>& { return neighbors[i]; }, 
        [](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) {
    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(
    std::size_t num_dims, 
    Index_ num_cells, 
    const Input_* data, 
    const knncolle::Builder<Index_, Input_, Distance_, Matrix_>& knn_method,
    const BuildSnnGraphOptions& options) 
{
    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