choose_classic_markers.hpp

Go to the documentation of this file.

#ifndef SINGLEPP_CHOOSE_CLASSIC_MARKERS_HPP
#define SINGLEPP_CHOOSE_CLASSIC_MARKERS_HPP
 
#include "defs.hpp"
 
#include "tatami/tatami.hpp"
 
#include "Markers.hpp"
 
#include <vector>
#include <cmath>
#include <stdexcept>
#include <algorithm>
#include <set>
#include <cstddef>
 
namespace singlepp {
 
inline std::size_t number_of_classic_markers(std::size_t num_labels) {
    return std::round(500.0 * std::pow(2.0/3.0, std::log(static_cast<double>(num_labels)) / std::log(2.0)));
}
 
struct ChooseClassicMarkersOptions {
    int number = -1;
 
    int num_threads = 1;
};
 
template<typename Value_, typename Index_, typename Label_>
Markers<Index_> choose_classic_markers(
    const std::vector<const tatami::Matrix<Value_, Index_>*>& representatives,
    const std::vector<const Label_*>& labels,
    const ChooseClassicMarkersOptions& options)
{
    auto nrefs = representatives.size();
    if (nrefs != labels.size()) {
        throw std::runtime_error("'representatives' and 'labels' should have the same length");
    }
    if (nrefs == 0) {
        throw std::runtime_error("'representatives' should contain at least one entry");
    }
    auto ngenes = representatives.front()->nrow();
 
    std::size_t nlabels = 0;
    for (decltype(nrefs) r = 0; r < nrefs; ++r) {
        const auto& current = *representatives[r];
 
        auto nrows = current.nrow();
        if (nrows != ngenes) {
            throw std::runtime_error("all entries of 'representatives' should have the same number of rows");
        }
 
        auto ncols = current.ncol();
        if (ncols) {
            auto curlab = labels[r];
            nlabels = std::max(nlabels, static_cast<std::size_t>(*std::max_element(curlab, curlab + ncols)) + 1);
        }
    }
 
    // Generating mappings.
    std::vector<std::vector<std::pair<bool, Index_> > > labels_to_index(nrefs);
    for (decltype(nrefs) r = 0; r < nrefs; ++r) {
        auto& current = labels_to_index[r];
        current.resize(nlabels);
        auto ncols = representatives[r]->ncol();
        auto curlab = labels[r];
 
        for (decltype(ncols) c = 0; c < ncols; ++c) {
            auto& info = current[curlab[c]];
            if (info.first) {
                throw std::runtime_error("each label should correspond to no more than one column in each reference");
            }
            info.first = 1;
            info.second = c;
        }
    }
 
    Markers<Index_> output(nlabels);
    for (auto& x : output) {
        x.resize(nlabels);
    }
 
    int actual_number = options.number;
    if (actual_number < 0) {
        actual_number = number_of_classic_markers(nlabels);
    } 
    if (actual_number > static_cast<int>(ngenes)) {
        actual_number = ngenes;
    }
 
    // Generating pairs for compute; this sacrifices some memory for convenience.
    std::vector<std::pair<Label_, Label_> > pairs;
    {
        std::set<std::pair<Label_, Label_> > pairs0;
        for (decltype(nrefs) r = 0; r < nrefs; ++r) {
            auto ncols = representatives[r]->ncol();
            auto curlab = labels[r];
            for (decltype(ncols) c1 = 0; c1 < ncols; ++c1) {
                for (decltype(c1) c2 = 0; c2 < c1; ++c2) {
                    pairs0.emplace(curlab[c1], curlab[c2]);
                }
            }
        }
        pairs.insert(pairs.end(), pairs0.begin(), pairs0.end()); // already sorted by the std::set.
    }
 
    auto npairs = pairs.size();
    tatami::parallelize([&](int, decltype(npairs) start, decltype(npairs) len) {
        std::vector<std::pair<Value_, Index_> > sorter(ngenes);
        std::vector<Value_> rbuffer(ngenes), lbuffer(ngenes);
        std::vector<std::shared_ptr<tatami::MyopicDenseExtractor<Value_, Index_> > > rextractors(nrefs), lextractors(nrefs);
 
        for (decltype(npairs) p = start, end = start + len; p < end; ++p) {
            auto curleft = pairs[p].first;
            auto curright = pairs[p].second;
 
            for (decltype(ngenes) g = 0; g < ngenes; ++g) {
                sorter[g].first = 0;
                sorter[g].second = g;
            }
 
            for (decltype(nrefs) i = 0; i < nrefs; ++i) {
                const auto& curavail = labels_to_index[i];
                auto lcol = curavail[curleft];
                auto rcol = curavail[curright];
                if (!lcol.first || !rcol.first) {
                    continue;                            
                }
 
                // Initialize extractors as needed.
                auto& lext = lextractors[i];
                if (!lext) {
                    lext = representatives[i]->dense_column();
                }
                auto lptr = lext->fetch(lcol.second, lbuffer.data());
 
                auto& rext = rextractors[i];
                if (!rext) {
                    rext = representatives[i]->dense_column();
                }
                auto rptr = rext->fetch(rcol.second, rbuffer.data());
 
                for (decltype(ngenes) g = 0; g < ngenes; ++g) {
                    sorter[g].first += lptr[g] - rptr[g]; 
                }
            }
 
            // At flip = 0, we're looking for genes upregulated in right over left,
            // as we're sorting on left - right in increasing order. At flip = 1,
            // we reverse the signs to we sort on right - left.
            for (int flip = 0; flip < 2; ++flip) {
                if (flip) {
                    for (auto& s : sorter) {
                        s.first *= -1;
                    }
                }
 
                // partial sort is guaranteed to be stable due to the second index resolving ties.
                std::partial_sort(sorter.begin(), sorter.begin() + actual_number, sorter.end());
 
                std::vector<Index_> stuff;
                stuff.reserve(actual_number);
                for (int g = 0; g < actual_number && sorter[g].first < 0; ++g) { // only negative values (i.e., positive log-fold changes for our comparison).
                    stuff.push_back(sorter[g].second); 
                }
 
                if (flip) {
                    output[curleft][curright] = std::move(stuff);
                } else {
                    output[curright][curleft] = std::move(stuff);
                }
            }
        }
    }, npairs, options.num_threads);
 
    return output;
}
 
// Overload for the non-const case.
template<typename Value_, typename Index_, typename Label_>
Markers<Index_> choose_classic_markers(
    const std::vector<tatami::Matrix<Value_, Index_>*>& representatives,
    const std::vector<const Label_*>& labels,
    const ChooseClassicMarkersOptions& options)
{
    std::vector<const tatami::Matrix<Value_, Index_>*> rep2(representatives.begin(), representatives.end());
    return choose_classic_markers(rep2, labels, options);
}
template<typename Value_, typename Index_, typename Label_>
Markers<Index_> choose_classic_markers(const tatami::Matrix<Value_, Index_>& representative, const Label_* labels, const ChooseClassicMarkersOptions& options) {
    return choose_classic_markers(std::vector<const tatami::Matrix<Value_, Index_>*>{ &representative }, std::vector<const Label_*>{ labels }, options);
}
 
}
 
#endif