random.h

Go to the documentation of this file.

#ifndef FPMAS_RANDOM_H
#define FPMAS_RANDOM_H
 
#include "generator.h"
#include "distribution.h"
#include <numeric>
 
namespace fpmas { namespace random {
    extern const std::mt19937::result_type default_seed;
 
    extern fpmas::random::mt19937::result_type seed;
 
    extern DistributedGenerator<> rd_choices;
 
    template<typename T, typename Generator_t>
        std::vector<T> local_choices(
                const std::vector<T>& local_items, std::size_t n,
                Generator_t& gen) {
            random::UniformIntDistribution<std::size_t> local_random_item(
                    0, local_items.size()-1);
            std::vector<T> random_items(n);
            for(std::size_t i = 0; i < n; i++)
                random_items[i] = local_items[local_random_item(gen)];
            return random_items;
        }
 
    template<typename T>
        std::vector<T> local_choices(
                const std::vector<T>& local_items, std::size_t n) {
            return local_choices(local_items, n, rd_choices);
        }
 
    template<typename T>
        std::vector<std::vector<T>> distributed_choices(
                api::communication::MpiCommunicator& comm,
                const std::vector<T>& local_items, std::size_t n) {
            communication::TypedMpi<std::size_t> int_mpi(comm);
            std::vector<std::size_t> item_counts
                = int_mpi.allGather(local_items.size());
 
            std::unordered_map<int, std::size_t> requests;
 
            random::DiscreteDistribution<std::size_t> rd_process(item_counts);
            for(std::size_t i = 0; i < n; i++)
                requests[rd_process(rd_choices)]++;
            requests = int_mpi.allToAll(requests);
 
            std::unordered_map<int, std::vector<T>> random_items;
            for(auto item : requests)
                random_items[item.first]
                    = local_choices(local_items, item.second, rd_choices);
 
            communication::TypedMpi<std::vector<T>> item_mpi(comm);
 
            std::vector<std::vector<T>> results(comm.getSize());
            for(auto item : item_mpi.allToAll(random_items))
                results[item.first] = item.second;
 
            return results;
        }
 
    template<typename T>
        std::vector<T> split_choices(
                api::communication::MpiCommunicator& comm,
                const std::vector<T>& local_items, std::size_t n) {
            // A deterministic generator that produces the same sequence on all
            // processes
            static mt19937_64 gen(seed);
 
            communication::TypedMpi<std::size_t> int_mpi(comm);
            std::vector<std::size_t> item_counts
                = int_mpi.allGather(local_items.size());
 
            random::DiscreteDistribution<std::size_t> rd_process(item_counts);
 
            std::unordered_map<int, std::size_t> requests;
 
            // Determines how many items will be selected from this process, so
            // that the sum of local_items_count on all processes is equal to
            // n.
            std::size_t local_items_count = 0;
            for(std::size_t i = 0; i < n; i++) {
                // Same result on all processes
                int process = rd_process(gen);
                if(process == comm.getRank())
                    local_items_count++;
            }
 
            return local_choices(local_items, local_items_count);
        }
    template<typename T, typename Generator_t>
        std::vector<T> local_sample(
                const std::vector<T>& local_items, std::size_t n,
                Generator_t& gen) {
            std::vector<T> sample(std::min(local_items.size(), n));
 
            for(std::size_t i = 0; i < sample.size(); i++)
                sample[i] = local_items[i];
 
 
            for(std::size_t i = sample.size(); i < local_items.size(); i++) {
                UniformIntDistribution<std::size_t> rd_index(0, i);
                std::size_t k = rd_index(gen);
                if(k < sample.size())
                    sample[k] = local_items[i];
            }
            return sample;
        }
 
    template<typename T>
        std::vector<T> local_sample(
                const std::vector<T>& local_items, std::size_t n
                ) {
            return local_sample(local_items, n, rd_choices);
        }
 
    
    template<typename Index_t, typename Generator_t>
        std::vector<Index_t> sample_indexes(
                Index_t begin, Index_t end,
                std::size_t n,
                Generator_t& gen
                ) {
            // (process, local_index)
            std::vector<Index_t> result_indexes(n);
 
            /* Initializes the indexes reservoir */
            Index_t p = begin;
            for(std::size_t i = 0; i < n; i++) {
                result_indexes[i] = p++;
            }
            /**************************************/
 
            std::size_t i = n;
            while(p != end) {
                UniformIntDistribution<std::size_t> random_k(0, i);
                std::size_t k = random_k(gen);
                if(k < n)
                    result_indexes[k] = p;
                ++i;
                ++p;
            }
 
            return result_indexes;
        }
 
    
    template<typename K>
    class Index {
        private:
            const std::map<K, std::size_t>* item_counts;
            typename std::map<K, std::size_t>::const_iterator it;
            std::size_t _offset;
 
        public:
            Index() = default;
 
            Index(const std::map<K, std::size_t>* item_counts)
                : item_counts(item_counts) {
                }
            Index(
                    const std::map<K, std::size_t>* item_counts,
                    typename std::map<K, std::size_t>::const_iterator it
                    )
                : item_counts(item_counts), it(it), _offset(0) {
                }
 
            Index(
                    const std::map<K, std::size_t>* item_counts,
                    K key,
                    std::size_t offset
                    ) :
                item_counts(item_counts),
                it(this->item_counts->find(key)),
                _offset(offset) {
                }
 
            K key() const {
                return it->first;
            };
            std::size_t offset() const {
                return _offset;
            }
 
            static Index begin(const std::map<K, std::size_t>* item_counts);
 
            static Index end(const std::map<K, std::size_t>* item_counts);
 
            Index& operator++();
 
            Index operator++(int);
 
            Index operator+(std::size_t n) const;
            
            static std::size_t distance(
                    const Index& i1, const Index& i2);
 
            bool operator==(const Index& index) const;
            bool operator!=(const Index& index) const;
 
            bool operator<(const Index& index) const;
    };
 
    template<typename K>
        Index<K> Index<K>::begin(const std::map<K, std::size_t>* item_counts) { 
            Index<K> _begin(item_counts, item_counts->begin());
            while(_begin.it != item_counts->end() && _begin.it->second == 0)
                ++_begin.it;
            return _begin;
        }
 
    template<typename K>
        Index<K> Index<K>::end(const std::map<K, std::size_t>* item_counts) {
            Index<K> _end(item_counts, item_counts->end());
            return _end;
        }
 
    template<typename K>
        Index<K>& Index<K>::operator++() {
            if(it != item_counts->end()) {
                if(_offset >= it->second-1) {
                    _offset=0;
                    ++it;
                    // Skips processes without items
                    while(it != item_counts->end() && it->second == 0)
                        ++it;
                } else {
                    ++_offset;
                }
            }
            return *this;
        }
 
    template<typename K>
        Index<K> Index<K>::operator++(int) {
            Index<K> index = *this;
            ++*this;
            return index;
        }
 
    template<typename K>
    Index<K> Index<K>::operator+(std::size_t n) const {
        Index<K> index = *this;
        // index.id->second <=> index.item_counts[index.key()]
        if(index._offset+n < index.it->second) {
            // Stay in current process, nothing to do
            index._offset+=n;
        } else {
            // Need to go to next process
            n -= index.it->second - index.offset();
            index._offset = 0;
            ++index.it;
            while(n >= index.it->second) {
                n-= index.it->second;
                ++index.it;
            }
            index._offset = n;
        }
        return index;
    }
    
    template<typename K>
    std::size_t Index<K>::distance(
            const Index<K>& i1, const Index<K>& i2) {
        if(i1.it != i1.item_counts->end() && i2.it != i2.item_counts->end()) {
            if(i1.key() == i2.key()) {
                return i2.offset() - i1.offset();
            } else {
                std::size_t result = i1.it->second - i1.offset();
                auto it = i1.it;
                ++it;
                while(it->first != i2.it->first) {
                    result+=it->second;
                    ++it;
                }
                result+=i2.offset();
                return result;
            }
        }
        if(i1.it == i1.item_counts->end() && i2.it == i2.item_counts->end())
            return 0;
        // Since i2 reachable from i1, only i2 can be end (if i1 is end and not
        // i2, unexpected result)
 
        std::size_t result = i1.it->second - i1.offset();
        auto it = i1.it;
        ++it;
        while(it != i1.item_counts->end()) {
            result+=it->second;
            ++it;
        }
        return result;
    }
 
    template<typename K>
        bool Index<K>::operator==(const Index<K>& index) const {
            if(it == item_counts->end() && index.it == index.item_counts->end())
                return true;
            if(it != item_counts->end() && index.it != index.item_counts->end())
                return key() == index.key() && offset() == index.offset();
            // One is end and other is not
            return false;
        }
 
    template<typename K>
        bool Index<K>::operator!=(const Index<K>& index) const {
            return !(*this==index);
        }
 
    template<typename K>
        bool Index<K>::operator<(const Index<K>& index) const {
            if(this->it == this->item_counts->end()
                    && index.it != index.item_counts->end())
                // This is end, index is not
                return false;
            if(index.it == index.item_counts->end()
                    && this->it != this->item_counts->end())
                // Index is end, this is not
                return true;
            if(index.it == index.item_counts->end()
                    && this->it == this->item_counts->end())
                // Both end (ie equal)
                return false;
            // Index and this are not end
            if(this->key() == index.key())
                return this->offset() < index.offset();
            return this->key() < index.key();
        }
 
    typedef Index<int> DistributedIndex;
 
    template<typename T>
        std::vector<std::vector<T>> distributed_sample(
                api::communication::MpiCommunicator &comm,
                const std::vector<T> &local_items, std::size_t n
                ) {
 
            communication::TypedMpi<std::size_t> int_mpi(comm);
            std::vector<std::size_t> item_counts
                = int_mpi.allGather(local_items.size());
            std::map<int, std::size_t> item_counts_map;
            for(std::size_t i = 0; i < item_counts.size(); i++)
                item_counts_map[i] = item_counts[i];
 
            // index = (process, local_offset)
            std::vector<DistributedIndex> result_indexes
                = sample_indexes(
                        DistributedIndex::begin(&item_counts_map),
                        DistributedIndex::end(&item_counts_map),
                        std::min(
                            n,
                            std::accumulate(item_counts.begin(), item_counts.end(), 0ul)
                            ), rd_choices);
            
            // Requests required items
            communication::TypedMpi<std::vector<std::size_t>> request_mpi(comm);
            std::unordered_map<int, std::vector<std::size_t>> requests;
            for(auto& item : result_indexes)
                requests[item.key()].push_back(item.offset());
            requests = request_mpi.allToAll(requests);
 
            // Exchange items
            communication::TypedMpi<std::vector<T>> item_mpi(comm);
            std::unordered_map<int, std::vector<T>> items;
            for(auto item : requests)
                for(auto index : item.second)
                    items[item.first].push_back(local_items[index]);
            items = item_mpi.allToAll(items);
 
            // Result transformation
            std::vector<std::vector<T>> results(comm.getSize());
            for(auto item : items)
                results[item.first] = item.second;
            return results;
        }
 
    template<typename T>
        std::vector<T> split_sample(
                api::communication::MpiCommunicator &comm,
                const std::vector<T> &local_items, std::size_t n
                ) {
            // A deterministic generator that produces the same sequence on all
            // processes
            static mt19937_64 gen(seed);
 
            communication::TypedMpi<std::size_t> int_mpi(comm);
            std::vector<std::size_t> item_counts
                = int_mpi.allGather(local_items.size());
            std::map<int, std::size_t> item_counts_map;
            for(std::size_t i = 0; i < item_counts.size(); i++)
                item_counts_map[i] = item_counts[i];
 
            // index = (process, local_offset)
            std::vector<DistributedIndex> result_indexes
                = sample_indexes(
                        DistributedIndex::begin(&item_counts_map),
                        DistributedIndex::end(&item_counts_map),
                        std::min(
                            n,
                            std::accumulate(item_counts.begin(), item_counts.end(), 0ul)
                            ), gen);
 
            // Because gen is sequential, the same result_indexes list is built
            // on all processes
            
            std::vector<T> result;
            for(auto& item : result_indexes)
                if(item.key() == comm.getRank())
                    result.push_back(local_items[item.offset()]);
            return result;
        }
}}
#endif