Matlab: Sum corresponding values if index is within a range

Question

I have been going crazy trying to figure a way to speed this up. Right now my current code talks ~200 sec looping over 77000 events. I was hoping someone might be able to help me speed this up because I have to do about 500 of these.

Problem: I have arrays (both 200000x1) that correspond to Energy and Position of a hit over 77000 events. I have the range of each event separated into two arrays, event_start and event_end. First thing I do is look for the position in a specific range, then I put the correspond energy in its own array. To get what I need out of this information, I loop through each event and its corresponding start/end to sum up all the energies from each it hit. My code is below:

    indx_pos = find(pos>0.7 & pos<2.0);
    energy = HitEnergy(indx_pos);

    for i=1:n_events
        Etotal(i) = sum(energy(find(indx_pos>=event_start(i) …
        & indx_pos<=event_end(i))));    
    end

Sample input & output:

% Sample input 
% pos and energy same length

n_events = 3;
event_start = [1 3 7]';
event_end   = [2 6 8]';

pos = [0.75 0.8 2.1 3.6 1.9 0.5 21.0 3.1]';
HitEnergy = [0.002 0.004 0.01 0.0005 0.08  0.1 1.7 0.007]';


%  Sample Output
Etotal = 0.0060
         0.0800
              0

Are the start/end ranges overlapping? If not, it's easy to do this with `accumarray`. Also, you don't need the call to `find` inside your loop. — Jonas, Apr 01 '15 at 14:28
The start/end ranges are not overlapping. I've never used accumarray but I'll look into that. Thanks. — UtahnFloridian, Apr 01 '15 at 17:26

Divakar · Accepted Answer · 2015-04-01T19:36:41.583

2

Approach #1: Generic case

One approach with bsxfun and matrix-multiplication -

mask = bsxfun(@ge,indx_pos,event_start.') & bsxfun(@le,indx_pos,event_end.')
Etotal = energy.'*mask

This could be a bit memory-hungry if indx_pos has lots of elements in it.

Approach #2: Non-overlapping start/end ranges case

One can use accumarray for this special case like so -

%// Setup ID array for use in accumarray later on
loc(numel(pos))=0; %// Fast pre-allocation scheme
valids = event_end+1<=numel(pos);
loc(event_end(valids)+1) = -1*(1:sum(valids));
loc(event_start) = loc(event_start)+(1:numel(event_end));
id = cumsum(loc);

%// Set elements as zeros in HitEnergy that do not satisfy the criteria:
%// pos>0.7 & pos<2.0
HitEnergy_select = (pos>0.7 & pos<2.0).*HitEnergy(:);

%// Discard elments in HitEnergy_select & id that have IDs as zeros  
HitEnergy_select = HitEnergy_select(id~=0);
id = id(id~=0);

%// Accumulate summations as done inside the loop in the original code 
Etotal = accumarray(id(:),HitEnergy_select);

edited Apr 01 '15 at 19:36

answered Apr 01 '15 at 15:01

Divakar

218,885
19
262
358

@ABD_worley Check out newly added approach too! – Divakar Apr 01 '15 at 18:10
Clearly there was a possibility to simplify `accumarray` for the non-overlapping case. I just couldn't think of it! – knedlsepp Apr 01 '15 at 18:27
@ABD_worley So, did you get any speedup with the second approach? Kinda curious about it! – Divakar Apr 01 '15 at 18:37
Hm. Interestingly I don't see that much of a difference, even though I thought my approach is more computationally intensive. – knedlsepp Apr 01 '15 at 18:42
1

@Divakar The second approach was crazy fast! 1.5 sec compared to my 200 sec. Thanks. – UtahnFloridian Apr 01 '15 at 19:12
@ABD_worley Wow, that's massive! Appreciate the update!! – Divakar Apr 01 '15 at 19:13

score 2 · Answer 2 · edited May 23 '17 at 11:52

The problem is that for every event you are searching the entire vector indx_pos. Constrain your search inside the loop to only the range from event_start(i) to event_end(i):

for i = 1:n_events
    I = event_start(i):event_end(i);
    posIIsWithinRange = pos(I)>0.7 & pos(I)<2.0;
    Etotal(i) = sum(HitEnergy(I(posIIsWithinRange)));
end

You could also use a vectorized version based on run length decoding and vectorizing the notion of colon. (Download the functions coloncatrld and runLengthDecode.)

I = coloncatrld(event_start, event_end);
energy = HitEnergy(I);
eventNum = runLengthDecode(event_end - event_start+1);
posIIsWithinRange = pos(I)>0.7 & pos(I)<2.0;
Etotal = accumarray(eventNum(posIIsWithinRange), energy(posIIsWithinRange), [n_events,1]);

This is similar to Divakar's Approach #2 with the addition that it should work for overlapping ranges too.

These two method you presented were about 2 sec too compared to my 200 sec. You guys are awesome. Thank you so much! — UtahnFloridian, Apr 01 '15 at 19:16

Matlab: Sum corresponding values if index is within a range

2 Answers2

Linked