Simple GA v15b

//**************************************************************
// File:GA_v15.sce
// Authors: Gerd Doeben-Henisch
// Version Start: January-18, 2010
// Version Last: May-23, 2011
//******************************************************************
// Function: Implements a simple GA according to Goldberg (1989)
//
//*****************************************************************
// Apri-l5, 2011 : Adding a mutation operator
//
//***************************************************
// April-8, 2011 : Adding an automatic generation of a population with random
values
//
//****************************************************
// April-21, 2011: Adding a frequency counter for binary strings
//
//**************************************************
//
// April-22-27, 2011: Adding an evaluation of the quality of sets of solutions
//
//**************************************************
//
// April-28, 2011: Computing the probability of certain patterns
//
//*************************************************
//
// April-29, 2011: GA without fitness
//
//*************************************************
//
// May-01, 2011: Virtual fitness function for comparison
//
//*************************************************
//
// May-05, 2011: Conditioned Mutation
//
//*************************************************
//
// May-12, 2011: Fitness ratios, convergence times
//
//*************************************************
//
// May-18, 2011: Purified version of ga()
// and : more converting functions vec22dec()
//
//*************************************************
//
// May-22, 2011: Improving usage of parameter p in case 
// of crossoverPrep with helpfunctions: the flag is fixed to 'l+7' and
// freed from p
// Corrections to crossoverPrep(); mismatch with while - conditions
//
//*************************************************
//
// May-23, 2011: Added a counter for the number of genomes in a POP which have
// reached 100\%
//
//*************************************************
// Helpful scilab-functions from the scilab libraries
//
// nancumsum ? This function returns the cumulative sum of the values of a
matrix
// tabul ? frequency of values of a matrix or vector
// variance ? variance of the values of a vector or matrix
// strange ? range
// stdevf ? standard deviation
// bin2dec := translate a binary string of '1','0' into a decimal number
//
// SW-STRUCTURE
//
// All data are organized in a dynamical table. Left hand the actual genes, in
the  middle
// supporting parameters and to the right intermediate modifications of the
genes to the left.
//
// POP FORMAT
// l := length of strings
// p := number of cells between string <1...l> and <l+p+1, ..., l+p+l>
// Pos 1-l := String
// Pos l+1 := Decimal Value of binary string
// Pos l+2 := Fitness for l+1 and l+p
// Pos l+3 := Percentage of overall fitness
// Pos l+4 := Expected count according to fitness
// Pos l+5 := Realized count
// Pos l+6 := 2nd Decimal value of a compund fitness
// Pos l+7 := crossover flag
// p := 7

// Variable for overall Fitness

FITNESS_ALL = 0

// Variable for average Fitness

AFITNESS = 0


//***************************************************
//Function to generate automatically a population with random values
// Input:
// l := length of strings
// n := size of population
// Output:
// A population according to the above defined format
//

function[POPX]=popgen(n,l)
  
  //Generate a matrix with only '0's
  
  POPX = zeros(n,l)
  MaxCells = n*l
  
  // Distribute randomly '1's
  
  for i=1:n
    for j=1:l
      if( rand() < 0.5 ) then POPX(i,j) = 0, else POPX(i,j) = 1, end
    end
    end
         
       endfunction
       
       
       

//***************************************************
//Function to generate random numbers between [1,up]
// 
// n := number of choices

function[Xs]=randInt(n, up)
  
  Xs = []
  
  for i=1:n
    x=0
    while(x==0)
      x = round(up*rand())
      end
    Xs(i) = x
      end
 
endfunction



//***************************************************
//Function to compute descriptive parameters of a population
// Input:
// l := length of strings
// n := size of population
// Output:
// Pc := Percentage of possible space (PSpace)
// Pb := Probability of gitting a certain individual of PSpace

function[Pc, Pb]=popDescr(n,l)
  
  // Cardinality of PSpace
  
  PSPC = 2^l
  
  // Percentage Pc
  
  Pc = n/(PSPC/100)
  
  // Probability Pb
  
  Pb = n/2^l
         
endfunction




//***************************************************
//Function to compute the frequencies of the elements of a population
//
// Input:
// l := length of strings
// n := size of population
// N := number of events
// show := if '1' then display the POPX-matrix
//
// Output:
// FX := Vector of freuencies ordered along the decimal values of a string (1st
column)
// FXDif := Vector of Differences compared to theoretical expectation in % (2nd
column)
// MEAN := The mean value of the differences
// STD := The standard deviation of the differences
//
// Idea: Generate a sequence of populations, count the occurences, and check the
degree of agreement between 
// empirical values and theoretica lexpectations

function[DIST2,STD, MEAN, FREQ,STD1, MEAN1, FREQ1,FX,
POPX]=frequCheck(l,n,N,show)
  
  
  // Install the counters
  
  r= (2^l)
  c= 3
  
  FX = zeros(r,c)
  
  // Generate an internal index for display
  
  for j=1:r
    FX(j,1)=j-1
    end
  
  // Generate the populations
  
  for i=1:N
  
  [POPX]=popgen(n,l)
  
  
  if show==1 then printf("POPX first\n"), disp(POPX), end
  
  //Compute decimal value ('0' is placed at index '1' !!!)
  // Count occurences 
  
  for j=1:n
    v=POPX(j,:),  d = vec2dec(v,l), [POPX(j,l+1)]= d, FX(d+1,2)=FX(d+1,2)+1 
  end//n 
  
  if show==1 then printf("POPX second\n"), disp(POPX), end
  
end //N

for j=1:r
  
    FX(j,3)=(FX(j,2)/n)/N //Normal frequency
    FX(j,4)=(FX(j,2)/n)/(((1/2^l)*N)/100) // Frequency as percentage
  end//m 
  
  
  disp(FX)
  
  FREQ1=tabul(FX(:,3))
  MEAN1=mean(FREQ1(:,1)) 
  STD1= stdevf(FREQ1(:,1), FREQ1(:,2))
  FREQ=tabul(FX(:,4))
  MEAN=mean(FREQ(:,1)) 
  STD= stdevf(FREQ(:,1), FREQ(:,2))
  MAX=max(FREQ(:,1))
  MIN=min(FREQ(:,1))
  DIST=MAX-MIN
  DIST2=DIST/2
  
  
  printf("Number of Events n * N = %d\n",n*N)
  
  // Show Plot
  
 // clf(), xdel,
 // plot2d(FREQ(:,1), FREQ(:,2))
         
endfunction


//***************************************************
//Function to compute the frequencies of the elements of a population
// and using a fitness function to compute 'virtual' fitness to compare
// the populations with others by fitness. But the function does not make use of
these values.
//
// Input:
// l := length of strings
// n := size of population
// N := number of events
// show := if '1' then display the POPX-matrix
//
// Output:
// FX := Vector of freuencies ordered along the decimal values of a string (1st
column)
// FXDif := Vector of Differences compared to theoretical expectation in % (2nd
column)
// MEAN := The mean value of the differences
// STD := The standard deviation of the differences
//
// Idea: Generate a sequence of populations, count the occurences, and check the
degree of agreement between 
// empirical values and theoretica lexpectations

function[FITNESS_ALL_PERC1,DIST2,STD, MEAN, FREQ,STD1, MEAN1, FREQ1,FX,
POPX]=frequCheck1(l,n,N,show)
  
  
  // Install the counters
  
  r= (2^l)
  c= 3
  
  FX = zeros(r,c)
  FITNESS_ALLLOG1=[]
  
  // Generate an internal index for display
  
  for j=1:r
    FX(j,1)=j-1
    end
  
  // Generate the populations
  
  for i=1:N
  
  [POPX]=popgen(n,l)
  
  
  if show==1 then printf("POPX first\n"), disp(POPX), end
  
  //Compute decimal value ('0' is placed at index '1' !!!)
  // Count occurences 
  
  for j=1:n
    v=POPX(j,:),  
    d = vec2dec(v,l), [POPX(j,l+1)]= d, //decimal values
    FX(d+1,2)=FX(d+1,2)+1, //occurences
    [POPX(j,l+2)]= fitness1(POPX(j,l+1)) //fitness values 
  end//n
  
  [FITNESS_ALL]=fitness(POPX,l)
  FITNESS_ALLLOG1(i)=FITNESS_ALL
  
  if show==1 then printf("POPX second\n"), disp(POPX), end
  
end //N

for j=1:r
  
    FX(j,3)=(FX(j,2)/n)/N //Normal frequency
    FX(j,4)=(FX(j,2)/n)/(((1/2^l)*N)/100) // Frequency as percentage
  end//m 
  
  
  disp(FX)
  
  FREQ1=tabul(FX(:,3))
  MEAN1=mean(FREQ1(:,1)) 
  STD1= stdevf(FREQ1(:,1), FREQ1(:,2))
  FREQ=tabul(FX(:,4))
  MEAN=mean(FREQ(:,1)) 
  STD= stdevf(FREQ(:,1), FREQ(:,2))
  MAX=max(FREQ(:,1))
  MIN=min(FREQ(:,1))
  DIST=MAX-MIN
  DIST2=DIST/2
  
  
  printf("Number of Events n * N = %d\n",n*N)
  
// Compute the maximal value and the percentage of success

[MAXFIT]=maxfit02(l,n)
FITNESS_ALL_PERC1=[]

for i=1:N, FITNESS_ALL_PERC1(i) = FITNESS_ALLLOG1(i)/(MAXFIT/100), end

// Show graphical results

  if show==2 then
     clf(), xdel, 
     plot2d([1:1:N], [FITNESS_ALL_PERC1]),
   end
         
endfunction

//***************************************************
//Function to repeat frequency  checks
//
// Input:
// l := length of strings
// n := size of population
// N := Number of events
// K := how often steps should be applied
//
// Output:
// FX := Vector of freuencies ordered along the decimal values of a string (1st
column)
// FREQ := Vector of frequencies and their occurences
// MEAN := The mean value of the differences
// STD := The standard deviation of the difference
// DIST2 := Half distances between MAX and MIN


  function[MEANX, DIST2X, MEAN1X, STD1X,DIST20,
MEAN0]=freqCheckRep(l,n,N,K,show)

  
  // Counter for MEANX, DIST2X...
  
  MEAN1X=[]
  STD1X =[]
  
  MEANX=[]
  DIST2X =[]
  
 for j=1:K
     
    [DIST2,STD, MEAN, FREQ,STD1, MEAN1, FREQ1,FX, POPX]=frequCheck(l,n,N,show) 
    
    
      MEAN1X(j,1) = MEAN1
      STD1X(j,1) = STD1
      MEANX(j,1) = MEAN
      DIST2X(j,1) = DIST2
          
        end //j
        
        
        MEAN10 = mean(MEAN1X(:,1))
        STD10 = mean(STD1X(:,1))
        MEAN0 = mean(MEANX(:,1))
        DIST20 = mean(DIST2(:,1))
        
        
       
  // Show Plot
  
 // clf(), xdel,
//  plot2d([DIST2SUM(:,1)], [DIST2SUM(:,4)])
         
endfunction


//***************************************************
//Function to repeat frequency  checks with gassimple0()
//
// Input:
// l := length of strings
// n := size of population
// N := Number of events
// K := how often steps should be applied
//
// Output:
// FX := Vector of freuencies ordered along the decimal values of a string (1st
column)
// FREQ := Vector of frequencies and their occurences
// MEAN := The mean value of the differences
// STD := The standard deviation of the difference
// DIST2 := Half distances between MAX and MIN


  function[MEANX, DIST2X, MEAN1X, STD1X,DIST20,
MEAN0]=freqCheckRep0(POP,l,p,n,N,K, MThreshold,show)

  
  // Counter for MEANX, DIST2X...
  
  MCount = 0
  MEAN1X=[]
  STD1X =[]
  
  MEANX=[]
  DIST2X =[]
  XDIM=[]
  
  for j=1:K
    
    XDIM(j,1)=j*N
    [MCount,DIST2,STD, MEAN, FREQ,STD1, MEAN1,
FREQ1,POP,FX]=gasimple0(POP,l,p,n,N, j,MThreshold,MCount,show)
    
    
     disp(MCount)
    
      MEAN1X(j,1) = MEAN1
      STD1X(j,1) = STD1
      MEANX(j,1) = MEAN
      DIST2X(j,1) = DIST2
          
        end //j
        
        
        MEAN10 = mean(MEAN1X(:,1))
        STD10 = mean(STD1X(:,1))
        MEAN0 = mean(MEANX(:,1))
        DIST20 = mean(DIST2(:,1))
        
        
       
       // Show Plot
       
  
 if show==2 then clf(), xdel, plot2d([XDIM(:,1)], [DIST2X(:,1)]), end
         
endfunction




//***************************************************
//Function to compute the changing quality of solution sets
//
// Input:
// l := length of strings
// n := size of population
// N := Number of events
// K := how often should be incremented
//
// Output:
// MEAN=:= Mean value of events
// DIST20 := Mean value of 1/2-Distance of events between MAX and MIN
// DIST2X := Series of DIST0s
// MEANX := Series of  MEANOs
// FX := Summary of all values


  function[MDX]=solutionSetApproach(l,n,N,K,show)

  
  // Counter
  
  MDX=zeros(K,6)
  
  for j=1:K
    
    disp(j)
     
    [DIST2,STD, MEAN, FREQ,STD1, MEAN1, FREQ1,FX, POPX]=frequCheck(l,n,j*N,show)
 
    
        
    MDX(j,1) = j*N
    MDX(j,2) = MEAN
    MDX(j,3) = DIST2
    MDX(j,4) = STD
    MDX(j,5) = MEAN1
    MDX(j,6) = STD1 
    
          
   end //j
        
  // Show Plot
  
  clf(), xdel,
  
  x=[MDX(:,1)], 
  y1=[MDX(:,3)], subplot(3,1,1)
  plot2d(x,y1), xtitle=("DIST2[%]")
  y2=[MDX(:,4)],subplot(3,1,2)
  plot2d(x,y2), xtitle=("STD[%]")
  y3=[MDX(:,6)], subplot(3,1,3)
  plot2d(x,y3), xtitle=("STD")
         
endfunction





//***************************************************
//Function to integrate multiple frequency simulations
//
// Input:
// l := length of strings
// n := size of population
// step := Inrement for number of events
// K := how often steps should be applied
//
// Output:
// DIST2SUM := Summary of distances DIST2 with regard to number of events 


  function[DIST2SUM]=distSum(l,n,step,K,show)

  
  // Counter for DIST2
  
  DIST2SUM=zeros(K,4)
  
 for j=1:K
      disp(j)
    
    [DIST2,STD, MEAN, FREQ,FX, POPX]=frequCheck(l,n,j*step,show) 
    
    
    DIST2SUM(j,1) = j*step
    DIST2SUM(j,2) = MEAN
    DIST2SUM(j,3) = DIST2
    DIST2SUM(j,4) = DIST2SUM(j,3)/(MEAN/100)
    
   end //j

  
  // Show Plot
  
  clf(), xdel,
  plot2d([DIST2SUM(:,1)], [DIST2SUM(:,4)])
         
endfunction

//***************************************************
//Function to compute the overall fitness of a matrix POP
// Sum up all the l+2-th positions
// l := length of fitness strings

function[FITNESS_ALL]=fitness(POP,l)
  
 FITNESS_ALL=0;
[r,c]=size(POP);
for j=1:r
  FITNESS_ALL=FITNESS_ALL+POP(j,l+2)
end

endfunction



//***************************************************
//Function to compute the maximal individual fitness of a matrix POP
// Compute all the l+2-th positions
// l := length of fitness strings

function[MFITNESS]=maxfitness(POP,l)
  
  // Extracts column l+2 and selects the maximal values
  
MFITNESS = max(POP([:],l+2))

endfunction


//***************************************************
//Function to compute the relative fitness of a string
// along with the average fitness 
// Compute all the l+3-th positions
// l := length of fitness strings

function[POP,AFITNESS]=rfitness(POP,l, FITNESS_ALL)

[r,c]=size(POP);
for j=1:r
  POP(j,l+3)=POP(j,l+2)/FITNESS_ALL
end

AFITNESS=FITNESS_ALL/r

endfunction


//***************************************************
//Function to compute the fitness ratio r=f1/f0
//
// l := length of fitness strings
// f1 := fitness of all individuals with a allele value '1' at a position
// f0 := fitness of all individuals with a allele value '0' at a position 
// r = f1/f0
//
// Assumption: POP has the relativ fitness available at l+3


function[POP,ratio,f1,f0]=fratio10(POP,l,show)

[r,c]=size(POP);
f1=0;
f0=0;

for i=1:r
  for j=1:l
    if (POP(i,j)== 1) then f1=f1+POP(i,l+3),
    else f0=f0+POP(i,l+3)
    end
    end //j
end//i

ratio=f1/f0

if show == 2 then disp(POP),  disp(f1), disp(f0), disp(ratio), end

endfunction

//*************************************************************
// Computing the worst case convergence time tcw
//
// p := number of parameters between strings left and right (actually  p=5)
// l := length of strings ( actually l=5)
// n := number of elements in POP (actually n=4)
// POP := a predefined population (can be done automatically)
// 
// f1 := fitness of all individuals with a allele value '1' at a position
// f0 := fitness of all individuals with a allele value '0' at a position 
// r = f1/f0
// tcw := worst case convergence time
//
// Assumption: POP has the relativ fitness available at l+3



function[POP,tcw,tcabin,ratio,f1,f0]=ftcw(POP,l,p,n,show)
  
  for i=1:n, v=POP(i,:),  
    d= vec2dec(v,l),  //decimal value
    [POP(i,l+1)]= d,  //decimal value
  end
  
  for i=1:n,[POP(i,l+2)]= fitness1(POP(i,l+1))
  end
  
  [FITNESS_ALL]=fitness(POP,l)
  
  [MFITNESS]=maxfitness(POP,l)
  
  [POP,AFITNESS]=rfitness(POP,l,FITNESS_ALL)
  
  [POP,ratio,f1,f0]=fratio10(POP,l,show)
  
  tcw = log((n-1)^2)/log(ratio)
  tcabin = log(n-1)/log(ratio)
  
endfunction




//***************************************************
//Function to compute the relative count of a string
// along with the realized new count by 'rounding up' 
// Compute all the l+4 and l+5-th positions
//
// Problem: if the sum of the new proposals is  bigger than n
// Solution: 
// l := length of fitness strings
// n :=  number of individuals in POP


function[POP]=newMember(POP,l,n)

[r,c]=size(POP);

for j=1:r
  POP(j,l+4)=POP(j,l+3)*n
  POP(j,l+5)=round(POP(j,l+4))
end

//Check for sum of proposals

r2 = 0
for j=1:r, r2=r2+POP(j,l+5), end

while r2>n
  z=round(rand()*r)
  if z==0 then z=1,end //if
  POP(z,l+5) = POP(z,l+5)-1
  r2=0, for j=1:r, r2=r2+POP(j,l+5), end//for
  end//while r2

endfunction



//***************************************************
// Dummy function to replace newMember in simplegas0
//
// along with the realized new count by 'rounding up' 
// Compute all the l+4 and l+5-th positions
// l := length of fitness strings
// n :=  number of individuals in POP

function[POP]=newMember0(POP,l,n)

[r,c]=size(POP);

for j=1:r
  POP(j,l+4)=1
  POP(j,l+5)=1
end


endfunction

//***************************************************
// Make a copy of a string from row_old to row_new: L --> R
// where the copy starts at position l+p+1
// l := length of string
// j := row_old
// k := row_new
// p := number of parameters (P=5)

function[POP]=strcpyLR(POP,l,p,j,k)
  
  show=0
[r,c]=size(POP);

// Testing the boundaries

if (j < 1) then if show == 1 then printf('ERROR: Position of old String outside
of Matrix! r = %d, j = %d', r,j),end, j=1,
  elseif  (j > r) then  if show == 1 then printf('ERROR: Position of old String
outside of Matrix!, r = %d, j = %d', r,j),end, j=r,
  elseif (k < 1) then  if show == 1 then printf('ERROR: Position of new String
outside of Matrix!, r = %d, k = %d', r,k),end, k=1
  elseif (k > r) then if show == 1 then  printf('ERROR: Position of new String
outside of Matrix!, r = %d, k = %d', r,k), end, k=r
  end
  
  // Making a copy
  
for i=1:l
    POP(k,l+p+i) = POP(j,i)
    end

endfunction

//***************************************************
// Make a copy of a string from row_old to row_new: R --> L
// where the copy starts at position l+p+1
// l := length of string
// j := row_old
// k := row_new
// p := number of parameters (p=5)

function[POP]=strcpyRL(POP,l,p,j,k)
  
  show=0
[r,c]=size(POP);

// Testing the boundaries

if (j < 1) then if show == 1 then  printf('ERROR: Position of old String outside
of Matrix! r = %d, j = %d', r,j),end, j=1
  elseif  (j > r) then if show == 1 then  printf('ERROR: Position of old String
outside of Matrix!, r = %d, j = %d', r,j),end, j=r
  elseif (k < 1) then  if show == 1 then printf('ERROR: Position of new String
outside of Matrix!, r = %d, k = %d', r,k), end, k=1
  elseif (k > r) then if show == 1 then  printf('ERROR: Position of new String
outside of Matrix!, r = %d, k = %d', r,k),end, k=r
  end
  
  // Making a copy
  
for i=1:l
    POP(j,i) = POP(k,l+p+i)
    end

endfunction

//***************************************************
// Make a copy of all strings from the old rows to the new ones
// where the new copies will start at position l+p+1
// l := length of string
// j := row_old
// k := row_new
// p := number of parameters
// r2 := memory of position for new strings

function[POP]= newPop(POP,l,p,n)

[r,c]=size(POP);
r2 = 1
  
  for j=1:r
    
    if POP(j,l+7) > 0 then r3 = POP(j,l+7)-1,
      
      for k=r2:(r2+r3), //printf('\n k = %d\n', k)
        [POP]=strcpyLR(POP,l,p,j,k)
      end
      
      r2=k+1
      
      end
    end

endfunction

//***************************************************
// Make a copy of all strings from the old rows to the new ones
// where the new copies will start at position l+p+1
//
// Problem: If the sum of new copies is greater than n then this causes a fault!
// Solution: At least check the sum before operation
//
// l := length of string
// j := row_old
// k := row_new
// p := number of parameters
// r2 := memory of position for new strings

function[POP]= newPop2(POP,l,p,n)

[r,c]=size(POP);

//Check for sum of proposals

r2 = 0

for j=1:r, r2=r2+POP(j,l+5), end

if r2 >n then error('In newPop2() sum of proposals >n!!!'),end
  
  
r2 = 1
  
for j=1:r
    
    if POP(j,l+5) > 0 then r3 = POP(j,l+5)-1,
      
      for k=r2:(r2+r3), //printf('\n k = %d\n', k)
        [POP]=strcpyLR(POP,l,p,j,k)
      end
      
      r2=k+1
      
      end
    end

endfunction


//***************************************************
// Prepare the crossover operations within a population POP 
// by randomly selecting the new strings and copy them onto the 
// POP base area  [1,l] 
//
// General Assumption: the number of members n is even!!!
//
// The strings are assumed to be in the area starting at l+p+1
// The parameter at l+p has been changed to a flag '1' := not yet mated
// j := target position in the base area for copy action


function[POP]= crossoverPrep(POP,l,p,n)

[r,c]=size(POP);

// Set a flag at position l+7 with '1'
  
  for j=1:r
     POP(j,l+7)=1
   end
   
// Select randomly the strings for transfer   

   j = 1 //:= Baseline for filling up with strings
   
 while(j < n+1)
   
   S = 1 //:= Flag for searching a string to be copied
   
   while(S == 1)
         
         k1 = round(n * rand())
         if k1==0 then k1=1
         end //if
   
   if POP(k1,l+7) <> 0 then POP=strcpyRL(POP,l,p,j,k1)
     //mshow(POP,n,l+p+l)
     POP(k1,l+7) = 0,
      S=0
    else S=1
    end //if
    
   end //while == S
    
       //      printf('k1 = %d\n',k1)
       //   printf('j = %d\n',j)
       
       j = j+1
       
   // mshow(POP,n,l+p+l)
 
 end // while == j 
   
    
endfunction

//***************************************************
// Apply a crossover operation onto two adjacent strings at intersection x
// The strings are assumed to be randomly paired in the area at position 1
//
// p := number of parameters between strings left and right (actually  p=5)
// l := length of strings ( actually l=5)
// n := number of elements in POP (actually n=4)

function[POP]= crossover(POP,l,p,n)
  
  // Take a first row
  r=1 
  while(r < n)
    
   // printf('r = %d\n',j)
   
   // Select randomly a point in this row
   
   x = round(l * rand()), if(x == 0) then x=1,end,
   
        // printf('x = %d\n',x)
         
         for i=x:l, m=POP(r+1,i)
           POP(r+1,i) = POP(r,i)
           POP(r,i) = m
         end
         
         r=r+2
    end

endfunction


//***************************************************
// Apply a mutation operation onto one randomly selected string 
// Select randomly a point in the string and replace it's
// value by it's inverse '1' --> '0', '0' --> '1'

function[POP]= mutation(POP,l,p,n)
  
  // Find a point in a string
  
     c = round(l * rand())
         if c == 0 then c=1  // This does not change anything
         end
         
 // Find a row in the table
         
     r = round(n * rand())
         if r == 0 then r=1  // This does not change anything
         end
         
 // replace the value
         
         if(POP(r,c) == 1) then POP(r,c) = 0,
         else POP(r,c) = 1
           end
           
 printf('mutation point at row = %d col = %d\n',r,c)         

endfunction

//***************************************************
// Print the content of a matrix
//
// d := depth of matrix 'downwards'
// w := width of  matrix from left to right



function[M]= mshow(M,d,w)
  
  for j=1:d,// printf('\n j= %3.1d : ',j)
    for i=1:w, //printf(' %3.1d ',M(j,i)) 
      end
      //printf('\n')
    end

endfunction

//***************************************************
// Translate strings with binaries '0', '1' as decimal numbers
//
// v = vector of binaries from a POP-matrix (left part)
// l := length of string
// D := decimal computed out of binaries


function[D]= vec2dec(v,l)
  
  str=string(v(1:l))
  for i=2:l, str(1)=str(1)+str(i)
  end
  D=bin2dec(str(1))
  
endfunction

//***************************************************
// Translate strings with binaries '0', '1' as decimal numbers
//
// v = vector of binaries from a POP-matrix (left part)
// l := length of string
// D := decimal computed out of binaries


function[D]= vec22dec(v,l1,l2)
  
  
  str=string(v(l1:l2))
  for i=2:l2-l1+1, str(1)=str(1)+str(i)
  end
  D=bin2dec(str(1))
  
endfunction

//***************************************************
// Translate strings with binaries '0', '1' and n-many compartments as decimal
numbers
//
// v = vector of binaries from a POP-matrix (left part)
// l := length of string
// b := number of bins for one number
// D := Array of decimals


function[D]= vec222dec(v,l,b,show)
  
  D=[]
  k=1
  
  for j=1:b:l+1-b
    if show==1 then printf('j = %d ',j),end
      
  str=string(v(j:j+b-1))
  for i=2:b, str(1)=str(1)+str(i)
  end //i
  
  if show==1 then printf('k = %d ',k),end
  
  D(k)=bin2dec(str(1))
  
  
  if show==1 then printf('D(k) = %d\n',D(k)),end
  
  k=k+1
  
  end //j
  
endfunction

//***************************************************
// maxfit01 builds a string of length l and computes the decimal value
//
// l := length of string
// MF := decimal computed out of binaries


function[MF]= maxfit01(l)
  
  FitString=""
  v=ones(1,l)
  str=string(v)
  for i=1:l, FitString = FitString + str(i)
  end
  MF=bin2dec(FitString)
  
endfunction

//***************************************************
// maxfi02 computes with maxfit01 the decimal value of a string with length l
and
// then compares the maximal value according to the fitness function y=x^2 
// multiplied by the number n of strings in a population 
//
// l := length of string
// MF := decimal computed out of binaries


function[MAXFIT]= maxfit02(l,n)
  
  MAXFIT=n*(maxfit01(l)^2)
  
endfunction

//***************************************************
// Simple fitness-function f=(x^2)
//
// D := decimal computed out of binaries


function[F]= fitness1(D)
  
  F=D^2
  
endfunction

//***************************************************
// Fitness-function f for a world W which is realized as a table t:D --> R
// IF (d,r) in t then f(d1,d2)=100, ELSE f(d1,d2)=0
//
// d1,d2 := decimal values of population POP
// W := a world organized as a table
// n := number of elements in W
// F := fitness value


function[F]= fitness2(d1,d2,W)
  
  [r,c]=size(W)
  upper =10
  goal = upper*2
  
  i=1
  while i < r+1, 
    if W(i,1) == d1 & W(i,2) == d2 then F=goal, i=r+1,
    else F=round(upper * rand()), i=i+1,
    end //if
    end//while
  
endfunction

//***************************************************
// Fitness-function f for a world W which is realized as a table t:D --> R
// IF (d,r) in t then f(d1,d2)=goal, ELSE f(d1,d2)=random(upper)
//
// D := Array of pairs of decimal values
// W := a world organized as a table
// F := fitness value


function[F,goal]= fitness3(D,W)
  
  [r,c]=size(W)
  [r2,c2]=size(D)
  F2=[]
  
  upper =10  // range of numbers for non-goals [0,upper]
  goal = upper*2  
  
  j=1 //Index for D
  for i=1:r //Look to every pair in W
    if W(i,1) == D(j) & W(i,2) == D(j+1) then F2(i)=goal,
    else F2(i)=round(upper * rand()),
    end //if
  j=j+2
end // for i

F=sum(F2)
  
endfunction


//*************************************************************
// Function to compute only the fitness of a population
//
// Input:
// p := number of parameters between strings left and right (actually  p=5)
// l := length of strings
// n := number of elements in POP
// Output
// POP := a population
// MFITNESS := Maximal Fitness of each individual

function[POP,MFITNESS, FITNESS_ALL]=popFit(l,p,n)
  
  [POP]=popgen(n,l+p)
  
  MFITNESS=0
 
  
  for j=1:n, v=POP(j,:),  [POP(j,l+1)]= vec2dec(v,l)
  end
  
  for i=1:n,[POP(i,l+2)]= fitness1(POP(i,l+1))
  end
  
  [FITNESS_ALL]=fitness(POP,l)
  
  [MFITNESS]=maxfitness(POP,l)
  
endfunction


//*************************************************************
// Function to compute  the maximal fitness of an individuum 
// for several populations
//
// Input:
// p := number of parameters between strings left and right (actually  p=5)
// l := length of strings
// n := number of elements in POP
// run := number of choices
// Output
// POP := a population
// MFITNESS := Maximal Fitness of each individual
// MFITNESSX := Array with maximal fitness

function[POP,MEAN, MFITNESSX]=popFitX(l,p,n,run)
  
  MAXFITUP = ((2^l)^2)
  
  for k=1:run
      [POP]=popgen(n,l+p)
  
  for j=1:n, v=POP(j,:),  [POP(j,l+1)]= vec2dec(v,l)
  end
  
  for i=1:n,[POP(i,l+2)]= fitness1(POP(i,l+1))
  end
  
  [MFITNESSX(k)]=maxfitness(POP,l)
  
end

for k=1:run
  MFITNESSX(k)=MFITNESSX(k)/(MAXFITUP/100)
end
MEAN = mean(MFITNESSX)

// Show graphical results

//clf(), xdel, 

//plot2d([1:1:run], MFITNESSX)
  
endfunction



//*************************************************************
// Function to compute  the maximal fitness of an individuum 
// for several populations
// for increasing n
//
// Input:
// p := number of parameters between strings left and right (actually  p=5)
// l := length of strings
// n := number of elements in POP
// run := number of choices
// Output
// POP := a population
// MFITNESS := Maximal Fitness of each individual
// MFITNESSX := Array with maximal fitness

function[MEAN, MEANX]=popFitXN(l,p,m,run)
  
  for i=1:m
    
  [POP,MEAN, MFITNESSX]=popFitX(l,p,i,run)
  
  MEANX(i) = MEAN
  
end

MEAN = mean(MEANX)

// Show graphical results

clf(), xdel, 

plot2d([1:1:m], MEANX)
  
endfunction


//*************************************************************
// Function to compute  the maximal fitness of a population 
// for several generations
//
// Input:
// p := number of parameters between strings left and right (actually  p=5)
// l := length of strings
// n := number of elements in POP
// run := number of generations
// Output
// POP := a population
// FITNESS_ALL := Maximal Fitness of a population
// FITNESS_ALLX := Array with maximal fitness
// MEAN := Mean value of array

function[STD,MEAN, FITNESS_ALLX]=popMaxFitX(l,p,n,run)
  
  MAXFITUP = ((2^l)^2)*n
  
  for k=1:run
      [POP]=popgen(n,l+p)
  
  for j=1:n, v=POP(j,:),  [POP(j,l+1)]= vec2dec(v,l)
  end
  
  for i=1:n,[POP(i,l+2)]= fitness1(POP(i,l+1))
  end
 
  [FITNESS_ALLX(k)]=fitness(POP,l)
  
   end

for k=1:run
  FITNESS_ALLX(k)=FITNESS_ALLX(k)/(MAXFITUP/100)
end
MEAN = mean(FITNESS_ALLX)
mf=tabul(FITNESS_ALLX)
STD=stdevf(mf([:],1), mf([:],2))

// Show graphical results

//clf(), xdel, 

//plot2d([1:1:run], MFITNESSX)
  
endfunction



//*************************************************************
// All Functions Unified
//
// p := number of parameters between strings left and right (actually  p=5)
// l := length of strings ( actually l=5)
// n := number of elements in POP (actually n=4)
// POP := a predefined population (can be done automatically)
// run := number of cycles
// MThreshold := number of cycles which have to be waited until the mutation
operators will be applied

function[POP,FITNESS_ALLLOG]=gasimple(POP,l,p,n,run, MThreshold)
  
  MCount = 0 
  
  FITNESS_ALLLOG=[]
  
for cyc = 1:run  
  
  [M]= mshow(POP,n,l+p)
  
  for j=1:n, v=POP(j,:),  [POP(j,l+1)]= vec2dec(v,l)
  end
  
  for i=1:n,[POP(i,l+2)]= fitness1(POP(i,l+1))
  end
  
  [FITNESS_ALL]=fitness(POP,l)
  FITNESS_ALLLOG(cyc)=FITNESS_ALL
  
  [MFITNESS]=maxfitness(POP,l)
  
  [POP,AFITNESS]=rfitness(POP,l,FITNESS_ALL)
  
  [POP]=newMember(POP,l,n)
  [POP]=newPop(POP,l,p,n)
  [POP]= crossoverPrep(POP,l,p,n)
  [POP]= crossover(POP,l,p,n)
  
  printf('Mutationcount = %d\n', MCount)
  
  if(MCount > MThreshold) then [M]= mshow(POP,n,l+p), 
    [POP]= mutation(POP,l,p,n), 
    [M]= mshow(POP,n,l+p), 
    MCount=0,
  end
  
  MCount = MCount + 1
  
end //for == run

// Compute the maximal value and the percentage of success

[MAXFIT]=maxfit02(l,n)
FITNESS_ALL_PERC=[]

for i=1:run, FITNESS_ALL_PERC(i) = FITNESS_ALLLOG(i)/(MAXFIT/100), end

// Show graphical results

clf(), xdel, 
//plot2d([1:1:run], FITNESS_ALLLOG) 
plot2d([1:1:run], FITNESS_ALL_PERC) 



endfunction

//*************************************************************
// A system with GA without any additional parameters
//
// p := number of parameters between strings left and right (actually  p=5)
// l := length of strings ( actually l=5)
// n := number of elements in POP (actually n=4)
// POP := a predefined population (can be done automatically)
// N := number of cycles
// MThreshold := number of cycles which have to be waited until the mutation
operators will be applied

function[FITNESS_ALL_PERC,POP]=ga(POP,l,p,n,N, MThreshold,show)
  
  MCount = 0
  FITNESS_ALLLOG=[]
  
for cyc = 1:N 
  
  for j=1:n, v=POP(j,:),  
    d= vec2dec(v,l),  //decimal value
    [POP(j,l+1)]= d,  //decimal value
  end
  
  for i=1:n,[POP(i,l+2)]= fitness1(POP(i,l+1))
  end
  
  [FITNESS_ALL]=fitness(POP,l)
  FITNESS_ALLLOG(cyc)=FITNESS_ALL
  
  [MFITNESS]=maxfitness(POP,l)
  
  [POP,AFITNESS]=rfitness(POP,l,FITNESS_ALL)
  
  [POP]=newMember(POP,l,n)
  [POP]=newPop(POP,l,p,n)
  [POP]= crossoverPrep(POP,l,p,n)
  [POP]= crossover(POP,l,p,n)
  
  if show==1 then
    printf('Mutationcount = %d\n', MCount),
    end
    
  
  if(MCount > MThreshold) then 
    [POP]= mutation(POP,l,p,n), 
     if show==1 then  disp(POP),  end
      MCount=0,
  end
  
  MCount = MCount + 1
  
end //for == N
  
  printf("Number of Events n * N = %d\n",n*N)
  
  // Compute the maximal value and the percentage of success 

[MAXFIT]=maxfit02(l,n)
FITNESS_ALL_PERC=[]

for i=1:N, FITNESS_ALL_PERC(i) = FITNESS_ALLLOG(i)/(MAXFIT/100), end

// Show graphical results

if show==2 then
clf(), xdel,  
plot2d([1:1:N], FITNESS_ALL_PERC),
end


endfunction



//*************************************************************
// A system with GA
//
// p := number of parameters between strings left and right (actually  p=5)
// l := length of strings ( actually l=5)
// n := number of elements in POP (actually n=4)
// POP := a predefined population (can be done automatically)
// N := number of cycles
// MThreshold := number of cycles which have to be waited until the mutation
operators will be applied

function[FITNESS_ALL_PERC,DIST2,STD, MEAN, FREQ,STD1, MEAN1, FREQ1,FX,
POP]=ga0(POP,l,p,n,N, MThreshold,show)
  
 // Install the counters
  
  r= (2^l)
  c= 6
  
  FX = zeros(r,c)
  MCount = 0 
  FITNESS_ALLLOG=[]
  
   // Generate an internal index for display
  
  for j=1:r
    FX(j,1)=j-1
    end
 
  
for cyc = 1:N 
  
  for j=1:n, v=POP(j,:),  
    d= vec2dec(v,l),  //decimal value
    [POP(j,l+1)]= d,  //decimal value
    FX(d+1,2)=FX(d+1,2)+1, //occurences
  end
  
  for i=1:n,[POP(i,l+2)]= fitness1(POP(i,l+1))
  end
  
  [FITNESS_ALL]=fitness(POP,l)
  FITNESS_ALLLOG(cyc)=FITNESS_ALL
  
  [MFITNESS]=maxfitness(POP,l)
  
  [POP,AFITNESS]=rfitness(POP,l,FITNESS_ALL)
  
  [POP]=newMember(POP,l,n)
  [POP]=newPop(POP,l,p,n)
  [POP]= crossoverPrep(POP,l,p,n)
  [POP]= crossover(POP,l,p,n)
  
  if show==1 then
    printf('Mutationcount = %d\n', MCount),
    end
    
  
  if(MCount > MThreshold) then 
    [POP]= mutation(POP,l,p,n), 
     if show==1 then  disp(POP),  end
      MCount=0,
  end
  
  MCount = MCount + 1
  
end //for == N

for j=1:r
  
    FX(j,3)=(FX(j,2)/n)/N //Normal frequency
    FX(j,4)=(FX(j,2)/n)/(((1/2^l)*N)/100) // Frequency as percentage
  end//m 

  
  disp(FX)
  
  FREQ1=tabul(FX(:,3))
  MEAN1=mean(FREQ1(:,1)) 
  STD1= stdevf(FREQ1(:,1), FREQ1(:,2))
  FREQ=tabul(FX(:,4))
  MEAN=mean(FREQ(:,1)) 
  STD= stdevf(FREQ(:,1), FREQ(:,2))
  MAX=max(FREQ(:,1))
  MIN=min(FREQ(:,1))
  DIST=MAX-MIN
  DIST2=DIST/2
  
  
  printf("Number of Events n * N = %d\n",n*N)
// Compute the maximal value and the percentage of success

[MAXFIT]=maxfit02(l,n)
FITNESS_ALL_PERC=[]

for i=1:N, FITNESS_ALL_PERC(i) = FITNESS_ALLLOG(i)/(MAXFIT/100), end

// Show graphical results

if show==2 then
clf(), xdel,  
plot2d([1:1:N], FITNESS_ALL_PERC),
end



endfunction



//*************************************************************
// GA algorithm without fitness
//
// p := number of parameters between strings left and right (actually  p=5)
// l := length of strings ( actually l=5)
// n := number of elements in POP (actually n=4)
// POP := a predefined population (can be done automatically)
// N := number of events
// K := Number of repetitions from calling function
// MThreshold := number of events which have to be waited until the mutation
operators will be applied

function[MCount,DIST2,STD, MEAN, FREQ,STD1, MEAN1,
FREQ1,POP,FX]=gasimple0(POP,l,p,n,N, K,MThreshold,MCount,show)
  
  
   // Install the counters
  
  r= (2^l)
  c= 5
 
  
  FX = zeros(r,c)
  
  // Generate an internal index for display
  
  for j=1:r
    FX(j,1)=j-1
  end
  
 
  
  for cyc = 1:N*K //Start Cyle for events
    

    if show ==1 then  disp(cyc), end
    
      for j=1:n
    v=POP(j,:),  d = vec2dec(v,l), [POP(j,l+1)]= d, FX(d+1,2)=FX(d+1,2)+1 
  end//n
  
   if show ==1 then  printf("\nFX after Count\n\n"), disp(FX), end
    
  [POP]=newMember0(POP,l,n) //Prepares copies by inserting '1' at position l+4
and l+5
  [POP]=newPop(POP,l,p,n) //Copies old strings to the right
  [POP]= crossoverPrep(POP,l,p,n)
  [POP]= crossover(POP,l,p,n)
  
 if show ==1 then printf('Mutationcount = %d\n', MCount),disp(POP), end
  
  if(MCount > MThreshold) then   [POP]= mutation(POP,l,p,n), 
    
    if show ==2 then printf('Mutationcount (MCount, cyc)= (%d,%d)\n', MCount,
cyc), end
    
    MCount = 0, end
  
  MCount = MCount + 1
  
  if show==1 then printf("POP after Operators\n"), disp(POP), end
 
  
end //cyc=N

for j=1:r
  
    FX(j,3)=(FX(j,2)/n)/(N*K) //Normal frequency
    FX(j,4)=(FX(j,2)/n)/(((1/2^l)*(N*K))/100) // Frequency as percentage
  end//m 
  
  
  disp(FX)
  
  FREQ1=tabul(FX(:,3))
  MEAN1=mean(FREQ1(:,1)) 
  STD1= stdevf(FREQ1(:,1), FREQ1(:,2))
  FREQ=tabul(FX(:,4))
  MEAN=mean(FREQ(:,1)) 
  STD= stdevf(FREQ(:,1), FREQ(:,2))
  MAX=max(FREQ(:,1))
  MIN=min(FREQ(:,1))
  DIST=MAX-MIN
  DIST2=DIST/2
  
  
  printf("Number of Events n * N * K = %d\n",n*N*K)
  
endfunction
  

//*************************************************************
// GA algorithm with fitness and conditioned mutation
//
// Works like 'normal' GA but uses mutation only if the actual maximal fitness
// stays 'constant' over MT-many cycles and is below te theoretical maximum
//
// p := number of parameters between strings left and right (actually  p=5)
// l := length of strings ( actually l=5)
// n := number of elements in POP (actually n=4)
// POP := a predefined population (can be done automatically)
// N := number of events
// MT := mutation trigger; number of values to monitor for to trigger mutation

function[FITNESS_ALL_PERC,DIST2,STD, MEAN, FREQ,STD1, MEAN1, FREQ1,FX,
POP]=ga1(POP,l,p,n,N, MT,show)
  
  
  //Theoretical Maximum for Fitness according to y=x^2
  
  TMaxFit = ((2^l)-1)^2*n
  
 // Install the counters
  
  r= (2^l)
  c= 6
  
  FX = zeros(r,c)

  MCount = 0 
  FITNESS_ALLLOG=[]
  MTriggerEstimator=[]
  
   // Generate an internal index for display
  
  for j=1:r
    FX(j,1)=j-1
    end
 
  
for cyc = 1:N 
  
  for j=1:n, v=POP(j,:),  
    d= vec2dec(v,l),  //decimal value
    [POP(j,l+1)]= d,  //decimal value
    FX(d+1,2)=FX(d+1,2)+1, //occurences
  end
  
  for i=1:n,[POP(i,l+2)]= fitness1(POP(i,l+1))
  end
  
  [FITNESS_ALL]=fitness(POP,l)
  FITNESS_ALLLOG(cyc)=FITNESS_ALL
  
  [MFITNESS]=maxfitness(POP,l)
  
  if show ==1 then disp(FITNESS_ALL), end
    
    MTriggerEstimator(cyc)=FITNESS_ALL
    
    
  if show ==1 then disp(MTriggerEstimator(cyc)), end
  
  [POP,AFITNESS]=rfitness(POP,l,FITNESS_ALL)
  
  [POP]=newMember(POP,l,n)
  [POP]=newPop(POP,l,p,n)
  [POP]= crossoverPrep(POP,l,p,n)
  [POP]= crossover(POP,l,p,n)
  
   if show ==1 then disp(MCount),   end
    
  
  if(MCount > MT) then 
    
     if show ==1 then disp(MTriggerEstimator(cyc)),
disp(MTriggerEstimator(cyc-MT)), disp(TMaxFit), end
    if (MTriggerEstimator(cyc) == MTriggerEstimator(cyc-MT) &
MTriggerEstimator(cyc) < TMaxFit) then [POP]= mutation(POP,l,p,n), 
      if show==1 then  disp(POP),  end
      end
      MCount=0,
  end
  
  MCount = MCount + 1
  
end //for == N

for j=1:r
  
    FX(j,3)=(FX(j,2)/n)/N //Normal frequency
    FX(j,4)=(FX(j,2)/n)/(((1/2^l)*N)/100) // Frequency as percentage
  end//m 

  
  disp(FX)
  
  FREQ1=tabul(FX(:,3))
  MEAN1=mean(FREQ1(:,1)) 
  STD1= stdevf(FREQ1(:,1), FREQ1(:,2))
  FREQ=tabul(FX(:,4))
  MEAN=mean(FREQ(:,1)) 
  STD= stdevf(FREQ(:,1), FREQ(:,2))
  MAX=max(FREQ(:,1))
  MIN=min(FREQ(:,1))
  DIST=MAX-MIN
  DIST2=DIST/2
  
  
  printf("Number of Events n * N = %d\n",n*N)
// Compute the maximal value and the percentage of success

[MAXFIT]=maxfit02(l,n)
FITNESS_ALL_PERC=[]

for i=1:N, FITNESS_ALL_PERC(i) = FITNESS_ALLLOG(i)/(MAXFIT/100), end

// Show graphical results

if show==2 then
clf(), xdel,  
plot2d([1:1:N], FITNESS_ALL_PERC),
end



endfunction

//*************************************************************
// A system with GA for a world with a table
//
// A world W organized as a table  maps a set D into a set R
// l := length of strings
// p := number of cells between string <1...l> and <l+p+1, ..., l+p+l>
// Pos 1-l := String
// Pos l+1 := Decimal Value of binary string
// Pos l+2 := Fitness for l+1 and l+p
// Pos l+3 := Percentage of overall fitness
// Pos l+4 := Expected count according to fitness
// Pos l+5 := Realized count
// Pos l+6 := 2nd Decimal value of a compund fitness
// Pos l+p := flag for crossover
// p := 7
// POP := a predefined population (can be done automatically)
// n := number of elements in POP
// N := number of cycles
// MT := number of cycles which have to be waited until the mutation operators
will be applied

function[FITNESS_ALL_PERC,POP]=gaw1(W,POP,l,p,n,N, MT,l1,show)
  
  MCount = 0
  FITNESS_ALLLOG=[]
  
  for cyc = 1:N 
    
   if show == 1 then disp('cyc ='),disp(cyc),disp(' '),end
  
  for j=1:n, v=POP(j,:),  
    d1= vec22dec(v,1,l1),  //decimal value1
    d2= vec22dec(v,l1+1,l),  //decimal value2
    [POP(j,l+1)]= d1,  //decimal value 
    [POP(j,l+6)]= d2,  //decimal value 
  end
  
  
  if show == 1 then printf('New Numbers in POP '), disp(POP), end
  
  for i=1:n,[POP(i,l+2)]= fitness2(POP(i,l+1),POP(i,l+6),W)
  end
  
  
  [FITNESS_ALL]=fitness(POP,l)
  
  if show == 1 then printf('FITNESS_ALL = %f\n',FITNESS_ALL), end
  
  FITNESS_ALLLOG(cyc)=FITNESS_ALL
  
  [MFITNESS]=maxfitness(POP,l)
  
  
  [POP,AFITNESS]=rfitness(POP,l,FITNESS_ALL)
  
  if show == 1 then disp(POP), end
  
  [POP]=newMember(POP,l,n)
  [POP]=newPop2(POP,l,p,n)
  
  if show == 1 then disp(POP), end
  
  [POP]= crossoverPrep(POP,l,p,n)
  [POP]= crossover(POP,l,p,n)
  
    if show == 1 then disp(POP),printf('MCount = %d\n',MCount), end

  if(MCount > MT) then 
    [POP]= mutation(POP,l,p,n), printf('Mutation at cycle = %d\n',cyc),
     if show==1 then  disp(POP),  end
      MCount=0,
  end
  
  MCount = MCount + 1
  
end //for == N
  
  printf("Number of Events n * N = %d\n",n*N)
  
  // Compute the maximal value and the percentage of success 
  
MAXFIT=20*n   //The value '20' is set in the function fitness2()
FITNESS_ALL_PERC=[]

for i=1:N, FITNESS_ALL_PERC(i) = FITNESS_ALLLOG(i)/(MAXFIT/100), end

// Show graphical results

if show==2 | show == 1 then
clf(), xdel,  
plot2d([1:1:N], FITNESS_ALL_PERC),
end


endfunction


//*************************************************************
// A system with GA for a world with a table and with 
// multi-compartment genomes
//
// A world W organized as a table  maps a set D into a set R
// l := length of strings
// p := number of cells between string <1...l> and <l+p+1, ..., l+p+l>
// Pos 1-l := String
// Pos l+1 := Decimal Value of binary string
// Pos l+2 := Fitness for l+1 and l+p
// Pos l+3 := Percentage of overall fitness
// Pos l+4 := Expected count according to fitness
// Pos l+5 := Realized count
// Pos l+6 := 2nd Decimal value of a compund fitness
// Pos l+7 := flag for crossover
// p := 7
// POP := a predefined population (can be done automatically)
// n := number of elements in POP
// N := number of cycles
// MT := number of cycles which have to be waited until the mutation operators
will be applied
// V := version

function[FITNESS_ALL_PERC,MAXGENOMECOUNT,POP,V]=gaw2(W,POP,l,n,N, MT,b,show,V)
  
  p = 7
  MCount = 0
  FITNESS_ALLLOG=[]
  
  for cyc = 1:N 
    
   if show == 1 then disp('cyc ='),disp(cyc),disp(' '),end

    
  for i=1:n , v=POP(i,:), [D]= vec222dec(v,l,b,show),  [POP(i,l+2),goal]=
fitness3(D,W),  end //for i
  
  if show == 1 then printf('New Numbers in POP '), disp(POP), end
  
  
  [FITNESS_ALL]=fitness(POP,l)
  
  if show == 1 then printf('FITNESS_ALL = %f\n',FITNESS_ALL), end
  
  FITNESS_ALLLOG(cyc)=FITNESS_ALL
  
  [MFITNESS]=maxfitness(POP,l)
  
  
  [POP,AFITNESS]=rfitness(POP,l,FITNESS_ALL)
  
 // if show == 1 then disp(POP), end
  
  [POP]=newMember(POP,l,n)
  
  if show == 1 then disp(POP), end  
  
  [POP]=newPop2(POP,l,p,n)
  
  if show == 1 then disp(POP), end
  
  [POP]= crossoverPrep(POP,l,p,n)
  [POP]= crossover(POP,l,p,n)
  
    if show == 1 then disp(POP),printf('MCount = %d\n',MCount), end

  if(MCount > MT) then 
    [POP]= mutation(POP,l,p,n), printf('Mutation at cycle = %d\n',cyc),
     if show==1 then  disp(POP),  end
      MCount=0,
  end
  
  MCount = MCount + 1
  
end //for == N
  
  printf("Number of Events n * N = %d\n",n*N)
  
  // Compute the maximal value and the percentage of success 
  
  [r,c]=size(POP)
  [r2,c2]=size(W)
  
  MAXFIT=goal*r2*r //The value 'goal' is set in the function fitness3()
  MAXFITGENOME = goal*r2  // To count the complete genomes in a population
FITNESS_ALL_PERC=[]

MAXGENOMECOUNT = 0
for i=1:r, 
  if POP(i,l+2) == MAXFITGENOME then MAXGENOMECOUNT =MAXGENOMECOUNT+1, end,
  end
for i=1:N, FITNESS_ALL_PERC(i,1)=i, FITNESS_ALL_PERC(i,2) =
FITNESS_ALLLOG(i)/(MAXFIT/100), end

// Show graphical results

if show==2 | show == 1 then
  clf(), xdel, 
  scf(V), 
plot2d([1:1:N], FITNESS_ALL_PERC(:,2)),
end


endfunction



//*************************************************************
// Computing the worst case convergence time tcw
//
// p := number of parameters between strings left and right (actually  p=5)
// l := length of strings ( actually l=5)
// n := number of elements in POP (actually n=4)
// POP := a predefined population (can be done automatically)
// 
// f1 := fitness of all individuals with a allele value '1' at a position
// f0 := fitness of all individuals with a allele value '0' at a position 
// r = f1/f0
// tcw := worst case convergence time
//
// Assumption: POP has the relativ fitness available at l+3



function[POP,tcw,tcabin,ratio,f1,f0]=ftcw2(W,POP,l,p,n,l1,show)
  
  
   for j=1:n, v=POP(j,:),  
    d1= vec22dec(v,1,l1),  //decimal value1
    d2= vec22dec(v,l1+1,l),  //decimal value2
    [POP(j,l+1)]= d1,  //decimal value 
    [POP(j,l+6)]= d2,  //decimal value 
  end
  
  for i=1:n,[POP(i,l+2)]= fitness2(POP(i,l+1),POP(i,l+6),W)
  end
  
  [FITNESS_ALL]=fitness(POP,l)
  
  [MFITNESS]=maxfitness(POP,l)
  
  [POP,AFITNESS]=rfitness(POP,l,FITNESS_ALL)
  
  [POP,ratio,f1,f0]=fratio10(POP,l,show)
  
  tcw = log((n-1)^2)/log(ratio)
  tcabin = log(n-1)/log(ratio)
  
endfunction




//***********************************************************************
// Examples
//
POP = [0 1 1 0 1 13 169 0 0 0; 1 1 0 0 0 24 576 0 0 0; 0 1 0 0 0  8 64 0 0 0; 1
0 0 1 1 19 361 0 0 0;]

W=[1 3; 2 1; 3 2]

W2=[1 3; 2 1; 3 2; 4 4; 5 7; 6 5]



Gerd Doeben-Henisch 2012-03-31