I'm looking to eliminate the following for loop:

I'm looking to eliminate the following for loop:

for j=1:96
 rn(ijk(j))=rn(ijk(j))+mp(j);
end

Background: rn is a 200x3x3 matrix, ijk is a 1800x1 matrix and np is a 12x8 matrix. ijk describes the mapping between mp and (the much larger) rn. For example, when j=1, mp(j)=0.84 and ijk(1)=601. This loop is embedded in a much larger loop.

There has got to be a way to eliminate the for loop presented here. Any suggestions?

Answer 1

[Sean de Wolski]

Perhaps:

idx = 1:96;
rn(ijk(idx))=rn(ijk(idx))+mp(idx);

[John]

Thanks for the prompy reply, Sean. I tried this before, but it didn't give me the same matrix. Actually, I tried

rn(ijk(1:96))=rn(ijk(1:96))+mp(1:96);

And then I tried your suggestion.

They did not produce the same matrix rn as:

for j=1:96

rn(ijk(j))=rn(ijk(j))+mp(j);

end

Maybe I'm way off base, but I figured what you said (and what I tried ealier) should work?

[Alex]

I just played around with a similar setup and I'm getting the same answer using either method.

One way that it could be causing different answers is if the data types are different (ex. doubles vs cells).

You can further debug this to see what is going on by reducing the problem. Start with a smaller idx (even 1) to see if you are getting the same results both ways. If 1 works, make idx bigger until you find the problem.

[Sean de Wolski]

Could you try each method in its entirety in two different functions and then compare the results after? Perhaps you're performing the same operation twice or not preallocating the matrices the same way (i.e. different sized results)

Answer 2

[John]

Here was the problem: my vector ijk pointed to same index in rn at multiple times. For example, consider

mp=[0.1;0.2;0.3]; ijk=[1,2,1]; rn=zeros(4,1);

The vectorized function gives: idx = 1:3; rn(ijk(idx))=rn(ijk(idx))+mp(idx);

rn =

    0.3000
    0.2000
         0
         0

The loop expression gives:

for j=1:3 rn(ijk(j))=rn(ijk(j))+mp(j); end

rn =

    0.4000
    0.2000
         0
         0

I'm still not sure how to vectorize the code for this special case, but the answers provided by Sean and Alex are correct to the question.

[Sean de Wolski]

See me new Answer!

Answer 3

[Sean de Wolski]

%Data:
mp=[0.1;0.2;0.3];
ijk=[1,2,1];
rnsz = [4 1]; %the size you want rn to be:

%Engine:
rn = accumarray(ijk',mp',rnsz); % Build rn; note transpose of subs/vals

[John]

Winner! Awesome!

Wish I had known this command earlier.

Thanks again.