home *** CD-ROM | disk | FTP | other *** search
/ Collection of Hack-Phreak Scene Programs / cleanhpvac.zip / cleanhpvac / TIERRA40.ZIP / BEAGLE / DIVERSE.C < prev    next >
C/C++ Source or Header  |  1992-09-11  |  25KB  |  805 lines
  1.  
  2. /* diverse.c  15-8-92  code for computing the diversity index */
  3. /*** Diverse: Version 4.0  Copyright (c) 1990, 1991, 1992  Tom Ray ***/
  4.  
  5. /* Thresholding decisions:
  6.  
  7. 1) When the threshold is set at one cell, the behavior is the same
  8.    as if thresholds were not used.
  9. 2) Size classes are counted as above thereshold only when one of their
  10.    genotypes is above threshold.
  11. 3) The age of a class is based on when it first appeared, not when it
  12.    crossed the threshold.
  13. 4) When a class crosses a threshold, its entire population is counted,
  14.    not just the above threshold part of the population.
  15. */
  16.  
  17. #include <stdio.h>
  18. #include <stdlib.h>
  19. #include <string.h>
  20. #include <math.h>
  21. #include <time.h>
  22. #include <ctype.h>
  23. #ifdef __TURBOC__
  24. #include <alloc.h>
  25. #endif /* __TURBOC__ */
  26.  
  27. #ifdef __STDC__           /* ANSI prototyping */
  28. #define P_(A) A
  29. #define const const
  30. #else                 /* non-ANSI prototyping */
  31. #define P_(A) ()
  32. #define const
  33. #endif
  34.  
  35. typedef unsigned int   Uint;
  36. typedef unsigned long  Ulong;
  37. typedef long int       I32s;
  38.  
  39. struct siz {
  40.     long   Time;
  41.     long   NumCell;
  42.     long   NumSize;
  43.     float  SizeDiv;
  44.     long   AgeSize;
  45.     } ;
  46.  
  47. struct gen {
  48.     long   NumGeno;
  49.     float  GenoDiv;
  50.     long   AgeGeno;
  51.     } ;
  52.  
  53. struct cum {
  54.     double  Time;
  55.     double  NumCell;
  56.     double  NumSize;
  57.     double  SizeDiv;
  58.     double  AgeSize;
  59.     double  NumGeno;
  60.     double  GenoDiv;
  61.     double  AgeGeno;
  62.     } cu, icu;
  63.  
  64. struct last_out {
  65.     Ulong  time;   /* += itime % 1000000 */
  66.     Ulong  ctime;  /* count of millions of cpu cycles */
  67.     Ulong  itime;  /* time interval since last birth or death */
  68.     char   bd;
  69.     Uint   size;
  70.     char   label[4];
  71.     } lo;
  72.  
  73. struct gene_dat {
  74.     struct tnode  *t;
  75.     long          index;
  76.     } ;
  77.  
  78. struct pop_dat {
  79.     long  pop;  /* current population of this size or genotype */
  80.     long  popt; /* current threshold population of this size or genotype */
  81.     long  age;  /* current age of this size or genotype */
  82.     } ;
  83.  
  84. struct tnode {
  85.     long            size;  /* genome size */
  86.     int             ntg;   /* number of threshold genotypes */
  87.     struct pop_dat  sd;    /* pop_dat for this size class */
  88.     int             gsize; /* allocated size of *g array */
  89.     struct pop_dat  *g;    /* array of pop_dat structures */
  90.     struct tnode    *l;    /* left sub-tree */
  91.     struct tnode    *r;    /* right sub-tree */
  92.     } ;
  93.  
  94. struct Totals {
  95.     I32s    TotPopTh;
  96.     I32s    TotPopTo;
  97.     I32s    TotSizes;
  98.     I32s    TotGenos;
  99.     double  time;
  100.     } to;
  101.  
  102. struct DivDat {
  103.     float   SizeDiv;
  104.     float   GenoDiv;
  105.     I32s    NumSize;
  106.     I32s    NumGeno;
  107.     I32s    NumCellSTh;
  108.     I32s    NumCellSTo;
  109.     I32s    NumCellGTh;
  110.     I32s    NumCellGTo;
  111.     double  tot_size_age;
  112.     double  tot_geno_age;
  113.     } dd;
  114.  
  115. struct range {
  116.     double  n;  /* minimum value */
  117.     double  x;  /* maximum value */
  118.     I32s    m;  /* count of values */
  119.     int     f;  /* first time through = 1 */
  120.     } timer, sdivr, gdivr, nsr, ngr, ncr, asar, agar;
  121.  
  122. void main P_((int  argc, char  **argv));
  123. void doranges P_(());
  124. void doaverages P_(());
  125. void breakrange P_(());
  126. struct tnode * AddTree P_((struct tnode  *p));
  127. void divCountTree P_((struct tnode  *t));
  128. void minmax2 P_((double  v, struct range  *r));
  129. int Lbl2Int P_((char  s[]));
  130. int t_read P_((char  data[], struct last_out  *lo, int  *first,
  131.      int  *genotypes));
  132.  
  133. int     genotypes, thr;
  134. double  AgeSize, AgeGeno;
  135. char    infile[13], ofile[13], bifile[13], bofile[13], data[81];
  136. FILE    *inf, *ouf, *our, *iouf;
  137.  
  138. void main (argc, argv)
  139. int  argc;
  140. char  **argv;
  141. {   unsigned long  size, fsize = 0;
  142.     I32s           mtime = 0, otime = 0, iamtime = 0, iaotime = 0, AvgFreq;
  143.     int            binum = 1, bonum = 1;
  144.     int            first = 1, first2 = 1, first3 = 1, format = 0, bdrecs = 0;
  145.     char           c, buf[80];
  146.     struct tnode   *troot = NULL;
  147.     struct gen     ge;
  148.     struct siz     si;
  149.     double         tmp;
  150. /* for debugging:
  151.     Ulong          dtime = 0, a = 0, b = 1;
  152.  */
  153.  
  154.     timer.n = timer.x = sdivr.n = sdivr.x = gdivr.n = gdivr.x = 0.;
  155.     nsr.n = nsr.x = ngr.n = ngr.x = ncr.n = ncr.x = 0.;
  156.     asar.n = asar.x = agar.n = agar.x = 0.;
  157.     timer.f = sdivr.f = gdivr.f = nsr.f = ngr.f = ncr.f = asar.f = agar.f = 1;
  158.     timer.m = sdivr.m = gdivr.m = nsr.m = ngr.m = ncr.m = asar.m = agar.m = 0;
  159.     to.TotPopTh = to.TotPopTo = to.TotSizes = to.TotGenos = 0;
  160.     to.time = 0.;
  161.  
  162.     printf ("diverse tool: to accept default values, hit Enter\n\n");
  163.  
  164.     sprintf(infile,"break.1");
  165.     printf("input file (default: %s) = ", infile);
  166.     gets(buf);
  167.     if (strlen(buf))
  168.         sprintf(infile, "%s", buf);
  169.  
  170.     inf = fopen(infile,"r");
  171.     if(inf == NULL)
  172.     {   printf("file %s not found\n", infile);
  173.         return ;
  174.     }
  175.     fclose(inf);
  176.  
  177.     thr = 2;
  178.     printf("threshold (default: %d) = ", thr);
  179.     gets(buf);
  180.     if (strlen(buf))
  181.     {   sscanf(buf,"%d", &thr);
  182.     }
  183.  
  184.     AvgFreq = 1000000L;
  185.     printf("Frequency for Average Output (default: %ld) = ", AvgFreq);
  186.     gets(buf);
  187.     if (strlen(buf))
  188.     {   sscanf(buf,"%ld", &AvgFreq);
  189.     }
  190.  
  191.     bdrecs = 0;
  192.     if (bdrecs)
  193.         sprintf(buf,"y");
  194.     else
  195.         sprintf(buf,"n");
  196.     printf("output birth death records: y or n (default: %s) = ", buf);
  197.     gets(buf);
  198.     if (strlen(buf))
  199.     {   if(!strcmp(buf,"n"))
  200.             bdrecs = 0;
  201.         else bdrecs = 1;
  202.     }
  203.     if (bdrecs)
  204.     {   sprintf(buf,"divdat.1");
  205.         printf("output file (default: %s) = ", buf);
  206.         gets(buf);
  207.         if (strlen(buf))
  208.             sscanf(buf,"%[^.]", ofile);
  209.         else
  210.             sprintf(ofile,"divdat");
  211.         strcpy(bofile,ofile);
  212.         strcat(ofile,".1");
  213.  
  214.         sprintf(buf,"binary"); format = 0;
  215.         printf("output format: binary or ascii (default: %s) = ", buf);
  216.         gets(buf);
  217.         if (strlen(buf))
  218.         {   if(!strcmp(buf,"binary"))
  219.                 format = 0;
  220.             else format = 1;
  221.         }
  222.  
  223.         size = 4096;
  224.         printf("break size (default: %lu) = ", size);
  225.         gets(buf);
  226.         if (strlen(buf))
  227.         {   sscanf(buf,"%lu", &size);
  228.         }
  229.         size *= 1024;
  230.  
  231.         our = fopen("brkrange","w");
  232.     }
  233.  
  234.     sscanf(infile,"%[^.].%d", bifile, &binum);
  235.     inf = fopen(infile,"r");
  236.     if(inf == NULL)
  237.     {   printf("\n Input File %s NOT found. \n", infile);
  238.         exit(1);
  239.     }
  240.     lo.time = lo.ctime = lo.size = 0;
  241.     cu.Time    = 0.; cu.NumCell = 0.; cu.NumSize = 0.; cu.SizeDiv = 0.;
  242.     cu.AgeSize = 0.; cu.NumGeno = 0.; cu.GenoDiv = 0.; cu.AgeGeno = 0.;
  243.     icu.Time    = 0.; icu.NumCell = 0.; icu.NumSize = 0.; icu.SizeDiv = 0.;
  244.     icu.AgeSize = 0.; icu.NumGeno = 0.; icu.GenoDiv = 0.; icu.AgeGeno = 0.;
  245.     for(;;)
  246.     {   if(fgets(data,80,inf) == NULL)
  247.         {   binum++;
  248.             sprintf(infile,"%s.%d", bifile, binum);
  249.             fclose(inf);
  250.             inf = fopen(infile,"r");
  251.             if(inf == NULL)
  252.                 break ;
  253.             if(fgets(data,80,inf) == NULL)
  254.                 break ;
  255.         }
  256.         t_read(data, &lo, &first, &genotypes);
  257.  
  258. /* begin: for debugging purposes */
  259. /*      if (to.time > dtime)
  260.             a = b;
  261.  */
  262. /* end: for debugging purposes */
  263.  
  264.         if(first2)
  265.         {   first2 = 0;
  266. #ifdef __TURBOC__
  267.             printf("\nTime = %4ld million  Coreleft = %6lu\r",
  268.                 0L, coreleft());
  269. #else  /* __TURBOC__ */
  270.             printf("\nTime = %4d million\r", first2);
  271. #endif /* __TURBOC__ */
  272.             fflush(stdout);
  273.             if (bdrecs)
  274.             {   if(format)
  275.                 {   ouf = fopen(ofile,"w");
  276.                     if(genotypes)
  277.                         fprintf(ouf,"11 %.1lf\n", to.time);
  278.                     else
  279.                         fprintf(ouf,"10 %.1lf\n", to.time);
  280.                 }
  281.                 else
  282.                 {   ouf = fopen(ofile,"wb");
  283.                     c = (char) format;
  284.                     fwrite(&c,sizeof(char),1,ouf);
  285.                     c = (char) genotypes;
  286.                     fwrite(&c,sizeof(char),1,ouf);
  287.                     fwrite(&to.time,sizeof(double),1,ouf);
  288.                 }
  289.             }
  290.             if (AvgFreq)
  291.             {   iouf = fopen("averages","w");
  292.                 fprintf(iouf,
  293.                     "       Time NumCell NumSize SizeDiv    AgeSize");
  294.                 if(genotypes)
  295.                     fprintf(iouf, " NumGeno GenoDiv    AgeGeno");
  296.                 fprintf(iouf,"\n\n");
  297.             }
  298.         }
  299.         troot = AddTree(troot);
  300.         dd.NumCellSTh = 0;
  301.         dd.NumCellSTo = 0;
  302.         dd.NumSize = 0;
  303.         dd.SizeDiv = 0;
  304.         dd.tot_size_age = 0.;
  305.         if(genotypes)
  306.         {   dd.NumCellGTh = 0;
  307.             dd.NumCellGTo = 0;
  308.             dd.NumGeno = 0;
  309.             dd.GenoDiv = 0;
  310.             dd.tot_geno_age = 0.;
  311.         }
  312.  
  313.         divCountTree(troot);
  314.  
  315. #ifdef ERROR
  316.         if(dd.NumCellSTh != to.TotPopTh)
  317.             printf("dd.NumCellSTh = %ld  to.TotPopTh = %ld\n",
  318.                 dd.NumCellSTh, to.TotPopTh);
  319.         if(dd.NumCellSTo != to.TotPopTo)
  320.             printf("dd.NumCellSTo = %ld  to.TotPopTo = %ld\n",
  321.                 dd.NumCellSTo, to.TotPopTo);
  322.         if(dd.NumSize != to.TotSizes)
  323.             printf("dd.NumSize = %ld  to.TotSizes = %ld\n",
  324.                 dd.NumSize, to.TotSizes);
  325.         if(genotypes)
  326.         {   if(dd.NumGeno != to.TotGenos)
  327.                 printf("dd.NumGeno = %ld  to.TotGenos = %ld\n",
  328.                     dd.NumGeno, to.TotGenos);
  329.             if(dd.NumCellGTh != dd.NumCellSTh)
  330.                 printf("dd.NumCellGTh = %ld  dd.NumCellSTh = %ld\n",
  331.                     dd.NumCellGTh, dd.NumCellSTh);
  332.             if(dd.NumCellGTo != dd.NumCellSTo)
  333.                 printf("dd.NumCellGTo = %ld  dd.NumCellSTo = %ld\n",
  334.                     dd.NumCellGTo, dd.NumCellSTo);
  335.         }
  336. #endif /* ERROR */
  337.  
  338.         if ((!genotypes && to.TotSizes) || (genotypes && to.TotGenos))
  339.         {   AgeSize = dd.tot_size_age / (double) to.TotSizes;
  340.             if (genotypes)
  341.                 AgeGeno = dd.tot_geno_age / (double) to.TotGenos;
  342.             if (bdrecs)
  343.             {   if (format)
  344.                 {   if (genotypes)
  345.                         fsize += 1 +
  346.                             fprintf(ouf,"%lx %ld %ld %g %lx %ld %g %lx\n",
  347.                         lo.itime, to.TotPopTh, dd.NumSize, dd.SizeDiv,
  348.                         (Ulong) AgeSize, dd.NumGeno, dd.GenoDiv,
  349.                         (Ulong) AgeGeno);
  350.                     else fsize += 1 + fprintf(ouf,"%lx %ld %ld %g %lx\n",
  351.                         lo.itime, to.TotPopTh, dd.NumSize, dd.SizeDiv,
  352.                         (Ulong) AgeSize);
  353.                 }
  354.                 else
  355.                 {   si.Time = lo.itime;
  356.                     si.NumCell = to.TotPopTh;
  357.                     si.NumSize = dd.NumSize;
  358.                     si.SizeDiv = dd.SizeDiv;
  359.                     si.AgeSize = (Ulong) AgeSize;
  360.                     fwrite(&si, sizeof(struct siz), 1, ouf);
  361.                     fsize += sizeof(struct siz);
  362.                     if(genotypes)
  363.                     {   ge.NumGeno = dd.NumGeno;
  364.                         ge.GenoDiv = dd.GenoDiv;
  365.                         ge.AgeGeno = (Ulong) AgeGeno;
  366.                         fwrite(&ge, sizeof(struct gen), 1, ouf);
  367.                         fsize += sizeof(struct gen);
  368.                     }
  369.                 }
  370.             }
  371.         }
  372.         tmp = (double) lo.itime;
  373.         cu.Time    += tmp; icu.Time    += tmp;
  374.         tmp = (double) lo.itime * to.TotPopTh;
  375.         cu.NumCell += tmp; icu.NumCell += tmp;
  376.         tmp = (double) lo.itime * dd.NumSize;
  377.         cu.NumSize += tmp; icu.NumSize += tmp;
  378.         tmp = (double) lo.itime * dd.SizeDiv;
  379.         cu.SizeDiv += tmp; icu.SizeDiv += tmp;
  380.         tmp = (double) lo.itime * AgeSize;
  381.         cu.AgeSize += tmp; icu.AgeSize += tmp;
  382.         if (genotypes)
  383.         {   tmp = (double) lo.itime * dd.NumGeno;
  384.             cu.NumGeno += tmp; icu.NumGeno += tmp;
  385.             tmp = (double) lo.itime * dd.GenoDiv;
  386.             cu.GenoDiv += tmp; icu.GenoDiv += tmp;
  387.             tmp = (double) lo.itime * AgeGeno;
  388.             cu.AgeGeno += tmp; icu.AgeGeno += tmp;
  389.         }
  390.         if (bdrecs && (fsize > size || first3))
  391.         {   first3 = 0;
  392.             breakrange();
  393.         }
  394.         if (bdrecs && fsize > size)
  395.         {   fsize = 0;
  396.             bonum++;
  397.             fclose(ouf);
  398.             sprintf(ofile,"%s.%d", bofile, bonum);
  399.             if(format)
  400.             {   ouf = fopen(ofile,"w");
  401.                 if(genotypes)
  402.                     fprintf(ouf,"11 %.1lf\n", to.time);
  403.                 else
  404.                     fprintf(ouf,"10 %.1lf\n", to.time);
  405.             }
  406.             else
  407.             {   ouf = fopen(ofile,"wb");
  408.                 c = (char) format;
  409.                 fwrite(&c,sizeof(char),1,ouf);
  410.                 c = (char) genotypes;
  411.                 fwrite(&c,sizeof(char),1,ouf);
  412.                 fwrite(&to.time,sizeof(double),1,ouf);
  413.             }
  414.         }
  415.  
  416.         to.time += (double) lo.itime;
  417.         mtime = (I32s) (to.time / 1000000L);
  418.         if (mtime > otime)
  419.         {
  420. #ifdef __TURBOC__
  421.             printf("Time = %4ld million  Coreleft = %6lu\r",
  422.                 mtime, coreleft());
  423. #else  /* __TURBOC__ */
  424.             printf("Time = %4ld million\r", mtime);
  425. #endif /* __TURBOC__ */
  426.             fflush(stdout);
  427.             otime = mtime;
  428.         }
  429.         minmax2((double) to.time,&timer);
  430.         minmax2((double) to.TotPopTh,&ncr);
  431.         minmax2((double) dd.NumSize,&nsr);
  432.         minmax2((double) dd.SizeDiv,&sdivr);
  433.         minmax2((double) AgeSize,&asar);
  434.         if(genotypes)
  435.         {   minmax2((double) dd.GenoDiv,&gdivr);
  436.             minmax2((double) dd.NumGeno,&ngr);
  437.             minmax2((double) AgeGeno,&agar);
  438.         }
  439.         if (AvgFreq)
  440.         {   iamtime = (I32s) (to.time / AvgFreq);
  441.             if (iamtime > iaotime)
  442.             {   doaverages();
  443.                 icu.Time    = 0.; icu.NumCell = 0.; icu.NumSize = 0.;
  444.                 icu.SizeDiv = 0.; icu.AgeSize = 0.; icu.NumGeno = 0.;
  445.                 icu.GenoDiv = 0.; icu.AgeGeno = 0.;
  446.                 iaotime = iamtime;
  447.             }
  448.         }
  449.     }
  450.     doranges();
  451.     fclose(our);
  452.     if (bdrecs)
  453.     {   breakrange();
  454.         fclose(ouf);
  455.     }
  456.     if (AvgFreq)
  457.         fclose(iouf);
  458. }
  459.  
  460. void doaverages()
  461. {   fprintf(iouf,"%11.0lf %7.1lf %7.1lf %7.2lf %10.0lf",
  462.         to.time, icu.NumCell / icu.Time, icu.NumSize / icu.Time,
  463.         icu.SizeDiv / icu.Time, icu.AgeSize / icu.Time);
  464.     if(genotypes)
  465.     {   fprintf(iouf," %7.1lf %7.2lf %10.0lf",
  466.             icu.NumGeno / icu.Time, icu.GenoDiv / icu.Time,
  467.             icu.AgeGeno / icu.Time);
  468.     }
  469.     fprintf(iouf,"\n");
  470. }
  471.  
  472. void breakrange()
  473. {   if(genotypes)
  474.         fprintf(our,"%12s %.0lf %.0lf %.0lf %g %ld %.0lf %g %ld\n",
  475.             ofile, timer.x, ncr.x, nsr.x,
  476.             sdivr.x, (Ulong) asar.x, ngr.x,
  477.             gdivr.x, (Ulong) agar.x);
  478.     else
  479.         fprintf(our,"%12s %.0lf %.0lf %.0lf %g %ld\n",
  480.             ofile, timer.x, ncr.x, nsr.x,
  481.             sdivr.x, (Ulong) asar.x);
  482. }
  483.  
  484. void doranges()
  485. {   FILE  *ouf;
  486.  
  487.     ouf = fopen("divrange","w");
  488.     fprintf(ouf, "           Minimum      Maximum      Average Number\n\n");
  489.     fprintf(ouf,"Time    %10.0lf %12.0lf              %ld\n",
  490.         timer.n, timer.x, timer.m);
  491.     fprintf(ouf,"NumCell %10.0lf %12.0lf %12.1lf %ld\n",
  492.         ncr.n,   ncr.x,   cu.NumCell / cu.Time, ncr.m);
  493.     fprintf(ouf,"NumSize %10.0lf %12.0lf %12.1lf %ld\n",
  494.         nsr.n,   nsr.x,   cu.NumSize / cu.Time, nsr.m);
  495.     fprintf(ouf,"SizeDiv %10lf %12lf %12lf %ld\n",
  496.         sdivr.n, sdivr.x, cu.SizeDiv / cu.Time, sdivr.m);
  497.     fprintf(ouf,"AgeSize %10.0lf %12.0lf %12.1lf %ld\n",
  498.         asar.n,  asar.x,  cu.AgeSize / cu.Time, asar.m);
  499.     if(genotypes)
  500.     {   fprintf(ouf,"NumGeno %10.0lf %12.0lf %12.1lf %ld\n",
  501.             ngr.n,   ngr.x,   cu.NumGeno / cu.Time, ngr.m);
  502.         fprintf(ouf,"GenoDiv %10lf %12lf %12lf %ld\n",
  503.             gdivr.n, gdivr.x, cu.GenoDiv / cu.Time, gdivr.m);
  504.         fprintf(ouf,"AgeGeno %10.0lf %12.0lf %12.1lf %ld\n",
  505.             agar.n,  agar.x,  cu.AgeGeno / cu.Time, agar.m);
  506.     }
  507.     fclose(ouf);
  508. }
  509.  
  510. struct tnode * AddTree(p)
  511. struct tnode  *p;
  512. {   int  i, j, osize;
  513.     struct pop_dat  *pd;
  514.  
  515.     if(p == NULL) /* this is a new size class */
  516.     {   if(lo.bd == 'd')
  517.         {   printf("new node is a death, exiting\n");
  518.             doranges();
  519.             exit(0);
  520.         }
  521.         p = (struct tnode  *) calloc(1,sizeof(struct tnode));
  522.         if(p == NULL)
  523.         {   printf("calloc failure, exiting\n");
  524.             doranges();
  525.             exit(0);
  526.         }
  527.         p->size = lo.size;         /* initialize this size */
  528.         p->sd.pop = 1;
  529.         p->sd.age = (-lo.itime);
  530.         p->sd.popt = 0;
  531.         to.TotPopTo++;
  532.         if(genotypes)               /* allocate genotype array */
  533.         {   i = Lbl2Int(lo.label); /* & initialize first genotype */
  534.             p->gsize = i + 1;
  535.             p->g = (struct pop_dat  *)
  536.                 calloc(p->gsize,sizeof(struct pop_dat));
  537.             if(p->g == NULL)
  538.             {   printf("calloc failure, exiting\n");
  539.                 doranges();
  540.                 exit(0);
  541.             }
  542.             pd = p->g + i;
  543.             pd->age = (-lo.itime);
  544.             pd->pop = 1;
  545.             pd->popt = 0;
  546.             if (pd->pop == thr)
  547.             {   to.TotGenos++;
  548.                 if (!p->ntg)
  549.                     to.TotSizes++;
  550.                 p->ntg++;
  551.                 to.TotPopTh++;
  552.             }
  553.         }
  554.         else if (p->sd.pop == thr)
  555.         {   to.TotSizes++;
  556.             to.TotPopTh++;
  557.             p->ntg++;
  558.         }
  559.         if (p->ntg)
  560.         {   p->sd.popt++;
  561.             pd->popt++;
  562.         }
  563.         p->l = p->r = NULL;
  564.     }
  565.     else if(lo.size < p->size) /* if not this size */
  566.         p->l = AddTree(p->l);
  567.     else if(lo.size > p->size)
  568.         p->r = AddTree(p->r);
  569.     else                        /* this size, but not new */
  570.     {   if(genotypes) /* reallocating gsize array */
  571.         {   i = Lbl2Int(lo.label);
  572.             if(i >= p->gsize)
  573.             {   osize = p->gsize;
  574.                 p->gsize = i + 5;
  575.                 if(p->g == NULL)
  576.                     p->g = (struct pop_dat  *)
  577.                         calloc(p->gsize, sizeof(struct pop_dat));
  578.                 else
  579.                     p->g = (struct pop_dat  *)
  580.                         realloc(p->g, p->gsize * sizeof(struct pop_dat));
  581.                 if(p->g == NULL)
  582.                 {   printf("realloc failure, exiting\n");
  583.                     doranges();
  584.                     exit(0);
  585.                 }
  586.                 for(j = osize; j < p->gsize; j++)
  587.                 {   pd = p->g + j;
  588.                     pd->pop = pd->age = 0;
  589.                 }
  590.             }
  591.         }
  592.         if(lo.bd == 'b') /* this is a birth */
  593.         {   if(!p->sd.pop)
  594.                 p->sd.age = (-lo.itime);
  595.             p->sd.pop++;
  596.             to.TotPopTo++;
  597.             if(genotypes)
  598.             {   pd = p->g + i;
  599.                 if(!pd->pop)
  600.                     pd->age = (-lo.itime);
  601.                 pd->pop++;
  602.                 if (pd->pop == thr)
  603.                 {   to.TotGenos++;
  604.                     if (!p->ntg)
  605.                         to.TotSizes++;
  606.                     to.TotPopTh += pd->pop;
  607.                     pd->popt = pd->pop;
  608.                     p->sd.popt += pd->pop;
  609.                     p->ntg++;
  610.                 }
  611.                 else if (pd->pop > thr)
  612.                 {   to.TotPopTh++;
  613.                     pd->popt++;
  614.                     p->sd.popt++;
  615.                 }
  616.             }
  617.             else if (p->sd.pop == thr)
  618.             {   to.TotSizes++;
  619.                 to.TotPopTh += p->sd.pop;
  620.                 p->sd.popt = p->sd.pop;
  621.                 p->ntg++;
  622.             }
  623.             else if (p->sd.pop > thr)
  624.             {   to.TotPopTh++;
  625.                 p->sd.popt++;
  626.             }
  627.         }
  628.         else /* this is a death */
  629.         {   p->sd.pop--;
  630.             to.TotPopTo--;
  631.             if(!p->sd.pop)
  632.                 p->sd.age = 0;
  633.             if(genotypes)
  634.             {   pd = p->g + i;
  635.                 pd->pop--;
  636.                 if(!pd->pop)
  637.                     pd->age = 0;
  638.                 if (pd->pop == thr - 1)
  639.                 {   to.TotGenos--;
  640.                     p->ntg--;
  641.                     if (!p->ntg)
  642.                         to.TotSizes--;
  643.                     to.TotPopTh -= thr;
  644.                     p->sd.popt -= thr;
  645.                     pd->popt = 0;
  646.                 }
  647.                 else if (pd->pop >= thr)
  648.                 {   to.TotPopTh--;
  649.                     pd->popt--;
  650.                     p->sd.popt--;
  651.                 }
  652.             }
  653.             else if (p->sd.pop == thr - 1)
  654.             {   to.TotSizes--;
  655.                 to.TotPopTh -= pd->pop;
  656.                 p->sd.popt = 0;
  657.                 p->ntg--;
  658.             }
  659.             else if (p->sd.pop >= thr)
  660.             {   to.TotPopTh--;
  661.                 p->sd.popt--;
  662.             }
  663.         }
  664.     }
  665.     return p;
  666. }
  667.  
  668. void divCountTree(t)
  669. struct tnode  *t;
  670. {   struct pop_dat  *pd;
  671.     int i;
  672.  
  673.     if(t != NULL)
  674.     {   if(t->sd.pop > 0)
  675.         {   t->sd.age += lo.itime;
  676.             dd.NumCellSTo += t->sd.pop;
  677.             if (t->ntg)
  678.             {   dd.NumSize++;
  679.                 dd.tot_size_age += (double) t->sd.age;
  680.                 dd.NumCellSTh += t->sd.popt;
  681.                 dd.SizeDiv -= ((float) t->sd.popt / (float) to.TotPopTh)
  682.                     * (log((float) t->sd.popt / (float) to.TotPopTh));
  683.             }
  684.             if(genotypes) for(i = 0; i < t->gsize; i++)
  685.             {   pd = t->g + i;
  686.                 if(pd->pop > 0)
  687.                 {   pd->age += lo.itime;
  688.                     dd.NumCellGTo += pd->pop;
  689.                     if (pd->pop >= thr)
  690.                     {   dd.GenoDiv -= ((float) pd->popt / (float) to.TotPopTh)
  691.                             * (log((float) pd->popt / (float) to.TotPopTh));
  692.                         dd.NumGeno++;
  693.                         dd.tot_geno_age += (double) pd->age;
  694.                         dd.NumCellGTh += pd->popt;
  695.                     }
  696.                 }
  697.             }
  698.         }
  699.         divCountTree(t->l);
  700.         divCountTree(t->r);
  701.     }
  702. }
  703.  
  704. void minmax2(v, r)
  705. double  v;
  706. struct range  *r;
  707. {   if(r->f)
  708.     {   r->f = 0;
  709.         r->n = r->x = v;
  710.     }
  711.     if(v < r->n) r->n = v; if(v > r->x) r->x = v;
  712.     r->m++;
  713. }
  714.  
  715. int Lbl2Int(s)
  716. char  s[];
  717. {   if(s[0] == '-') return 0;
  718.     return 1 + (s[2]- 'a') + (26 * (s[1] - 'a')) + (676 * (s[0] - 'a'));
  719. }
  720.  
  721. int t_read(data, lo, first, genotypes)
  722. char  data[];
  723. struct last_out  *lo;
  724. int  *first;
  725. int  *genotypes;
  726. {   struct last_out  ti;
  727.     int    nargs;
  728.     char   v2[9], v3[9], v4[9];
  729.  
  730.     sscanf(data,"%s", v2);
  731.     if(!strcmp(v2,"num_sp")) return 0;
  732.     nargs = sscanf(data,"%lx%s%s%s", &ti.time, v2, v3, v4);
  733.     lo->itime = ti.time;
  734.     if(*first)
  735.     {   *first = 0;
  736.         if(nargs == 4) *genotypes = 1;
  737.         else *genotypes = 0;
  738.         lo->time += ti.time;     /* assumes lo structure initialized to zero */
  739.         if(lo->time >= 1000000L)
  740.         {   lo->time %= 1000000L;
  741.             lo->ctime++;
  742.         }
  743.         lo->bd = v2[0];
  744.         sscanf(v3,"%u", &lo->size);
  745.         if(*genotypes) strcpy(lo->label,v4);
  746.         else strcpy(lo->label,"");
  747.     }
  748.     else
  749.     {   lo->time += ti.time;
  750.         if(lo->time >= 1000000L)
  751.         {   lo->time %= 1000000L;
  752.             lo->ctime++;
  753.         }
  754.         if(*genotypes) switch(nargs)
  755.         {   case 1: break;
  756.             case 2:
  757.             {   if(isdigit(v2[0]))
  758.                 {   sscanf(v2,"%u", &lo->size); break; }
  759.                 else
  760.                 {   if(strlen(v2) == 1)
  761.                     {   lo->bd = v2[0]; break; }
  762.                     else
  763.                     {   strcpy(lo->label,v2); break; }
  764.                 }
  765.             }
  766.             case 3:
  767.             {   if(isdigit(v2[0]))
  768.                 {   sscanf(v2,"%u", &lo->size);
  769.                     strcpy(lo->label,v3);
  770.                 }
  771.                 else
  772.                 {   lo->bd = v2[0];
  773.                     if(isdigit(v3[0]))
  774.                         sscanf(v3,"%u", &lo->size);
  775.                     else
  776.                         strcpy(lo->label,v3);
  777.                 }
  778.                 break;
  779.             }
  780.             case 4:
  781.             {   lo->bd = v2[0];
  782.                 sscanf(v3,"%u", &lo->size);
  783.                 strcpy(lo->label,v4);
  784.                 break;
  785.             }
  786.         }
  787.         else switch(nargs)
  788.         {   case 1: break;
  789.             case 2:
  790.             {   if(isdigit(v2[0]))
  791.                     sscanf(v2,"%u", &lo->size);
  792.                 else
  793.                     lo->bd = v2[0];
  794.                 break;
  795.             }
  796.             case 3:
  797.             {   lo->bd = v2[0];
  798.                 sscanf(v3,"%u", &lo->size);
  799.                 break;
  800.             }
  801.         }
  802.     }
  803.     return 1;
  804. }
  805.