BuildSnnGraph.hpp

Go to the documentation of this file.

#ifndef SCRAN_BUILDSNNGRAPH_HPP
#define SCRAN_BUILDSNNGRAPH_HPP
 
#include "../utils/macros.hpp"
 
#if __has_include("igraph.h")
#include "igraph.h"
#include "igraph_utils.hpp"
#endif
 
#include <vector>
#include <algorithm>
#include <memory>
 
#include "knncolle/knncolle.hpp"
 
namespace scran {
 
class BuildSnnGraph {
public:
    enum Scheme { RANKED, NUMBER, JACCARD };
 
    struct Defaults {
        static constexpr int neighbors = 10;
 
        static constexpr Scheme weighting_scheme = RANKED;
 
        static constexpr bool approximate = false;
 
        static constexpr int num_threads = 1;
    };
private:
    int num_neighbors = Defaults::neighbors;
    Scheme weight_scheme = Defaults::weighting_scheme;
    bool approximate = Defaults::approximate;
    int nthreads = Defaults::num_threads;
 
public:
    BuildSnnGraph& set_neighbors(int k = Defaults::neighbors) {
        num_neighbors = k;
        return *this;
    }
 
    BuildSnnGraph& set_approximate(bool a = Defaults::approximate) {
        approximate = a;
        return *this;
    }
 
    BuildSnnGraph& set_weighting_scheme(Scheme w = Defaults::weighting_scheme) {
        weight_scheme = w;
        return *this;
    }
 
    BuildSnnGraph& set_num_threads(int n = Defaults::num_threads) {
        nthreads = n;
        return *this;
    }
 
public:
    struct Results {
        size_t ncells;
 
#if __has_include("igraph.h")
        std::vector<igraph_integer_t> edges;
#else
        std::vector<int> edges;
#endif
 
#if __has_include("igraph.h")
        std::vector<igraph_real_t> weights;
#else
        std::vector<double> weights;
#endif
 
#if __has_include("igraph.h")
        igraph::Graph to_igraph() const {
            igraph_vector_int_t edge_view;
            igraph_vector_int_view(&edge_view, edges.data(), edges.size());
            return igraph::Graph(&edge_view, ncells, /* directed = */ false);
        }
#endif
    };
 
private:
    template<class Indices_, class Function_>
    Results run(const std::vector<Indices_>& indices, Function_ get_index) const {
        size_t ncells = indices.size();
 
        // Not parallel-frendly, so we don't construct this with the neighbor search
        std::vector<std::vector<std::pair<int, int> > > hosts(ncells);
        for (size_t i = 0; i < ncells; ++i) {
            hosts[i].push_back(std::make_pair(0, i)); // each point is its own 0-th nearest neighbor
            const auto& current = indices[i];
            int counter = 1;
            for (auto x : current) {
                hosts[get_index(x)].push_back(std::make_pair(counter, i));
                ++counter;
            }
        }
 
        // Constructing the shared neighbor graph.
        std::vector<std::vector<int> > edge_stores(ncells);
        std::vector<std::vector<double> > weight_stores(ncells);
 
#ifndef SCRAN_CUSTOM_PARALLEL
        #pragma omp parallel num_threads(nthreads)
        {
#else
        SCRAN_CUSTOM_PARALLEL([&](size_t, size_t start, size_t length) -> void {
#endif
 
            std::vector<int> current_score(ncells);
            std::vector<int> current_added;
            current_added.reserve(ncells);
 
#ifndef SCRAN_CUSTOM_PARALLEL
            #pragma omp for
            for (size_t j = 0; j < ncells; ++j) {
#else
            for (size_t j = start, end = start + length; j < end; ++j) {
#endif
 
                const auto& current_neighbors = indices[j];
                int nneighbors = current_neighbors.size();
 
                for (int i = 0; i <= nneighbors; ++i) {
                    // First iteration treats 'j' as the zero-th neighbor.
                    // Remaining iterations go through the neighbors of 'j'.
                    const int 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'.
                    for (const auto& h : hosts[cur_neighbor]) {
                        auto othernode = h.second;
 
                        if (othernode < j) { // avoid duplicates from symmetry in the SNN calculations.
                            int& existing_other = current_score[othernode];
                            if (weight_scheme == RANKED) {
                                // Recording the lowest combined rank per neighbor.
                                int currank = h.first + i;
                                if (existing_other == 0) { 
                                    existing_other = currank;
                                    current_added.push_back(othernode);
                                } else if (existing_other > currank) {
                                    existing_other = currank;
                                }
                            } else {
                                // 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() * 2);
 
                for (auto othernode : current_added) {
                    int& otherscore = current_score[othernode];
                    double finalscore;
                    if (weight_scheme == RANKED) {
                        finalscore = static_cast<double>(nneighbors) - 0.5 * static_cast<double>(otherscore);
                    } else {
                        finalscore = otherscore;
                        if (weight_scheme == 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();
 
#ifndef SCRAN_CUSTOM_PARALLEL
            }
        } 
#else
            }
        }, ncells, nthreads);
#endif
 
        // Collating the total number of edges.
        size_t nedges = 0;
        for (const auto& w : weight_stores) {
            nedges += w.size();
        }
 
        Results output;
        output.ncells = ncells;
 
        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.reserve(nedges * 2);
        for (const auto& e : edge_stores) {
            output.edges.insert(output.edges.end(), e.begin(), e.end());
        }
 
        return output;
    }
 
public:
    Results run(size_t ndims, size_t ncells, const double* mat) const {
        std::unique_ptr<knncolle::Base<> > ptr;
        if (!approximate) {
            ptr.reset(new knncolle::VpTreeEuclidean<>(ndims, ncells, mat));
        } else {
            ptr.reset(new knncolle::AnnoyEuclidean<>(ndims, ncells, mat));
        }
        return run(ptr.get());
    }
 
    template<class Algorithm>
    Results run(const Algorithm* search) const {
        // Collecting neighbors.
        size_t ncells = search->nobs();
        std::vector<std::vector<int> > indices(ncells);
 
#ifndef SCRAN_CUSTOM_PARALLEL
        #pragma omp parallel for num_threads(nthreads)
        for (size_t i = 0; i < ncells; ++i) {
#else
        SCRAN_CUSTOM_PARALLEL([&](size_t, size_t start, size_t length) -> void {
        for (size_t i = start, end = start + length; i < end; ++i) {
#endif
            auto neighbors = search->find_nearest_neighbors(i, num_neighbors);
            auto& current = indices[i];
            for (auto x : neighbors) {
                current.push_back(x.first);
            }
        }
#ifdef SCRAN_CUSTOM_PARALLEL
        }, ncells, nthreads);
#endif
 
        return run(indices);
    }
 
    template<typename Index_, typename Distance_>
    Results run(const knncolle::NeighborList<Index_, Distance_>& neighbors) const {
        return run(neighbors, [](const std::pair<Index_, Distance_>& x) -> Index_ { return x.first; });
    }
 
    template<class Indices_>
    Results run(const std::vector<Indices_>& indices) const {
        typedef typename std::remove_cv<typename std::remove_reference<decltype(std::declval<Indices_>()[0])>::type>::type Element;
        return run(indices, [](Element i) -> Element { return i; });
    }
};
 
}
 
#endif