# Local heapsort in C++

Suppose we need to sort a sequence and we know that every sequence component is within $d$ steps from its correct position. In such a case we can use a local heapsort: we take $d$ as an argument and set $w = d + 1$ (window width). We load the heap with $w$ first sequence components after which we pop the heap and add the next component to the heap (a sliding window). Finally, we dump the leftovers to the end of the range to be sorted:

local_heapsort.hpp

// Created Apr 5, 2017 by Rodion "rodde" Efremov

#ifndef CODERODDE_UTIL_LOCAL_HEAPSORT
#define CODERODDE_UTIL_LOCAL_HEAPSORT

#include <algorithm>
#include <iterator>
#include <queue>
#include <vector>

namespace coderodde {
namespace util {

template<typename RandomIt, typename Cmp>
void local_heapsort(RandomIt begin,
RandomIt end,
Cmp cmp,
size_t max_distance)
{
using value_type =
typename std::iterator_traits<RandomIt>::value_type;

size_t range_distance = std::distance(begin, end);
max_distance = std::min(max_distance, range_distance);
size_t window_width = max_distance + 1;

auto window_end = begin;
std::priority_queue<value_type, std::vector<value_type>, Cmp> heap;

auto current_iterator = begin;
auto next_iterator = begin;

for (auto iter = current_iterator; iter != next_iterator; ++iter)
{
heap.push(*iter);
}

while (next_iterator != end)
{
*current_iterator = heap.top();
heap.pop();
heap.push(*next_iterator);
++current_iterator;
++next_iterator;
}

while (!heap.empty())
{
*current_iterator = heap.top();
heap.pop();
++current_iterator;
}
}

} // End of namespace coderodde::util
} // End of namespace coderodde

#endif // CODERODDE_UTIL_LOCAL_HEAPSORT


main.cpp

#include "local_heapsort.hpp"
#include <algorithm>
#include <chrono>
#include <iostream>
#include <iterator>
#include <random>
#include <vector>

using namespace std;
using namespace std::chrono;
using coderodde::util::local_heapsort;

static void small_demo()
{
cout << "Sorting: ";
vector<int> vec = { 3, 2, 1, 5, 4, 8, 7, 6 };
copy(vec.begin(), vec.end(), ostream_iterator<int>(cout, " "));
cout << endl;
local_heapsort(vec.begin(), vec.end(), std::greater<int>(), 2);
cout << "Sorted:  ";
copy(vec.begin(), vec.end(), ostream_iterator<int>(cout, " "));
cout << endl;
}

static const size_t ARRAY_LENGTH = 50 * 1000 * 1000;
static const size_t MAX_DISTANCE = 10;

static vector<int> get_random_int_array_with_max_distance(size_t array_length,
size_t max_distance)
{
vector<int> ret(array_length);

for (int i = 0; i < array_length; ++i)
{
ret[i] = i;
}

// Shuffle small chunks in such that any array component is at most
// 'max_distance' position from its eventual position:
size_t chunk_size = max_distance + 1;
random_device rd;
mt19937 g(rd());

for (size_t i = 0; i < array_length; i += chunk_size)
{
shuffle(&ret[i], &ret[min(i + chunk_size, array_length)], g);
}

return ret;
}

static void benchmark()
{
vector<int> array1 = get_random_int_array_with_max_distance(ARRAY_LENGTH,
MAX_DISTANCE);
vector<int> array2(array1);

auto tp1 = high_resolution_clock::now();
// For fairness, we pass the same comparator to std::sort
//as to local_heapsort.
sort(array1.begin(), array1.end(), std::less<int>());
auto tp2 = high_resolution_clock::now();

cout << "std::sort in "
<< duration_cast<milliseconds>(tp2 - tp1).count()
<< " milliseconds."
<< endl;

tp1 = high_resolution_clock::now();
local_heapsort(array2.begin(),
array2.end(),
std::greater<int>(),
MAX_DISTANCE);
tp2 = high_resolution_clock::now();

cout << "local_heapsort in "
<< duration_cast<milliseconds>(tp2 - tp1).count()
<< " milliseconds."
<< endl;

cout << "Algorithms agree: "
<< boolalpha
<< equal(array1.begin(),
array1.end(),
array2.begin())
<< endl;
}

int main() {
small_demo();
benchmark();
return 0;
}


Performance figures

My demo program gives the following output:

Sorting: 3 2 1 5 4 8 7 6
Sorted:  1 2 3 4 5 6 7 8
std::sort in 1845 milliseconds.
local_heapsort in 1359 milliseconds.
Algorithms agree: true


Critique request

Please tell me anything that comes to mind. Especially, I wish to know how to make my code more idiomatic; say, is there opportunity to replace some explicit loops with facilities from <algorithm>?

• Is the get_random_... sufficient for generating input for the function? So that it would be possible to generate lots of different input for the sorting function. – Incomputable Apr 5 '17 at 20:11
• For once, I can see what you are trying to achieve - pity it doesn't show up in the code. – greybeard Apr 6 '17 at 2:15
• @coderodde, I've got some sweet things: benchmarks! Posted it as different answer, and added my reasoning about results into the post. – Incomputable Apr 7 '17 at 22:13

Bug:

I found this post, which says that if Cmp is a function pointer, priority_queue will default construct it to nullptr, which is straight undefined behavior. So, passing cmp into priority_queue becomes a must.

I believe that idiomatic code is the code that conforms C++ concepts, so I will try to cover those as well.

Formal specification:

I believe the formal specification is the best thing to check conformance.

From the names of template type parameters I inferred following C++ concepts:

• RandomAccessIterator
• Compare = BinaryPredicate + semantics
• CopyConstructible (may be other requirements on element type)

Code has good conformance, except Compare. If you'll follow long chain of links on cppreference, at the end it becomes clear that BinaryPredicate is ought to be only CopyConstructible. The priority_queue will try to default construct it, which is a violation of concept. The fix would be just to pass the cmp into constructor. Fortunately, there is a constructor that takes just Compare and container defaulted to default constructed container (not that there is abundance of constructors).

Shortcut:

Every time iterator is created and std::advance() is used, it is possible to just create and initialize using std::next(). Instead of:

auto iterator = another_iterator;


could be:

auto iterator = std::next(another_iterator, distance);


Default comparison:

It is possible to use default type for the type parameter. Unfortunately in C++11 declaration will be pretty long:

    template<typename RandomIt,
typename Cmp =
typename std::less<std::iterator_traits<RandomIt>::value_type>>
void local_heapsort(RandomIt begin,
RandomIt end,
size_t max_distance,
Cmp cmp = {})


Though it is possible to make specialization in std until switching to newer C++ version. In C++14 it would be just typename Cmp = std::less<>.

Minor feature:

I believe the container type should be also a template parameter, but a last, defaulted to vector, one. In some cases people could deploy sophisticated allocators.

Smaller things:

Usually begin and end are named first and last.

The algorithm could check the size of the first window, so if the range would have small length, the code wouldn't invoke undefined behavior. I believe the decrease maintenance cost will outweigh decrease in performance (which probably won't be visible). If there would be some function to get/check the window length or guidelines documented somewhere, may be check could be removed.

You could try to detach the sliding window into own function. I tried it, but the result was not-good-not-bad, which is kind of weird. It also has strong coupling with my transform iterator, which is not great. So, better design is certainly needed.

The main.cpp is a total carnage, but the static keyword on each function got my special attention. It doesn't look like it is needed. Is it an artifact of writing primarily Java? :)

• static - I thought those functions are not to be visible outside. – coderodde Apr 6 '17 at 16:01
• @coderodde, yes, but I don't think they are worth writing on test files. It is used so rarely that I've never seen them yet. – Incomputable Apr 6 '17 at 16:07

@Incomputable already said most of what I was going to say. So this is not as much an answer as it is curious fact.

Seeing your problem formulation, I remembered one thing from my CS classes, many years ago:

Bubble Sort (and it's friend, Shaker Sort) is fast!... When the data is almost sorted.

And this is exactly the case you have, so I implemented Shaker Sort (which is simply Bubble sort but you bubble both ways) and compared it to your Local Heapsort:

template<typename RandomIt, typename Cmp>
void shakersort(RandomIt begin, RandomIt end, Cmp cmp){
auto changed = true;
while(changed){
changed = false;
auto prev = begin;
auto next = std::next(prev);

while(next != end){
if(cmp(*prev, *next)){
using std::swap;
swap(*next, *prev);
changed=true;
}
next = std::next(next);
prev = std::next(prev);
}

if(changed){
changed = false;
prev = std::prev(prev);
next = std::prev(next);
while(prev != begin){
if(cmp(*prev, *next)){
using std::swap;
swap(*next, *prev);
changed=true;
}
prev = std::prev(prev);
next = std::prev(next);
}
if(cmp(*prev, *next)){
using std::swap;
swap(*next, *prev);
changed=true;
}
}
}
}


Please excuse the code duplication, I'm not asking for a review ;)

And when I run it on my machine:

~ \$ g++ -march=native -O3  --std=c++14 main.cpp && ./a.out
Sorting: 3 2 1 5 4 8 7 6
Sorted:  1 2 3 4 5 6 7 8
Shakersort in 859 milliseconds.
local_heapsort in 987 milliseconds.
Algorithms agree: true


It turns out that although your algorithm is faster than std::sort, Shaker Sort is even faster in this case. This is also the reason that you use Shaker Sort in the final stages of Quick Sort when the ranges get small.

• Out of curiosity: does -DNDEBUG make any difference for std::sort – Maikel Apr 6 '17 at 6:42
• Like, "noooooo!" – coderodde Apr 6 '17 at 15:29
• @coderodde If it's any comfort, your local_heapsort did beat regular bubblesort by 2% or so :) – Emily L. Apr 6 '17 at 20:48

I'm posting a different answer because this one is a benchmarking and comparsion of the presented code, whereas other one is purely about the coding style and nitpicks

Setup:

I won't dig much deep into the setup, the only thing I want to say is that I use libc++, which turns out to be crucial in terms of performance of std::sort. I was watching youtube and listening music as if benchmark was not running, to introduce some load on the system.

I also added Emily's code into the mix.

Results:

The first thing that came to mind was to benchmark the code on increasing sizes and randomly generated window length. It uses simple caching technique to speed up benchmarking. I increased the size with step of 10.

class generator
{
shino::random_int_generator<> gen;
public:
using input_type = std::size_t;

generator()
{}
generator(generator&& other) = default;

std::tuple<std::vector<int>, std::size_t> operator()(input_type input)
{
static std::vector<int> v;
static std::size_t window_length = 1;
if (v.size() == input)
{
return std::make_tuple(v, window_length);
}
v.resize(input);
std::uniform_int_distribution<std::size_t> dist(1, input - 1);
window_length = dist(gen.get_engine());
for (int i = 0; i < v.size(); ++i)
{
v[i] = i;
}

for (std::size_t i = 0; i < v.size(); i += window_length)
{
std::shuffle(&v[i], &v[std::min(i + window_length, v.size())], gen.get_engine());
}

return std::make_tuple(v, window_length);
}
};


And the results were ... Confusing?

As it is possible to see, shaker sort was performing terrible, but local heapsort was slightly outperforming std:sort. In the beginning, I thought one of the sorting algorithms were not doing it correctly, but after debugging all of them produced the same output for given input. So, does it mean that local heapsort is the fastest?

Yes and no, as it turns out ...

If we'll look carefully at what Emily said:

Bubble Sort (and it's friend, Shaker Sort) is fast!... When the data is almost sorted.

(Emphasize mine)

So, how much is that almost? It turns out that there is no definition and no way to calculate sortedness. Nevertheless, the maximum distance actually affects sortedness a lot. In the case above, the maximum distance from sorted position was varying wildly, so I decided to write another generator to check the theory. The window length was increasing by 1, and the vector size was fixed to 50'000 elements.

class different_winsize_generator
{
shino::random_int_generator<> gen;
std::size_t vsize;
public:
using input_type = std::size_t;

different_winsize_generator(std::size_t size):
vsize(size)
{}

different_winsize_generator(different_winsize_generator&& other) = default;

std::tuple<std::vector<int>, std::size_t> operator()(input_type input)
{
static std::vector<int> v(vsize);
static std::size_t window_length = input;
if (window_length == input) //return cached value
{
return std::make_tuple(v, window_length);
}
window_length = input;
std::uniform_int_distribution<std::size_t> dist(1, input - 1);
for (int i = 0; i < v.size(); ++i)
{
v[i] = i;
}

for (std::size_t i = 0; i < v.size(); i += window_length)
{
std::shuffle(&v[i], &v[std::min(i + window_length, v.size())], gen.get_engine());
}

return std::make_tuple(v, window_length);
}
};


And the results were just what I thought!

Here is the last one performing sort on 500'000 integers (excluded from listing below) and window size = [5; 10]

It makes things much more clear.

Conclusion:

Shaker sort should be used on literally sorted (window length <= ~10) data, otherwise probably even insertion sort will be faster. std::sort, not surprisingly to me, lost again. I believe the secret of local heapsort remaining steady was std::priority_queue. The situation might change if the data won't fit into L3 cache, but I didn't have time to benchmark that, since it will take quite a lot of time (1M integers were being very slow on shaker sort)

Full code with debugging (ran separately in debug build), pretty cryptic though:

#include "src/benchmark_v2.hpp"
#include "src/random_engine.hpp"

#include <queue>
#include <iostream>
#include <algorithm>
#include <cassert>
#include <functional>

//author: coderodde http://codereview.stackexchange.com/users/58360/coderodde
template<typename RandomIt, typename Cmp>
void local_heapsort(RandomIt begin,
RandomIt end,
Cmp cmp,
size_t max_distance)
{
using value_type =
typename std::iterator_traits<RandomIt>::value_type;

size_t range_distance = std::distance(begin, end);
max_distance = std::min(max_distance, range_distance);
size_t window_width = max_distance + 1;

auto window_end = begin;
std::priority_queue<value_type, std::vector<value_type>, Cmp> heap(cmp);

auto current_iterator = begin;
auto next_iterator = begin;

if (begin + max_distance <= end)
{
for (auto iter = current_iterator; iter != next_iterator; ++iter)
{
heap.push(*iter);
}
}
else
{
std::cout << "something went wrong with generator... \n";
}

while (next_iterator != end)
{
*current_iterator = heap.top();
heap.pop();
heap.push(*next_iterator);
++current_iterator;
++next_iterator;
}

while (!heap.empty())
{
*current_iterator = heap.top();
heap.pop();
++current_iterator;
}
}

//author: Emily L. http://codereview.stackexchange.com/users/36120/emily-l
template<typename RandomIt, typename Cmp>
void shakersort(RandomIt begin, RandomIt end, Cmp cmp){
auto changed = true;
while(changed){
changed = false;
auto prev = begin;
auto next = std::next(prev);

while(next != end){
if(cmp(*prev, *next)){
using std::swap;
swap(*next, *prev);
changed=true;
}
next = std::next(next);
prev = std::next(prev);
}

if(changed){
changed = false;
prev = std::prev(prev);
next = std::prev(next);
while(prev != begin){
if(cmp(*prev, *next)){
using std::swap;
swap(*next, *prev);
changed=true;
}
prev = std::prev(prev);
next = std::prev(next);
}
if(cmp(*prev, *next)){
using std::swap;
swap(*next, *prev);
changed=true;
}
}
}
}

class generator
{
shino::random_int_generator<> gen;
public:
using input_type = std::size_t;

generator()
{}
generator(generator&& other) = default;

std::tuple<std::vector<int>, std::size_t> operator()(input_type input)
{
static std::vector<int> v;
static std::size_t window_length = 1;
if (v.size() == input)
{
return std::make_tuple(v, window_length);
}
v.resize(input);
std::uniform_int_distribution<std::size_t> dist(1, input - 1);
window_length = dist(gen.get_engine());
for (int i = 0; i < v.size(); ++i)
{
v[i] = i;
}

for (std::size_t i = 0; i < v.size(); i += window_length)
{
std::shuffle(&v[i], &v[std::min(i + window_length, v.size())], gen.get_engine());
}

return std::make_tuple(v, window_length);
}
};

class different_winsize_generator
{
shino::random_int_generator<> gen;
std::size_t vsize;
public:
using input_type = std::size_t;

different_winsize_generator(std::size_t size):
vsize(size)
{}

different_winsize_generator(different_winsize_generator&& other) = default;

std::tuple<std::vector<int>, std::size_t> operator()(input_type input)
{
static std::vector<int> v(vsize);
static std::size_t window_length = input;
if (window_length == input) //return cached value
{
return std::make_tuple(v, window_length);
}
window_length = input;
std::uniform_int_distribution<std::size_t> dist(1, input - 1);
for (int i = 0; i < v.size(); ++i)
{
v[i] = i;
}

for (std::size_t i = 0; i < v.size(); i += window_length)
{
std::shuffle(&v[i], &v[std::min(i + window_length, v.size())], gen.get_engine());
}

return std::make_tuple(v, window_length);
}
};

int main()
{
auto shaker = [](std::vector<int>& v, std::size_t)
{
shakersort(v.begin(), v.end(), std::greater<>{});
std::cout << v.front() << '\n';
};

auto localheapsort = [](std::vector<int>& v, std::size_t max_distance)
{
local_heapsort(v.begin(), v.end(), std::greater<>{}, max_distance);
std::cout << v.front() << '\n'; //produce some side effect
};

auto stdsort = [](std::vector<int>& v, std::size_t)
{
std::sort(v.begin(), v.end());
std::cout << v.front() << '\n'; //ditto
};

auto size_benchmark = shino::benchmarker(generator{}, shaker, localheapsort, stdsort);
auto windowsize_benchmark = shino::benchmarker(different_winsize_generator(50'000),
shaker, localheapsort, stdsort);

for (std::size_t i = 1000; i <= 10'000; i+= 10)
{
size_benchmark.time(i, 3);
//        v1 = gen(i);
//        v2 = gen(i);
//        v3 = gen(i);
//        assert(v1 == v2);
//        assert(v2 == v3);
//
//        std::apply(localheapsort, v1);
//        std::apply(stdsort, v2);
//        std::apply(shaker, v3);
//        assert(v1 == v2);
//        assert(v2 == v3);

std::cout << "current size: " << i << '\n';
}

size_benchmark.save_as<std::chrono::microseconds>("benchmarks/local-heap-sort/size-benchmark/benchmarks.txt",
{"benchmarks/local-heap-sort/size-benchmark/shakersort.txt",
"benchmarks/local-heap-sort/size-benchmark/localheapsort.txt",
"benchmarks/local-heap-sort/size-benchmark/stdsort.txt"},
"vector<int> size", "microseconds");

for (std::size_t i = 5; i <= 1000; ++i)
{
windowsize_benchmark.time(i, 5);
}

windowsize_benchmark.save_as<std::chrono::microseconds>("benchmarks/local-heap-sort/window-size-benchmark/benchmarks.txt",
{"benchmarks/local-heap-sort/window-size-benchmark/shakersort.txt",
"benchmarks/local-heap-sort/window-size-benchmark/localheapsort.txt",
"benchmarks/local-heap-sort/window-size-benchmark/stdsort.txt"},
"window length", "microseconds");
}

• Good, comprehensive testing. Nicely done. – Emily L. Apr 8 '17 at 10:54
• @EmilyL., thanks! It's really exciting for me to get compliment from such experienced programmer as you. – Incomputable Apr 8 '17 at 18:16
• On the last graph with window length = 5, it's interesting to see std::sort win over local_heapsort only for this window size (and smaller I assume). I know many implementations of std::sort use shakersort (or similar) when a partition becomes small (5-10) elements and that might explain why it does well here. The partition steps hardly does any work and whatever pivot strategy is working ideally so the whole thing basically reduces to shakersort with overhead for splitting the partitions. – Emily L. Apr 8 '17 at 20:31