Number of imbracted loops as a parameter

Hello to the community,

I'll have to use imbricated loops but their number will vary.

For instance:

if parameter = 3

For i1 = 1 : N
    For i2 = 1 : N
          For i3 = 1 : N
            FORMULA;
          end;
     end;

end;

This is the case where I have to use 3 imbricated loops. But it might happen that I'll have to use 4, 5, 6 or 7. I'd like to avoid the 'if' statement of course.

Any idea ?

Thanks a lot!

(sorry for my bad english) Alicia

[Image Analyst]

Perhaps your English is better than mine. I had no idea what "imbricated" meant until I Googled it. Apparently it means "arrange (scales, sepals, plates, etc.) so that they overlap like roof tiles." (I usually call for loops like that "nested".)

Answer 1

[Alicia]

Hello Guillaume,

Thanks for your very quick and precious answer !

I don't know if I can avoid loops: I need to calculate something like this

( Sum from i=0 to N of a(i)*P(i) )^m

That's why I thought about doing m loops from 0 to N to get all the combinations : if m = 2 and N = 1

= a(0)a(0)P(0)P(0) + a(0) a(1) P(0) P(1) + a(1) a(0) P(1) P(0) + a(1) a(1) P(1) P(1)

Do you think I could something else ?

About your solution: if m=7 then do I have to set :

    maxindices = [N N N N N N N]
    curindices = [0 0 0 0 0 0 0]

General form:

    maxindices = ones(1,m)*N;
     curindices = zeros(1,m);

And then the a(curindice)*p(curindice) for the formula part And finally,

curindices(curloop+1:end) = 0 % because The loop is in fact from 0 to N

Am I right ?

[Guillaume]

Right, first, continuing with my solution, yes, maxindices is as you've shown. You could declare it as:

   maxindices = ones(1, M) * N;
   %or assuming N == numel(a) == numel(p)
   maxindices = ones(1, M) * numel(a);

Matrix indices start at one in Matlab, a(0) is not valid, so curindices needs to be all ones:

curindices = ones(1, M);

As for the formula part, you need to calculate the product of a(curindices) and p(curindices):

curprod = prod(a(curindices)) * prod(p(curindices));

You can then add curprod to a running sum.

curloop needs to stay as it is. It tracks which index is to be incremented (which of the i(n) of your original example). As mentioned, you can't count from 0 to N, so you'll have to count from 1 to N+1, so just change maxindices to be N+1.

---------

Now, what you are calculating looks similar to a kronecker tensor product. Assuming that N == number of elements in a and p, you could thus write your code as:

aprod = 1;
pprod = 1;
for mindex = 1:M
  aprod = kron(aprod, a);
  pprod = kron(pprod, p);
end
result = sum(aprod .* pprod);

[Alicia]

Indeed! And it is correct, I tested the code with a simple example.

But coul you explain your code ? I read your link but I don't really understand why do you start with aprod and pprod with value 1 ?

Thanks a lot!

[Guillaume]

1 is just the identity value for kron. You could do also start with

aprod = a;
pprod = p;

but then you need to do the loop one less time, i.e.:

for mindex=2:M
   aprod = kron(aprod, a);
   ...

The first version has the advantage of covering the degenerate case where M = 0. If it's never going to happen, then you can use the 2nd version.

And if M is never less than two, you can even do:

aprod = kron(a, a);
pprod = kron(p, p);
for mindex = 3:M
  ...

To close off, it would have been nice if you'd credited me for the answer rather than yourselve.

Answer 2

[Guillaume]

First of all, consider if you can operate on matrices instead of scalars. That is, avoid the loops altogether. It all depends on what FORMULA does, but most matlab operations can be performed on matrices.

If you can't avoid loops, here's what you can do:

    maxindices = [N N N N N]; %upper bounds of all the in, as many as indices.
    curindices = [1 1 1 1 1]; %current loop indices, this is your i1, i2, i3, ...
    while true      
        FORMULA(curindices); %curindices(1) is i1, curindices(2) is i2, ...
        %now all there is to do is increment the correct loop.
        %this will be the last loop that has not reached its max.
        curloop = find(curindices < maxindices, 1, 'last');
        if isempty(curloop) 
            %all indices have reached their maximum. All done, break the while loop.
            break;
        else
            curindices(curloop) = curindices(curloop) + 1; %increment current loop.
            curindices(curloop+1:end) = 1; %reset all indices after current loop.
        end
    end