Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
"not a number" message

Subject: "not a number" message

From: Nick Stucky-Mack

Date: 3 Jan, 2009 02:59:02

Message: 1 of 4

Below is an excerpt of my code. I keep getting the error message "??? Attempted to access relAbundCol(NaN,:); index must be a positive integer or logical." The weird thing is that this only shows up after it has run successfully multiple times (the code below is part of a larger loop). I can't seem to figure out what is causing it to come up as not a number. Help!
Relevant excerpt of the code is below.


b = 0.25; %(*birth rate*)
s = 10^-2; %(*decrease in death rate for more fit cells*)
uA = 10^-3; %(*mutation rate for advantageous mutants*)
uN = 10^-2; %(*mutation rate for neutral mutants*)
topTime = 600; %(*number of time steps for each run*)
posMut = 30; %(*length of our mutation vector*)
NumRuns = 10000; %(*number of runs*)

for q = 1:NumRuns,
 abundance = zeros(posMut);
 abundance(1,1) = 1;
    time = 0;

relAbundCol = cumsum(abundance);

for i=2:posMut,
   relAbundCol(posMut,i) = relAbundCol(posMut,i) + relAbundCol(posMut,i-1);
end

relAbundOne = relAbundCol(posMut,:)/sum(sum(abundance));

Subject: "not a number" message

From: Roger Stafford

Date: 3 Jan, 2009 04:11:02

Message: 2 of 4

"Nick Stucky-Mack" <njstucky@fas.harvard.edu> wrote in message <gjmk9m$5ke$1@fred.mathworks.com>...
> Below is an excerpt of my code. I keep getting the error message "??? Attempted to access relAbundCol(NaN,:); index must be a positive integer or logical." The weird thing is that this only shows up after it has run successfully multiple times (the code below is part of a larger loop). I can't seem to figure out what is causing it to come up as not a number. Help!
> Relevant excerpt of the code is below.
>
> b = 0.25; %(*birth rate*)
> s = 10^-2; %(*decrease in death rate for more fit cells*)
> uA = 10^-3; %(*mutation rate for advantageous mutants*)
> uN = 10^-2; %(*mutation rate for neutral mutants*)
> topTime = 600; %(*number of time steps for each run*)
> posMut = 30; %(*length of our mutation vector*)
> NumRuns = 10000; %(*number of runs*)
>
> for q = 1:NumRuns,
> abundance = zeros(posMut);
> abundance(1,1) = 1;
> time = 0;
>
> relAbundCol = cumsum(abundance);
>
> for i=2:posMut,
> relAbundCol(posMut,i) = relAbundCol(posMut,i) + relAbundCol(posMut,i-1);
> end
>
> relAbundOne = relAbundCol(posMut,:)/sum(sum(abundance));

  I would advise you to show more of the significant parts of your code if you want to receive useful answers rather than wild guesses. There is very little that can be deduced from what you show here.

Roger Stafford

Subject: "not a number" message

From: Nick Stucky-Mack

Date: 3 Jan, 2009 04:36:01

Message: 3 of 4

Alright, thanks Roger! I've included the entirety of the code here.


%(*this is the code for the branching tree algorithm for cancer mutations*)

%(*this is where we define our parameters*)

b = 0.25; %(*birth rate*)
s = 10^-2; %(*decrease in death rate for more fit cells*)
uA = 10^-3; %(*mutation rate for advantageous mutants*)
uN = 10^-2; %(*mutation rate for neutral mutants*)
topTime = 600; %(*number of time steps for each run*)
posMut = 30; %(*length of our mutation vector*)
NumRuns = 10000; %(*number of runs*)



deaths = zeros(posMut);
bplusdeaths = zeros(posMut);
for j=1:posMut,
   deaths(1,j) = b*(1 - s);
end

for i=2:posMut,
   for j=1:posMut,
      deaths(i,j)=deaths(i - 1,1)*(1 - s);
end
end

for j=1:posMut,
   for i=1:posMut,
      bplusdeaths(i,j) = b + deaths(i,j);
end
end




totalAbundance = zeros(posMut);
averageAbundance = zeros(posMut);
goodRunsCounter = 0;

%(*here is where we loop through each of the runs*)

for q = 1:NumRuns,
 abundance = zeros(posMut);
 abundance(1,1) = 1;
    time = 0;
 



%(*find gamma*)
while(time<topTime)
   gamma = 0;
   for i=1:posMut,
      gamma = gamma + dot(abundance(i,:),bplusdeaths(i,:));
   end
if(gamma == 0)
break;
end


%(*Exponential distribution*)
deltaT = -log(rand)/(gamma);



handOfGodCol = rand(1);
  
relAbundCol = cumsum(abundance);

for i=2:posMut,
   relAbundCol(posMut,i) = relAbundCol(posMut,i) + relAbundCol(posMut,i-1);
end

relAbundOne = relAbundCol(posMut,:)/sum(sum(abundance));
  

j=1;
chosenMutantN = 1;


%(*this chooses which neutral mutant we are choosing to have something happen to*)

while(relAbundOne(j)<=handOfGodCol)
chosenMutantN = j + 1;
j = j + 1;
end
  


handOfGodRow = rand(1);


relAbundRow = zeros(posMut);
relAbundRow(1,chosenMutantN) = abundance(1,chosenMutantN);
for i=2:posMut,
   relAbundRow(i,chosenMutantN) = abundance(i,chosenMutantN) + relAbundRow(i-1,chosenMutantN);
end

relAbundTwo = relAbundRow(:,chosenMutantN)/sum(abundance(:,chosenMutantN));


j=1;
chosenMutantA = 1;


%(*this chooses which advantageous mutant we are choosing to have something happen*)

while(relAbundTwo(j)<=handOfGodRow)
chosenMutantA = j + 1;
j = j + 1;
end



handOfGodAgain = rand(1);
  
%(*this chooses what happens to that chosen mutant*)

total = abundance(chosenMutantA,chosenMutantN)*bplusdeaths(chosenMutantA,1);
  birth = abundance(chosenMutantA,chosenMutantN)*b*(1 - uA - uN)/total;
  mutationN = abundance(chosenMutantA,chosenMutantN)*b*(1 - uA)/total;
  mutationA = abundance(chosenMutantA,chosenMutantN)*b/total;
  
if(handOfGodAgain < birth)
abundance(chosenMutantA,chosenMutantN) = abundance(chosenMutantA,chosenMutantN) + 1;
elseif(handOfGodAgain < mutationN)
abundance(chosenMutantA,chosenMutantN + 1) = abundance(chosenMutantA,chosenMutantN + 1) + 1;
elseif(handOfGodAgain < mutationA)
              abundance(chosenMutantA + 1,chosenMutantN) = abundance(chosenMutantA + 1,chosenMutantN) + 1;
              else
abundance(chosenMutantA,chosenMutantN) = abundance(chosenMutantA,chosenMutantN) - 1;
end

  
  time = time + deltaT;
end



 %(*this is where we find the average outputs of the runs*)

if(gamma ~= 0)
  goodRunsCounter = goodRunsCounter + 1;
  disp('we just got a good run!');
  disp(q);
  for v = 1:posMut,
     for w=1:posMut,
   totalAbundance(v,w) = totalAbundance(v,w) + abundance(v,w);
     end
  end
end
end
disp('Total Abundance Matrix = ');
disp(totalAbundance);
disp('Number of Runs with surviving cells = ');
disp(goodRunsCounter);
for g=1:posMut,
   for h=1:posMut,
if(goodRunsCounter~=0)
averageAbundance(g,h)=totalAbundance(g,h)/goodRunsCounter;
end
   end
end
disp('Average Abundance Matrix = ');
disp(averageAbundance);


"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gjmogm$evm$1@fred.mathworks.com>...
> "Nick Stucky-Mack" <njstucky@fas.harvard.edu> wrote in message <gjmk9m$5ke$1@fred.mathworks.com>...
> > Below is an excerpt of my code. I keep getting the error message "??? Attempted to access relAbundCol(NaN,:); index must be a positive integer or logical." The weird thing is that this only shows up after it has run successfully multiple times (the code below is part of a larger loop). I can't seem to figure out what is causing it to come up as not a number. Help!
> > Relevant excerpt of the code is below.
> >
> > b = 0.25; %(*birth rate*)
> > s = 10^-2; %(*decrease in death rate for more fit cells*)
> > uA = 10^-3; %(*mutation rate for advantageous mutants*)
> > uN = 10^-2; %(*mutation rate for neutral mutants*)
> > topTime = 600; %(*number of time steps for each run*)
> > posMut = 30; %(*length of our mutation vector*)
> > NumRuns = 10000; %(*number of runs*)
> >
> > for q = 1:NumRuns,
> > abundance = zeros(posMut);
> > abundance(1,1) = 1;
> > time = 0;
> >
> > relAbundCol = cumsum(abundance);
> >
> > for i=2:posMut,
> > relAbundCol(posMut,i) = relAbundCol(posMut,i) + relAbundCol(posMut,i-1);
> > end
> >
> > relAbundOne = relAbundCol(posMut,:)/sum(sum(abundance));
>
> I would advise you to show more of the significant parts of your code if you want to receive useful answers rather than wild guesses. There is very little that can be deduced from what you show here.
>
> Roger Stafford
>

Subject: "not a number" message

From: Roger Stafford

Date: 3 Jan, 2009 07:09:07

Message: 4 of 4

"Nick Stucky-Mack" <njstucky@fas.harvard.edu> wrote in message <gjmpvh$6hd$1@fred.mathworks.com>...
> Alright, thanks Roger! I've included the entirety of the code here.
> .......

A few comments, Nick:

  I notice that in your first trip on each run through the

 while(time<topTime)

loop, every column except one in 'abundance' is initially all zeros, and yet you divide by the sum of one of these columns in

 relAbundTwo = ....
  relAbundRow(:,chosenMutantN)/sum(abundance(:,chosenMutantN));

If 'chosenMutantN' is anything but 1 at this point you could be in for trouble from the resultant infinity. One of the ways a NaN is produced is by dividing an infinity by another infinity or by subtracting two infinities. Another way is zero divided by zero.

  There is a division by 'total' which could be zero that first trip since it comes from an element of 'abundance'

 total = abundance(chosenMutantA,chosenMutantN)*...
          bplusdeaths(chosenMutantA,1);

and when you then calculate 'birth' you have a sure-fire NaN from the zero divided by zero situation there.

  It would be wise to make use of your debugger facilities to catch stuff like this yourself as it happens so you can analyze it for errors or unwise coding. After all, if you are going to be an effective programmer, it is absolutely essential to know how to debug your own code. No-one else can do as good a job as you yourself. I speak from long experience.

  You need much more robust code that is protected against certain kinds of unplanned consequences. For example in

 while(relAbundOne(j)<=handOfGodCol)
  chosenMutantN = j + 1;
  j = j + 1;
 end

you have no guard against j going beyond the end of the 'relAbundOne' array limit in case all its values are too small.

Roger Stafford

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us