elimination of conscutive regions (generalization: ones with zeros between)

I need effectively eliminate (by zeroing) the consecutive "1's" between "-1's" and start/end of column at each column of matrix A, which now can be separated by any number of zeroes. The number of consecutive "1's" between "-1's" and start/end of column is > N. This is a non-trivial generalization of my previous Question.

Again, typical size(A) = [100000,1000].

See example:

A =
     1    -1     0
     0     1     1
     0     1     1
     1     1     0
     0     0     1
     1    -1     0
    -1     1     1
    -1     0    -1
     1     1     1
     0     1    -1

For N = 2 the expected result is

Aclean =
     0    -1     0
     0     0     0
     0     0     0
     0     0     0
     0     0     0
     0    -1     0
    -1     0     0
    -1     0    -1
     1     0     1
     0     0    -1

For N = 3 the expected result is

Aclean =
    1    -1     0
    0     1     0
    0     1     0
    1     1     0
    0     0     0
    1    -1     0
   -1     1     0
   -1     0    -1
    1     1     1
    0     1    -1

[Image Analyst]

I don't understand the question, or examples. The "number of consecutive "1's" between "-1's" is what it is. It's determined by A in advance and presumably A was created without regard to any N.

And for the outputs, you are changing the bottom of the middle column even though those 1's are not between ANY -1's. What's the correct rule?

[Michal]

I need to eliminate (zeroed) all "1's" which are located between "-1's" and start/end at each column of A separately. The number of eliminated "1's" between "-1's" is greater than user defined N.

So in my example, in 1st column of A are three "1's" and any number of "0's", separated by two "-1's" and then one "1's". For N = 2 I will eliminate only first 3 > (N=2) "1's". For N = 3 I do not eliminate any "1's" from 1st columns of A. Similarly for all columns.

[Matt J]

which now can be separated by any number of zeroes

This is the definition of non-consecutive. Why, then, are you saying the eliminated 1's are to be consecutive?

[Michal]

Well OK, any sequence of "1's" and "0's" separated by start/end and "-1's".

Answer 1

[Bruno Luong]

A = [1    -1     0;
     0     1     1;
     0     1     1;
     1     1     0;
     0     0     1;
     1    -1     0;
    -1     1     1;
    -1     0    -1;
     1     1     1;
     0     1    -1];

N = 3;
[m,n] = size(A);
Aclean = A;
for j=1:n
    Aj = [-1; A(:,j); -1];
    i = find(Aj == -1);
    c = histc(find(Aj==1),i);
    b = c <= N;
    im = i(b);
    ip = i([false; b(1:end-1)]);
    a = accumarray(im,1,[m+2,1])-accumarray(ip,1,[m+2,1]);
    mask = cumsum(a);
    mask(i) = 1;
    Aclean(:,j) = Aclean(:,j).*mask(2:end-1);
end
Aclean

[Bruno Luong]

Cleaner version

[m,n] = size(A);
Aclean = A;
for j=1:n
    Aj = A(:,j);
    im1 = find(Aj == -1);
    i1 = find(Aj==1);
    [c,loc] = histc(i1,[-Inf; im1; Inf]);
    b = c > N;
    Aclean(i1(b(loc)),j) = 0;
end

[Michal]

Is there any serious reason to still use old "histc" instead "histcounts"? Both functions has same performance and "histc" will be removed in future releases.

[Bruno Luong]

No, I'm just used to HISTC.

Answer 2

[Michal]

I am still looking for better (faster) solution. Any idea how to improve so far best solution:

A = [1    -1     0;
     0     1     1;
     0     1     1;
     1     1     0;
     0     0     1;
     1    -1     0;
    -1     1     1;
    -1     0    -1;
     1     1     1;
     0     1    -1];
N = 3;
sep = A==-1;
sep(1,:) = true;
idx = cumsum(sep(:));
sep(1,:) = A(1,:)==-1;
num = accumarray(idx, A(:)==1);
iff = num <= N;
Aclean  = reshape(sep(:)|iff(idx), size(A)) .* A;
Aclean

Big test matrix:

A = double(rand(100000,1000)>.4) - double(rand(100000,1000)>.65);

Bruno's code (N = 5): Elapsed time is 4.848747 seconds.

My code (N = 5): Elapsed time is 3.257089 seconds.

[jonas]

Looks quite efficient to me. The reshape is quite costly (about 20-25% of total time on my machine), so if you can remove that somehow you may increase the performance further.

[Bruno Luong]

@jonas: your statement is misleading: RESHAPE by itself costs almost nothing. The time is the operation | and .*

But yeah the code is quite good. There is only one non-optimize step, when step to be restore

sep(1,:) = A(1,:)==-1;

Under this statement MATLAB makes one more copying of big array. If user can afford extra memory, he could do this

sep = A==-1;
idx = sep;
idx(1,:) = true;
idx = cumsum(idx(:));
iff = accumarray(idx, A(:)==1) <= N;
Aclean = reshape(sep(:)|iff(idx), size(A)) .* A;

[Michal]

Could be possible to create mask term in my code

mask = reshape(sep(:)|iff(idx), size(A))

or yours mask with zeros (ones in your code) only on positions where are actually "1's" at matrix A??? I am using mask term for additional processing. Any idea how to do it effectively?

Of course it is a simple task but it add some extra overhead and I need fast as possible solution, because this function is called many times for different matrices A.

One possible solution in my code is:

mask = ~xor(mask, A)|mask

but, may be there exist any faster method?

[jonas]

@Bruno:

It is not misleading. The reshape (implying the entire operation minus the element-wise multiplication) is about 15-20% (perhaps 25% was a bit exaggerated). You are right in that the | operation is the costly part of this operation. However, you can tone down the rhetoric, misleading sounds almost as if it was intentional. I did split the .*A from reshape prior to running the code...

[Bruno Luong]

It is misleading, at least to me, when I read your comment, as if you said the time is causes by RESHAPE. But clearly you wants to tell the "statement with reshape is quite costly".

I don't think I use the word "misleading" in the sense other than what it means, no more no less.

[Bruno Luong]

@Michal I don't understand your last comment. What is "mask term"?

[jonas]

What I implied was that he could save plenty of time by maintaining the 2D structure, if that is at all possible. Thank you for your valuable input though.

[Michal]

@Bruno The mask term is at last computing command line in my code.

   mask.*A

where mask is RESHAPE() function discussed above.

[Bruno Luong]

Is this is what you want?

sep = A==-1;
idx = sep;
idx(1,:) = true;
idx = cumsum(idx(:));
is1 = A==1;
iff = accumarray(idx, is1(:)) <= N;
mask = reshape(sep(:)|iff(idx), size(A));
Aclean = mask .* A;
mask = mask & is1; % mask of 1s that are not turned off

[Michal]

No, the only way I so far found is

mask = ~xor(mask, A)|mask

as I mentioned in my previous comment.

[Bruno Luong]

It's not the same thing than

~xor(Aclean, A)

[Michal]

Yes … nice simplification!

[Bruno Luong]

And what about

mask = Aclean==A

?

Answer 3

[Bruno Luong]

That's true, somehow it's a 1D scanning problem.

I think we are close to the limit of MATLAB can do, if faster speed is still needed, then one should go to MEX programming route instead of torturing MATLAB to squeeze out the last once of speed.

[Bruno Luong]

Sorry, it should be in the comment section