31a78c76bcee832020f25d297f55cceff971d77f
[afpopgen] / 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, "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                         default:
92                                   fprintf(stderr, "Usage: %s [-s seed] [-r run] [-y startyear] [-eh] [-p hex:freq,freq[:freq,freq...]\n",
93                                                   av[0]);
94                                   exit(EXIT_FAILURE);
95                                   break;
96                 }
97         }
98
99         var_init(seed);
100         
101         if (defaultpop) {
102                 hex = find_hex(map, 4788);
103                 hex->freq[0] = hex->freq[1] = 2000;
104                 hex->pop = hex->freq[0] + hex->freq[1];
105         }
106
107 #if 0
108         r = var_rhyper(1, 5, 5);
109         printf("rhg = %u\n", r);
110 #endif
111
112         /* dump out the coordinates of the hexes */
113 #if 0
114         printf("map init complete (%p)\n", map);
115         double hexvert[14];
116         int i;
117         printf("{\n");
118         for (hex = map; hex; hex = hex->next) {
119                 HL_vertices(hex->hex, hexvert);
120                 printf("\t{ \"hex\": %lld, \"coords\": [", hex->hex);
121                 for (i=0; i<5; i++) {
122                         printf("%f,%f, ", hexvert[i*2], hexvert[i*2+1]);
123                 }
124                 printf("%f,%f", hexvert[5*2], hexvert[5*2+1]);
125                 printf("] }\n");
126         }
127 #endif
128
129
130         printf("{");
131         printf(" \"run\": %d, \"seed\": %d, \"gens\": %d", run, seed, gens);
132         printf(",\n\"options\": { \"capfactor\": %f }", capfactor);
133         printf(",\n\"gendata\": [\n");
134
135         dumpgenjson(map, 0, year);
136
137         for (gen = 1; gen <= gens; gen++) {
138                 year += 25;
139                 nextgen(map, year, holocene + eemian > 0);
140                 printf(",\n");
141                 dumpgenjson(map, gen, year);
142         }
143         printf("\n]");
144
145         printf(",\n\"hexlist\": ");
146         hexlistjson(map);
147
148         printf("\n}\n");
149
150         return 0;
151 }
152
153 void dumpgenjson(struct afrhex *map, int gen, int year) {
154         struct afrhex *hex;
155         int count = 0;
156         unsigned int maxpop = 0;
157
158         printf("{ \"gen\": %d, \"year\": %d, \"pop\": %llu, \"hexes\": [\n",
159                         gen, year, total_pop(map)
160               );
161         for (hex = map; hex; hex = hex->next) {
162                 if (hex->pop == 0) continue;
163                 if (hex->pop > maxpop) {
164                         maxpop = hex->pop;
165                 }
166                 count++;
167                 printf("%s\t{ \"hex\": %lld, \"pop\": %u, \"freq\": [%u, %u], "
168                                 "\"climate\": \"%s\", \"cap\": %u, \"outpop\": [%u, %u], "
169                                 "\"inpop\": [%u, %u]"
170                                 " }",
171                                 count > 1 ? ",\n" : "",
172                                 hex->hex, hex->pop, hex->freq[0], hex->freq[1],
173                                 hex->climate->code, hex->climate->cap, 
174                                 hex->outmigrate[0], hex->outmigrate[1],
175                                 hex->inmigrate[0], hex->inmigrate[1]
176                       );
177         }
178         printf(" ],\n\t\"maxpop\": %u}", maxpop);
179 }
180
181 void hexlistjson(struct afrhex *map) {
182         struct afrhex *hex;
183         int count = 0;
184
185         printf("[");
186         for (hex = map; hex; hex = hex->next) {
187                 printf("%s%s%lld",
188                                 count > 0 ? "," : "",
189                                 count % 10 == 0 ? "\n" : "",
190                                 hex->hex
191                       );
192                 count++;
193         }
194         printf("]");
195 }
196
197