How can I apply a random number generator in a nested loop?

Hi all,

I'm new to matlab, but here is my dilemma: I want to select the cells with the highest values in a matrix (A) within a given range (determined by a radius). The selection of the highest value only happens if a given probability (which considers values t) is high enough. To do this, I am reading in a values (t), and applying this to a probability at every iteration ( for a max of 120 values corresponding to 120 iterations). My results, which is a plot of the cell selected every iteration, is a vertical line across the matrix (which does not make sense). I think it may be an issue with the for loop involving the random number generator. I would appreciate any suggestions concerning how to "fix" this, and I hope my description makes some kind of sense. Below is a shortened version of the code:

for i=1:maxtime  % iteration over time steps to the maximum 120 time steps

      xstart = xpos-range; % These are start positions
      ystart = ypos-range; 

      for t=1:120; %iteration over all values 1-120 
  current_t=temp_hourly(t,:); % this updates the value every iteration

prob_thermoreg= current_t-threshold_t; % determines the probability 

rn=rand; % generates random number 

if prob_thermoreg>rn;

%the logic below selects the cell with the highest value
maxtreecover = 0.0;  %find value larger than this
maxx_tree = xstart;
maxy_tree = ystart; 

for jj=1:20  %double the radius set at 10 cells
for ii=1:20  

          if A(ystart + ii,xstart + jj) > maxtreecover; % if the value of A at different ii and jj is greater than max 0.0, select this cell
              maxtreecover = A(ystart + ii,xstart + jj); 
              maxx_tree = xstart + jj;
              maxy_tree = ystart + ii;
          end
  end
  end
   end
      end

%The rest of the code updates the x and y coordinates

end

Answer 1

[Walter Roberson]

Your line

current_t=temp_hourly(t,:); % this updates the value every iteration

gives the impression that current_t will be a vector. If so then

prob_thermoreg= current_t-threshold_t; % determines the probability

would be a vector. Then when you did

if prob_thermoreg>rn;

with it being a vector that test would be the equivalent of

if all(prob_thermoreg>rn)

which would fail if any member of prob_thermoreg was less than the random number.

Could you confirm whether current_t will be a vector or a scalar?

[Stephanie Diaz]

Hi Walter,

Current_t will read in from "temp_hourly" as one value every iteration, so as a scalar. Temp_hourly is a 120X1 vector and so for iteration 1, value 1 is read in (essentially the first current_t value), until value 120. Current_t gets incorporated into the probability and the rn determines whether there will then be an action afterwards, for each iteration. I realize I do not have an "else" statement detailing what happens if the rn is greater than the probability, do you think this is also causing results to be strange?

[Walter Roberson]

" My results, which is a plot of the cell selected every iteration"

You do not show any code saving results each iteration, or plotting each iteration.

[Stephanie Diaz]

Hi Walter,

Below is the rest of the code that plots the results and updates variables. This section comes right after the last 5 consecutive "ends":

newmaxtreecover = maxtreecover;  % This is cell selected (with max value)
  newposx =  maxx_tree;  % This is the new x-coordinate of the cell
  newposy =  maxy_tree;  % This is the new y-coordinate of the cell

    xpossav(i) = xpos;  % save values of x-coordinate
    ypossav(i) = ypos;  % save values of y-coordinate
    timesav(i) = i;

    xpos = newposx % Update x-coordinate
    ypos = newposy % Update y-coordinate

    end

    figure %5
    plot(timesav,xpossav,timesav,ypossav), xlabel ('Time'), ylabel ('X, Y'), title ('')

    figure %5
    plot(xpossav,ypossav), xlabel ('X'), ylabel ('Y'), title ('')

      end

[Walter Roberson]

Suppose that for some reason all of the random tests failed. Then you would never updated maxx_tree or maxy_tree anywhere, so xpos and ypos would not change. But the code posted does not show what xpos or ypos or maxx_tree or maxy_tree are initialized to, which makes it difficult to predict what will be plotted.

If you had said "horizontal line" instead of "vertical line", then that is what I would expect from the first plot if the x and y never change; and in that same situation, the second plot would be plotting the constant x, y pair for each entry, which would leave nothing visible (if you added '-*' to the plot() then you would see a single marker.)

[Stephanie Diaz]

Hi Walter,

I'm not sure if this will be of much help, but below is the entire relevant code. Xpos and ypos are both initialized at 100. I think I understand what you're saying, but the probabilities are almost always greater than rn, which means maxy_tree and maxx_tree should update, correct?

       clear all
          close all

          A=dlmread('treecover_layer2.txt'); %tree cover map 
          B=dlmread('water_clip.txt'); %euclidean distance to water
          C=dlmread('evi_wet_l5.txt');%EVI wet season
          temp=dlmread('temp.txt'); %hourly temperature for 5 days in wet season
          temp_hourly=temp(:,2);%reads in only the column for temperature in the file (having 120 rows)
          threshold_t=30;
          xpos = 100;  % x-coordinate of initial position 
          ypos = 100;   % y-coordinate of initial position 
          maxtime = 120;  % number of time steps (24 hours times 5 days)
          range = 10;  % visual range 

          for i=1:maxtime  % iteration over time steps to the maximum 120 time steps
          % starting from xstart = xpos-range and ystart = ypos-range
              xstart = xpos-range; 
              ystart = ypos-range; 

              for t=1:120; %iteration over all t values 1-120 (5 days)
          current_t=temp_hourly(t,:); % updates current_t every iteration

          prob_thermoreg=1/(1+exp(-.05*(current_t-threshold_t))); % determmines probability given current_t

          rn=rand; %random number generator

          if prob_thermoreg>rn;
                    % The logic below will find the cell with the max value from the starting positions, xpos, ypos.
          maxtreecover = 0.0;  %find max larger than this
          maxx_tree = xstart;
          maxy_tree = ystart; 

          for jj=1:20  %double the radius 
          for ii=1:20  

                  if A(ystart + ii,xstart + jj) > maxtreecover; % if the value of A at different ii and jj is greater than max tree cover 
                      maxtreecover = A(ystart + ii,xstart + jj); %selection of the cell with max value
                      maxx_tree = xstart + jj;
                      maxy_tree = ystart + ii;
                  end
          end
          end
           end
              end

          newmaxtreecover = maxtreecover;  % This is cell selected
          newposx =  maxx_tree;  % This is the new x-coordinate 
          newposy =  maxy_tree;  % This is the new y-coordinate 

          xpossav(i) = xpos;  % save values of x-coordinate
          ypossav(i) = ypos;  % save values of y-coordinate
          timesav(i) = i;

          xpos = newposx % Update x-coordinate
          ypos = newposy % Update y-coordinate

          end

          image(A,'CDataMapping','scaled')

          figure %5
          plot(timesav,xpossav,timesav,ypossav), xlabel ('Time'), ylabel ('X, Y'), title ('')

          figure %5

      plot(xpossav,ypossav), xlabel ('X'), ylabel ('Y'), title ('')

[Walter Roberson]

Your code does not assign to maxx_tree or maxy_tree in the case that none of the random tests are satisfied. If your code does not crash on the line

            newposx =  maxx_tree;  % This is the new x-coordinate

then that could only be because one of the random tests was satisfied, contradicting your earlier statements.

Can you attach your code and your various .txt so we can test your code?

[Stephanie Diaz]

Here are the 2 text files necessary for the code and the actual code. I dont think I mentioned this, but eventually I will need to add various "else" statements if the probability is less than rn, and an alternative action is taken. I was hoping that expanding the code piece by piece would be easiest.

Thanks again for the help!

[Walter Roberson]

We need water_clip.txt and evi_wet_l5.txt

[Stephanie Diaz]

I called in those two files because I will need them later on once I expand the code, but for now I only need the two I attached. Should I not call in those water_clip and evi_wet_l5 if I'm not using them in the current code? I've attached them in the case that I do.

[Walter Roberson]

I am away from my computer for a couple of hours but will test when I return.

[Walter Roberson]

I do not see a vertical line as such. I see that most of the change in x and y occurs in the first step, and the rest in the second step, and then no change in x or y. Look at diff(xpossav) and diff(ypossav)

If you put in a disp() or fprintf() after the rand() you will see that the random test is satisfied on most (but not all) t values for every i.

If you look at your code, in each step you are looking for the highest point in A within a certain distance of the current point, and you move there, staying where you are if you are already at the highest point in that range. Then in the next iteration you repeat searching the same distance from where you are currently. Consider now what happens when you reach the highest point in that distance: there is no way to move off of it. So you go to the local peak and you get stuck there.

And that's what you are seeing: you move from the starting point, and you quickly get stuck at a local maxima. (You might even have gotten stuck at the global maxima.)

[Stephanie Diaz]

I see exactly what you are saying, thank you! I will add some sort of if statement to prevent this, perhaps increasing the distance for selecting the maximum point.