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

   1 #include <stdint.h>
   2 #include <stdlib.h>
   3 #include <string.h>
   4 #include <math.h>
   5 #include <stdio.h>
   6
   7 #ifndef LOCI
   8 #define LOCI 1
   9 #endif
  10
  11 #include "africa.h"
  12 #include "variate.h"
  13
  14 static struct climate climates[] = {
  15         { "Af", "Tropical Rain Forest", 3, 0, { 161,26,20 } },
  16         { "Am", "Tropical Monsoon", 45, 0, { 255,33,26 } },
  17         { "As", "Dry Savanna", 50, 0, { 255,163,150 } },
  18         { "Aw", "Wet Savanna", 125, 0, { 255,205, 199 } },
  19         { "Awf", "Wet Savanna forest", 100, 0, { 255,205, 199 } },
  20         { "BSh", "Hot Steppe", 12, 0, { 211,139,27 } },
  21         { "BSk", "Cold Steppe", 12, 0, { 207,169,81 } },
  22         { "BWk", "Cold Desert", 3, 0, { 254,245,91 } },
  23         { "BWh", "Hot Desert", 3, 0, { 255,197,13 } },
  24         { "Cfa", "Warm Temperate", 125, 0, { 14,43,13 } },
  25         { "Cfb", "Maritime", 50, 0, { 21,68,22 } },
  26         { "Csa", "Med", 50, 0, { 17,214,16 } },
  27         { "Csb", "Med", 50, 0, { 140,230,18 } },
  28         { "Cwa", "Warm Termperate", 125, 0, { 184,103,24 } },
  29         { "Cwb", "Maritime", 125, 0, { 154,100,25 } },
  30         { "Lf", "Lake", 0, 0, { 0,0,245 } },
  31         { {0} }
  32 };
  33
  34 struct alt_climate_search {
  35         int hex, year;
  36 };
  37
  38 static int alt_climate_compare(const void *key, const void *climate) {
  39         const struct climaterange *v;
  40         const struct alt_climate_search *k;
  41
  42         v = climate;
  43         k = key;
  44
  45         if (v->hex > k->hex) return -1; /* later hex */
  46         if (v->hex < k->hex) return 1; /* earlier hex */
  47         if (v->end < k->year) return 1; /* ends before year */
  48         if (v->start > k->year) return -1; /* hasn't started yet */
  49
  50         return 0;
  51 }
  52
  53 static int compare_cr(const void *av, const void *bv) {
  54         const struct climaterange *a, *b;
  55         a = av;
  56         b = bv;
  57         if (a->hex < b->hex) return -1;
  58         if (a->hex > b->hex) return 1;
  59         if (a->start < b->start) return -1;
  60         if (a->start > b->start) return 1;
  61         return 0;
  62 }
  63
  64 struct climate *find_climate(char *code) {
  65         int i;
  66         for (i=0; climates[i].code[0]; i++) {
  67                 if (!strcmp(code, climates[i].code)) {
  68                         return &climates[i];
  69                 }
  70         }
  71         return 0;
  72 }
  73
  74 struct afrhex *search_hex(struct afrhex *map, int64_t hex) {
  75         for (; map; map = map->next) {
  76                 printf("checking %lld for %lld\n", map->hex, hex);
  77                 if (map->hex == hex) {
  78                         printf("found\n");
  79                         return map;
  80                 }
  81         }
  82         return map;
  83 }
  84
  85 struct afrhex *find_hex(struct afrhex *map, int64_t hex) {
  86         for (; map; map = map->next) {
  87                 if (map->hex == hex) {
  88                         return map;
  89                 }
  90         }
  91         return map;
  92 }
  93
  94 uint64_t total_pop(struct afrhex *map) {
  95         uint64_t count = 0;
  96
  97         for (; map; map = map->next) {
  98                 count += map->freq[0];
  99                 count += map->freq[1];
 100         }
 101         return count/2;
 102 }
 103
 104 /* set up climate cap */
 105 /* base is people per 100 sq km */
 106 /* we divide by an additional 4.0 for the small pop model */
 107 #define SQRT3 1.73205080756887729352
 108 void climate_init(double width, double capfactor) {
 109         double area;
 110         int i;
 111
 112         area = width * width * SQRT3 / 2.0;
 113         for (i = 0; climates[i].code[0]; i++) {
 114                 climates[i].cap = climates[i].basecap * area / 100.0 * capfactor;
 115         }
 116
 117 }
 118
 119 /* width in km */
 120 struct afrhex *map_init(struct maphex *maphexes, double width, double capfactor) {
 121         int i, j;
 122         struct afrhex *hex, *prev, *map;
 123         hex = prev = 0;
 124
 125         climate_init(width, capfactor);
 126
 127         for (i=0; maphexes[i].hex; i++) {
 128                 hex = malloc(sizeof *hex);
 129                 hex->hex = maphexes[i].hex;
 130
 131                 for (j=0; j < 2; j++) {
 132                         hex->freq[j] = 0;
 133                         hex->outmigrate[j] = 0;
 134                         hex->inmigrate[j] = 0;
 135                 }
 136
 137                 hex->climate = find_climate(maphexes[i].climate);
 138                 hex->map = &maphexes[i];
 139                 hex->pop = 0;
 140                 hex->migrants = 0;
 141                 hex->attract = 0;
 142                 hex->next = prev;
 143                 prev = hex;
 144         }
 145
 146         map = hex;
 147
 148         /* pre-compute hex pointers to adjacent hexes */
 149         for (hex = map; hex; hex = hex->next) {
 150                 for (i = 0; i < hex->map->nadj; i++) {
 151                         hex->map->adjhex[i] = find_hex(map, hex->map->adjacent[i]);
 152                 }
 153         }
 154
 155         return map;
 156 }
 157 struct climaterange *find_alt_climate(int hex, int year) {
 158         static int sorted = 0;
 159         static size_t nmemb = -1;
 160         struct alt_climate_search s;
 161
 162 #if 1
 163         if (nmemb == -1) {
 164                 nmemb = 0;
 165                 while (altclimates[nmemb].code[0]) nmemb++;
 166         }
 167
 168         if (!sorted) {
 169                 qsort(altclimates, nmemb, sizeof(struct climaterange), compare_cr);
 170                 sorted = 1;
 171         }
 172
 173         s.hex = hex;
 174         s.year = year;
 175
 176         return bsearch(&s, altclimates, nmemb, sizeof(struct climaterange), alt_climate_compare);
 177
 178 #else
 179         int i;
 180         //fprintf(stderr, "scanning altclimate %d %d\n", hex, year);
 181         for (i=0; altclimates[i].code[0]; i++) {
 182                 if (altclimates[i].hex == hex && altclimates[i].start <= year && altclimates[i].end >= year) {
 183                         return &altclimates[i];
 184                 }
 185         }
 186         return 0;
 187 #endif
 188 }
 189
 190 /* TODO implement change */
 191 struct climate *climateyear(struct afrhex *hex, int year) {
 192         struct climate *c = 0;
 193         struct climaterange *alt;
 194
 195         alt = find_alt_climate(hex->hex, year);
 196
 197         if (alt) {
 198                 c = find_climate(alt->code);
 199         } else {
 200                 c = find_climate(hex->map->climate);
 201         }
 202 #if 0
 203         if (c != hex->climate || hex->hex == 187)
 204         fprintf(stderr, "%d got new climate %s in %d\n", hex->hex, c->code, year);
 205 #endif
 206
 207         return c;
 208 }
 209
 210 #define POP(x) ( ((x[0]) + (x[1])) / 2)
 211
 212 void nextgen(struct afrhex *hexlist, int year, unsigned int climate) {
 213         struct afrhex *hex;
 214
 215         for (hex = hexlist; hex; hex = hex->next) {
 216                 hex->outmigrate[0] = hex->outmigrate[1] = 0;
 217                 hex->inmigrate[0] = hex->inmigrate[1] = 0;
 218                 if (climate) {
 219                         hex->climate = climateyear(hex, year);
 220                 }
 221         }
 222
 223         /* growth and emmigration */
 224         for (hex = hexlist; hex; hex = hex->next) {
 225                 unsigned int pop;
 226                 unsigned int new[2];
 227
 228                 hex->pop = POP(hex->freq);
 229                 pop = (hex->freq[0] + hex->freq[1])/2;
 230                 pop = hex->pop;
 231                 if (pop == 0 ) {
 232                         hex->attract = hex->climate->cap;
 233                         continue; /* nothing else to do here */
 234                 }
 235
 236 #if 0
 237                 printf("hex %lld (%s) pop = [%u, %u] (%d/%d ind)\n", hex->hex, hex->climate->code, hex->freq[0], hex->freq[1], hex->pop, hex->climate->cap);
 238 #endif
 239                 /* TODO if too small, say under 25, either move everyone or die */
 240                 if (pop > hex->climate->cap) {
 241                         pop *= 0.6; /* should probably round down */
 242                 } else {
 243                         pop *= 1.15; /* should probably round up */
 244                 }
 245
 246                 /* pick pop * 2 since we're diploid */
 247                 var_rmultinom(pop * 2, hex->freq, new);
 248                 hex->freq[0] = new[0] ; hex->freq[1] = new[1];
 249                 hex->pop = POP(new);
 250
 251                 hex->migrants = hex->pop - hex->climate->cap;
 252                 if (hex->migrants > 0) {
 253                         var_rhyperv(hex->migrants * 2, new, hex->outmigrate);
 254                         /* TODO could just set it to a quarter of the cap or something */
 255                 }
 256 #if 0
 257                 printf("hex %lld new pop = [%u, %u] (%d ind) (-%u,-%u) (+%u,+%u)\n",
 258                                 hex->hex, hex->freq[0], hex->freq[1], hex->pop,
 259                                 hex->outmigrate[0], hex->outmigrate[1],
 260                                 hex->inmigrate[0], hex->inmigrate[1]
 261                                 );
 262 #endif
 263
 264                 hex->attract = (int)sqrt((double)hex->climate->cap);
 265                 if (hex->attract < hex->climate->cap - hex->pop) {
 266                         hex->attract = hex->climate->cap - hex->pop;
 267                 }
 268         }
 269
 270         /* now handle immigration */
 271         for (hex = hexlist; hex; hex = hex->next) {
 272                 unsigned int attr[18];
 273                 int i;
 274                 int tattr = 0;
 275                 struct afrhex *target;
 276
 277                 if (hex->migrants <= 0) continue;
 278                 /* gather attractiveness of migration targets */
 279 #if 0
 280                 printf("hex %lld migrating: ( ", hex->hex);
 281 #endif
 282                 for (i = 0; i < hex->map->nadj; i++) {
 283                         attr[i] = hex->map->adjhex[i]->attract;
 284                         tattr += attr[i];
 285 #if 0
 286                         printf("%lld:%u ", hex->map->adjhex[i]->hex, attr[i]);
 287 #endif
 288
 289                 }
 290 #if 0
 291                 printf(")\n");
 292 #endif
 293                 /* no where to go ? */
 294                 if (tattr == 0) {
 295                         hex->outmigrate[0] = 0;
 296                         hex->outmigrate[1] = 0;
 297                         hex->migrants = 0;
 298                         continue;
 299                 }
 300
 301                 i = var_pick(i, attr);
 302 #if 0
 303                 printf("choose %d = %lld\n", i+1, hex->map->adjhex[i]->hex);
 304 #endif
 305                 target = hex->map->adjhex[i];
 306                 target->inmigrate[0] += hex->outmigrate[0];
 307                 target->inmigrate[1] += hex->outmigrate[1];
 308         }
 309
 310         /* and finally sum up population so we're ready for the next generation */
 311         for (hex = hexlist; hex; hex = hex->next) {
 312                 hex->freq[0] += hex->inmigrate[0];
 313                 hex->freq[1] += hex->inmigrate[1];
 314                 hex->freq[0] -= hex->outmigrate[0];
 315                 hex->freq[1] -= hex->outmigrate[1];
 316                 hex->pop = (hex->freq[0] + hex->freq[1])/2;
 317 #if 0
 318                 if (hex->pop > 0) {
 319                         printf("hex %lld new pop = [%u, %u] (%d/%u ind) (-%u,-%u) (+%u,+%u)\n",
 320                                         hex->hex, hex->freq[0], hex->freq[1], hex->pop, hex->climate->cap,
 321                                         hex->outmigrate[0], hex->outmigrate[1],
 322                                         hex->inmigrate[0], hex->inmigrate[1]
 323                                         );
 324                 }
 325 #endif
 326         }
 327
 328 }