Subdividing intervals that contain the largest error values

I'm trying to accomplish the following:

Start out with 1 interval on the real line. This interval is processed in some function and out comes an array of errors (size < 10 usually). The maximum element of this array is also stored and the interval is put onto a queue sorted by this maximum error.

Then we proceed by taking the interval with the largest error from the queue and split this interval into 2 equal parts. We then do the same operations as above on both these intervals and then put them on the queue as well.

I will demonstrate how I solved this below. But high performance is essential in this application, so I would very much appreciate your input on how I can speed things up here. Is std::vector a bad choice for example?

Here is my Interval class.

Interval.h

#include <queue>

class Interval
{
public:
Interval(size_t dim, double a, double b);
void UpdateInterval(double, double);

std::vector<double> errors;      // vector of length dim
double maxError, xmin, xmax, center;
private:
friend bool operator<(const Interval& lhs, const Interval& rhs) { return lhs.maxError < rhs.maxError; }
};

Interval.cpp

#include "Interval.h"

Interval::Interval(size_t fdim, double a, double b)
{
xmin = a;
xmax = b;
center = (xmax + xmin) / 2;

errors.resize(fdim);
maxError = std::numeric_limits<double>::max();
}

void Interval::UpdateInterval(double a, double b)
{
xmin = a;
xmax = b;
center = (xmax + xmin) / 2;
}

In this simplified example I simply populate the vectors with random doubles. I don't care about the performance of this step since I will populate these vectors in a more sophisticated way.

Source.cpp

#include "Interval.h"

typedef std::priority_queue < Interval, std::vector<Interval>, std::less<Interval>> PriorityQueue;

static void SplitInterval(Interval &i1, Interval&i2)
{
size_t dim = i1.errors.size();
double a = i1.xmin;
double b = i1.xmax;
double center = (a + b) / 2.0;

i1.UpdateInterval(a, center);
i2 = Interval(dim, center, b);
}

int main()
{
{
PriorityQueue queue;

Interval i1(5, 0.0, 100.0);
Interval i2(5, 0.0, 100.0);

std::generate(i1.errors.begin(), i1.errors.end(), std::rand);
queue.push(std::move(i1));

for (int counter = 0; counter < 1000; counter++)
{
// A trick because top() returns a const T&
Interval i1 = std::move(const_cast<Interval&>(queue.top()));
queue.pop();

// Splits i1 in half and put the left halft in i1 and the right half in i2
SplitInterval(i1, i2);

// Populate the error arrays with random data
std::generate(i1.errors.begin(), i1.errors.end(), std::rand);
std::generate(i2.errors.begin(), i2.errors.end(), std::rand);

// Calculate the maximum error of each interval
i1.maxError = *max_element(i1.errors.begin(), i1.errors.end());
i2.maxError = *max_element(i2.errors.begin(), i2.errors.end());

// Push the interval onto the queue and it will be sorted according to it's maxError value
queue.push(std::move(i1));
queue.push(std::move(i2));
}
}
return 0;
}

I decided to push( value_type&& value ) instead of push( const value_type &value ) since it's much more efficient (~30% faster). We avoid doing a lot of copies. The drawback is the little trick we have to do for each top() operation since top() returns const value_type &. You have to be really delicate here not to ruin the order of the queue.

You might wonder why I need to have the error array as a class member since I never do anything with it? In my more realistic application I need to subtract the errors from a global error before splitting the intervals in two. So that is the reason.

Really appreciate you feedback on how I can optimize the performance of this in terms of CPU time. Tweak it? Rewrite it from the ground?

• Is there an upper bound on the size of errors? – T.C. May 20 '16 at 23:36
• Not in theory, but I have never encountered a situation where we have needed this vector to be greater than 5 elements – DoubleTrouble May 21 '16 at 18:48
• In that case, consider using a std::array<double, 5> or a small vector if you don't want a size upper bound (more recent versions of boost has one). – T.C. May 21 '16 at 18:55
• Using std::array wade this code run ~40% fast. That is a significant increase in performance! – DoubleTrouble May 21 '16 at 21:47

Trips to the heap are expensive; avoid them if possible. When errors is not too big, the cost of moving a vector is probably similar to the cost of just copying the data.

1) If the dimension of errors is known at compile time and you don't need intervals with different errors sizes to be the same type, then make it a template parameter, and make errors a std::array:

template<size_t Dim>
struct Interval {
// ...
std::array<double, Dim> errors;
// ...
};

2) If it is not known at compile time, but you know an upper bound (5, according to your comments) that is small, then use a std::array or a boost::static_vector or a boost::small_vector. In all three the storage is embedded inside the object itself so that you don't need a heap trip. With std::array you may need to track the actual size yourself. small_vector switches to heap allocation if the embedded storage is not large enough; static_vector throws an exception in that case.

3) If you do either of the above, then you don't need the icky std::move/const_cast combo (well, if you use small_vector and pass a large dimension, it might still be of benefit, but in the usual case you don't need it).

4) If you have some idea about the maximum number of elements in your queue, you can try to reserve() the needed capacity upfront to avoid reallocations. In all sane implementations of vector, this should work:

std::vector<Interval> buffer;
buffer.reserve(1001); // or whatever the maximum is
PriorityQueue queue(std::less<Interval>(), std::move(buffer));

You can also derive a class from PriorityQueue, add a reserve member function that just calls c.reserve(), and use that, since the adapters expose the underlying container as a protected data member named c.

• I'm glad you submitted a comment. I will award you the bounty if nothing else happens in this topic! I have a few followup questions though: Is reading/writing in the free store expensive as well, or only allocation? My tests indicate that moving is faster than copying even if the vector is small (5 elements). Is that surprising? Maybe I should do a test with a simple C-style dynamic array and see what kind of performance I get. – DoubleTrouble May 23 '16 at 19:20
• @DoubleTrouble Copying a vector means doing a fresh allocation - that's the expensive part. Copying the data shouldn't be too expensive. – T.C. May 23 '16 at 19:21
• Preallocating space for the the elements in the queue is not something I have considered, thanks for showing me that! – DoubleTrouble May 23 '16 at 19:22

Mmmmm

// A trick because top() returns a const T&
Interval i1 = std::move(const_cast<Interval&>(queue.top()));

A trick that is illegal.

Also I am not convinced this does not break the underlying priority queue. As the top value when popped is unwound out of a heap using the compare operation that inspects the internal members of the object. Its not until it is removed from the heap that it is removed from the underlying vector.

Really not much else to review.
I am not sure what performance things you expect us to find the code here has very little to do with performance.

• It's gross but legal, I think. The underlying object is not const, and the ordering uses <, which only looks at maxError, which is not affected by the move. – T.C. May 20 '16 at 20:30
• @T.C.: Are you sure maxError is not affected by move. There is no user defined version and thus the compiler generated move (Which I am not sure of the definition). – Martin York May 20 '16 at 22:16
• Yes, that just do a memberwise move. And moving a double is a copy. – T.C. May 20 '16 at 22:29
• @T.C.: OK may not be broken as it stands (through luck). But it is now highly vulnerable to maintenance changes. The easiest move implementation will use the swap technique which will leave change the maxError state. When that changes happens you will get no warning that your code is now suddenly broken. – Martin York May 20 '16 at 23:59
• Thanks for your comments Loki and T.C. To make this less dangerous I could just write my own move constructor and not move maxError. – DoubleTrouble May 21 '16 at 18:42

Let me offer a bit of criticism on Interval: use an initializer list (the impact of performance is likely tiny though). Also, could maxError be a constant? One possible rewrite would be this (I moved the member variables to the private side, but perhaps you don't care about this).

class Interval
{
public:
Interval(size_t dim, double a, double b)
: xmin(a), xmax(b), center((xmin + xmax) / 2), errors(dim), maxError(std::numeric_limits<double>::max())
{

}

void UpdateInterval(double a, double b)
{
xmin = a;
xmax = b;
center = (xmax + xmin) / 2;
}

private:
friend bool operator<(const Interval& lhs, const Interval& rhs) { return lhs.maxError < rhs.maxError; }

std::vector<double> errors;      // vector of length dim
double xmin;
double xmax;
double center;
const double maxError;
};

In this way you avoid the (in this case) useless default initializations of the member variables the compiler performs, and take advantage of constructors (e.g., for the vector).

This will not impact performance greatly. It's challenging to give a better suggestion, as I don't know what your actual use case is. But do take the engineering approach: if what you have doesn't meet your needs, try something else. Concretely, focus on data structures here. I'd bet it's faster to use a single vector, whose size is a multiple of two storing "pairs" (start1, end1) etc. instead of using Interval objects. (Perhaps you'd need a second vector for some "manual" bookkeeping, but you should get the idea).

• Even though this comment doesn't concern performance that much I still appreciate the feedback. Thanks! – DoubleTrouble May 23 '16 at 19:10

Your Interval class has a member center which stores (xmax + xmin) / 2. However, in your SplitInterval method you're not using the stored value. Instead, you're working out the same thing and using that to perform the split. Either, lose the member from Interval, or update SplitInterval to use the already calculated member value.

static void SplitInterval(Interval &i1, Interval&i2)
{
size_t dim = i1.errors.size();
double a = i1.xmin;
double b = i1.xmax;
/*********** This is already in i1.center  ************/
double center = (a + b) / 2.0;

i1.UpdateInterval(a, center);
i2 = Interval(dim, center, b);
}
• Hi, thanks for your comment. I agree with you completely. I need it as a member for other reasons, so I should definitely update SplitInterval – DoubleTrouble May 21 '16 at 18:38