// <pre>
// EVOGOD
// SciLab program to simulate the evolution of religion
// Version:              1.0
// Modification date:    1/28/2008
// Copyright 2008 James W. Dow
// This program is released under the GNU General Public License Version 2, June 1991. See http://www.gnu.org/licenses/gpl.txt

// SYSTEM PARAMETERS

ns = 140; ip1 = 10; le = 65; mc = 5; dcr = 0.01; dcu = 0.01; dcur = 0.01; dfir = 0.010; dfiu = -0.005; dfisu = -0.02; dp = 0; lr = 0; ur = 1; lv = 0;  uu = 1; dist = 0; popgrow = 0; greenbeard = 1; mult =1; cmr = 0.5; cmu = 0.5; icp = 0.2;

// MODEL

// Build the initial population matrix (a) of agents . We don't know the size yet. It is easier to build it and then compute the size. Each agent is a row.
nva = 6;  //  Number of variables (columns) for each agent vector
// The column variables for each agent are
          // 1. age (ag)
          // 2. initial probability of communicating a real experiences (ipr). A real experience is an experience that gives the holder knowledge that can enhance the fitness of others. This knowledge can be communicated.
          // 3. initial probability of communicating unreal experiences (ipu)  A number representing the probability that the agent will personally have an unreal experience. An unreal experience is something that can be communicated but contains no knowledge that can enhance fitness. The probability variables pr and pu represent biological evolution in the brain.
          // 4. developing probability of communicating a real experience (cpr). This changes during the simulation
          // 5. developing probability of communicating an unreal experience (cpu).
          // 6.fitness (fi) The number of offspring produced over a lifetime.

// INITIALIZATION of ages and number of agents
     // The number of agents in cohort k is ip1 - (k-1)*dp.
// Function for the number of agents in cohort k
function y = na(k), y = ip1 - (k-1) * dp, endfunction;

// Set up ages
// i indexes the cohort, j indexes the column of the a matrix.
// Do the first cohort. a(1) contains the age. The maximum number of cohorts is le, the initial life expectancy. However there are ilkley to be fewer if the population in each cohort decrements (dp) by a large amount
a = ones(na(1),nva);
// Do remaining cohorts
for i = 2:le;
     if na(i) < 0 then break; end;
         b = ones(na(i),nva) * i ;
         a = cat(1,a,b);
     end;

// Determine the size of the population
// Function for the number of rows in a matrix.
function y = sz(x), y=size(x,'r'); endfunction;
// Set initial population size
ips = sz(a);

// INITIAL COMMUNICATION PROBABILITIES
// Set up initial probalities of communicating a real experience, a(:,2) and, an unreal experience, a(:,3).
// Set up a random vector of  probabilities. If dist = 0 use a uniform distribution. If dist = 1 use a binomial distribution.
select dist ;
     case 0 then;   // Uniform distribution
          // Function to set up a random vector of uniformly distributed probabilities between l and u.
          function y = rv(l,u), y = (u - l) * rand(ips,1) + l, endfunction;
          a(:,2) = rv(lr,ur);
          a(:,3) = rv(lv,uu);
     case 1 then;   // Binomial distribution
          a(:,2) = grand(ips,1,'bin',ip1,ur) / ip1;
          a(:,3) = grand(ips,1,'bin',ip1,uu) / ip1;
     end;
// Learned communication probabilities. Start with icp and add according to inherited fitness.
a(:,4:5) = ones(ips,2)*icp;
// INITIAL FITNESS
a(:,6) = ones(ips, 1);

// Set up recording for the steps.
rps = ones(ns,1);		// Population size
mcp = zeros(ns,2); // Mean of the probability of commnicating real and unreal information

// STEP THROUGH THE SIMULATION
// Note that in this simulation it is possible for an agent to communicate with himself. This does not seem to be a conceptual problem and occurs in reality.

// Functions to provide a set of x random indices to agents in a. greenbeard determines the distribution from which the random indices are selected.

select greenbeard;

     case 0 then ; // The probabilities of selection are distributed uniformly over the population
          function y=ria(x), y1=grand(1,'prm',(1:sz(a))'), y=y1(1:x), endfunction;

     case 1 then ; // The probabilities of selection are distributed according to the current probabilities, cpu, of generating unreal communications. The higher the cpu of an agent, the greater the chance of it being selected.
          function y=ria(x);
          ps = size(a,'r')
          // Normalize the cpu so that they sum to 1
          y1 = a(:,5) / sum(a(:,5))
          // Form a cumulative distribution
          y2 = cumsum(y1)
          y = []    // Initialize y as a blank vector.
          for i = 1:x
               y3 = y2 < rand() // A T-F vector for lower values
               y =[y, sum(y3)+1] // Counts agents lower in the cumulative distribution and adds the count to the vector, thus creating the index of the area within which the random number lies.
               end
               endfunction;
     end;

// START STEPPING
for s = 1:ns;   // Do the step
     // Calculate population size for this step
     ps = sz(a);
     nrc = 0;  // monitor total real communications
     nuc = 0;  // monitor unreal communications
     for j = 1:ps; // Do the agents

     // Communicate REAL experiences
     // The number of agents communicated with (nca) will depend on the communnicators probablity of communicating real information. It will be a proportion of mc
          nca = round(a(j,4) * mc) ; // Get the number of agents to receive real communications from this one.
          // Select a set of agents to receive real communications. Note that i is a vector of indices.
          i = ria(nca); // select a subset of agents in a at random. i is a column vector. It can be empty
          // Have them learn to communicate real experiences
          if ~ (i == []) then;
               // Change the probability of the receivers communicating real information. This will depend on the inherited probability of communicating unreal real information
               a(i,4) = min(ones(size(i,'r'),1),a(i,4) + a(i,2)*dcr);
               // Add to their fitness.
               if mult == 0 then;
                    a(i,6) = a(i,6) + (cmr * dfir); // Multiply the amount by a constant, cmr
                    else; // Multiply by the receivers learned probability of communicating real information a(i,2)
                    a(i,6) = a(i,6) + (a(i,4) * dfir);
                    end;
               end;

          // Communicate UNREAL experiences. The number of agents contacted will depend on the probability of communicating unreal information
          nca = a(j,5) * mc ; // Calculate number
          nuc = nuc + nca;  // Add to unreal communications
          // Select a set of agents randomly
          i = ria(nca); // note that i is a vector and possibly null if there are less than .5 agents to contact
          // Are there agents to contact
          if ~ (i == []) then;
               // Change the probability of the receivers communicating unreal information
               a(i,5) = min(ones(size(i,'r'),1),a(i,5) + a(i,3)*dcu);
               // Change the probability of the receiver's communicating real information
               a(i,4) = min(ones(size(i,'r'),1),a(i,4) + a(i,2)*dcur);
               //Change the fitness of the receiver. dfiu will normally be negative
               if mult == 0 then;
                    a(i,6) = a(i,6) + (cmu * dfiu); // Multiply by a constant
                    else;
                    a(i,6) = a(i,6) + (a(i,5) * dfiu); //Multiply by the recivers learned probability of making unreal communications
                    end;
               end;
      // Make the sender pay a cost for sending unreal information. This is the cost of costly signaling, and can occur even if there are no receivers
          a(j,6) = a(j,6) + dfisu; // Note dfisu is normally negative
          end; // Agents are done. j is free

// EVOLVE THE POPULATION OF AGENTS
	// Advance the ages
     a(:,1) = a(:,1) + 1;

	//Find agents who will die
	ida = []; //Start with a blank index vector
     for i = 1:sz(a);
          if a(i,1) > le then; ida = [ida ; i]; // Add the index to ida
               end;
          end;
     // Index of dead agents (ida) is ready
     // Add offspring of dead agents
     // Calculate the number of offspring of each dead agents
     no = round(a(ida,6)); // vector of the number of offspring of agents indexed by ida
     // Form new offspring with age 1 and fitness 1 (not inherited)
     // Maintain the probabilities of communication of the parents
     for i = 1:sz(ida); // Scan dead agents
          j = ida(i) ; // index of the dead agent in a
          for k=1:no(i); // Scan number of offspring
               new = [1, a(j,2), a(j,3), icp, icp, 1];
               // Catenate new agent at the end of a
               a($+1,:) = new;
               end;
          end;
     a(ida,:) = [];     // Remove the dead ones from a.

// POPULATION GROWTH CONTROL
// Function to  randomly select x agents selected from a uniform distribution.

     function y=riau(x), y1=grand(1,'prm',(1:sz(a))'), y=y1(1:x), endfunction;

     select popgrow ;
          case 0 then; // Population can grow or decline. Do nothing
          case 1 then;// Control the increase or the decline
               dps = sz(a) - ips; // Calculate the change in population size
               if dps > 0 then ; // There is an increase in population
                    // Remove dps agents at random
                    i = riau(dps); // Get a uniform set at random
                    a(i,:) = [] ; // Kill them off
               elseif dps < 0 then ; // There is a decrease in population
                    i = riau(-dps) ; // Pick a set of agents at random
                    // Have them reproduce
                    // Maintain the probabilities of communication of the parents
                    new = [ones(-dps,1), a(i,2:3), a(i,2:3), ones(-dps,1)];
                    // Catenate new agents at the end of a
                    a = [a ; new];
                    end; // End if
          // There is neither an increase of a decrease
          end; // End of population control

// save results for this step

     ps = sz(a);    // Calculate population size
     rps(s,1) = ps; // Record the population size
     mcp(s,:) = [mean(a(:,2)), mean(a(:,3))]; // record means of pr and pu
// immediate display
     printf("step %i population %i \tmean pr %f mean pu  %f \n",s,ps,mcp(s,1),mcp(s,2));

     end; // Move to next step
// End of simulation


// PLOT RESULTS
xbasc(0);
scf(0);    // Open first plot window.
plot2d(1:ns,rps);	// Population size per step
ax = get("current_axes");
ax.x_label.text = "Step";
ax.x_label.font_size = 3;
ax.y_label.font_size = 3;
ax.sub_ticks = [3,1];
xbasc(1);
scf(1);   // Open second plot window.
plot2d(1:ns,mcp);   // Probability of real and unreal communications per step
ax = get("current_axes");
ax.x_label.text = "Step";
ax.x_label.font_size = 3;
ax.y_label.font_size = 3;
ax.sub_ticks = [3,1];
// Plot histograms of pr and pu
xbasc(2); scf(2);   // Open third plot window
title("Final Distribution of ipr", 'fontsize', 4);
x = nfreq(int(a(:,2)*10));bar (x(:,1),x(:,2))
xbasc(3); scf(3);   // Open fourth plot window
title("Final Distribution of ipu", 'fontsize', 4);
x = nfreq(int(a(:,3)*10));bar (x(:,1),x(:,2));
// </pre>