multiple functions in code to define different pathways

Hi, i have the following code, at the moment it gives me results like the one shown in picture 1 as it searches along one path, however i would like it do as shown in picture 2(sorry about the poor drawing) by searching along different paths starting from different. I know how initiate different points such as indicated in line 12 of the code but i dont know how to allow the code take into consideration these different points i think the solution is to have different functions for "front". i have attached the matrix HC7.

Lx = 51; % Defines the size of the lattice . 
Ly = 51;
HCEx=1200000000;
pt =HC7;
IM = zeros(Ly,Lx); % Initialization of the lattice
% IM(i,j)=2
nstop=0;
t =1;
tend =1; % Number of "time steps".
k = 1; % Counting index for front propagation . 
q = 1; % Counting index for image writing .
% Initial step. The invader starts from a "point" or a channel 
%IM(Ly,50)=1; % Point IP
IM(Ly,10)=1; % Channel IP
% Finds the indices of the largest element in a given row in pt and invades (=1) the corresponding element in IM
%Defining the Trap spill point
while tend==1 && nstop==0  
% Check if there are trapped clusters
IM=bwlabel(IM-1,4); % Labels defending clusters , by setting invader sites = 0.
IM=IM+1; % Keeps the labeling , and puts the invader sites = 1.
for i = 1:Ly
s = sort(IM(i,1:end)); % Sorts the ith row in ascending order. 
if s(1)==1 % Tests if the ith row contains an invader site .
for j = 1:Lx
bin = front(i-1,j, IM,Lx,Ly); % If the ith row contains an invader site , all the elements in the ith row are being tested in the function front to see if they are part of the front .
if bin==1
% pfront is the front matrix where the first column
% contains the probability and the second and third
% columns contain the corresponding indices for this front site.
pfront(k,1) = pt(i-1,j); 
pfront(k,2) = i-1; 
pfront(k,3) = j;
k = k+1;
end
end
end
end
% Probability  y-coordinate xcoordinate.
% Finds the row index of the smallest probability in pfront and invades (=1) the corresponding element in IM
e = find(pfront(:,1)==min(pfront(:,1))); 
IM(pfront(e,2),pfront(e,3)) = 1;
clear pfront 
k = 1;
% Put invader sites = 1, and defending sites = 0.
IM(IM~= 1)=0;
% The following lines saves the image for each iteration .
imagesc(IM),colorbar,drawnow
% The following lines saves the image for each iteration .
end
%Computes the total volume

function [bin] = front(a,b,IM,Lx,Ly)
% The site is part of the front if one of the neighboring sites is an invader.
if IM(a,b)~=IM(1,1) % Takes trapped clusters into account If the site % tested is not part of the large defending cluster it should not be counted as a front site .
bin=0;
elseif (b+1)==(Lx+1)
if (IM(a+1,b)==1)||(IM(a-1,b)==1)||(IM(a,b-1)==1) 
    bin = 1;
else
bin = 0;
end
elseif (b-1)==(1-1)
if (IM(a+1,b)==1)||(IM(a-1,b)==1)||(IM(a,b+1)==1)
    bin = 1;
else
bin = 0;
end
elseif (a+1)==(Ly+1)
if (IM(a-1,b)==1)||(IM(a,b+1)==1)||(IM(a,b-1)==1)
    bin = 1;
else
bin = 0;
end
elseif (a-1)==(1-1)
if (IM(a+1,b)==1)||(IM(a,b+1)==1)||(IM(a,b-1)==1) 
    bin = 1;
else
bin = 0;
end
else
if (IM(a+1,b)==1)||(IM(a-1,b)==1)||(IM(a,b+1)==1)||(IM(a,b-1)==1) 
    bin = 1;
else
bin = 0;
end
end
end

Answer 1

[Ben Frankel]

There are some bounds-checking bugs in your code:

1. front(1, 1, IM, Lx, Ly) will crash.
2. front(1, Lx, IM, Lx, Ly) will crash.
3. front(Ly, 1, IM, Lx, Ly) will crash.
4. front(Ly, Lx, IM, Lx, Ly) will crash.

In other words, "front" crashes on the corners of the grid.

Here is an alternative implementation of "front" that bounds-checks properly:

function [bin] = front(a, b, IM, Lx, Ly)
bin = 0;
y_coords = [a+1 a-1 a   a  ];
x_coords = [b   b   b+1 b-1];
for ii = 1:4
    y = y_coords(ii);
    x = x_coords(ii);
    if y == 0 || y == Ly + 1 || x == 0 || x == Lx + 1
        continue
    end
    if IM(y, x) == 1
        bin = 1;
    end
end
end

Furthermore,

1. The line "bin = front(i-1,j, IM,Lx,Ly);" will crash whenever i == 1. You should make sure that i > 1 first.
2. You keep clearing pfront, so once the program is done (no cells left to invade) it will crash. You should use "pfront = [];" to empty the array instead and then check if the array is empty before using it.

I also have a few suggestions / pointers:

1. "e = find(pfront(:,1)==min(pfront(:,1)))" would be better written as "[~, e] = min(pfront(:, 1))" because the min function actually has two return values: the minimum value, and its index.
2. nstop and tend seem like they are never used in the code you posted (tend is always 1, and nstop is always 0).
3. Instead of "s = sort(IM(i,1:end)); if s(1)==1" you can write "if any(IM(i, :) == 1)".

Also, you have a bug in your algorithm. You keep adding cells that have already been invaded to pfront, preventing other cells from being invaded (because the already-invaded-cells have smaller probability so "[~, e] = min(pfront(:, 1))" will just select a cell that has already been invaded).

What you can do is iterate over the entire grid, and add each cell to your list of candidates (pfront) if it is A) part of the main defending cluster (check its IM label) and B) part of the front. B will be taken care of by the "front" function, so we just have to take care of A.

To see if a cell is part of the main defending cluster, we can verify that the defending cluster that it IS part of touches the edge of the grid. If we surround the grid in a border of defenders before calling bwlabel, it will consider all defending clusters that are adjacent to the border to be the same cluster, so all non-trapped clusters will have the same label. So if we do this:

IM = bwlabel(padarray(IM-1, [1 1], 1), 4) + 1;  % Add border, then use bwlabel
not_trapped = IM(1, 1);  % Sample label from border
IM = IM(2:end-1, 2:end-1);  % Remove border

Then to verify A we just check "IM(i, j) == not_trapped".

With all that out of the way, we can make some modifications to get the behavior you're looking for. If you want multiple sources to grow simultaneously, you need to distinguish your pfront values by the connected component of invaders it is a part of. You could use an additional bwlabel to get the connected components of invaders, and add another component to each pfront entry that represents the source a cell is on the front of.

However, this introduces an additional complication. What if a cell is on the front of two different connected components at the same time? If we iterate as before, the "front" function will only tell us that a cell is on some front, but it won't tell us which front (or fronts) it's on, now that there are multiple fronts. We can modify the "front" function to return a list of fronts that a given cell is a part of:

function [bin] = front(a, b, source, Lx, Ly)
bin = [];
y_coords = [a+1 a-1 a   a  ];
x_coords = [b   b   b+1 b-1];
for ii = 1:4
    y = y_coords(ii);
    x = x_coords(ii);
    if y == 0 || y == Ly + 1 || x == 0 || x == Lx + 1
        continue
    end
    if source(y, x)
        bin = [bin, source(y, x)];
    end
end
end

And that "source" argument isn't IM anymore. Instead, source should be the result of a second bwlabel on invaders: "source = bwlabel(IM == 1)". Now adding a cell to the fronts that it is a part of looks like this:

if IM(i, j) == not_trapped
    fronts = front(i, j, source, Lx, Ly);
    for n = 1:numel(fronts)
        pfront(k, 1) = fronts(n);
        pfront(k, 2) = pt(i, j);
        pfront(k, 3) = i;
        pfront(k, 4) = j;
        k = k+1;
    end
end

And finally, the invasion step now looks like this (update the front of each source separately):

fronts = unique(pfront(:, 1));
for n = 1:numel(fronts)
    candidates = pfront(pfront(:, 1) == fronts(n), 2:end);
    [~, e] = min(candidates(:, 1));
    IM(candidates(e, 2), candidates(e, 3)) = 1;
end

[Nkenna Amadi]

thank you very much for your comprehensive replies, ive been able to fix most of what you just mentioned, im just not sure how to implement your correction on the bounds checking, like where exactly does that go ? do i have to delete/replace any other subsequent codes with ? a bit more explanation on that would be really appreciated. With respect to the behaviour i expect, yes the invader should be able to spread anywhere adjacent to the existing invader, however my main point is i want to be able have multiple invaders(or source points if you would like) at specific locations i define and for the simulation to operate on all these invaders. Once again thank you very much for your detailed reply.

[Ben Frankel]

Thanks for the clarification, that makes sense. The bounds-checking snippet is for the function "front". It considers the four adjacent cells (X) to the cell you're looking at (O at position (a, b)):

   X 
  XOX
   X

It skips the cells that are out of bounds:

if y == 0 || y == Ly + 1 || x == 0 || x == Lx + 1
    continue
end

Then it checks if the X is an invader:

if IM(y, x) == 1
    bin = 1;
end

I'll update my answer based on your clarification.

EDIT: Actually, can you explain why you check IM(a, b) ~= IM(1, 1) in "front" and why you use bwlabel? If you're just spreading from invaders, you shouldn't have to worry about connected components.

[Nkenna Amadi]

For the bwlabel, the invader needs to be connected to the source point at the bottom which was defined by IM(Ly,1:Lx)=1; Essentially a common day example would be a shower head being the source. At the moment it just comes out from one hole, what i need is like again a case like shower head with multiple holes in which the invader(water) can come out from. so in this case if i define these points as IM(Ly,10)=1, IM(Ly,20)=1, IM(Ly,30)=1 etc i want the simulation to run considering all these points as outlets where subsequent invasions need to be connected to these defined points. The use of IM(a, b) ~= IM(1, 1) is there to take into account something called trapping. Essentially it trapped clusters into account If the site tested is not part of the large defending cluster it should not be counted as a front site.

   X 
  XOX if all adjecent to O have already been invaded O cannot be 
   X  invaded cus it is trapped

[Nkenna Amadi]

also could you show me the exact placement for the fronts bound check, is that at the beginning or the end of the function

[Ben Frankel]

A potential problem with comparing to IM(1, 1) is that it assumes that the upper left corner is always part of the "large defending cluster". If an invader spreads from the very bottom to the very top of the grid, it will cut the left side of the defending cluster off from the right side. In such a situation, how do you define which defending cluster is the large one?

Correct me if I'm wrong, but this is my take on what you're looking for: Every iteration you create a list of candidates for invasion, and a cell is a candidate for invasion if and only if it is adjacent to another invaded cell AND its connected component is not entirely surrounded by invaded cells.

If that's the case, your approach almost works, but not quite, for the reason outlined in the first paragraph. To fix that issue, you can add a border of defenders to your matrix with "IM = padarray(IM, [1 1])". I'll update my answer.

[Nkenna Amadi]

yes you are correct about the iteration process. looking forward to the update

[Nkenna Amadi]

hi, the answer you just provided i dont think it works. also could you show me the exact placement for the fronts bound check, is that at the beginning or the end of the function. i also dont understand what you mean iterate over the entire grid as i presumed i was already doing that. may be when you have access to matlab and can show me your updated code things would be easier. sorry for my poor understanding, im relatively new to matlab.

[Ben Frankel]

The code provided for front is the entire function :)

You aren't iterating over the entire grid, as you only visit a row i-1 if the row i has an invader in it. Your looping structure looks like this:

for i = 1:Ly 
    if any(IM(i, :) == 1)
        for j = 1:Lx
            % Visit the cell (i-1, j)

And I'm suggesting that it should look like this:

for i = 1:Ly
    for j = 1:Lx
        % Visit the cell (i, j)

[Nkenna Amadi]

i tried implementing your corrections like you stated, however it doesnt seem to work for some reason

[Ben Frankel]

What behavior do you see? I have access to MATLAB now, and it seems to work for me.

[Nkenna Amadi]

it just keeps crashing and giving error signs, can i see your final compiled code that you ran ?

[Ben Frankel]

Sure, but it's nothing I haven't already said (except the line "IM(IM ~= 1) = 2" which is just so defenders are all drawn in the same color).

while true
    IM = bwlabel(padarray(IM - 1, [1 1], 1), 4) + 1;
    not_trapped = IM(1, 1);
    IM = IM(2:end-1, 2:end-1);

      k = 1;
      pfront = [];
      for i = 1:Ly
          for j = 1:Lx
              if front(i, j, IM, Lx, Ly) && IM(i, j) == not_trapped
                  pfront(k, 1) = pt(i, j);
                  pfront(k, 2) = i;
                  pfront(k, 3) = j;
                  k = k+1;
              end
          end
      end

      if isempty(pfront)
          break;
      end

      [~, e] = min(pfront(:, 1));
      IM(pfront(e, 2), pfront(e, 3)) = 1;

      IM(IM ~= 1) = 2;
      imagesc(IM), colorbar, drawnow
  end

[Nkenna Amadi]

thank you so much for your patience and help. The code works fine now. this is what you have done i presume, and the result is as shown in picture 1. Youve fixed the boundary problem and front issues i had. However my expected behaviour that im really wanting is that if i add IM(Ly,20)=1; % Point IP just below where we have IM(Ly,10)=1; % Point IP the simulation still only considers IM(Ly,10)=1; % Point IP as shown in picture 2. What i really want is as shown in picture 3(i have manually updated the results for IM by adding ones, for illustration purposes) where it considers both IM(Ly,10)=1 and IM(Ly,20)=1; % Point IP. So in essence it is meant to grow from both source points. I probably didnt explain like this before thats probably why.

Lx = 51; % Defines the size of the lattice . 
Ly = 51;
HCEx=1200000000;
pt =HC7;
IM = zeros(Ly,Lx); % Initialization of the lattice
% IM(i,j)=2
nstop=0;
t =1;
tend =1; % Number of "time steps".
k = 1; % Counting index for front propagation . 
q = 1; % Counting index for image writing .
% Initial step. The invader starts from a "point" or a channel 
IM(Ly,10)=1; % Point IP

% Finds the indices of the largest element in a given row in pt and invades (=1) the corresponding element in IM
%Defining the Trap spill point
while true
    IM = bwlabel(padarray(IM - 1, [1 1], 1), 4) + 1;
    not_trapped = IM(1, 1);
    IM = IM(2:end-1, 2:end-1);
      k = 1;
      pfront = [];
      for i = 1:Ly
          for j = 1:Lx
              if front(i, j, IM, Lx, Ly) && IM(i, j) == not_trapped
                  pfront(k, 1) = pt(i, j);
                  pfront(k, 2) = i;
                  pfront(k, 3) = j;
                  k = k+1;
              end
          end
      end
      if isempty(pfront)
          break;
      end
      [~, e] = min(pfront(:, 1));
      IM(pfront(e, 2), pfront(e, 3)) = 1;
      IM(IM ~= 1) = 2;
      imagesc(IM), colorbar, drawnow
end

function [bin] = front(a, b, IM, Lx, Ly)
bin = 0;
y_coords = [a+1 a-1 a   a  ];
x_coords = [b   b   b+1 b-1];
for ii = 1:4
    y = y_coords(ii);
    x = x_coords(ii);
    if y == 0 || y == Ly + 1 || x == 0 || x == Lx + 1
        continue
    end
    if IM(y, x) == 1
        bin = 1;
    end
end
end

[Ben Frankel]

You want them to grow simultaneously, as in, two invasions per draw? The way it is now, only the one on the left grows because it always has the lower pt values.

If they are to both grow simultaneously, you need to have separate pfronts for each connected component of invaders. You could use an additional bwlabel to get the connected components of invaders, and add another dimension to pfront that represents the connected component a cell is on the front of.

However, this introduces an additional complication. What if a cell is on the front of two different connected components at the same time? If we iterate as before, the "front" function will only tell us that a cell is on a front, but it won't tell us WHICH front it's on (now that there are multiple fronts). We can modify the "front" function to return a list of fronts that a given cell is a part of.

Also, if two sources combine, should they be considered the same source thereafter? (I'm going to assume yes, and update my answer)

[Nkenna Amadi]

Yes if they combine they should be considered same source thereafter. The issue with regards to different components at the same time may not be an issue later on as i develop the code, as i am thinking of adding a rate of growth later on so for example source 2 would grow 2x faster than source one.

[Nkenna Amadi]

amazing ! thank you so much ! i had an abritary way of adding a rate at which each source goes but do you have a suggestion on the best way to introduce a rate at which it source grows and then if the source combine the new rate would be a function of the rate of the two sources so it would be faster. I know you've done alot for me already !

[Ben Frankel]

To do that you'd need some form of memory between iterations (how do you know that the source you're looking at is the result of two sources combining after the previous iteration?). I don't think there's an easy way to do it, but it's definitely possible.

Here's my suggestion:

1. Before your main loop begins, do a preliminary pass over the grid to find all initial sources. To remember them, create an array of "representative invader cells", i.e., sample one invader cell from each source. You might also be able to store each source's corresponding rate in the same array, along with some other metadata.
2. At the beginning of each iteration, after creating the "source" matrix, you need to figure out which sources have merged. You know that two sources have merged if their representative invader cells have the same "source" value.
3. Based on whichever sources have merged, update your representative invader cells and their rates.

As for manifesting those different rates, I would do the following:

1. Keep track of the current "tick" that starts at 0 and increments at the start of every iteration.
2. Based on each source's metadata (rate and some track of how many times it has invaded) you can determine how many more times it has to invade this tick.
3. If that number is 0, skip the source.
4. If that number is 1, invade.
5. If that number is 2 or higher... that's non-trivial to implement. Make sure your rates are slower than the tick rate so this never happens.