pd.if.org Git - afpopgen/blob - sim/sim.c

   1 #define _POSIX_C_SOURCE 200112L
   2 #include <stdio.h>
   3 #include <stdlib.h>
   4 #include <string.h>
   5 #include <time.h>
   6 #include <errno.h>
   7 #include <unistd.h>
   8
   9 #include "africa.h"
  10 #include "variate.h"
  11 #include "hexagon.h"
  12 void dumpgenjson(struct afrhex *map, int gen, int year);
  13 void hexlistjson(struct afrhex *map);
  14
  15 static void setpop(struct afrhex *map, int hexid, int *pop) {
  16         struct afrhex *hex;
  17         int i;
  18         int sum = 0;
  19
  20         hex = find_hex(map, hexid);
  21         for (i=0; i < 2; i++) {
  22                 hex->freq[i] = pop[i];
  23                 sum += pop[i];
  24         }
  25         if (sum % 2) {
  26                 fprintf(stderr, "allele frequency for hex %d not even\n", hexid);
  27         }
  28         hex->pop = sum/2;
  29 }
  30
  31 static void setpopstr(struct afrhex *map, char *spec) {
  32         char *resume;
  33         char *start;
  34         long hex;
  35         int freq[2];
  36         int i;
  37
  38         errno = 0;
  39         hex = strtol(spec, &resume, 10);
  40         if (errno) goto error;
  41
  42         for (i = 0; i < 2; i++) {
  43                 if (*resume != ':' && *resume != ',') goto error;
  44                 start = resume+1; /* skip the colon or comma */
  45                 errno = 0;
  46                 freq[i] = strtol(start, &resume, 10);
  47                 if (errno) goto error;
  48         }
  49         setpop(map, hex, freq);
  50         return;
  51
  52 error:
  53         fprintf(stderr, "unable to parse popspec '%s'\n", spec);
  54         fprintf(stderr, "format is hex:freq,freq\n");
  55         exit(EXIT_FAILURE);
  56 }
  57
  58 int main (int ac, char *av[]) {
  59         struct afrhex *map;
  60         struct afrhex *hex;
  61         int gen = 0;
  62         int opt;
  63         struct timespec ts;
  64         /* defaults below */
  65         int year = -140000;
  66         unsigned int seed = 0;
  67         int run = 1;
  68         int eemian = 0;
  69         int holocene = 0;
  70         int defaultpop = 1;
  71         int gens = 1000;
  72         double capfactor = 0.25;
  73         char tag[16] = "";
  74
  75         map = map_init(maphexes, 100.0, capfactor);
  76
  77         clock_gettime(CLOCK_REALTIME, &ts);
  78         seed = ((unsigned long)(ts.tv_sec) ^ (unsigned long)ts.tv_nsec) & 0xffffffff;
  79
  80         while ((opt = getopt(ac, av, "f:r:s:ehp:y:g:m:t:d:")) != -1) {
  81                 switch (opt) {
  82                         case 's': seed = atoi(optarg); break;
  83                         case 'r': run = atoi(optarg); break;
  84                         case 'y': year = atoi(optarg); break;
  85                         case 'g': gens = atoi(optarg); break;
  86                         case 'm': capfactor = strtod(optarg, 0); climate_init(100.0, capfactor); break;
  87                         case 't': strcpy(tag, optarg);
  88                         case 'e': eemian = 1; break;
  89                         case 'h': holocene = 1; break;
  90                         case 'p': setpopstr(map, optarg); defaultpop = 0;break;
  91                         case 'f': if (freopen(optarg, "w", stdout) == NULL) {
  92                                           perror("sim");
  93                                           exit(EXIT_FAILURE);
  94                                   }
  95                                   break;
  96                         default:
  97                                   fprintf(stderr, "Usage: %s [-s seed] [-r run] [-y startyear] [-eh] [-p hex:freq,freq[:freq,freq...]\n",
  98                                                   av[0]);
  99                                   exit(EXIT_FAILURE);
 100                                   break;
 101                 }
 102         }
 103
 104         var_init(seed);
 105
 106         if (defaultpop) {
 107                 hex = find_hex(map, 4788);
 108                 hex->freq[0] = hex->freq[1] = 2000;
 109                 hex->pop = hex->freq[0] + hex->freq[1];
 110         }
 111
 112 #if 0
 113         r = var_rhyper(1, 5, 5);
 114         printf("rhg = %u\n", r);
 115 #endif
 116
 117         /* dump out the coordinates of the hexes */
 118 #if 0
 119         printf("map init complete (%p)\n", map);
 120         double hexvert[14];
 121         int i;
 122         printf("{\n");
 123         for (hex = map; hex; hex = hex->next) {
 124                 HL_vertices(hex->hex, hexvert);
 125                 printf("\t{ \"hex\": %lld, \"coords\": [", hex->hex);
 126                 for (i=0; i<5; i++) {
 127                         printf("%f,%f, ", hexvert[i*2], hexvert[i*2+1]);
 128                 }
 129                 printf("%f,%f", hexvert[5*2], hexvert[5*2+1]);
 130                 printf("] }\n");
 131         }
 132 #endif
 133
 134
 135         printf("{");
 136         printf(" \"run\": %d, \"seed\": %d, \"gens\": %d", run, seed, gens);
 137         printf(",\n\"options\": { \"capfactor\": %g }", capfactor);
 138         printf(",\n\"gendata\": [\n");
 139
 140         dumpgenjson(map, 0, year);
 141
 142         for (gen = 1; gen <= gens; gen++) {
 143                 year += 25;
 144                 nextgen(map, year, holocene + eemian > 0);
 145                 printf(",\n");
 146                 dumpgenjson(map, gen, year);
 147         }
 148         printf("\n]");
 149
 150         printf(",\n\"hexlist\": ");
 151         hexlistjson(map);
 152
 153         printf("\n}\n");
 154
 155         return 0;
 156 }
 157
 158 void dumpgenjson(struct afrhex *map, int gen, int year) {
 159         struct afrhex *hex;
 160         int count = 0;
 161         unsigned int maxpop = 0;
 162
 163         printf("{ \"gen\": %d, \"year\": %d, \"pop\": %llu, \"hexes\": [\n",
 164                         gen, year, total_pop(map)
 165               );
 166         for (hex = map; hex; hex = hex->next) {
 167                 if (hex->pop == 0) continue;
 168                 if (hex->pop > maxpop) {
 169                         maxpop = hex->pop;
 170                 }
 171                 count++;
 172                 printf("%s\t{ \"hex\": %lld, \"pop\": %u, \"freq\": [%u, %u], "
 173                                 "\"climate\": \"%s\", \"cap\": %u, \"outpop\": [%u, %u], "
 174                                 "\"inpop\": [%u, %u]"
 175                                 " }",
 176                                 count > 1 ? ",\n" : "",
 177                                 hex->hex, hex->pop, hex->freq[0], hex->freq[1],
 178                                 hex->climate->code, hex->climate->cap,
 179                                 hex->outmigrate[0], hex->outmigrate[1],
 180                                 hex->inmigrate[0], hex->inmigrate[1]
 181                       );
 182         }
 183         printf(" ],\n\t\"maxpop\": %u}", maxpop);
 184 }
 185
 186 void hexlistjson(struct afrhex *map) {
 187         struct afrhex *hex;
 188         int count = 0;
 189
 190         printf("[");
 191         for (hex = map; hex; hex = hex->next) {
 192                 printf("%s%s%lld",
 193                                 count > 0 ? "," : "",
 194                                 count % 10 == 0 ? "\n" : "",
 195                                 hex->hex
 196                       );
 197                 count++;
 198         }
 199         printf("]");
 200 }
 201
 202