Overlapping time-intervals WITHOUT for/while loops?

The best way to ask this question is via a clear example. Please look at the two timelines (e.g. time in seconds) A and B:

It is clear that the intervals for each timeline are:

 intervals_a = 
    0 1
    1 4
    4 7
    7 9

 intervals_b = 
    0 2
    2 3
    3 5
    5 8

Notice that the first a-interval overlaps the first b-interval. The second a-interval overlaps the first, second, and third b-intervals, and so on.

Ultimately, I need an output which shows the indices of a-intervals which overlap with b-intervals. In this case, we have:

 output = 
    1 1   % 1st a-interval overlaps 1st b-interval
    2 1   % 2nd a-interval overlaps 1st b-interval
    2 2   % 2nd a-interval overlaps 2nd b-interval
    2 3   % 2nd a-interval overlaps 3rd b-interval
    3 3   % etc...
    3 4
    4 4

The big challenge is: The solution cannot contain for/while loops ("why" is irrelevant). Can this be done efficiently using vector / matrix / array / sort, or other tools? Thanks in advance!

Answer 1

[Cedric]

EDIT as it's not a HW

There are several methods for doing this, in particular based on ARRAYFUN/BSXFUN that perform the FOR loop(s) that you mention internally, which usually increases efficiency and code conciseness. As you already developed a working code based on loops (that you could I guess easily rewrite using the two aforementioned functions), I give you an alternative solution based on CUMSUM that you can adapt to your needs:

 a = [0, 1, 4, 7, 9] ;                          % From your initial statement.
 b = [0, 2, 3, 5, 8] ;

 m = 1 + max([a,b]) ;
 intId = zeros(m, 2) ;
 intId(1+[a, m+b]) = 1 ;                        % (linear indexing)
 overlaps = unique(cumsum(intId, 1), 'rows') ;

Note that you might want to call UNIQUE with a 3rd argument 'stable' (if relevant). Executing this code produces:

 overlaps =
     1     1
     2     1
     2     2
     2     3
     3     3
     3     4
     4     4
     4     5
     5     5

Former set of hints

I suspect that it is a HW, so I can't give you a direct answer, but here are a few hints..

• The output doesn't have the size of the input (vectors a and b or intervals arrays), so the solution is probably not some direct, smart type of indexing or relational operation.
• The number of elements of the output is not a multiple of the number of elements of the inputs, so it is unlikely that the solution is based on a RESHAPE, and/or a REPMAT.
• ARRAYFUN and BSXFUN are usually a good way to loop over elements of arrays.
• ARRAYFUN can generate non-uniform output if needed.
• Relational operators are available as functions, e.g. GT for >.

[Hamad]

This exercise forms a small part of a high-frequency trading system I am building, based on lead/lag linkages across stock index futures markets. This is one of my postdoctoral projects, and is certainly not "a HW".

Doing it with nested while-loops is easy, but array tools will enable better parallelism and possibly move onto the GPU later. This is why I ask the question here.

I appreciate your hints, but it is not an answer. Therefore, in this case I cannot click "accept" or "vote" in favor of what you have written. You may wish to delete your answer so that others can have an opportunity to contribute. Many thanks for your guidance. Best wishes.

[Cedric]

Ok, my apologize for the mistake; you write " solution cannot contain for/while loops", which looks quite much like the HW statements that we often see here. I guess that it is the reason why you got no answer before I decided to post a few hints. I updated my original list of hints with an answer based on CUMSUM.

[Hamad]

Many thanks for your contribution. The code runs well for large arrays and faster than my nested while loops. However, there is a problem with your code. If you pass in another test set, suppose

 a = [2,6,9,10];
 b = [1,4,6,7];

then your program yields

 overlaps =

     0     0
     0     1
     1     1
     1     2
     2     3
     2     4
     3     4
     4     4

whereas the correct solution given by the nested loop is

 oo =

     1     1
     1     2
     2     3

Note that "overlaps" contains "oo", but this is not always true particularly for large arrays. Note also that "overlaps" in this case contains zeros and fours. Adding a leading 0 to A and B mitigates this problem, but still yields incorrect intervals. I have not had time to go through the code in detail, but I would say your solution comes VERY close to solving this puzzle. I thank you for that.

[Cedric]

The best I can do at this point is to explain how the code works, so you can tailor it to your needs.

The principle is to create vectors of time interval IDs with a step of one ID per second (as boundaries seem to be integers), and to extract unique configurations when we pair these vectors for a and b.

Taking

 >> a = [0, 1, 4, 7, 9] ;

we want to build the following vector of 10 interval IDs:

 intId_a = [1 2 2 2 3 3 3 4 4 5]
            0 1 2 3 4 5 6 7 8 9  <= that correspond to these times [s]

Note that we should works with column vectors, but I make examples using row vectors because it is simpler to display on the forum. One way to achieve this (intId_a) is to create a vector of zeros ans set to 1 elements that correspond to interval boundaries. For that, we can use a+1 as an index (as a starts at 0):

 >> intId_a = zeros(1, 10) ;
 >> intId_a(a+1) = 1
 intId_a =
     1     1     0     0     1     0     0     1     0     1

and compute a cumulative sum of its elements:

 >> intId_a = cumsum(intId_a)
 intId_a =
     1     2     2     2     3     3     3     4     4     5

Doing the same for b, we get:

 >> b = [0, 2, 3, 5, 8] ;
 >> intId_b = zeros(1, 10) ;
 >> intId_b(b+1) = 1 ;
 >> intId_b = cumsum(intId_b)
 intId_b =
     1     1     2     3     3     4     4     4     5     5

Now if we CAT them together:

 >> intId = [intId_a; intId_b]
 intId =
     1     2     2     2     3     3     3     4     4     5
     1     1     2     3     3     4     4     4     5     5

you realize that it is roughly your solution, but that we should eliminate (at least) columns that appear multiple times. We can do this using UNIQUE on the transpose of this array, setting it up to return unique rows:

 >> unique(intId.', 'rows')
 ans =
     1     1
     2     1
     2     2
     2     3
     3     3
     3     4
     4     4
     4     5
     5     5

You see that this solution is valid if a and b both start at the same instant and that this instant is t=0 (we use a+1 and b+1 for this reason). We also used the same length of 10 for both vectors intId_a and intId_b, which is the overall max time boundary+1. Finally, you see that the end of the intervals is given an ID (5) in this process.

The difference between this example and my solution is that in the latter I compute this 10 as m=1+max([a,b]) and I use linear indexing in a 2D array of zeros in place of CAT-ing intermediary variables, which hopefully makes it more efficient.

Now it's easy to tailor this to your needs, e.g. if the assumption that both a and b start at t=0 falls, or even worse, that they don't start at the same instant. One way to achieve this could be..

 o = -min([a, b]) + 1 ;     % Offset for translat. overall min time to 1.
 n = o + max([a,b]) ;       % Number of seconds in full range of times.
 intId = zeros(n, 2) ;
 intId(o+[a, n+b]) = 1 ;
 overlaps = unique(cumsum(intId, 1), 'rows', 'stable')

now this code will make it possible to work with negative times, or very large initial times, but it will leave ID=0 at the beginning if both time ranges start at different times. This is a good behavior as it tells that for an interval in a given vector there is no vis-a-vis in the other vector. Running it with e.g.

 a = [100, 101, 102, 105, 108] ;
 b = [103, 106, 109] ;

gives

 overlaps =
     1     0
     2     0
     3     0
     3     1
     4     1
     4     2
     5     2
     5     3

which makes sense as the intervals [100-101[, [101-102[, and a chunk of [102-103[ in a have no vis-a-vis in b. The same happens at the end of intervals and in some sense you don't want the smallest of the upper boundaries to be given an interval ID. So we have to filter these 'boundary cases' out.

There are three options for that: pre-processing the inputs to prevent boundary cases to create unwanted solutions, trying to complexify a bit the code that described above so the computation manages boundary cases, and last, post-processing the output to eliminate unwanted cases.

Here is an example how to code the first option, where we reduce a and b to their contiguous part and eliminate the upper boundaries so they don't add interval IDs.

 % - Pre-process a and b.
 lb = max(min(a), min(b)) ;      % Lower boundary.
 ub = min(max(a), max(b)) ;      % Upper boundary.
 app = [lb, a(a>lb & a<ub)] ;    % Pre-proc. version of a.
 bpp = [lb, b(b>lb & b<ub)] ;    % Pre-proc. version of b.

 % - Build array of overlaps.
 o = -min([app,bpp]) + 1 ;       % Offset to shift a and b so min ID is 1.
 n = o + max([app,bpp]) ;        % Number of seconds in full range of times.
 intId = zeros(n, 2) ;
 intId(o+[app,n+bpp]) = 1 ;
 overlaps = unique(cumsum(intId, 1), 'rows', 'stable')

This code outputs:

 overlaps =                     overlaps =
     1     1                        1     1
     2     1                        1     2
     2     2                        2     3
     2     3
     3     3
     3     4
     4     4

respectively for the first and second cases that you submitted.

Note that a feature (maybe unwanted) of this solution is that interval IDs start (at 1) at the first overlap, so intervals of b that cover times before the lowest boundary of a would not be counted.

If you fully understand what we've done here, it won't be complicated to finalize tailoring this approach to your needs.

[Hamad]

Many thanks for your answer!

[Cedric]

My pleasure!

Note that I added 'stable' in the call to UNIQUE (which prevents this function to order its output) and I am not sure that you need/want it.

[Cedric]

I just realized that with the preprocessing, you can get a little improvement by replacing the line

 o = -min([app,bpp]) + 1 ;       % Offset to shift a and b so min ID is 1.

with

 o = -lb + 1 ;                   % Offset to shift a and b so min ID is 1.

[Hamad]

Indeed. You may be interested to know that one application for this exercise is based on one of my recent PhD projects on Statistical Arbitrage. Thanks again and best wishes. Hamad.

[Cedric]

Thank you for the link, it is always interesting to see in which context what we help coding is finally used!

Best wishes,

Cedric

[Hamad]

Cedric, a further MATLAB challenge awaits you, if you feel up to it!

https://www.mathworks.co.uk/matlabcentral/answers/69413-how-many-times-does-each-number-appear

Best wishes, Hamad

[Image Analyst]

It's your homework problem. Don't you feel up to it? Why challenge Cedric to do it for you? At least put some effort into it and show some attempt at solving it, and then we can correct that or suggest improvements to it. Anyway, I did give you a hint. It's your move now.

Answer 2

[Nima Nikmehr]

Cedric's answer was intersing. Is this code usable for any number of matrices (a,b,c,d,...)?

[Cedric]

Yes, but you have to add the correct offsets when building intId, for example

intId(1+[a, m+b, 2*m+c, 3*m+d]) = 1 ;

The best way to understand the approach is to run my original example and to display intId before the call to CUMSUM and after. Running

a = [0, 1, 4, 7, 9] ;                          % From your initial statement.
b = [0, 2, 3, 5, 8] ;
m = 1 + max([a,b]) ;
intId = zeros(m, 2) ;
intId(1+[a, m+b]) = 1

we get

intId =
     1     1
     1     0
     0     1
     0     1
     1     0
     0     1
     0     0
     1     0
     0     1
     1     0
     

This shows that we built an array of zeros with 1s at locations where the interval changes. Then, by using CUMSUM along dimension 1, we get interval IDs (see the doc of CUMSUM to understand why if needed):

cumsum(intId, 1)

which outputs:

ans =
     1     1
     2     1
     2     2
     2     3
     3     3
     3     4
     3     4
     4     4
     4     5
     5     5

Now this output has duplicates (each row that has no interval change has the same values as the row above), so we pass it to UNIQUE (by row) to get only the rows with changes.

Back to the first block of code, we see that one way to place the 1s without computing the row/col indices of each 1 (which is easy to do but unnecessary), is to use vectors a and b as linear indices (again, seach for that topic online if you are not familiar). For this, however, we must add an offset to b to get values that index linearly the second column of intId, and this offset is the number of rows, m. If we needed a 3rd, 4th, etc column, for extra vectors c, d, etc., we would just have to add the correct offsets (2*m , 3*m, etc.), which is what the expression at the top of this comment does.

Now if vectors don't all start at 0 or have different lengths, it is no big deal but you have to undertsand the approach and update the code a bit, as illustrated in my comments to Hamad.

[Nima Nikmehr]

Thank you Cedric! Your code works. I appreciete your time on this question.