multiple functions in code to define different pathways

1 view (last 30 days)
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

Answers (1)

Ben Frankel
Ben Frankel on 26 Jul 2018
Edited: Ben Frankel on 27 Jul 2018
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
  17 Comments
Nkenna Amadi
Nkenna Amadi on 27 Jul 2018
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
Ben Frankel on 27 Jul 2018
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.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!