sparse_vector.h

Go to the documentation of this file.

#ifndef LARDATAOBJ_UTILITIES_SPARSE_VECTOR_H
#define LARDATAOBJ_UTILITIES_SPARSE_VECTOR_H


// C/C++ standard library
#include <cstddef> // std::ptrdiff_t
#include <stdexcept> // std::out_of_range()
#include <vector>
#include <ostream>
#include <iterator> // std::distance()
#include <algorithm> // std::lower_bound(), std::max()
#include <numeric> // std::accumulate

#       include <type_traits> // std::is_integral


namespace lar {

// -----------------------------------------------------------------------------
// ---  utility classes for sparse_vector
// ---

template <typename T>
class const_value_box {
        typedef const_value_box<T> this_t;
                public:
        typedef T value_type; 
        
        const_value_box() {}
        
        explicit const_value_box(value_type new_value): value(new_value) {}
        
        this_t& operator= (value_type) { return *this; }
        
        operator value_type() const { return value; }
        operator const value_type& () const { return value; }
        
                protected:
        value_type value; 
}; // class const_value_box


template <typename T>
class value_const_iterator:
        public std::iterator<std::random_access_iterator_tag, T>
{
        typedef std::iterator<std::random_access_iterator_tag, T> base_t; 
        typedef value_const_iterator<T> this_t; 
        
                public:
        
        using typename base_t::value_type;
        using typename base_t::difference_type;
        
        value_const_iterator(): value() {}
        
        value_const_iterator(value_type new_value): value(new_value) {}
        
        value_const_iterator(value_type new_value, difference_type offset):
                index(offset), value(new_value) {}
        
        value_type operator* () const { return value; }
        
        value_type operator[] (difference_type) const { return value; }
        
        this_t& operator++() { ++index; return *this; }
        
        this_t operator++(int) { this_t copy(*this); ++index; return copy; }
        
        this_t& operator--() { --index; return *this; }
        
        this_t operator--(int) { this_t copy(*this); --index; return copy; }
        
        this_t& operator+= (difference_type ofs) { index += ofs; return *this; }
        
        this_t& operator-= (difference_type ofs) { index -= ofs; return *this; }
        
        this_t operator+ (difference_type ofs) const
                { return this_t(value, index + ofs); }
        
        this_t operator- (difference_type ofs) const
                { return this_t(value, index - ofs); }
        
        difference_type operator- (const this_t& iter) const
                { return index - iter.index; }
        
        bool operator== (const this_t& as) const { return index == as.index; }
        bool operator!= (const this_t& as) const { return index != as.index; }
        
        bool operator< (const this_t& as) const { return index < as.index; }
        bool operator<= (const this_t& as) const { return index <= as.index; }
        bool operator> (const this_t& as) const { return index > as.index; }
        bool operator>= (const this_t& as) const { return index >= as.index; }
        
                protected:
        difference_type index{0}; 
        value_type value; 
}; // class value_const_iterator<>


template <typename T>
value_const_iterator<T> operator+ (
        typename value_const_iterator<T>::difference_type ofs,
        value_const_iterator<T>& iter
)
        { return { iter.value, iter.index + ofs }; }

/* // not ready yet...
template <typename T>
class value_iterator: public value_const_iterator<T> {
        using base_t = value_const_iterator<T>;
        using this_t = value_iterator<T>;
        using const_value_box<value_type> = value_box;
                public:
        
        using typename base_t::value_type;
        using typename base_t::reference;
        
        value_box operator* () const { return value_box(value); }
        value_box operator[] (difference_type) const { return value_box(value); }
        
}; // class value_iterator<>
*/



//------------------------------------------------------------------------------
//---  lar::range_t<SIZE>
//---
template <typename SIZE>
class range_t {
                public:
        typedef SIZE size_type; 
        typedef std::ptrdiff_t difference_type; 
        
        static_assert
                (std::is_integral<size_type>::value, "range_t needs an integral type");
        
        size_type offset; 
        size_type last; 
        
        range_t(): offset(0), last(0) {}
        
        range_t(size_type from, size_type to):
                offset(from), last(std::max(from, to)) {}
        
        
        void set(size_type from, size_type to)
                { offset = from; last = std::max(from, to); }
        
        size_type begin_index() const { return offset; }
        
        size_type end_index() const { return last; }
        
        size_type relative_index(size_type index) const { return index - offset; }
        
        size_type size() const { return last - offset; }
        
        void resize(size_type new_size) { last = offset + new_size; }
        
        void move_head(difference_type shift) { offset += shift; }
        
        void move_tail(difference_type shift) { last += shift; }
        
        bool empty() const { return last <= offset; }
        
        bool includes(size_type index) const
                { return (index >= offset) && (index < last); }
        
        bool includes(const range_t& r) const
                { return includes(r.begin_index()) && includes(r.end_index()); }
        
        bool overlap(const range_t& r) const
                { return (begin_index() < r.end_index()) && (end_index() > r.begin_index()); }
        
        bool separate(const range_t& r) const
                { return (begin_index() > r.end_index()) || (end_index() < r.begin_index()); }
        
        bool borders(size_type index) const
                { return (index >= offset) && (index <= last); }
        
        bool operator< (const range_t& than) const { return less(*this, than); }
        
        bool operator== (const range_t& as) const
                { return (offset == as.offset) && (last == as.last); }
        
        bool is_valid() const { return last >= offset; }
        
        static bool less(const range_t& a, const range_t& b)
                { return a.offset < b.offset; }
        static bool less(const range_t& a, size_type b)
                { return a.offset < b; }
        static bool less(size_type a, const range_t& b)
                { return a < b.offset; }
        
        typedef bool (*less_int_range)(size_type, const range_t& b);
}; // class range_t



// -----------------------------------------------------------------------------
// ---  lar::sparse_vector<T>
// ---

template <typename T>
class sparse_vector {
        typedef sparse_vector<T> this_t;
        // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        //  - - - public interface
                public:
        //  - - - types
        typedef T value_type; 
        typedef std::vector<value_type> vector_t;
        typedef typename vector_t::size_type size_type;                
        typedef typename vector_t::difference_type difference_type;
        typedef typename vector_t::pointer pointer;
        typedef typename vector_t::const_pointer const_pointer;
//      typedef typename vector_t::reference reference;
//      typedef typename vector_t::const_reference const_reference;
        
        // --- declarations only ---
        class reference;
        class const_reference;
        class iterator;
        class const_iterator;
        
        class datarange_t;
        // --- ----------------- ---
        
        typedef std::vector<datarange_t> range_list_t;
        typedef typename range_list_t::iterator range_iterator;
        typedef typename range_list_t::const_iterator range_const_iterator;
        
        
        // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
        //  - - - public methods
        sparse_vector(): nominal_size(0), ranges() {}
        
        
        sparse_vector(size_type new_size): nominal_size(0), ranges()
                { resize(new_size); }
        
        sparse_vector(const vector_t& from, size_type offset = 0):
                nominal_size(0), ranges()
                { add_range(offset, from.begin(), from.end()); }
        
        sparse_vector(sparse_vector const&) = default;
        
        sparse_vector(sparse_vector&& from)
                : nominal_size(from.nominal_size)
                , ranges(std::move(from.ranges))
                { from.nominal_size = 0; }
        
        sparse_vector& operator=(sparse_vector const&) = default;
        
        sparse_vector& operator=(sparse_vector&& from)
                {
                        ranges = std::move(from.ranges);
                        nominal_size = from.nominal_size;
                        from.nominal_size = 0;
                        return *this;
                } // operator= (sparse_vector&&)
        
        sparse_vector(vector_t&& from, size_type offset = 0):
                nominal_size(0), ranges()
                { add_range(offset, std::move(from)); }
        
        
        ~sparse_vector() = default;
        
        
        //  - - - STL-like interface
        void clear() { ranges.clear(); nominal_size = 0; }
        
        
        size_type size() const { return nominal_size; }
        
        bool empty() const { return size() == 0; }
        
        size_type capacity() const { return nominal_size; }
        
        void resize(size_type new_size);
        void resize(size_type new_size, value_type def_value);
        
        iterator begin();
        iterator end();
        const_iterator begin() const;
        const_iterator end() const;
        const_iterator cbegin() const { return begin(); }
        const_iterator cend() const { return end(); }
        
        value_type operator[] (size_type index) const;
        
        reference operator[] (size_type index);
        
        
        //  - - - special interface
        
        
        // --- element test
        bool is_void(size_type index) const;
        
        
        bool back_is_void() const
                { return ranges.empty() || (ranges.back().end_index() < size()); }
        
        
        size_type count() const;
        
        

        value_type& set_at(size_type index, value_type value);
        
        void unset_at(size_type index);
        
        void push_back(value_type value) { resize(size() + 1, value); }
        
        void push_back(value_type value, value_type thr)
                {
                        if (is_zero(value, thr)) resize(size() + 1);
                        else                     push_back(value);
                } // push_back()
        
        

        template <typename ITER>
        void assign(ITER first, ITER last) { clear(); append(first, last); }
        
        template <typename CONT>
        void assign(const CONT& new_data) { clear(); append(new_data); }
        
        void assign(vector_t&& new_data) { clear(); append(std::move(new_data)); }
        
        
        
        
        // --- range location
        
        const range_list_t& get_ranges() const { return ranges; }
        
        size_type n_ranges() const { return ranges.size(); }
        
        const datarange_t& range(size_t i) const { return ranges[i]; }
        
        range_const_iterator begin_range() const { return ranges.begin(); }
        
        range_const_iterator end_range() const { return ranges.end(); }
        
        

        range_const_iterator find_range_iterator(size_type index) const;
        range_iterator find_range_iterator(size_type index);
        

        const datarange_t& find_range(size_type index) const;
        datarange_t& find_range(size_type index);
        
        void make_void_around(size_type index);
        

        template <typename ITER>
        const datarange_t& add_range(size_type offset, ITER first, ITER last);
        
        template <typename CONT>
        const datarange_t& add_range(size_type offset, const CONT& new_data)
                { return add_range(offset, new_data.begin(), new_data.end()); }
        
        const datarange_t& add_range(size_type offset, vector_t&& new_data);
        

        template <typename ITER>
        const datarange_t& append(ITER first, ITER last)
                { return add_range(size(), first, last); }
        
        template <typename CONT>
        const datarange_t& append(const CONT& range_data)
                { return add_range(size(), range_data); }
        
        const datarange_t& append(vector_t&& range_data)
                { return add_range(size(), std::move(range_data)); }
        
        

        void void_range(range_iterator iRange) { ranges.erase(iRange); }
        void void_range(size_t iRange) { ranges.erase(ranges.begin() + iRange); }
        
        
        bool is_valid() const;
        
        
        void make_void(iterator first, iterator last);
        
        
        bool optimize() { return optimize(min_gap()); }
        bool optimize(size_t) { return false; /* no optimization implemented yet */ }
        
        
        
        static constexpr value_type value_zero{0}; 
        
        static value_type abs(value_type v) { return (v < value_zero)? -v: v; }
        
        static value_type is_zero(value_type v) { return v == value_zero; }
        
        static value_type is_zero(value_type v, value_type thr)
                { return abs(v - value_zero) <= thr; }
        
        static value_type is_equal(value_type a, value_type b)
          { return is_zero(abs(a - b)); }
        
        static value_type is_equal(value_type a, value_type b, value_type thr)
          { return is_zero(abs(a - b), thr); }
        
        
        static size_t expected_vector_size(size_t size);
        
        static size_t min_gap();
        
        static bool should_merge(
                const typename datarange_t::base_t& a,
                const typename datarange_t::base_t& b
                );
        
        // --------------------------------------------------------------------------
                protected:
        
        size_type nominal_size; 
        range_list_t ranges; 
        

        range_iterator find_next_range_iter(size_type index)
                { return find_next_range_iter(index, ranges.begin()); }
        range_const_iterator find_next_range_iter(size_type index) const
                { return find_next_range_iter(index, ranges.begin()); }


        range_iterator find_next_range_iter(size_type index, range_iterator rbegin);
        range_const_iterator find_next_range_iter
                (size_type index, range_const_iterator rbegin) const;
        
        size_type minimum_size() const
                { return ranges.empty()? 0: ranges.back().end_index(); }
        
        range_iterator insert_range(range_iterator iInsert, const datarange_t& data)
                { return data.empty()? iInsert: ranges.insert(iInsert, data); }
        range_iterator insert_range(range_iterator iInsert, datarange_t&& data)
                { return data.empty()? iInsert: ranges.insert(iInsert, std::move(data)); }
        
        range_iterator eat_range_head(range_iterator iRange, size_t index);
        
        datarange_t& merge_ranges(range_iterator iRange);
        
        size_type fix_size();
        
}; // class sparse_vector<>


} // namespace lar

template <typename T>
std::ostream& operator<< (std::ostream& out, const lar::sparse_vector<T>& v);



// -----------------------------------------------------------------------------
// --- sparse_vector::datarange_t definition
// ---

template <typename T>
class lar::sparse_vector<T>::datarange_t: public range_t<size_type> {
                public:
        typedef range_t<size_type> base_t; 
        
        typedef typename vector_t::iterator iterator;
        typedef typename vector_t::const_iterator const_iterator;
        
        datarange_t(): base_t(), values() {}
        
        datarange_t(const base_t& range): base_t(range), values(range.size()) {}
        
        template <typename ITER>
        datarange_t(size_type offset, ITER first, ITER last):
                base_t(offset, offset + std::distance(first, last)),
                values(first, last)
                {}

        datarange_t(size_type offset, vector_t&& data):
                base_t(offset, offset + data.size()), values(data)
                {}
        
        
        iterator get_iterator(size_type index)
                { return values.begin() + index - base_t::begin_index(); }
        const_iterator get_iterator(size_type index) const
                { return values.begin() + index - base_t::begin_index(); }
        
        iterator begin() { return values.begin(); }
        iterator end() { return values.end(); }
        const_iterator begin() const { return values.begin(); }
        const_iterator end() const { return values.end(); }
        const_iterator cbegin() const { return values.cbegin(); }
        const_iterator cend() const { return values.cend(); }
        
        void resize(size_t new_size)
                { values.resize(new_size); fit_size_from_data(); }
        void resize(size_t new_size, value_type def_value)
                { values.resize(new_size, def_value); fit_size_from_data(); }
        
        value_type& operator[] (size_type index)
                { return values[base_t::relative_index(index)]; }
        const value_type& operator[] (size_type index) const
                { return values[base_t::relative_index(index)]; }
        
        const vector_t& data() const { return values; }
//      vector_t& data() { return values; }
        
        template <typename ITER>
        datarange_t& extend(size_type index, ITER first, ITER last);
        
        
        void move_head(size_type to_index, value_type def_value = value_zero);
        
        void move_tail(size_type to_index, value_type def_value = value_zero)
                { resize(base_t::relative_index(to_index), def_value); }
        
        
                protected:
        vector_t values; 
        
        void fit_size_from_data() { base_t::resize(values.size()); }
        
}; // class datarange_t



// -----------------------------------------------------------------------------
// --- sparse_vector iterators definition
// ---

template <typename T>
class lar::sparse_vector<T>::const_reference {
                protected:
        const value_type* ptr;
                public:
        const_reference(const value_type* pValue = 0): ptr(pValue) {}
        const_reference(const value_type& value): const_reference(&value) {}
        
        explicit operator value_type() const { return ptr? *ptr: value_zero; }
        operator const value_type&() const
                { return ptr? *ptr: value_zero; }
}; // lar::sparse_vector<T>::const_reference


template <typename T>
class lar::sparse_vector<T>::reference: public const_reference {
        friend class iterator;
        
        // This object "disappears" when assigned to: either the assignment is not
        // possible, and then a segmentation fault will occur, or the return value
        // is an actual C++ reference to the assigned value.
        // The same is true when explicitly converting it to a reference.
                public:
        reference(value_type* pValue = 0): const_reference(pValue) {}
        reference(value_type& value): reference(&value) {}
        
        reference& operator=(const reference&) = default;
        value_type& operator=(value_type v)
                { return const_cast<value_type&>(*const_reference::ptr) = v; }
        
        operator const_reference() const
                { return const_reference(const_reference::ptr); }
        explicit operator value_type&()
                { return const_cast<value_type&>(*const_reference::ptr); }
        
                protected:
        explicit reference(const const_reference& from): const_reference(from) {}
        
}; // lar::sparse_vector<T>::reference


// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
template <typename T>
class lar::sparse_vector<T>::const_iterator {
        //
        // This iterator fulfils the traits of an immutable forward iterator.
        //
        
                protected:
        friend class container_t;
        
        typedef sparse_vector<T> container_t;
        typedef typename container_t::size_type size_type;
        typedef typename container_t::range_list_t::const_iterator ranges_const_iterator;
        
        
                public:
        // so far, only forward interface is implemented;
        // backward is tricky...
        typedef std::forward_iterator_tag iterator_category;
        
        struct special {
                class begin {};
                class end {};
        };
        
        // types
        typedef typename container_t::value_type value_type;
        typedef typename container_t::difference_type difference_type;
        typedef typename container_t::pointer pointer;
        typedef typename container_t::reference reference;
        // note that this is not out little box, it's the normal C++ constant
        // reference from the underlying vector
        typedef typename vector_t::const_reference const_reference;
        
        
        const_iterator(): cont(nullptr), index(0), currentRange() {}
        
        const_iterator(const container_t& c, size_type offset):
                cont(&c), index(std::min(offset, c.size())), currentRange()
                { refresh_state(); }
        
        const_iterator(const container_t& c, const typename special::begin):
                cont(&c), index(0), currentRange(c.get_ranges().begin())
                {}
        
        const_iterator(const container_t& c, const typename special::end):
                cont(&c), index(c.size()), currentRange(c.get_ranges().end())
                {}
        
        const_reference operator[] (size_type offset) const
                { return (*cont)[index + offset]; }
        
        const_reference operator*() const;
        
        const_iterator& operator++();
        const_iterator operator++(int)
                { const_iterator copy(*this); const_iterator::operator++(); return copy; }
        
        
        const_iterator& operator+= (difference_type delta);
        const_iterator& operator-= (difference_type delta);
        const_iterator operator+ (difference_type delta) const;
        const_iterator operator- (difference_type delta) const;
        
        difference_type operator- (const const_iterator& iter) const;
        
        bool operator== (const const_iterator& as) const
                { return (cont == as.cont) && (index == as.index); }
        bool operator!= (const const_iterator& as) const
                { return (cont != as.cont) || (index != as.index); }
        bool operator< (const const_iterator& than) const
                { return (cont == than.cont) && (index < than.index); }
        bool operator> (const const_iterator& than) const
                { return (cont == than.cont) && (index > than.index); }
        bool operator<= (const const_iterator& than) const
                { return (cont == than.cont) && (index <= than.index); }
        bool operator>= (const const_iterator& than) const
                { return (cont == than.cont) && (index >= than.index); }
        
        range_const_iterator get_current_range() const { return currentRange; }
        
                protected:
        //
        // Logic of the internals:
        // - cont provides the list of ranges
        // - index is the absolute index in the container
        // - currentRange is a pointer to the "current" range, cached for speed
        // 
        // An iterator past-the-end has the index equal to the size of the
        // container it points to. Operations on it are undefined, except that
        // going backward by one decrement goes back to the container
        // (if not empty) TODO currently backward iteration is not supported
        // 
        // When the iterator is pointing to an element within a range,
        // currentRange points to that range.
        // When the iterator is pointing into the void of the vector, currentRange
        // points to the range next to that void area, or end() if there is none.
        // When the iterator is past-the-end, currentRange is not defined.
        // 
        // Conversely, when the index is equal (or greater) to the size of the
        // container, the iterator is not dereferenciable, the iterator points
        // past-the-end of the container.
        // When currentRange points to a valid range and the index
        // is before that range, the iterator is pointing to the void.
        // When currentRange points to a valid range and the index
        // is inside that range, the iterator is pointing to an item in that
        // range.
        // When currentRange points to a valid range, the index can't point beyond
        // that range.
        // When currentRange points to the past-to-last range, the iterator is
        // pointing to the void (past the last range).
        //
        const container_t* cont; 
        size_type index; 
        ranges_const_iterator currentRange; 
        
        void refresh_state();
        
}; // class lar::sparse_vector<T>::const_iterator


template <typename T>
class lar::sparse_vector<T>::iterator: public const_iterator {
        typedef typename const_iterator::container_t container_t;
        friend typename const_iterator::container_t;
        
                public:
        typedef typename const_iterator::reference reference;
        typedef typename const_iterator::const_reference const_reference;
        typedef typename const_iterator::special special;
        
        iterator(): const_iterator() {}
        
        iterator(container_t& c, size_type offset = 0):
                const_iterator(c, offset) {}
        
        iterator(const container_t& c, const typename special::begin _):
                const_iterator(c, _) {}
        
        iterator(const container_t& c, const typename special::end _):
                const_iterator(c, _) {}
        
        reference operator[] (size_type offset) const
                { return (*const_iterator::cont)[const_iterator::index + offset]; }
        
        iterator& operator++() { const_iterator::operator++(); return *this; }
        iterator operator++(int _) { return const_iterator::operator++(_); }
        
        iterator& operator+= (difference_type delta)
                { const_iterator::operator+=(delta); return *this; }
        iterator& operator-= (difference_type delta)
                { const_iterator::operator+=(delta); return *this; }
        iterator operator+ (difference_type delta) const
                { return const_iterator::operator+(delta); }
        iterator operator- (difference_type delta) const
                { return const_iterator::operator-(delta); }
        
        
        reference operator*() const
                { return reference(const_iterator::operator*()); }
        
//      /// Returns the current range internal value; use it at your own risk!!
//      range_iterator get_current_range() const
//              { return const_iterator::currentRange; }
        
                protected:
        
        // also worth conversion
        iterator(const_iterator from): const_iterator(from) {}
        
}; // class lar::sparse_vector<T>::iterator



// -----------------------------------------------------------------------------
// ---  implemetation  ---------------------------------------------------------
// -----------------------------------------------------------------------------

//------------------------------------------------------------------------------
//--- sparse_vector implementation

template <typename T>
constexpr typename lar::sparse_vector<T>::value_type
        lar::sparse_vector<T>::value_zero;


template <typename T> 
void lar::sparse_vector<T>::resize(size_type new_size) {
        if (new_size >= size()) {
                nominal_size = new_size;
                return;
        }
        
        // truncating...
        range_iterator iLastRange = find_next_range_iter(new_size);
        ranges.erase(iLastRange, ranges.end());
        if (!ranges.empty()) {
                range_iterator iLastRange = ranges.end() - 1;
                if (new_size == iLastRange->begin_index())
                        ranges.erase(iLastRange);
                else if (new_size < iLastRange->end_index())
                        iLastRange->resize(new_size - iLastRange->begin_index());
        } // if we have ranges
        
        // formally resize
        nominal_size = new_size;
} // lar::sparse_vector<T>::resize()


template <typename T> 
void lar::sparse_vector<T>::resize(size_type new_size, value_type def_value) {
        if (new_size == size()) return;
        if (new_size > size()) {
                
                if (back_is_void()) // add a new range
                        append(vector_t(new_size - size(), def_value));
                else // extend the last range
                        ranges.back().resize(new_size - ranges.back().begin_index(), def_value);
                
                nominal_size = new_size;
                return;
        }
        // truncating is the same whether there is a default value or not
        resize(new_size);
} // lar::sparse_vector<T>::resize()


template <typename T>
inline typename lar::sparse_vector<T>::iterator lar::sparse_vector<T>::begin()
        { return iterator(*this, typename iterator::special::begin()); }

template <typename T>
inline typename lar::sparse_vector<T>::iterator lar::sparse_vector<T>::end()
        { return iterator(*this, typename iterator::special::end()); }

template <typename T>
inline typename lar::sparse_vector<T>::const_iterator
        lar::sparse_vector<T>::begin() const
        { return const_iterator(*this, typename const_iterator::special::begin()); }

template <typename T>
inline typename lar::sparse_vector<T>::const_iterator
        lar::sparse_vector<T>::end() const
        { return const_iterator(*this, typename const_iterator::special::end()); }

template <typename T>
typename lar::sparse_vector<T>::value_type lar::sparse_vector<T>::operator[]
        (size_type index) const
{
        // first range not including the index
        range_const_iterator iNextRange = find_next_range_iter(index);
        
        // if not even the first range includes the index, we are in the void
        // (or in some negative index space, if index signedness allowed);
        if (iNextRange == ranges.begin()) return value_zero;
        
        // otherwise, let's take the previous range;
        // if it includes the index, we return its value;
        // or it precedes it, and our index is in the void: return zero
        const datarange_t& range(*--iNextRange);
        return (index < range.end_index())? range[index]: value_zero;
} // lar::sparse_vector<T>::operator[]


template <typename T>
typename lar::sparse_vector<T>::reference lar::sparse_vector<T>::operator[]
        (size_type index)
{
        // first range not including the index
        range_iterator iNextRange = find_next_range_iter(index);
        
        // if not even the first range includes the index, we are in the void
        // (or in some negative index space, if index signedness allowed);
        if (iNextRange == ranges.begin()) return reference();
        
        // otherwise, let's take the previous range;
        // if it includes the index, we return its value;
        // or it precedes it, and our index is in the void: return zero
        datarange_t& range(*--iNextRange);
        return (index < range.end_index())? reference(range[index]): reference();
} // lar::sparse_vector<T>::operator[]


template <typename T>
bool lar::sparse_vector<T>::is_void(size_type index) const {
        if (ranges.empty() || (index >= size()))
                throw std::out_of_range("empty sparse vector");
        // range after the index:
        range_const_iterator iNextRange = find_next_range_iter(index);
        return ((iNextRange == ranges.begin()) 
                || ((--iNextRange)->end_index() <= index));
} // lar::sparse_vector<T>::is_void()


template <typename T>
inline typename lar::sparse_vector<T>::size_type lar::sparse_vector<T>::count()
        const
{
        return std::accumulate(begin_range(), end_range(), size_type(0),
                [](size_type s, const datarange_t& rng) { return s + rng.size(); }
                );
} // count()


template <typename T>
typename lar::sparse_vector<T>::value_type& lar::sparse_vector<T>::set_at
        (size_type index, value_type value)
{
        // first range not including the index
        range_iterator iNextRange = find_next_range_iter(index);
        
        // if not even the first range includes the index, we are in the void
        // (or in some negative index space, if index signedness allowed);
        if (iNextRange != ranges.begin()) {
                // otherwise, let's take the previous range;
                // if it includes the index, we set the existing value;
                // or it precedes it, and our index is in the void, add the value as a
                // range
                datarange_t& range(*--iNextRange);
                if (index < range.end_index()) return range[index] = value;
        }
        // so we are in the void; add the value as a new range
        return const_cast<datarange_t&>(add_range(index, { value }))[index];
} // lar::sparse_vector<T>::set_at()


template <typename T>
void lar::sparse_vector<T>::unset_at(size_type index) {
        // first range not including the index
        range_iterator iNextRange = find_next_range_iter(index);
        
        // if not even the first range includes the index, we are in the void
        // (or in some negative index space, if index signedness allowed);
        if (iNextRange == ranges.begin()) return;
        
        // otherwise, let's take the previous range;
        // or it precedes the index, and our index is in the void
        datarange_t& range(*--iNextRange);
        if (index >= range.end_index()) return; // void already
        
        // it includes the index:
        // - if it's a one-element range, remove the range
        if (range.size() == 1) ranges.erase(iNextRange);
        // - if it's on a border, resize the range
        else if (index == range.begin_index())
                range.move_head(index + 1);
        else if (index == range.end_index() - 1)
                range.move_tail(index);
        // - if it's inside, break the range in two
        else if (index < range.end_index()) {
                // we are going to put a hole in a range;
                // this effectively creates two ranges;
                // first create the rightmost range, and add it to the list
                ranges.emplace(++iNextRange, index + 1,
                        range.begin() + range.relative_index(index + 1),
                        range.end()
                        );
                // then cut the existing one
                range.move_tail(index);
        }
} // lar::sparse_vector<T>::unset_at()


template <typename T>
typename lar::sparse_vector<T>::range_const_iterator
        lar::sparse_vector<T>::find_range_iterator(size_type index) const
{
        if (ranges.empty()) throw std::out_of_range("empty sparse vector");
        // range after the index:
        range_const_iterator iNextRange = find_next_range_iter(index);
        return ((iNextRange == ranges.begin()) 
                || (index >= (--iNextRange)->end_index()))?
                ranges.end(): iNextRange;
} // lar::sparse_vector<T>::find_range_iterator() const


template <typename T>
inline typename lar::sparse_vector<T>::range_iterator 
        lar::sparse_vector<T>::find_range_iterator(size_type index)
{
        return ranges.begin() + (
                (const_cast<const this_t*>(this)->find_range_iterator(index))
                - ranges.begin());
} // lar::sparse_vector<T>::find_range_iterator()


template <typename T>
const typename lar::sparse_vector<T>::datarange_t&
        lar::sparse_vector<T>::find_range(size_type index) const
{
        if (ranges.empty()) throw std::out_of_range("empty sparse vector");
        // range on the index:
        range_const_iterator iNextRange = find_range_iterator(index);
        if (iNextRange == ranges.end())
                throw std::out_of_range("index in no range of the sparse vector");
        return *iNextRange;
} // lar::sparse_vector<T>::find_range() const

template <typename T>
inline typename lar::sparse_vector<T>::datarange_t&
        lar::sparse_vector<T>::find_range(size_type index)
{
        return const_cast<datarange_t&>
                (const_cast<const this_t*>(this)->find_range(index));
} // lar::sparse_vector<T>::find_range()


template <typename T>
void lar::sparse_vector<T>::make_void_around(size_type index) {
        if (ranges.empty() || (index >= size()))
                throw std::out_of_range("empty sparse vector");
        // range after the index:
        range_iterator iNextRange = find_next_range_iter(index);
        if ((iNextRange == ranges.begin()) 
                || ((--iNextRange)->end_index() <= index))
        {
                return;
        }
        ranges.erase(iNextRange);
} // lar::sparse_vector<T>::make_void_around()


template <typename T> template <typename ITER>
const typename lar::sparse_vector<T>::datarange_t& lar::sparse_vector<T>::add_range
  (size_type offset, ITER first, ITER last)
{
        // insert the new range before the existing range which starts after offset
        range_iterator iInsert = std::upper_bound(
                ranges.begin(), ranges.end(), offset,
                typename datarange_t::less_int_range(datarange_t::less)
                );
        
        // is there a range before this, which includes the offset?
        if ((iInsert != ranges.begin()) && (iInsert-1)->borders(offset)) {
                // then we should extend it
                (--iInsert)->extend(offset, first, last);
        }
        else {
                // no range before the insertion one includes the offset of the new range;
                // ... we need to add it as a new range
                iInsert = insert_range(iInsert, { offset, first, last });
        }
        return merge_ranges(iInsert);
} // lar::sparse_vector<T>::add_range<ITER>()


template <typename T>
const typename lar::sparse_vector<T>::datarange_t& lar::sparse_vector<T>::add_range
  (size_type offset, vector_t&& new_data)
{
        // insert the new range before the existing range which starts after offset
        range_iterator iInsert = std::upper_bound(
                ranges.begin(), ranges.end(), offset,
                typename datarange_t::less_int_range(datarange_t::less)
                );
        
        // is there a range before this, which includes the offset?
        if ((iInsert != ranges.begin()) && (iInsert-1)->borders(offset)) {
                // then we should extend it
                (--iInsert)->extend(offset, new_data.begin(), new_data.end());
        }
        else {
                // no range before the insertion one includes the offset of the new range;
                // ... we need to add it as a new range;
                // it is not really clear to me [GP] why I need a std::move here, since
                // new_data is a rvalue already; in doubt, I have painted all this kind
                // of constructs with move()s, just in case
                iInsert = insert_range(iInsert, { offset, std::move(new_data) });
        }
        return merge_ranges(iInsert);
} // lar::sparse_vector<T>::add_range(vector)


template <typename T> 
void lar::sparse_vector<T>::make_void(iterator first, iterator last) {
        // iterators have in "currentRange" either the range they point into, 
        // or the range next to the void they point to
        
        if ((first.cont != this) || (last.cont != this)) {
                throw std::runtime_error
                        ("lar::sparse_vector::make_void(): iterators from alien container");
        }
        // if first is after next, no range is identified
        if (first >= last) return;
        
        range_iterator
                first_range
                        = ranges.begin() + (first.get_current_range() - ranges.begin()),
                last_range = ranges.begin() + (last.get_current_range() - ranges.begin());
        
        //--- "first": void up to this iterator ---
        // if first in the last void region, there is nothing to erase
        if (first_range == ranges.end()) return;
        
        // if first is in the middle of a valid range instead, we have to resize it
        if (first.index > first_range->begin_index()) {
                if (first_range == last_range) {
                        // we are going to erase a subset of a range;
                        // this effectively creates two ranges;
                        // first create the rightmost (last) range, and add it to the list
                        last_range = ranges.emplace(++last_range, last.index,
                                first_range->begin() + first_range->relative_index(last.index),
                                first_range->end()
                                );
                        // then cut the existing one
                        first_range->move_tail(first.index);
                        return;
                }
                first_range->move_tail(first.index);
                ++first_range; // from next range on, start voiding
        }
        
        //--- "last"
        // if the index is inside a range, remove up to it
        if ((last_range != ranges.end()) && (last.index > last_range->begin_index()))
                eat_range_head(last_range, last.index);
        
        // finally, remove entirely the ranges in between
        ranges.erase(first_range, last_range);
} // lar::sparse_vector<T>::make_void()


template <typename T>
bool lar::sparse_vector<T>::is_valid() const {
        // a sparse vector with no non-null elements can't be detected invalid
        if (ranges.empty()) return true;
        
        range_const_iterator iNext = ranges.begin(), rend = ranges.end();
        while (iNext != rend) {
                range_const_iterator iRange = iNext++;
                if (iRange->empty()) return false;
                if (iNext != rend) {
                        if (!(*iRange < *iNext)) return false;
                        if (!iRange->separate(*iNext)) return false;
                }
        } // while
        if (nominal_size < ranges.back().end_index()) return false;
        return true;
} // lar::sparse_vector<T>::is_valid()



// --- private methods

template <typename T> 
typename lar::sparse_vector<T>::range_iterator
        lar::sparse_vector<T>::find_next_range_iter
        (size_type index, range_iterator rbegin)
{
        // this range has the offset (first index) above the index argument:
        return std::upper_bound(
                rbegin, ranges.end(), index,
                typename datarange_t::less_int_range(datarange_t::less)
                );
} // lar::sparse_vector<T>::find_next_range_iter()

template <typename T> 
typename lar::sparse_vector<T>::range_const_iterator
        lar::sparse_vector<T>::find_next_range_iter
        (size_type index, range_const_iterator rbegin) const
{
        // this range has the offset (first index) above the index argument:
        return std::upper_bound(
                rbegin, ranges.end(), index,
                typename datarange_t::less_int_range(datarange_t::less)
                );
} // lar::sparse_vector<T>::find_next_range_iter() const


template <typename T> 
typename lar::sparse_vector<T>::datarange_t& lar::sparse_vector<T>::merge_ranges
  (range_iterator iRange)
{
        range_iterator iNext = iRange + 1;
        while (iNext != ranges.end()) {
                if (!iRange->borders(iNext->begin_index())) break;
                // iRange content dominates, but if iNext has data beyond it,
                // then we copy it
                if (iNext->end_index() > iRange->end_index()) {
                        iRange->extend
                                (iRange->end_index(), iNext->get_iterator(iRange->end_index()), iNext->end());
                }
                iNext = ranges.erase(iNext);
        } // while
        fix_size();
        return *iRange;
} // lar::sparse_vector<T>::merge_ranges()


template <typename T> 
typename lar::sparse_vector<T>::range_iterator lar::sparse_vector<T>::eat_range_head
        (range_iterator iRange, size_t index)
{
        if (index <= iRange->begin_index()) return iRange;
        if (index >= iRange->end_index()) return ranges.erase(iRange);
        iRange->move_head(index);
        return iRange;
} // lar::sparse_vector<T>::eat_range_head()


template <typename T>
typename lar::sparse_vector<T>::size_type lar::sparse_vector<T>::fix_size() {
        if (!ranges.empty())
                nominal_size = std::max(nominal_size, ranges.back().end_index());
        return nominal_size;
} // lar::sparse_vector<T>::fix_size()


// --- static methods

template <typename T>
inline size_t lar::sparse_vector<T>::expected_vector_size(size_t size) {
        // apparently, a chunk of heap memory takes at least 32 bytes;
        // that means that a vector of 1 or 5 32-bit integers takes the same
        // space; the overhead appears to be 8 bytes, which can be allocated
        return sizeof(vector_t)
                + std::max(size_t(32), (alignof(datarange_t)*size + 8));
} // lar::sparse_vector<T>::expected_vector_size()


template <typename T>
inline size_t lar::sparse_vector<T>::min_gap() {
        // we assume here that there is no additional overhead by alignment;
        // the gap adds the space of another datarange_t, including the vector,
        // its data and overhead from heap (apparently, 8 bytes);
        // 
        return (sizeof(datarange_t) + 8) / sizeof(value_type) + 1; // round up
} // lar::sparse_vector<T>::min_gap()


template <typename T>
inline bool lar::sparse_vector<T>::should_merge(
        const typename datarange_t::base_t& a,
        const typename datarange_t::base_t& b
        )
{
        size_type gap_size = (a < b)?
                b.begin_index() - a.begin_index() - a.size():
                a.begin_index() - b.begin_index() - b.size();
        return expected_vector_size(a.size() + b.size() + gap_size)
                <= expected_vector_size(a.size()) + expected_vector_size(b.size());
} // lar::sparse_vector<T>::should_merge()



// --- non-member functions
template <typename T>
std::ostream& operator<< (std::ostream& out, const lar::sparse_vector<T>& v) {
        
        out << "Sparse vector of size " << v.size() << " with "
                << v.get_ranges().size() << " ranges:";
        typename lar::sparse_vector<T>::range_const_iterator
                iRange = v.begin_range(), rend = v.end_range();
        while (iRange != rend) {
                out << "\n  [" << iRange->begin_index() << " - " << iRange->end_index()
                        << "] (" << iRange->size() << "):";
                typename lar::sparse_vector<T>::datarange_t::const_iterator
                        iValue = iRange->begin(), vend = iRange->end();
                while (iValue != vend) out << " " << (*(iValue++));
                ++iRange;
        } // for
        return out << std::endl;
} // operator<< (ostream, sparse_vector<T>)



// -----------------------------------------------------------------------------
// --- lar::sparse_vector<T>::datarange_t implementation
// ---
template <typename T> template <typename ITER>
typename lar::sparse_vector<T>::datarange_t& lar::sparse_vector<T>::datarange_t::extend
  (size_type index, ITER first, ITER last)
{
        size_type new_size = std::max(
                base_t::relative_index(index) + std::distance(first, last),
                base_t::size());
        base_t::resize(new_size);
        values.resize(new_size);
        std::copy(first, last, get_iterator(index));
        return *this;
} // lar::sparse_vector<T>::datarange_t::extend()


template <typename T>
void lar::sparse_vector<T>::datarange_t::move_head
        (size_type to_index, value_type def_value /* = value_zero */)
{
        difference_type delta = to_index - base_t::begin_index();
        if (delta == 0) return;
        base_t::move_head(delta);
        if (delta > 0) // shrink
                values.erase(values.begin(), values.begin() + delta);
        else { // extend
                values.insert(values.begin(),
                        value_const_iterator<value_type>(def_value),
                        value_const_iterator<value_type>(def_value) - delta
                        );
        }
} // lar::sparse_vector<T>::datarange_t::move_head()


// -----------------------------------------------------------------------------
// --- lar::sparse_vector<T>::const_iterator implementation
// ---
template <typename T>
typename lar::sparse_vector<T>::const_iterator&
lar::sparse_vector<T>::const_iterator::operator++() {
        // no container, not doing anything;
        // index beyond the end: stays there
        if (!cont || (index >= cont->size())) return *this;
        
        // increment the internal index
        ++index;
        
        // were we in a valid range?
        if (currentRange != cont->ranges.end()) {
                // if the new index is beyond the current range, pick the next
                if (currentRange->end_index() <= index) ++currentRange;
        }
        // if we have no valid range, we are forever in the void
        return *this;
} // lar::sparse_vector<T>::iterator::operator++()


template <typename T>
typename lar::sparse_vector<T>::const_iterator::const_reference
lar::sparse_vector<T>::const_iterator::operator*() const {
        // no container, no idea what to do
        if (!cont) throw std::out_of_range("iterator to no sparse vector");
        
        // index beyond the end: can't return any reference
        if (index >= cont->size()) return value_zero;
        
        // are we in a valid range? if not, we are past the last range
        if (currentRange == cont->ranges.end()) return value_zero;
        
        // if the index is beyond the current range, we are in the void
        if (index < currentRange->begin_index()) return value_zero;
        
        return (*currentRange)[index];
} // lar::sparse_vector<T>::const_iterator::operator*()


template <typename T>
typename lar::sparse_vector<T>::const_iterator&
        lar::sparse_vector<T>::const_iterator::operator+= (difference_type delta)
{
        if (delta == 1) return this->operator++();
        index += delta;
        if ((currentRange == cont->ranges.end())
                || !currentRange->includes(index)
                )
                refresh_state();
        return *this;
} // lar::sparse_vector<T>::const_iterator::operator+=()

template <typename T>
inline typename lar::sparse_vector<T>::const_iterator&
        lar::sparse_vector<T>::const_iterator::operator-= (difference_type delta)
        { return this->operator+= (-delta); }


template <typename T>
typename lar::sparse_vector<T>::const_iterator
        lar::sparse_vector<T>::const_iterator::operator+ (difference_type delta) const
{
        if ((currentRange == cont->ranges.end())
                || !currentRange->includes(index + delta)
                )
                return const_iterator(*cont, index + delta);
        const_iterator iter(*this);
        iter.index += delta;
        return iter;
} // lar::sparse_vector<T>::const_iterator::operator+()

template <typename T>
inline typename lar::sparse_vector<T>::const_iterator
        lar::sparse_vector<T>::const_iterator::operator- (difference_type delta) const
        { return this->operator+ (-delta); }


template <typename T>
inline typename lar::sparse_vector<T>::const_iterator::difference_type
        lar::sparse_vector<T>::const_iterator::operator-
        (const const_iterator& iter) const
{
        if (cont != iter.cont) {
                throw std::runtime_error("lar::sparse_vector::const_iterator:"
                        " difference with alien iterator");
        }
        return index -iter.index;
} // lar::sparse_vector<T>::const_iterator::operator-(const_iterator)


template <typename T>
void lar::sparse_vector<T>::const_iterator::refresh_state() {
        // update the currentRange
        // currentRange is the range including the current item, or next to it
        if (cont) {
                // the following is actually the range after the index:
                currentRange = cont->find_next_range_iter(index);
                // if the index is inside the previous index, then we want to move there:
                if (currentRange != cont->ranges.begin()) {
                        if ((currentRange - 1)->end_index() > index) --currentRange;
                }
        }
        else {
                currentRange = {};
        }
} // lar::sparse_vector<T>::const_iterator::refresh_state()


// -----------------------------------------------------------------------------
// --- lar::sparse_vector<T>::iterator implementation
// ---
//
// nothing new so far
//


#endif // LARDATAOBJ_UTILITIES_SPARSE_VECTOR_H