Using boolean matrix to govern While loop

Hi there,

I'm working on a radiation transport MC in which I adjust two variables in incremental steps to produce a table of outputs. The following code does what I want it to:

for a = 1:5 % variable 1
  for b = 1:5 % variable 2
    for n = 1:10^5 % monte carlo

        while x(a,b) < 10 % the tricky part!
          x(a,b) = x(a,b) + rand * (a*2) * (b*.1);
        end

      end
    end
  end

handwave everything's initialized and the physics actually make sense.

So as I said, the above nested For loops work...but now I want to do all 25 combinations of variables at the same time using matrices. I hope to cut down on the number of random variables I generate, and also I just think it's a cooler way to do business. The problem is the "while" loop; each neutron "n" may scatter a different number of times before reaching arbitrary distance "10".

So, now I have something like:

a = meshgrid(1:5)*2;
b = meshgrid(6:10)*.1;

for n = 1:10^5
  s = ones(5); % keep track of "live" combinations of variables
  x = zeros(5); % distance traveled in each simulation

    while sum(s(:)) ~= 0
      x = x + s.* rand .* a .* b; % only if associated "s" is still 1
      s = s .* (x<10); % if any simulations exceeded 10, kill them
    end
  end

Again, pretend that the physics make sense and that I'm storing useful data from the result of each "n". This code also works, but it's much slower than the nested for loops. I think that's because each time there are a bunch of "multiply by zero" calculations that are slowing me down.

So finally, my question! Does anyone know a way to not calculate the zeros? I sort of want to do like a "forEach s(:,:), while s(:,:)==true, calculate all x(:,:)'s together". Can anyone help me translate that into Matlab-speak?

Thanks very much,

Bryan

Answer 1

[Walter Roberson]

numA = 5;
numB = 5;

[a, b] = ndgrid(1:numA, 1:numB);
a = a(:) .* 2;
b = b(:) .* 0.1;
ab = a .* b;

numAB = length(ab);

for n = 1 : 10^5

    s = true(numAB, 1);
    x = zeros(numAB, 1);

    while any(s)
      x(s) = x(s) + rand(nnz(s),1) .* ab(s);
      s(s) = x(s) < 10;
    end

end

x = reshape(x, numA, numB);

[Bryan]

That's it! Thank you very much, sir!

A few breadcrumbs for future travelers:

• use "any(s)", not "any(x)". Just a typo.
• "rand(nnz(s),1)" not required unless a different number is actually desired for each simulation. Just using "rand" works.
• Matlab still understands x(s) if they're equal-sized matrices. He may have "vectorized" it for speed?

Works great though. Thanks again!

Bryan

[Walter Roberson]

Opps, fixed the any(s) typo.

Your original code used a different rand for each iteration, so that's what I used.

x(s) is using logical indexing; having them be the same size is typical.