Algorithms solutions: September 2016

Due to the fact that the Fibonacci Heap was created for when the use of decrease key is greater relative to other operations, Fibonacci Heaps in theory are useful for really big and dense graphs, while a D-ary Heap may be more fitting a for a sparse graph. There running times reflect this:

Disktra's algorithm with D-ary Heap: $Ο (E \log_{D} (V))$
Disktra's algorithm with Fibonacci Heap : $Ο (V \log_{2} (V) + E)$

In practice thought, a D-ary Heap performs better than a Fibonacci heap in most cases. Also a D-ary -heap with D greater than 2 may perform better because it produces less cache misses and virtual memory page faults.

Here is my implementation (maybe later I'll add another post with some benchmarks):

Fibonacci Heap:

#ifndef INDEX_FIB_HEAP_FIB_HEAP_H
#define INDEX_FIB_HEAP_FIB_HEAP_H

#include <algorithm>
#include <stdexcept>
#include <functional>

namespace mds {

    template<typename P, typename  Func = std::less<P>>
    class index_fib_heap {

        struct node {

            unsigned index;
            P prty;

            unsigned degree = 0;
            bool mark       = false;
            node* parent    = nullptr;
            node* childs    = nullptr;
            node* left      = nullptr;
            node* right     = nullptr;


            node(node* n) :
                prty(n->prty),
                index(n->index),
                degree(n->degree),
                mark(n->mark) {
            
            }

            node(unsigned pindex, const P& priority)
                : index(pindex),
                prty(priority) { }

            node(unsigned pindex, P&& priority)
                : index(pindex),
                prty(std::move(priority)) { }

            node* add_brother(node* n) {
                n->parent   = parent;
                n->left     = this;
                n->right    = right;
                right->left = n;
                right       = n;
                return this;
            }
            node* remove_self() {
                left->right = right;
                right->left = left;
                return this;
            }

            node* add_child(node* n) {
                if (childs == nullptr) {
                    childs    = n->to_single_list();
                    n->parent = this;
                }
                else {
                    childs->add_brother(n);
                }
                n->mark = false;
                ++degree;
                return this;
            }

            node* to_single_list() {
                left  = this;
                right = this;
                return this;
            }
            void destroy() {
                auto n = childs;
                for (auto i = 0u; i < degree; ++i) {
                    auto next = n->right;
                    n->destroy();
                    n = next;
                }
                delete this;
            }

        };

        node* copy_node(node* n) {
            auto new_n = new node(n);
            index_table[n->index] = new_n;
            return new_n;
        }

        node* deep_copy(node* root) {
            auto copy = copy_node(root);
            if (root->childs == nullptr) {
                return copy;
            }

            auto n = root->childs->right;
            auto child_copy = deep_copy(root->childs)
                              ->to_single_list();

            copy->childs = child_copy;
            child_copy->parent = copy;
            for (auto i = 1u; i < root->degree; ++i) {
                child_copy->add_brother(deep_copy(n));
                n = n->right;
            }
            return copy;
        }

        void consolidate() {
            unsigned bound = static_cast<unsigned>(
                            std::log(N) / LOG_OF_GOLDEN) + 1;
            auto degree_array = new node*[bound]();
            for (auto n = min; root_size > 0; --root_size) {
                auto parent = n;
                auto d = parent->degree;
                n = n->right;
                while (degree_array[d] != nullptr) {
                    auto child = degree_array[d];
                    if (cmp(child->prty, parent->prty)) {
                        std::swap(child, parent);
                    }
                    parent->add_child(child->remove_self());
                    degree_array[d++] = nullptr;
                }
                degree_array[d] = parent;
            }

            auto d = 0u;
            while (degree_array[d++] == nullptr);

            min = degree_array[d - 1]->to_single_list();
            for (; d < bound; ++d) {
                if (degree_array[d] != nullptr) {
                    add_to_root(degree_array[d]);
                }
            }
            delete[] degree_array;
            ++root_size;
        }

        node* remove_min() {
            if (empty()) {
                throw std::underflow_error("underflow");
            }

            auto deleted = min;
            add_childs_to_root(deleted);
            deleted->remove_self();
            --root_size;

            if (deleted == deleted->right) {
                min = nullptr;
            }
            else {
                min = min->right;
                consolidate();
            }
            --N;
            return deleted;
        }

        void add_to_root(node* n) {
            min->add_brother(n);
            if (cmp(n->prty, min->prty)) {
                min = n;
            }
            ++root_size;
        }

        void add_childs_to_root(node* n) {
            auto c = n->childs;
            auto d = n->degree;
            for (auto i = 0u; i < d; ++i) {
                auto next = c->right;
                add_to_root(c);
                c = next;
            }
        }

        template<typename PP>
        void decrease_priority(node *n, PP&& new_p) {
            if (cmp(n->prty, new_p)) {
                throw std::runtime_error("key is greater");
            }

            n->prty = std::forward<PP>(new_p);
            auto p = n->parent;

            if (p != nullptr 
                  && cmp(n->prty, p->prty)) {
                cut(n, p);
                cascading_cut(p);
            }
            if (cmp(n->prty, min->prty)) {
                min = n;
            }
        }
        void cut(node* child, node* parent) {
            child->remove_self();
            child->mark = false;
            --parent->degree;
            if (parent->degree == 0) {
                parent->childs = nullptr;
            }
            else if (parent->childs == child) {
                parent->childs = child->right;
            }

            min->add_brother(child);
            ++root_size;
        }

        void cascading_cut(node* n) {
            while (n->parent != nullptr && n->mark) {

                node* p = n->parent;
                cut(n, p);
                n = p;
            }
            n->mark = n->mark || n->parent;
        }

        void check_index(unsigned index) const {
            if (index >= max_size) {
                std::out_of_range("index out of range");
            }
        }

        static constexpr  double LOG_OF_GOLDEN = 0.4812118;

        unsigned max_size;
        node**   index_table;
        node*    min;
        unsigned N;
        unsigned root_size;
        Func     cmp;

    public:

        index_fib_heap(unsigned pmax_size, Func pcmp)
            : max_size(pmax_size),
            index_table(new node*[pmax_size]()),
            min(nullptr),
            N(0),
            root_size(0),
            cmp(pcmp) { }

        index_fib_heap(unsigned pmax_size)
            : index_fib_heap(pmax_size, Func()) { }


        index_fib_heap(index_fib_heap<P, Func> && fib) noexcept :
            min(fib.min),
            N(fib.N),
            root_size(fib.root_size),
            cmp(fib.cmp),
            index_table(fib.index_table),
            max_size(fib.max_size)
        {  
            fib.N = fib.max_size = fib.root_size = 0;
            fib.index_table = nullptr;
            fib.min = nullptr;
        }

        index_fib_heap(const index_fib_heap<P, Func> & fib) noexcept 
            : N(fib.N),
              root_size(fib.root_size),
              cmp(fib.cmp),
              max_size(fib.max_size),
              index_table(new node*[fib.max_size]()) 
        {
            if (fib.empty()) { 
                return; 
            }
            min    = deep_copy(fib.min)->to_single_list();
            auto n = fib.min->right;
            for (auto i = 1u; i < root_size; ++i) {
                min->add_brother(deep_copy(n));
                n = n->right;
            }
        }

        
        ~index_fib_heap() {
            auto n = min;
            for (auto i = 0u; i < root_size; ++i) {
                auto next = n->right;
                n->destroy();
                n = next;
            }
            delete[] index_table;
        }


        void pop() {
            auto min_node = remove_min();
            index_table[min_node->index] = nullptr;
            delete min_node;
        }

        std::pair<unsigned, P> top() const {
            return { min->index, min->prty };
        }

        template<typename PP>
        bool insert(unsigned index, PP&& prty) {
            check_index(index);
            auto &new_node = index_table[index];
            if (new_node != nullptr) {
                return false;
            }
            new_node = new node(index,
                std::forward<PP>(prty));
            if (min == nullptr) {
                min = new_node->to_single_list();
                ++root_size;
            }
            else {
                add_to_root(new_node);
            }
            ++N;
            return true;
        }

        template<typename PP>
        void decrease(unsigned index, PP&& prty) {
            check_index(index);
            auto node = index_table[index];
            if (node != nullptr) {
                decrease_priority(node,
                    std::forward<PP>(prty));
            }
        }

        bool empty() const {
            return min == nullptr;
        }

        unsigned size() const {
            return N;
        }

        bool contains(unsigned index) const {
            return index_table[index] != nullptr;
        }

        friend void swap(index_fib_heap<P, Func>& first, index_fib_heap<P, Func>& second) noexcept {
            using std::swap;
            swap(first.N, second.N);
            swap(first.max_size, second.max_size);
            swap(first.index_table, second.index_table);
            swap(first.min, second.min);
            swap(first.root_size, second.root_size);
            swap(first.cmp, second.cmp);
        }

        index_fib_heap<P, Func>& operator=(index_fib_heap<P, Func> other)  noexcept {
            swap(*this, other);
            return *this;
        }

    };
}
#endif //INDEX_FIB_HEAP_FIB_HEAP_H

D-ary Heap

#ifndef D_HEAP_D_HEAP_H
#define D_HEAP_D_HEAP_H
#include <vector>
#include <iostream>
#include <algorithm>
#include <climits>

namespace mds {
   
    template<class T, class F = std::less<T>>
    class DHeap
    {
        const static std::size_t defualt_D = 4;
    public:
   
        DHeap(std::size_t d, std::size_t pmax_size)
            : DHeap(d, pmax_size, F()) { }



        DHeap(std::size_t d, std::size_t pmax_N, const F& func)
            : D(d), cmp(func), 
              N(0), max_N(pmax_N), 
              keys(new T[pmax_N]), 
              key_pos(new std::size_t[pmax_N]), 
              pq(new std::size_t[pmax_N]) { 
            
            for (auto i = 0u; i < max_N; ++i) {
                key_pos[i] = SIZE_MAX;
            }
        }

        DHeap(std::size_t pmax_N, const F& func)
            : DHeap(defualt_D, pmax_N, func) { }

        DHeap(std::size_t pmax_N)
            : DHeap(pmax_N, F()) { }

        DHeap(DHeap<T, F> && heap) noexcept
            : cmp(heap.cmp), 
              keys(heap.keys), 
              key_pos(heap.key_pos), 
              pq(heap.pq),
              D(heap.D),
              N(heap.N),
              max_N(heap.max_N)
        {
            heap.keys = nullptr;
            heap.key_pos = nullptr;
            heap.pq = nullptr;
            heap.N = heap.max_N = 0;
        }

        

        ~DHeap() {
            delete[] keys;
            delete[] key_pos;
            delete[] pq;
        }


        void pop() {
            if (empty()) {
                throw new std::underflow_error("underflow");
            }
            exchange(0, --N);
            key_pos[pq[N]] = SIZE_MAX;
            heapify(0);
        } 

        bool insert(std::size_t index, const T& element) {
            return genericInsert(index, element);
        }
        
        bool insert(std::size_t index, T&& element) {
            return genericInsert(index, std::move(element));
        }

        void decrease(std::size_t index, const T& key) {
            decrease_template(index, key);
        }

        void decrease(std::size_t index, T&& key) {
            decrease_template(index, std::move(key));
        }

        bool empty() const {
            return N == 0;
        }

        std::pair<unsigned, T> top() const {
            if (empty()) {
                throw new std::underflow_error("underflow");
            }
            return { pq[0], keys[pq[0]] };
        }

    private:

        
        F           cmp;
        T           *keys;
        std::size_t *key_pos;
        std::size_t *pq;
        std::size_t D;
        std::size_t N;
        std::size_t max_N;

        void heapify(std::size_t index) {
            std::size_t first = 0;
            while ((first = firstChild(index)) < N) {
                unsigned to = std::min(N, firstChild(index + 1));
                unsigned min = find_min(first, to);

                if (!compare(min, index)) {
                    break; 
                }
                exchange(min, index);
                index = min;
            }
        }
        
        void swim(std::size_t index) {
            std::size_t p = 0;
            while (index > 0 
                   && compare(index, p = parent(index))) {
                exchange(p, index);
                index = p;
            }
        }

        bool compare(std::size_t i, std::size_t j) const {
            return cmp(keys[pq[i]], keys[pq[j]]);
        }

        void exchange(std::size_t i, std::size_t j) {

            using std::swap;
            swap(key_pos[pq[i]], key_pos[pq[j]]);
            swap(pq[i], pq[j]);
        }


        unsigned firstChild(std::size_t index) const {
            return (D * index) + 1;
        }
        std::size_t parent(std::size_t index) const {
            return (index - 1) / D;
        }
        std::size_t find_min(std::size_t first, std::size_t to) const {
            auto min = first;
            for (auto i = first + 1; i < to; ++i) {
                min = compare(i, min) ? i : min;
            }
            return min;
        }
        template<typename U>
        bool genericInsert(std::size_t index, U&& element) {
            check(index);
            if (contains(index)) {
                return false;
            }

            keys[index]    = std::forward<U>(element);
            key_pos[index] = N;
            pq[N]          = index;
            
            swim(N++);
            return true;
        }

        template<typename TT>
        void decrease_template(std::size_t index, TT&& key) {
            check(index);
            if (cmp(keys[index], key)) {
                throw new std::invalid_argument("invalid argument");
            }
            keys[index] = std::forward<TT>(key);
            swim(key_pos[index]);
        }

        template< typename E>
        friend std::ostream & operator<<(std::ostream &os, const DHeap<T, E>& p)
        {
            for (auto = 0u; i < p.N; ++i) {
                os << p.table[i] << ",";
            }
            return os;
        }

        void check(std::size_t index) {
            if (index >= max_N) {
                throw std::out_of_range("index out of range");
            }
        }
        bool contains(std::size_t i) {          
            return key_pos[i] != SIZE_MAX;
        }
    };

}

#endif /*D_HEAP_D_HEAP_H*/

Graph:

#ifndef GRAPH_GRAPH_H
#define GRAPH_GRAPH_H
#include <cstddef>
#include <vector>


namespace mds {

    template <typename W, typename V = std::size_t>
    class edge {
        
        V v;
        V w;
        W wght;

    public:
       

        edge(const V& pfrom, const V& pto, const W& pweight)
            : v(pfrom), w(pto), wght(pweight) { }

        edge(V&& pfrom, V&& pto, W&& pweight)
            : v(std::move(pfrom)), 
              w(std::move(pto)), 
              wght(std::move(pweight)) { }

        V from() const { return v; }
        V to() const { return w; }
        const W& weight() const { return wght; }

        friend std::ostream & operator<<(std::ostream &os, const edge<W, V>& e)
        { 
            return os << e.from() << " --("
                      << e.weight()
                      << ")--> "
                      << e.to();
        }
    };

    template<typename W>
    class graph {

    public:
        using adj_structure = std::vector<edge<W>>;

        graph(std::size_t V)
            : adj_list(std::vector<adj_structure>(V, std::vector<edge<W>>())) { }
        
        graph(): graph(0) { }

        graph(graph<W>&& g) 
            : edges(g.edges), 
            adj_list(std::move(g.adj_list)) 
        {
            g.adj_list.clear(); 
            g.edges = 0;
        }

        graph(std::size_t V, std::initializer_list<edge<W>> list) 
            : graph(V) 
        {  
            for (const auto &e : list) {
                add_edge(e);
            }
        }

        graph(const graph<W>& g)
            : adj_list(g.adj_list), edge(g.edges)
        { }

        void add_edge(const edge<W>& e) {
            check(e.from());
            check(e.to());
            adj_list[e.from()].push_back(e);
            ++edges;
        }

        void add_edge(edge<W>&& e) {
            check(e.from());
            check(e.to());
            adj_list[e.from()].push_back(std::move(e));
            ++edges;
        }

        const adj_structure& adj(std::size_t v) const {
            return adj_list[v];
        }

        std::size_t E() const {
            return edges;
        }

        std::size_t V() const {
            return adj_list.size();
        }

        graph<W>& operator=(graph<W> other) {
            swap(*this, other);
            return *this;
        }
        
       
        friend void swap(graph<W>& first, graph<W>& second) noexcept {
            using std::swap;

            swap(first.adj_list, second.adj_list);
            swap(first.edges, second.edges);
        }

    private:
        std::vector<adj_structure>  adj_list;
        std::size_t edges;

        void check(std::size_t v) const {
            if (v >= V()) {
                throw std::out_of_range("vertex out of range");
            }
        }
    };

}
#endif // GRAPH_GRAPH_H

Dijkstra's algorithm:

#ifndef DIJKSTRA_ALGORITHM_DIJKSTRA_ALGORITHM_H
#define DIJKSTRA_ALGORITHM_DIJKSTRA_ALGORITHM_H
#include "Graph.h"
#include <limits>
#include "index_fib_heap.h"
#include <type_traits>
#include <functional>
#include <deque>
#include <climits>

namespace mds {

 
    struct identity {
        template<typename U>
        constexpr auto operator()(U&& v) const noexcept
            -> decltype(std::forward<U>(v)) {
            return std::forward<U>(v);
        }
    };   

    template<typename W, 
             typename F = identity, 
             typename IPQ = mds::index_fib_heap<
                                typename std::remove_reference<
                                    decltype(F()(std::declval<W>()))
                                >::type
                            >
    >
    class dijkstra_algorithm {
    
        using Weight = typename std::remove_reference<
            decltype(F()(std::declval<W>()))
        >::type;
        
       F               map_w;
       IPQ             pq;
       const graph<W>  *g;
       const edge<W>   **edge_to;
       Weight          *distances;
       
       std::pair<std::size_t, Weight> next() {
           if (pq.empty()) {
               return { SIZE_MAX, { } };
           }
           auto min = pq.top();
           pq.pop();
           for (const auto &e : g->adj(min.first)) {
               relax(&e);
           }
           return min;
       }

       void relax(const edge<W>* edge) {
           auto to    = edge->to();
           auto new_d = distances[edge->from()]
                        + map_w(edge->weight());

           if (new_d < distances[to]) {
               distances[to] = new_d;
               edge_to[to]   = edge;
               if (!pq.insert(to, new_d)) {
                   pq.decrease(to, new_d);
               }    
           }
       }

       void check(std::size_t v) const {
           if (g == nullptr || v >= g->V()) {
               throw std::out_of_range("vertex out of range");
           }
       }

    public:

       dijkstra_algorithm(std::size_t source, const graph<W>* pg, F pmap_w, IPQ ppq) 
           : g(pg), 
             map_w(pmap_w), 
             edge_to(new const edge<W>*[pg->V()]()),
             distances(new Weight[pg->V()]),
             pq(std::move(ppq)) 
       {       
           for (std::size_t i = 0; i < pg->V(); ++i) {
               distances[i] = std::numeric_limits<Weight>::max();
           }
           pq.insert(source, distances[source] = { });   
       }

       dijkstra_algorithm(std::size_t source, const graph<W>* pg, F pmap_w)
           : dijkstra_algorithm(source, pg, pmap_w, { pg->V() }) { }

       dijkstra_algorithm(std::size_t source, const graph<W>* pg)
           : dijkstra_algorithm(source, pg, { }) { }
   
       dijkstra_algorithm(dijkstra_algorithm<W, F, IPQ> && dijkstra) noexcept
          : map_w(std::move(dijkstra.map_w)), 
           pq(std::move(dijkstra.pq)), 
           g(dijkstra.g), 
           edge_to(dijkstra.edge_to), 
           distances(dijkstra.distances) 
       { 
           dijkstra.g         = nullptr;
           dijkstra.edge_to   = nullptr;
           dijkstra.distances = nullptr;
       }

       dijkstra_algorithm(const dijkstra_algorithm<W, F, IPQ>& dijkstra)
           : map_w(dijkstra.map_w),
             pq(dijkstra.pq),
             g(dijkstra.g) 
       {
           edge_to   = new const edge<W>*[g->V()];
           distances = new Weight[g->V()];
           for (std::size_t i = 0; i < g->V(); ++i) {
               edge_to[i]   = dijkstra.edge_to[i];
               distances[i] = dijkstra.distances[i];
           }
       }

       ~dijkstra_algorithm() {       
           delete[] edge_to;
           delete[] distances;
       }

       Weight distance_to(std::size_t v) const {
           check(v);
           return distances[v];
       }

       bool has_path(std::size_t v) const {
          return distance_to(v) < std::numeric_limits<Weight>::max();
       }
       
       std::deque<const edge<W>*> path(std::size_t v) const {
           check(v);
           
           std::deque<const edge<W>*> path;
           auto e = edge_to[v];
           while (e != nullptr) {
               path.push_front(e);
               e = edge_to[e->from()];
           }
           return path;
       }

       friend void swap(dijkstra_algorithm<W, F, IPQ>& first, 
                        dijkstra_algorithm<W, F, IPQ>& second) noexcept {
           using std::swap;
           swap(first.map_w, second.map_w);
           swap(first.pq, second.pq);
           swap(first.g, second.g);
           swap(first.edge_to, second.edge_to);
           swap(first.distances, second.distances);
       
       }

       dijkstra_algorithm<W, F, IPQ>& operator=(dijkstra_algorithm<W, F, IPQ> other)  noexcept {
           swap(*this, other);
           return *this;
       }

       class dijkstra_iterator
           : public std::iterator<std::forward_iterator_tag,
           std::remove_cv_t<std::size_t>,
           std::ptrdiff_t,
           std::size_t*,
           std::size_t& >
       {

          dijkstra_algorithm<W, F, IPQ>*  dijkstra;
          std::pair<std::size_t, Weight> current;
         
       public:

           explicit dijkstra_iterator(dijkstra_algorithm<W, F, IPQ>*  pdijkstra)
               : dijkstra(pdijkstra) { 
               current = pdijkstra->next();
          
           }
           
           explicit dijkstra_iterator(std::size_t pcurrent)
               : current({ pcurrent, { } }) { }

           dijkstra_iterator& operator++ () {
               current = dijkstra->next();
               return *this;
           }
           dijkstra_iterator operator++ (int) {
               dijkstra_iterator tmp(*this);
               current = dijkstra->next();
               return tmp;
           }

           bool operator == (const dijkstra_iterator& rhs) const {
               return current.first == rhs.current.first;
           }

           bool operator != (const dijkstra_iterator& rhs) const {
               return current.first != rhs.current.first;
           }

           const std::pair<std::size_t, Weight>& operator* () const {
               return current;
           }

           const std::pair<std::size_t, Weight>&  operator-> () const {
               return current;
           }
       };

       using iterator = dijkstra_iterator;
       iterator begin() {
           return iterator(this);
       }

       iterator end()  {
           return iterator(SIZE_MAX);
       }
    };
}

#endif //DIJKSTRA_ALGORITHM_DIJKSTRA_ALGORITHM_H

Algorithms solutions

Monday, September 12, 2016

Dijkstra's algorithm for the SSP problem with a Fibonacci Heap and a D-ary Heap.

Fibonacci Heap:

D-ary Heap

Graph:

Dijkstra's algorithm:

Contributors

Blog Archive