Speed Up Finding Number of Elements Between Values Matlab

user_4700841 0
Viewed 69

I've created what is a fairly simple MATLAB script to simulate the behaviour discussed in this question over on Maths SE.

clearvars;
samples = 1000;
x = 256;
r=exprnd(1/20e6,1,samples); % Generate exponentially distributed randoms.
endTime = sum(r);
quickMean=sum(r(1:x))/x; % Quick calc the mean and median.
quickMedian=0.693 * quickMean;
p = cumsum(r); % Convert event deltas into timestamps
bitstream = false(1,samples);
time = 0;
lastTime = 0;
for i = 1:samples
    lastTime = time;
    time = time + quickMedian;
    if (numel(p(p < time & p > lastTime)) > 0)
        bitstream(i) = true;
    end
    if (time > p(end))
        break
    end
end
ratio = sum(bitstream)/samples;

The script seems to work, however, if I use a large number of samples (say a million), which would be beneficial, it really crawls.

I'm assuming that the problematic statement is this one:

p(p < time & p > lastTime)

Is there a more efficient way to check if any elements in an array fall between two values?

matlab
2 Answers

Let’s examine that whole expression:

numel(p(p < time & p > lastTime)) > 0

We can separate that out for clarity:

I = p < time & p > lastTime;
tmp = p(I);
n = numel(tmp);
n > 0

Here, the creation of tmp is pretty expensive: it looks at where I is true, and it copies those elements over to a new array. But the only thing you do with this array is seeing how many elements it has. Logically n will be equal to the number of true elements in I. And you don’t really need this number, you just need to know if it’s larger than 0. That is, you want to know if any of the elements in I is true. You can do so with any:

any(p < time & p > lastTime)
user_73287820

As I mentioned in my comment, we can use the fact that p is monotonically increasing and ignore values less than lastTime. If we find the last value for which p < time, only the values to the right can be greater than* time on the next iteration (lastTime).

clearvars;
samples = 1000;
x = 256;
r=exprnd(1/20e6,1,samples); % Generate exponentially distributed randoms.
endTime = sum(r);
quickMean=sum(r(1:x))/x; % Quick calc the mean and median.
quickMedian=0.693 * quickMean;
p = cumsum(r); % Convert event deltas into timestamps
bitstream = false(1,samples);
time = 0;
lastTime = 0;

% code is the same up to here ---

lastTimeIdx = 1; % index of (last value < lastTime) + 1
for i = 1:samples
    lastTime = time;
    time = time + quickMedian;

    valsInRange = p(lastTimeIdx:end) < time; % p > lastTime & p < time
    timeIdx = find(valsInRange, 1, 'last'); % returns [] or index

    if ~isempty(timeIdx)
        bitstream(i) = true;
        lastTimeIdx = lastTimeIdx + timeIdx; % update start of next search
    end
    if (time > p(end))
        break
    end
end
ratio = sum(bitstream)/samples;

*Actually, this is "greater than or equal to", but since the values of p are unique, they are the same thing.

Okay, I just tried histc in Octave. I'm embarrassed to say that it's ridiculously fast. Like 4 orders of magnitude faster. Here's the code I used, but histc is deprecated in MATLAB, and the binning for histcounts is different, so you may have to play with it a bit.

bitstream_hist = histc(p, [0:samples]*quickMedian) > 0;
bitstream_hist = bitstream_hist(1:samples);

One million samples finishes in a fraction of a second. Sorry I didn't think of this sooner.

user_13770970
Related