Vectorization fail with sub2ind?

I have the for loop code

 for m=1:numel(N,1)
    y  =  N(m,1);  % y matrix index
    x  =  N(m,2);  % x matrix index
    dy =  N(m,4);  % codfficient
    dx =  N(m,3);  % coefficient
    wts(y,x) = wts(y,x)+dx.*dy;
 end

I attempted to vectorize this as follows

 ind = sub2ind(size(wts),N(:,2),N(:,1));
 wts(ind) = wts(ind)+N(:,3).*N(:,4);

Unfortunately, for some reason the result is totally wrong.

This is a followup to a previous question : http://www.mathworks.com/matlabcentral/answers/119845-using-a-vector-of-indices-to-add-values-to-a-matrix

Someone later said that repeated values of "(x,y) / ind" will be erroneous. I don't know what that means, but there are many repeated values of the sub2ind index in my code (lots of addition).

How can this be fixed?

[per isakson]

The documentation, numel, Number of array elements, says that numel takes one argument. However,

    >> numel( magic(6),1 )
    ans =
         1

did you intend

    >> size( magic(6),1 )
    ans =
         6

What is the initial value of wts?

[Anand]

What do you mean wrong result? The vectorized code looks right to me.

[per isakson]

Why not show an executable example based on a for-loop?

[Christopher]

I have built the following script which I think is illustrative:

clear all
xygrd=10;                 % grid height, width
wts=zeros(xygrd,xygrd);   % preallocate
wts2=zeros(xygrd,xygrd);  % preallocate

% create points randomly allocated to the matrix xygrd*xygrd by randomly
% picking x,y indices
N1 = floor(rand(1000,2)*xygrd)+1; % x,y indices
N2 = rand(1000,2);                % point values (dx,dy)
N=[N1 N2];                        % x,y,dx,dy

for m=1:size(N,1)
    y  =  N(m,1);  % y matrix index
    x  =  N(m,2);  % x matrix index
    dy =  N(m,4);  % dy
    dx =  N(m,3);  % dx
    wts(y,x) = wts(y,x)+dx.*dy;
end

% vectorization
ind = sub2ind(size(wts2),N(:,2),N(:,1));
wts2(ind) = wts2(ind)+N(:,3).*N(:,4);

pcolor(wts-wts2); % should all be zeros

matrix wts (from the for loop) is completely different from wts2. I need a vectorized code for the for loop.

Thanks for thinking about this with me.

Answer 1

[Roger Stafford]

Here's a simple example to show you how repeated values of 'ind' will produce different results from your for-loop:

 wts = zeros(3,1);
 ind = [2;2;2];
 wts(ind) = wts(ind) + [3;7;19];
 wts
 ans = 
       0
      19
       0

In your vectorization you would have been expecting [0;29;0] as in a for-loop, but instead only the last addition is performed, the other two results having been overwritten.

To do what you want, use matlab's 'accumarray' function:

wts = wts + accumarray([N(:,1),N(:,2)],N(:,3).*N(:,4),size(wts));

It will perform all the additions, even with repetitions.

(Note: You appear to have N(:,2) and N(:,1) in reversed order in your call to 'sub2ind' which would produce further errors.)

[Christopher]

Thanks Roger, It appears you have found the correct method!

% diff vectorization
wts3 = wts3 + accumarray([N(:,1),N(:,2)],N(:,3).*N(:,4),size(wts));

is can be added to the script posted earlier, and it reproduces the for-loop.