How to optimize/parallelize the following clustering/joining algorithm:

Question

I have relatively-small algorithm that takes up ~60% of the total run-time of my scientific code (57 lines of 3600), so I would like to find a way to optimize what I'm doing and make the code order-independent so that I can apply a cilk_for parallel strcture.

Here's what it does, verbally: I have an std::vector of pointers to custom objects called Segment (vector<Segment*> newSegment). Each Segment contains a std::vector of integers (mesh indices). In this function, I would like to find any Segment that overlaps with any another, with overlap being defined as the member indices overlapping on the number line. If they do overlap, I would like to join them together (insert the A.indices into B.indices) and delete one (delete A).

ex. 1: A.indices={1,2,3} B.indices={4,5,6} do not overlap; do nothing

ex. 2: A.indices={1,2,4} B.indices={3,5,6} do overlap; A= deleted B.indices={1,2,3,4,5,6}

The overlaps are sparse, but existent.

Here's the current code:

main algorithm:

//make sure segments don't overlap
for (unsigned i = 0; i < newSegment.size(); ++i) {
    if (newSegment[i]->size() == 0) continue;
    for (unsigned j = i + 1; j < newSegment.size(); ++j) {
        if (newSegment[i]->size() == 0) continue;
        if (newSegment[j]->size() == 0) continue;
        int i1 = newSegment[i]->begin();
        int i2 = static_cast<int>(newSegment[i]->end());
        int j1 = newSegment[j]->begin();
        int j2 = static_cast<int>(newSegment[j]->end());
        int L1 = abs(i1 - i2); 
        int L2 = abs(j1 - j2); 
        int dist = max(i1,i2,j1,j2) - min(i1,i2,j1,j2);

        //if overlap, fold segments together
        //copy indices from shorter segment to taller segment
        if (dist <= L1 + L2) {
            unsigned more, less;
            if (newSegment[i]->slope == newSegment[j]->slope) {
                if (value_max[i] > value_max[j]) {
                    more = i;
                    less = j;
                } else {
                    more = j;
                    less = i;
                }
            } else if (newSegment[i]->size() == 1) {
                more = j; less = i;
            } else if (newSegment[j]->size() == 1) {
                more = i; less = j;
            } else assert(1 == 0);
              while(!newSegment[less]->indices.empty()) {
                unsigned index = newSegment[less]->indices.back();
                newSegment[less]->indices.pop_back();
                newSegment[more]->indices.push_back(index);
            }
        }
    }

}//end overlap check

//delete empty segments
vector<unsigned> delList;
for (unsigned i = 0; i < newSegment.size(); ++i) {
    if (newSegment[i]->size() == 0) {                            //delete empty
        delList.push_back(i);
        continue;
    }
}
while (delList.size() > 0) {
    unsigned index = delList.back();
    delete newSegment.at(index);
    newSegment.erase(newSegment.begin() + index);
    delList.pop_back();
}

Relevant Segment object class definition and member functions:

class Segment{

    public:
    Segment();
    ~Segment();

    unsigned size();
    int begin();
    unsigned end();
    std::vector<int> indices;
    double slope;
};

int Segment::begin() {
    if (!is_sorted(indices.begin(),indices.end()))      std::sort(indices.begin(),indices.end());
    if (indices.size() == 0) return -1; 
    return indices[0];
}

unsigned Segment::end() {
    if (!is_sorted(indices.begin(),indices.end()))    std::sort(indices.begin(),indices.end());
    return indices.back();
}

unsigned Segment::size() {
    unsigned indSize = indices.size();
    if (indSize == 1) {
        if (indices[0] == -1) return 0;
    }   
    return indSize;
}

Ideas:

1. Since I don't care about the order of the Segment objects, they could be in an orderless container?
2. In my algorithm, I find overlap by looking at the first and last indices of each segment. I do an std::is_sorted (and then maybe a std::sort) when I fetch the indices because the list can change when more indices are inserted. Maybe I could put the indices in a std::set rather than std::vector to save the explicit sort-checking/sorting?

3. I'm pretty sure that by editing the indices as I go, this makes it order-dependent. Perhaps, I could break the code into the following organization using the concept of an undirected graph to make it order-independent:

   • edge discovery (without modifying indices)
   • join clusters of connected nodes (Segment objects that overlap) using a graph traversal
   • delete empty Segment objects

Questions

1. Are either of the ideas above worthwhile or negligible to performance?
2. How else can I optimize it?
3. How (if not the above) can I make the algorithm order-independent?

Answer 1

The is_sorted() function is probably expensive, and so you should avoid it. Why not sort everything in one go right at the beginning before entering the loops?

The best way to optimize your code is by inventing a new algorithm which avoids the nested loops of N, because that has a complexity of O(N^2) (see "big-Oh notation".) See Bart van Ingen Schenau's comment below on how to achieve this.

Answer 2

I came to identical algorithm than @BartVanIngenSchenau in this comment Basically sort the set of segments based on the min element of each segment. Then two adjacent element overlap if and only if Segment[i].max >= Segment[i+1].min

But I would like to add that sorting looks unnecessary at all and only keeping the max and the min element. Just update them when merging segments. (segment1+segment2).min = min(segment1.min,segment2.min) and (segment1+segment2).max = max(segment1.max,segment2.max) Moreover if the segment are sorted by min element you have (Segment[i]+Segment[i+1]).min = segment[i].min (but this last thing could be premature optimization.) I noted + the merge of two segments.

For cache locality, the best for merging might be to have a layout simimar to the following layout

ptr_to_2nd_segment
n_elt_of_1st_segment,
min_elt_of_1st_segment,
[
[other_elts_of_1st_segment,]
max_elt_of_1st_segment,]

ptr_to_3rd_segment
n_elt_of_2nd_segment,
min_elt_of_2nd_segment,
[
[other_elts_of_2nd_segment,]
max_elt_of_2nd_segment,]

...

merging two elements in this configuration would be quite simple, it would be just a matter of updating ptr to next element, adding the number of elements,shifting the second segment elements and swapping the max elements if necessary. That would let some junk after each merge (8 bytes on 32 bits architecture and 16 bytes on 64 bits architecture). Knowing if you can support such junk is application dependent (moreover, you could doing a kind of garbage collection between two iteration of the algorithm.)

For parallelizing, once the set of segments are sorted by min element, you can divide in n part the set of segments and doing the merge independently. Then only merge at the border of each parts. But as @MikeNakis says in this comment as merging is quite memory bound, they might not be well parallelize