How to perform logical AND on intervals of contiguous locations

I have the following problem. Let's say I have the arrays

x = logical([0, 1, 0,0, 1,1, 0,0,0, 1,1,1, 0])
y = logical([0, 1, 1,0, 1,1, 0,1,0, 1,0,1, 0])

Array X has three intervals of 1's with indices 2:2, 5:6, and 10:12. I want to apply an "interval AND" operation to X, based on Y, in the following sense: for each interval of ones in X, if any element in Y is zero in that interval, the whole interval is zeroed, i.e., Z = intervalAND(X,Y) should be the same as

z = logical([0, 1, 0,0, 1,1, 0,0,0, 0,0,0, 0])

Let me explain. Since all(Y(2:2)) = 1, it produces ones in Z(2:2). The same happens in the second interval (5:6): Both Y(5) and Y(6) are true, producing ones in Z. However, there is a zero in Y(10:12) which zeroes the whole interval Z(10:12).

I know how to do it with a for loop:

d = diff(x);
pos = find(d == 1);
neg = find(d == -1);
z = x;
for k = 1:length(neg)
    interval = pos(k)+1 : neg(k);
    if ~all(y(interval))
        z(interval) = false;
    end
end

However, I need to vectorize it to make it run faster (I am working with huge arrays). Does someone know how to compute Z without using a for/while loop?

[James Tursa]

How huge? Are you open to a mex routine?

[Arturo Camacho Lozano]

To make it simpler, I explained it using an array, but x is actually a column of a matrix X of size 2000 times 100 (approx.). The script produces many matrices like that one. To give you an idea of the problem, whatever is inside the for loop is being called 40 million times (in a 4 minutes song, which is the application). I would prefer not using a mex routine, because Matlab is being used only for prototyping.

[James Tursa]

Is the algorithm running on each column individually, or running across the entire matrix as a whole? I.e., does the 1's logic extend across columns?

[Arturo Camacho Lozano]

Each column is independant. The same column vector Y is applied to all columns, though.

Answer 1

[Bruno Luong]

x = logical([0, 1, 0,0, 1,1, 0,0,0, 1,1,1, 0])
y = logical([0, 1, 1,0, 1,1, 0,1,0, 1,0,1, 0])

code without loop or groupping, on my bench test about 3 time faster than Stephen's accumarray solution

i = find(diff([0 x 0]));
n = histc(find(~y), i);
j = [1;-1]*(n(1:2:end)==0);
if x(end)
    i(end)=[];
    j(end)=[];
end
z = logical(cumsum(accumarray(i(:),j(:),[length(x),1])));

[Bruno Luong]

Code adapted for array (working along columns)

X = rand(2000,100) > 0.2;
Y = rand(2000,100) > 0.2;

% work along columns
tic
[m,n] = size(X);
Z = false(size(X));
for c=1:n
    i = find(diff([0; X(:,c); 0]));
    n = histc(find(~Y(:,c)), i);
    n = n(:).'; % EDIT: bug correction when find(~Y(:,c)) is a scalar
    j = [1;-1]*(n(1:2:end)==0);
    if X(end,c)
        i(end)=[];
        j(end)=[];
    end
    Z(:,c) = cumsum(accumarray(i(:),j(:),[m,1]));
end
toc % Elapsed time is 0.047264 seconds.

If Y columns unchanged, one can do

find(~Y(:,c))

once outside the loop.

[Bruno Luong]

X = rand(2000,100) > 0.2;
Y = rand(2000,100) > 0.2;

Fully vectorized code for matrices

(from tic/toc: Elapsed time is 0.010189 seconds)

[m,n] = size(X);
z = false(1,n);
Xpad = [z; X; z];
Ypad = [z; Y; z];
[i,j] = find(diff(Xpad));
k = i+(j-1)*(m+2); % == sub2ind(size(Xpad),i,j);
cX = cumsum(Xpad,1);
cY = cumsum(Ypad,1);
k1 = k(1:2:end);
k0 = k(2:2:end);
b = cY(k0)-cY(k1) < cX(k0)-cX(k1);
d = repmat([1;-1],[1 size(b,1)]);
d(:,b) = 0;
d = accumarray([i(:),j(:)],d(:),[m+1 n]);
Z = cumsum(d(1:m,:),1);

[Arturo Camacho Lozano]

How can this last version be modified for a single column vector y?

[Bruno Luong]

Matrix X and single column y

clear
X = rand(2000,100) > 0.2;
y = rand(2000,1) > 0.2;

tic
[m,n] = size(X);
z = false(1,n);
Xpad = [z; X; z];
Ypad = [0; y; 0];
[i,j] = find(diff(Xpad));
k = i+(j-1)*(m+2); % == sub2ind(size(Xpad),i,j);
cX = cumsum(Xpad,1);
cY = cumsum(Ypad,1);
k1 = k(1:2:end);
k0 = k(2:2:end);
i1 = i(1:2:end);
i0 = i(2:2:end);
b = cY(i0)-cY(i1) < cX(k0)-cX(k1);
d = repmat([1;-1],[1 size(b,1)]);
d(:,b) = 0;
d = accumarray([i(:),j(:)],d(:),[m+1 n]);
Z = cumsum(d(1:m,:),1);
toc

[Bruno Luong]

Note: You might convert Z to logical

[Bruno Luong]

Even better

[m,n] = size(X);
z = false(1,n);
[i,j] = find(diff([z; X; z]));
cY = cumsum([0; y; 0],1);
i1 = i(1:2:end);
i0 = i(2:2:end);
b = cY(i0)-cY(i1) == i0-i1;
d = [b, -b]';
d = accumarray([i(:),j(:)],d(:),[m+1 n]);
Z = cumsum(d(1:m,:),1);

[Bruno Luong]

One more variation

[m,n] = size(X);
z = false(1,n);
[i,j] = find(diff([z; X; z]));
cY = cumsum([0; y; 0],1);
i1 = i(1:2:end);
i0 = i(2:2:end);
b = repelem(cY(i0)-cY(i1) == i0-i1, 2, 1) & i<=m;
d = repmat([1;-1], length(i1), 1);
Z = cumsum(accumarray([i(b),j(b)], d(b), [m n]), 1);

[Arturo Camacho Lozano]

What is the advantage (if any) of the latter? (The code looks more complex and verbose.)

[Bruno Luong]

Faster. It does not create unecessary elements to accumulate then removed.

The one before is still OK if you prefer readable code.

[Arturo Camacho Lozano]

Thanks four answer. This implementation works way faster than other proposed implementations that use "splitapply". Besides, it was easy to port to the target language (Python).

Answer 2

[Mohammad Sami]

You can group based on the values of x.

gid = cumsum(x ~= circshift(x,1));
if(gid(1) == 0)
gid = gid + 1;
end
a = splitapply(@min,y,gid);
z = a(gid);

[Arturo Camacho Lozano]

Great! It needs a small adjusment, though. The second argument to splitapply should be x&y:

a = splitapply(@min, x&y, gid);

Let me clarify. If we change the fourth element of Y to 1, i.e.

x = logical([0, 1, 0,0, 1,1, 0,0,0, 1,1,1, 0])

y = logical([0, 1, 1,1, 1,1, 0,1,0, 1,0,1, 0])

the output should be the same as before:

z = logical([0, 1, 0,0, 1,1, 0,0,0, 0,0,0, 0])

because there was already a 0 in that position of X.

(Please edit the answer to accept it).

Answer 3

[Matt J]

Using group1s from

https://www.mathworks.com/matlabcentral/fileexchange/78008-group1s

>> xg=group1s(x)+1;
>> yg=splitapply(@all,y,xg);
>> z=yg(xg)

z =

  1×13 logical array

   0   1   0   0   1   1   0   0   0   0   0   0   0