Loops not working correctly

I have three different MATLAB files: randomStep.m, nearInfectious.m and singleSimulation.m all of which contain functions. The randomStep.m and nearInfectious.m have functions called randomStep and nearInfectious which are being called in the singleSimulation.m file.

The randomStep function takes in x and y coordinates (as vectors) of people, and makes them move to different parts of the screen, not exceeding xMax and yMax or going into negative values. The nearInfectious function has a population size (which will depend on the number of x and y coordinate inputs as vecotors) with 1 person starting out infected and the rest susceptible. I am storing the inital states (i.e. susceptible "s", infectious "i" or recovered "r") in an array called stateTable, which is 1 x populationSize. If anyone is within a radius of the infected person, they become infected and stateTable updates. This returns a logical array (1 x population size) where 1 indicates whether the people (person 1 to person populationSize) are within a radius of an infected person. 1 means they are within a radius of any infected person and 0 means they are not. Anyone who is infectious stays infectious.

In singleSimulation, I call the randomStep and nearInfectious functions with some input parameters. I want to return the time series for the susceptible, infectious and recovered people (the number of susceptible, infectious and recovered people each unit of time i.e. day 1, day 2). The recover probability is 0.2, and obviously that only occurs if someone is infected. So far, my singleSimulation.m file works fine except it isn't incrementing the rSeries values properly, because sometimes they decrease with each time step which cannot happen. How can I make it so that the "r" values don't decrease with each time step? Any help is greatly appreciated. Thank you.

Here is my randomStep.m code:

function [x,y] = randomStep(x,y,xMax,yMax,speed)        
    
    N=normrnd(0,1,1,length(x));
    for i=1:length(x)
        if(x(i)+N(i)*speed)>xMax
            x(i)=xMax-mod(x(i),xMax);
            
        elseif(x(i)+N(i)*speed)<0
            x(i)=abs(x(i)+N(i)*speed);
            
        else
            x(i)=x(i)+N(i)*speed;
           
        end       
    end
    
    N=normrnd(0,1,1,length(y));
    for j=1:length(y)
        if(y(j)+N(j)*speed)>yMax
            y(j)=yMax-mod(y(j),yMax);
            
        elseif(y(j)+N(j)*speed)<0
            y(j)=abs(y(j)+N(j)*speed);
            
        else
            y(j)=y(j)+N(j)*speed;           
        end        
    end
    disp(x)
    disp(y)
end

Here is my nearInfectious.m code:

function near = nearInfectious(x,y,states,radius) 
a=[];
b=[];

    for j=1:length(x)
        if states(j) == "i"
            a=[a x(j)];
            b=[b y(j)];         
        end   
    end 
disp(a);
disp(b);

near=[];


   for k=1:length(x)
       near_infected=false;
       x_coord=x(k);
       y_coord=y(k);
       
       for p=1:length(a)
          x_infected=a(p);
          y_infected=b(p);
          
          dist = sqrt( (x_coord-x_infected)^2 + (y_coord-y_infected)^2 );
          if dist <= radius
              near_infected=true;             
          end
       end
       near=[near,near_infected];     
   end
disp(near);
end

Here is my singleSimulation.m code:

[sSeries, iSeries, rSeries]=singleSimulation2(10,400,400,5,20,20,0.2)

function [sSeries, iSeries, rSeries]=singleSimulation2(populationSize, xMax, yMax, timeSteps, speed, radius, recoveryProb)

stateTable=["s","i","s","s","s","s","s","s","s","s"]
sSeries=[9,0,0,0,0];
iSeries=[1,0,0,0,0];
rSeries=[0,0,0,0,0];


for d=2:length(sSeries)

[x,y]=randomStep([103.2, 203.1, 300, 303.3, 244, 92, 37, 222, 299, 111],[45.2, 0.1, 20, 99, 108, 268, 399, 254, 77, 53],400,400,20);
%near=nearInfectious([103.2, 203.1, 300,303.3,244,92,37,222,299,111],[45.2, 0.1, 20,99,108,268,399,254,77,53],["s","i","s","s","s","s","s","s","s","s"],20);
near=nearInfectious(x,y,stateTable,150);


 for i=2:length(near)
     if(near(i)==true)
        
        stateTable(i)="i";
        
        if rand() <= recoveryProb && stateTable(i) =="i"
            stateTable(i)="r";            
        end
        
     end
     
 end
 
         
  disp(stateTable)    
  sSeries(d)=sum(stateTable=="s");
  iSeries(d)=sum(stateTable=="i");
  rSeries(d)=sum(stateTable=="r"); 

end
end

Answer 1

[Walter Roberson]

Your r count can decrease, because your code does not prevent r people from being re-infected.

[Khadeja Khan]

I tried to get r to not decrase by setting near to 0 whenever stateTable =="r" but that did not work. Any suggestions?

[Walter Roberson]

What is your revised code?

Also, why do you only permit someone to recover if they are near someone infectious ?

[Khadeja Khan]

If they are near someone infectious, that means they become infectious and they can recover with a certian probability.

[Walter Roberson]

If they are already infectious, then why is it that they need to be near someone else infectious in order to recover?

Your code

        if rand() <= recoveryProb && stateTable(i) =="i"
            stateTable(i)="r";            
        end

is inside the near(i)==true test, when it should not be

[Khadeja Khan]

Yup so I've moved it outside the near(i)==true test but still inside the second for loop and I am still incurring the same issue:

[sSeries, iSeries, rSeries]=singleSimulation2(10,400,400,5,20,20,0.2)

function [sSeries, iSeries, rSeries]=singleSimulation2(populationSize, xMax, yMax, timeSteps, speed, radius, recoveryProb)

stateTable=["s","i","s","s","s","s","s","s","s","s"]

sSeries=[9,0,0,0,0];

iSeries=[1,0,0,0,0];

rSeries=[0,0,0,0,0];

for d=2:length(sSeries)

[x,y]=randomStep([103.2, 203.1, 300, 303.3, 244, 92, 37, 222, 299, 111],[45.2, 0.1, 20, 99, 108, 268, 399, 254, 77, 53],400,400,20);

%near=nearInfectious([103.2, 203.1, 300,303.3,244,92,37,222,299,111],[45.2, 0.1, 20,99,108,268,399,254,77,53],["s","i","s","s","s","s","s","s","s","s"],20);

near=nearInfectious(x,y,stateTable,150);

for i=2:length(near)

if(near(i)==true)

stateTable(i)="i";

end

if rand() <= recoveryProb && stateTable(i) =="i"

stateTable(i)="r";

end

end

disp(stateTable)

sSeries(d)=sum(stateTable=="s");

iSeries(d)=sum(stateTable=="i");

rSeries(d)=sum(stateTable=="r");

end

end