Subdividing intervals that contain the largest error values

Question

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?

Answer 1

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.

Answer 2

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.

Answer 3

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).

Answer 4

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);
}