home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1997 July / macformat52.iso / mac / Shareware Plus / Educational / LEE 2.1 / Source / populati.c < prev    next >
Encoding:
C/C++ Source or Header  |  1996-08-07  |  17.8 KB  |  901 lines
  1. /* populati.c
  2.  *                         Copyright (1992)
  3.  *
  4.  *         Jeff Elman.  University of California, San Diego
  5.  *          Rik Belew.  University of California, San Diego
  6.  *      Stefano Nolfi.  Institute of Psychology, Rome.
  7.  *    Filippo Menczer.  University of California, San Diego
  8.  *        Greg Linden.  University of California, San Diego
  9.  *
  10.  *        This software may be redistributed without charge;
  11.  *                 this notice should be preserved.
  12.  */
  13.  
  14. #include "defs.h"
  15.  
  16.  
  17. /*
  18.  * main loop, let individuals evolve generation by generation
  19.  */
  20.  
  21. generati()
  22.  
  23. {
  24.     struct  indiv *ap;
  25.     char         buf1[40];
  26.     char         buf2[40];
  27.     int         flag;
  28.     int         k;
  29.     int         population;
  30.     register int     i;
  31.     register float    *input;
  32.     unsigned long         tempo, tempo0;
  33.     char   resfile[64];            /* the results file name      */
  34.     float reinforce;
  35.  
  36.  
  37.  
  38.     if (verbose>0) {
  39.         /*
  40.          * open the results file: the name is the seed
  41.          * and it contains:
  42.          * -- generation
  43.          * -- pop_size
  44.          * -- whatever else is printed by save_dat(),
  45.          *    depending on the experiment...
  46.          */
  47.         sprintf(resfile, "%lu.dat", seed);
  48.         resfp = fopen(resfile, "a");
  49.     }
  50.  
  51.  
  52.     /*
  53.      * go thru all generations requested
  54.      */
  55.     tempo0 = (unsigned long)time(NULL);
  56.     for (gen=startgeneration; gen < generations; gen++)
  57.     {
  58.  
  59.  
  60.         updateGenInteractive();
  61.  
  62.  
  63.         /*
  64.          * we open .err file
  65.          *
  66.          * ATTENTION: this file contains errors in a
  67.          * format different fron the original, due to
  68.          * having swapped the loops over indiv's and
  69.          * lifecycles!!!
  70.          */
  71.         if ((ecount > 0) && (errsig) && (verbose>0))
  72.         {
  73.             sprintf(buf1,"%d",gen);
  74.             genfilen(buf2,buf1,"err");
  75.             errfp = fopen(buf2, "w+");
  76.             if (errfp == NULL)
  77.             {
  78.                 printf("ERROR: Can't open .err file");
  79.                 exit(-1);
  80.             }
  81.         }
  82.  
  83.         if (verbose > 0) 
  84.         {
  85.             tempo = (unsigned long)time(NULL);
  86.             printf("\nRunning generation: %d,\t",gen);
  87.             printf("Pop: %d,\t",pop_size); 
  88.             printf("Time: %ld\n",tempo-tempo0);
  89.  
  90.             /*
  91.              * save data
  92.              */
  93.             fprintf(resfp, "%d\t%d", gen, pop_size);
  94.             save_dat(resfp);
  95.             fprintf(resfp, "\n");
  96.             fflush(resfp);
  97.  
  98.         }
  99.  
  100.             for (lifecycle=0; lifecycle < sweeps; lifecycle++)
  101.           {
  102.  
  103.                    if (verbose >= 2)
  104.                         printf("\ngen %d, time %d\n",gen,lifecycle);
  105.  
  106.  
  107.                /*
  108.                 * for each time, we let the individuals live;
  109.             * this is done stochastically, as 
  110.             * each individual is picked at random
  111.             * out of the current population
  112.                 */
  113.            population = pop_size;
  114.            for (k=0; k<population; k++)
  115.            {
  116.               ap = get_indiv(mrand(pop_size));
  117.               ap->lifecycle++;
  118.  
  119.  
  120.               EventsInteractive();
  121.  
  122.  
  123.             /*
  124.              * sense the world,
  125.              * activate the network
  126.              * and make an action
  127.              */
  128.             sense_world(ap);
  129.             activate_net(ap);
  130.             reinforce = act_in_world(ap);
  131.  
  132.  
  133.                         /*
  134.                          *  we print cycle by cycle informations
  135.                          */
  136.                         if ((printout) && (verbose>0))
  137.                            print_out(ap,lifecycle);
  138.  
  139.  
  140.             /*
  141.              * PREDICTION LEARNING:
  142.              * 1. we copy activations from input
  143.              * units to teaching inputs, in order
  144.              * for learning to compute the correct
  145.              * prediction error with the new
  146.              * sensory input (after the move)
  147.              * 2. we change weights
  148.              */
  149.             if (learn == 2)
  150.             {
  151.                 input = *ap->layerp;
  152.                 for (i=0; i<*layer_descp; i++) amoebat[i] = *input++;
  153.                 update_net(ap, amoebat);
  154.                 ap->error += global_error;
  155.                 if (errsig && (ecount > 0) && (((lifecycle + 1) % ecount) == 0) && (lifecycle > 0) && (verbose > 0))
  156.                 {
  157.                     fprintf(errfp,"%.2f ",ap->error);
  158.                     ap->error = 0.0;
  159.                 }
  160.             }
  161.  
  162.  
  163.             /*
  164.              * REINFORCEMENT LEARNING:
  165.              * we use change in energy as reinforcement signal:
  166.              * if positive, output is used as teaching input; 
  167.              * otherwise, the opposite of the output is used
  168.              */
  169.             if (learn == 1) r_learn(ap, reinforce);
  170.  
  171.  
  172.             /*
  173.              * if energy is high, reproduce:
  174.              * this is for now deterministic and 
  175.              * determined by the energy threshold
  176.              * ALPHA, which at the beginning is
  177.              * twice the energy expected average
  178.              */
  179.             if (ap->energy >= ALPHA) reproduce(ap);
  180.  
  181.             /*
  182.              * if energy is low, die
  183.              */
  184.             if (ap->energy <= 0.0) die(ap);
  185.              
  186.              
  187.             if (gDone)
  188.                 break;
  189.          
  190.            }
  191.  
  192.                    if ((ecount > 0) && (errsig) && (((lifecycle + 1) % ecount) == 0) && (verbose > 0))
  193.                      fprintf(errfp,"\n");
  194.  
  195.                   if (gDone)
  196.             break;
  197.  
  198.            /*
  199.             * replenish world food  
  200.             */
  201.            update_world();
  202.  
  203.               }
  204.  
  205.                if (verbose >= 2)
  206.                   printf("\n\nNext generation\n\n");
  207.  
  208.                if (((ecount > 0) || gDone) && (verbose>0))
  209.                    fclose(errfp);
  210.  
  211.         next_generation();
  212.  
  213.         if (gDone)
  214.             break;
  215.  
  216.            /*
  217.             * save this generation 
  218.             */
  219.         if ( gen == (generations - 1) || 
  220.            ((save_everyn > 0) && ((gen % save_everyn) == 0)))
  221.             
  222.             if (verbose>0)
  223.             {
  224.                 save_world(gen);
  225.                 save_all_indiv(gen);
  226.             }
  227.  
  228.            }
  229.  
  230.     if (verbose > 0)
  231.         fclose(resfp);
  232.  
  233. }
  234.  
  235.  
  236.  
  237.  
  238.  
  239. /*
  240.  * init first generation
  241.  */
  242.  
  243. init_pop(num_ind)
  244.  
  245.     int    num_ind;
  246. {
  247.     struct     indiv *ap;    
  248.     register int i;
  249.  
  250.     /*
  251.      * initialize the linked list of indiv structures;
  252.      * add_indiv() needs struct ptr, generation, id, parent id
  253.      * in the indiv struct (see defs.h) 
  254.      */
  255.     
  256.     a_head = (struct indiv *) malloc(sizeof(struct indiv));
  257.     if (a_head == NULL) {
  258.         fprintf(stderr, "a_head malloc err\n");
  259.         exit(-1);
  260.     }
  261.     ap = a_head;
  262.     ap->next = (struct indiv *)-1;
  263.     ap->last = (struct indiv *)0;
  264.     for (i = 0; i < num_ind; i++)
  265.     {
  266.         ap = add_indiv(ap, 0, i, 0);
  267.     }
  268.     end_head = ap;
  269.     end_head->id = -999;
  270. }
  271.  
  272.  
  273.  
  274. /*
  275.  * add an individual to the linked list of individual structures; take
  276.  * a valid (but empty) amoeba structure pointer, fill it out, and
  277.  * return the address of the next amoeba structure.  Other
  278.  * stuff (like weights) must get filled in by calling routine.
  279.  */
  280. struct indiv *
  281. add_indiv(ip, g, id, pid)
  282.  
  283.     struct  indiv *ip;
  284.     int    g;        /* generation number */
  285.     int    id;        /* id of this amoeba */
  286.     int    pid;        /* parent id */
  287. {
  288.     int i,j;
  289.     int size,count;
  290.  
  291.     /*
  292.      * we'd better not tromp on existing structures
  293.      */
  294.     if (ip->next != (struct indiv *)-1)
  295.            {
  296.         fprintf(stderr, "individual already exists\n");
  297.         exit(-1);
  298.            }
  299.         /*
  300.          * we allocate memory for the net
  301.          */
  302.         build_geno_net(ip);
  303.         build_pheno_net(ip);
  304.  
  305.        /*
  306.         * load actual indiv only for the initial generation
  307.         */
  308.        if (startgeneration == 0)
  309.        {
  310.     ip->generation = g;
  311.     ip->id         = id;
  312.     ip->pid        = pid;
  313.     ip->energy     = (float)mrand(ALPHA);
  314.     ip->gutsize    = gutsize;
  315.     ip->error      = 0.0;
  316.     ip->fit        = 0;
  317.     ip->lifecycle  = 0;
  318.     ip->worldx     = mrand(x_dim);
  319.     ip->worldy     = mrand(y_dim);
  320.     
  321.     switch (mrand(4)){
  322.         case 0:
  323.             ip->direction  = UP;
  324.             break;
  325.         case 1:
  326.             ip->direction  = LEFT;
  327.             break;
  328.         case 2:
  329.             ip->direction  = DOWN;
  330.             break;
  331.         case 3:
  332.             ip->direction  = RIGHT;
  333.             break;
  334.     }
  335.                 
  336.     for (i=0;i<types;i++) ip->gut[i] = 0;
  337.     count = 0;
  338.     for (i=0;i<nsensors;i++)
  339.     {
  340.  
  341.         ip->sensor_specs[i].system = s_type[i];
  342.  
  343.         size = sensor_size(ip->sensor_specs[i].system);
  344.         count += size;
  345.         for (j=0;j<*layer_descp;j++)
  346.             ip->sensor_specs[i].mask[j] = get_mask(count, size, j);    
  347.         for (j=0;j<types;j++)
  348.             ip->sensor_specs[i].complex[j] = 0;
  349.  
  350.         /*
  351.          * DEFAULT INITIALIZATION OF SENSOR COMPLEXES
  352.          * assuming    nsensors=const*types
  353.          *         MU_SENSOR_PROB=0.0 
  354.          * result is (e.g., COMPLEX_SIZE=1):
  355.          *    A(0),B(1),...,X(TYPES-1),A,B,...,X,...
  356.          */
  357.         for (j=0;j<COMPLEX_SIZE;j++)
  358.                     ++ip->sensor_specs[i].complex[i%types];
  359.  
  360.         /*
  361.          * alternatively we can use the following: 
  362.          * RANDOM INITIALIZATION OF SENSOR COMPLEXES
  363.          *    
  364.         for (j=0;j<COMPLEX_SIZE;j++)
  365.             ++ip->sensor_specs[i].complex[mrand(types)];
  366.          *
  367.          */
  368.  
  369.         ip->sensor_specs[i].orientation = s_orient[i];
  370.         ip->sensor_specs[i].range = s_range[i];
  371.  
  372.         /* location = ...; */
  373.     }
  374.     count = 0;
  375.     for (i=0;i<nmotors;i++)
  376.     {
  377.  
  378.         ip->motor_specs[i].system = m_type[i];
  379.  
  380.         size = motor_size(ip->motor_specs[i].system);
  381.         count += size;
  382.         for (j=0;j<*(layer_descp+nlayers-1);j++)
  383.             ip->motor_specs[i].mask[j] 
  384.                 = get_mask(count, size, j);
  385.         ip->motor_specs[i].power 
  386.             = motor_power(ip->motor_specs[i].system);
  387.         /* location, orientation = ...; */
  388.     }
  389.     /*
  390.      * update the world to insert the org. in its cell
  391.      */
  392.     ins_org(ip);
  393.         /*
  394.          * initialize net weights and biases
  395.          */
  396.     init_net(ip);
  397.        }
  398.  
  399.     /*
  400.      * get a new structure and stick it in 'next' element
  401.      */
  402.     ip->next = (struct indiv *)malloc(sizeof(struct indiv));
  403.     if (ip->next == NIL_POINTER)
  404.           {
  405.         fprintf(stderr, "malloc error on individual\n");
  406.         exit(-1);
  407.           }
  408.     /*
  409.      * put old guy's address the new one's structure so he can
  410.      * work backward
  411.      */
  412.     (ip->next)->last = ip;    
  413.     /*
  414.      * tells us the next guy not yet real
  415.      */    
  416.     (ip->next)->next = (struct indiv *)-1;    
  417.     return(ip->next);
  418. }
  419.  
  420.  
  421.  
  422. /*
  423.  * generate a children ip2 from a father ip1:
  424.  * genotype of father is used to generate offspring;
  425.  * for now, development is limited to a copy of the
  426.  * backup weights (genotype) onto the normal weights
  427.  * (phenotype) 
  428.  */
  429.  
  430. reproduce(ip1)
  431.  
  432.     struct  indiv *ip1;
  433.  
  434. {
  435.     struct  indiv *ip2;
  436.     float  **iiwp1;
  437.     float  **iibp1;
  438.     float  *iiw1;
  439.     float  *iib1;    
  440.     float  **wp2;
  441.     float  **bp2;
  442.         float  **iiwp2;
  443.         float  **iibp2;
  444.     float  *w2;
  445.     float  *b2;
  446.         float  *iiw2;
  447.         float  *iib2;
  448.     int    *l;
  449.     int    i,j,ii;
  450.  
  451.     /*
  452.      * get a new structure for the child 
  453.      */
  454.     ip2 = (struct indiv *)malloc(sizeof(struct indiv));
  455.     if (ip2 == NIL_POINTER)
  456.           {
  457.         fprintf(stderr, "malloc error on child\n");
  458.         exit(-1);
  459.           }
  460.     /*
  461.       * append child to end of population chain
  462.       */
  463.     (end_head->last)->next = ip2;
  464.     ip2->last = end_head->last;
  465.     ip2->next = end_head;
  466.     end_head->last = ip2;
  467.  
  468.     ip2->gutsize    = ip1->gutsize;
  469.     ip2->generation = ip1->generation + 1;
  470.     ip2->lifecycle  = 0;
  471.     ip2->id        = pop_size;
  472.     ip2->pid        = ip1->id;
  473.     ip2->error      = 0.0;
  474.     ip1->fit++;
  475.     ip2->fit        = 0;
  476.     ip1->energy    /= 2.0;
  477.     ip2->energy    = ip1->energy;
  478.     ip2->worldx     = mrand(x_dim);
  479.     ip2->worldy     = mrand(y_dim);
  480.     switch (mrand(4)){
  481.         case 0:
  482.             ip2->direction  = UP;
  483.             break;
  484.         case 1:
  485.             ip2->direction  = LEFT;
  486.             break;
  487.         case 2:
  488.             ip2->direction  = DOWN;
  489.             break;
  490.         case 3:
  491.             ip2->direction  = RIGHT;
  492.             break;
  493.     }
  494.         for (i=0;i<types;i++) ip2->gut[i] = 0;
  495.         for (i=0;i<nsensors;i++)
  496.         {
  497.                 ip2->sensor_specs[i].system = ip1->sensor_specs[i].system;
  498.                 for (j=0;j<*layer_descp;j++)
  499.                         ip2->sensor_specs[i].mask[j] = ip1->sensor_specs[i].mask[j];
  500.                 for (j=0;j<types;j++)
  501.                         ip2->sensor_specs[i].complex[j] = ip1->sensor_specs[i].complex[j];
  502.         ip2->sensor_specs[i].orientation = ip1->sensor_specs[i].orientation;
  503.         ip2->sensor_specs[i].range = ip1->sensor_specs[i].range;
  504.                 /* location = ...; */
  505.         }
  506.         for (i=0;i<nmotors;i++)
  507.         {
  508.                 ip2->motor_specs[i].system = ip1->motor_specs[i].system;
  509.                 for (j=0;j<*(layer_descp+nlayers-1);j++)
  510.                         ip2->motor_specs[i].mask[j] = ip1->motor_specs[i].mask[j];
  511.                 ip2->motor_specs[i].power = ip1->motor_specs[i].power;
  512.                 /* location, orientation = ...; */
  513.         }
  514.     /*
  515.      * update the world to insert the org. in its cell
  516.      */
  517.     ins_org(ip2);
  518.     /*
  519.      * update population size 
  520.      */ 
  521.     pop_size++;
  522.  
  523.     if (verbose >= 4)
  524.       printf("\n%d_%d generated %d_%d\n",ip1->generation,ip1->id,
  525.                                           ip2->generation,ip2->id);
  526.  
  527.     /*
  528.      * copy genotype from father to child  
  529.      * and mutate the child's genotype
  530.      */
  531.     build_geno_net(ip2);
  532.  
  533.     for (i=0,iibp1=ip1->iibiasp,iibp2=ip2->iibiasp,l=layer_descp; 
  534.         i<nlayers; i++,iibp1++,iibp2++,l++)
  535.  
  536.         for (ii=0,iib1=*iibp1,iib2=*iibp2; 
  537.             ii< *l; ii++,iib1++,iib2++)
  538.  
  539.             *iib2 = *iib1;
  540.  
  541.     for(i=0,iiwp1=ip1->iiweightp,iiwp2=ip2->iiweightp,l=layer_descp;
  542.         i<(nlayers-1); i++,iiwp1++,iiwp2++,l++)
  543.  
  544.         for (ii=0,iiw1=*iiwp1,iiw2=*iiwp2; 
  545.             ii< *l * *(l + 1); ii++,iiw1++,iiw2++)
  546.  
  547.             *iiw2 = *iiw1;
  548.  
  549.     mutate(ip2);
  550.  
  551.     /*
  552.      * develop phenotype from genotype:
  553.      * for now, the latter is made of the backup weights
  554.      * and development is simply a copy of these onto the 
  555.      * actual learnable weights;
  556.      * note that subsequent references to an individual 
  557.      * are to its phenotype! 
  558.      */
  559.     build_pheno_net(ip2);
  560.  
  561.     for (i=0,bp2=ip2->biasp,iibp2=ip2->iibiasp,l=layer_descp; 
  562.         i<nlayers; i++,bp2++,iibp2++,l++)
  563.  
  564.         for (ii=0,b2=*bp2,iib2=*iibp2; 
  565.             ii< *l; ii++,b2++,iib2++)
  566.  
  567.             *b2 = *iib2;
  568.  
  569.     for (i=0,wp2=ip2->weightp,iiwp2=ip2->iiweightp,l=layer_descp; 
  570.         i<(nlayers-1); i++,wp2++,iiwp2++,l++)
  571.  
  572.         for (ii=0,w2=*wp2,iiw2=*iiwp2; 
  573.             ii< *l * *(l + 1); ii++,w2++,iiw2++)
  574.  
  575.             *w2 = *iiw2;
  576. }
  577.  
  578.  
  579.  
  580. /*
  581.  * delete all of an indiv's data structures;
  582.  * delete him from linked list of amoeba structures;
  583.  * return pointer to the next guy
  584.  */
  585. struct indiv *
  586. delete_indiv(ip)
  587.  
  588.     struct    indiv *ip;
  589. {
  590.     struct indiv *next;
  591.  
  592.     if (ip->next == (struct indiv *) -1)
  593.     {
  594.         printf("error: deleting a non-individual\n");
  595.         exit(-1);
  596.     }
  597.     /*
  598.      * update world to delete org. from its cell
  599.      */
  600.     del_org(ip);
  601.     next = ip->next;
  602.     if (ip != a_head)
  603.     {
  604.         (ip->last)->next = ip->next;
  605.         (ip->next)->last = ip->last;
  606.     }
  607.     else
  608.     {
  609.         a_head = ip->next;
  610.         (ip->next)->last = (struct indiv *) 0;
  611.     }
  612.     destroy_network(ip);
  613.     free((char *) ip);
  614.     return(next);
  615. }
  616.  
  617.  
  618.  
  619. /*
  620.  * return the next indiv; if called with
  621.  * -1 initialize
  622.  */
  623.  
  624. struct indiv *
  625. getnext_indiv(flag)
  626.  
  627.     int    flag;
  628. {
  629.     static    struct indiv *ap;
  630.     if (flag == -1)
  631.              {
  632.            ap = a_head;
  633.            return(a_head);
  634.          }
  635.            else
  636.              {
  637.            if (ap->next == (struct indiv *)-1)
  638.                     {
  639.               return((struct indiv *) -1);
  640.             }
  641.                   else
  642.                     {
  643.               ap = ap->next;
  644.               return(ap);
  645.             }
  646.          }
  647. }
  648.  
  649.  
  650. /*
  651.  * return the amoeba whose id number is passed
  652.  */
  653.  
  654. struct indiv *
  655. get_indiv(id)
  656.  
  657.     int    id;
  658. {
  659.     struct    indiv *ap;
  660.  
  661.     ap = a_head;
  662.     while (ap->id != id)
  663.            {
  664.         if (ap->next == (struct indiv *) -1)
  665.                 {
  666.           printf("Unable to find indiv: %d\n", id);
  667.           exit(-1);
  668.         }
  669.             else
  670.           ap = ap->next; 
  671.        }
  672.     return(ap);
  673. }    
  674.  
  675.  
  676. /*
  677.  * mutate the genotype of an individual
  678.  * including sensor and motor systems
  679.  */
  680. mutate(ip)
  681.  
  682.     struct indiv *ip;
  683.  
  684. {
  685.  
  686.     float      **wp;
  687.     float      **bp;
  688.     float      *w;
  689.     float      *b;
  690.     int        *bldp;
  691.     int        m;
  692.     int        n1,n2,l;
  693.  
  694.  
  695.     for(m=0;m<mutations;m++)
  696.     {
  697.         l=mrand(nlayers-1);
  698.         n1 = mrand(*(layer_descp+l));
  699.         n2 = mrand(*(layer_descp+l+1));
  700.         wp = ip->iiweightp + l;
  701.         bldp = layer_descp + l;
  702.         w = *wp;
  703.  
  704.         if (verbose > 2)
  705.             printf("\nmutated layer %d from unit %d to %d",l,n1,n2);
  706.  
  707.         *(w + n1 + (n2 * *bldp)) += rans(mutations_range);
  708.         }
  709.  
  710.     if (mbiasestoo) 
  711.       for(m=0;m<bmutations;m++)
  712.     {
  713.         l=mrand(nlayers-1);
  714.         n1 = mrand(*(layer_descp+l+1));
  715.         bp = ip->iibiasp + l + 1;
  716.         b = *bp;
  717.  
  718.         if (verbose > 2)
  719.             printf("\nmutated bias of layer %d unit %d",l+1,n1);
  720.  
  721.         *(b + n1) += rans(mutations_range);
  722.     }
  723.  
  724.     mu_sensorymotor(ip);
  725. }
  726.  
  727.  
  728. /*
  729.  * ATTENTION: has been changed; no reproduction here,
  730.  *            moved into inner loop instead;
  731.  *    here we save the results for the current generation
  732.  */
  733. next_generation()
  734. {
  735.  
  736.     struct    indiv *ap;
  737.     struct    indiv *tap;
  738.     struct    indiv *best_ap;
  739.     int    flag;
  740.     float   v;
  741.     int     f;
  742.     int     maxbest;
  743.     char    genfile[64];
  744.     FILE    *fp;
  745.  
  746.     /*
  747.      * save energy info for each amoeba and store in .gen file;
  748.          */
  749.     if (verbose > 1)
  750.     {
  751.      sprintf(genfile, "%lu.gen", seed);
  752.      if ((fp=fopen(genfile, "a")) == NULL)
  753.      {
  754.         printf("error: opening file %s",genfile);
  755.         exit(-1);
  756.      }
  757.      flag = -1;
  758.      while ((ap = getnext_indiv(flag))->next != (struct indiv *) -1)
  759.      {
  760.         fprintf(fp,"%.2f ",ap->energy);
  761.         flag = 0;
  762.      }
  763.      fprintf(fp,"\n");
  764.      fclose(fp);
  765.     }
  766.  
  767.     /*
  768.      *  we save the best individuals
  769.      */
  770.     if (save_best > 0)
  771.     {
  772.      /*
  773.       * based on energy, sort to get the best amoebas 
  774.       * and save their indiv info in temporary amoeba struct.
  775.       */
  776.  
  777.     if (verbose >= 2)
  778.         printf("Best :");
  779.  
  780.     maxbest = (save_best < pop_size)? save_best : pop_size;
  781.     for(f=0;f<maxbest;f++)
  782.        {
  783.         flag = -1;
  784.         v = -99999.9;
  785.         while ((ap = getnext_indiv(flag))->next != (struct indiv *) -1)
  786.         {
  787.             if (ap->energy > v)
  788.             {
  789.                 v = ap->energy;
  790.                 best_ap = ap;
  791.             }
  792.             flag = 0;
  793.         }
  794.  
  795.         if (verbose >= 2)
  796.           printf("%d (%.2f) ",best_ap->id,best_ap->energy);
  797.  
  798.                 /*
  799.                  * we delete the best individual from
  800.          * the a_head chain
  801.          */
  802.                 if (best_ap == a_head)
  803.                 {
  804.                         a_head = best_ap->next;
  805.             a_head->last = (struct indiv *) 0;
  806.                 }
  807.                 else
  808.                 {
  809.                         (best_ap->last)->next = best_ap->next;
  810.                         (best_ap->next)->last = best_ap->last;
  811.                 }
  812.  
  813.         /*
  814.          * we add best_ap to t_head chain
  815.          */
  816.         if (f == 0)
  817.         {
  818.             t_head = best_ap;
  819.             tap = t_head;
  820.             tap->next = end_head;
  821.         }
  822.         else
  823.         {
  824.             tap->next = best_ap;
  825.             best_ap->last = tap;
  826.             tap = best_ap;
  827.             tap->next = end_head;
  828.         }
  829.        }   
  830.     
  831.     if (verbose >= 2) 
  832.         printf("\n");
  833.  
  834.  
  835.     if (verbose > 0)
  836.         save_best_indiv(maxbest);
  837.  
  838.  
  839.  
  840.  
  841.     /*
  842.       * re-append t_head to a_head 
  843.      */
  844.     if (end_head->last == (struct indiv *) 0)
  845.       {
  846.         a_head = t_head;
  847.         a_head->last = (struct indiv *) 0;
  848.       }
  849.     else
  850.       {
  851.         t_head->last = end_head->last;
  852.         (t_head->last)->next = t_head;
  853.       }
  854.         end_head->last = tap;
  855.     }
  856. }
  857.  
  858.  
  859.  
  860.  
  861. /*
  862.  * delete individual, both genotype and phenotype 
  863.  */
  864.  
  865. die(ap)
  866.  
  867.     struct indiv *ap;
  868.  
  869. {
  870.     int    i,j;
  871.  
  872.     /*
  873.      * update energy_reserve for checking energy conservation
  874.      * and drop gut content in current cell
  875.      */
  876.     if (ap->energy < 0.0) energy_reserve += ap->energy;
  877.     for (i=0;i<types;i++)
  878.         for (j=0;j<ap->gut[i];j++)
  879.             (void)ins_food(i,ap->worldx,ap->worldy,FALSE);
  880.  
  881.     /*
  882.      * update population size
  883.      */
  884.     pop_size--;
  885.     if (pop_size < 1)
  886.     {
  887.         printf("\nEXTINCTION OCCURRED!\n");
  888.         exit(1);
  889.     }
  890.  
  891.     /*
  892.      * remove the individual from memory;
  893.      * keep id's consecutive by changing the
  894.      * largest id to the id of the dead
  895.      */
  896.     (get_indiv(pop_size))->id = ap->id;
  897.     delete_indiv(ap);
  898. }
  899.  
  900.  
  901.