Advertisement
Guest User

Untitled

a guest
Feb 23rd, 2016
78
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. //
  2. // weasel.c - A variant of Dawkins' "Weasel" program for modeling drift vs selection.
  3. //
  4. // Author: keiths at www.theskepticalzone.com
  5. //
  6. // DOS/Windows modifications by petrushka at The Skeptical Zone.
  7. //
  8. // Compile under Linux using 'gcc -std=gnu99 -lm weasel.c -o weasel'.
  9. // Add a '#define WINDOWS' before compiling under Windows or DOS.
  10. //
  11.  
  12. #define VERSION "02232016"
  13.  
  14. #include <stdio.h>
  15. #include <stdlib.h>
  16. #include <string.h>
  17. #include <time.h>
  18. #include <unistd.h>
  19. #include <ctype.h>
  20. #include <math.h>
  21.  
  22. #ifdef WINDOWS
  23. #include <conio.h>
  24. #else
  25. #include <sys/select.h>
  26. #include <termios.h>
  27. #endif
  28.  
  29. // adjustable parameters
  30. #define POPULATION_SIZE 50000 // total population size
  31. #define NUM_SURVIVORS 1 // number of survivors per generation (must be less than POPULATION_SIZE)
  32. #define GENOME_LEN 28 // number of loci in each genome (each locus is one character)
  33. #define MUTATION_RATE 0.01 // probability that a locus is mutated per generation
  34. #define SELECTION_COEFF 5.0 // selection coefficient
  35. #define DISPLAY_INTERVAL 10 // display results every n generations, where n == DISPLAY_INTERVAL
  36. #define PAUSE_TIME 50 // pause time in milliseconds
  37. #define GENOMES_TO_DISPLAY 1 // in STEP_MODE, number of genomes to display at each step
  38. #define STEP_MODE 1 // if set, program will pause every n generations and display the genomes
  39. #define ENABLE_LATCHING 0 // if set, enable latching
  40. #define ENABLE_SELECTION 1 // if set, enable selection; otherwise, select survivors with no bias
  41.  
  42. #define CHARSET_SIZE 27 // 26 uppercase letters plus the space character
  43. #define NUM_LINES 300 // number of lines in the program's display "window"
  44. #define CHARS_PER_LINE 120 // for each line in the "window"
  45.  
  46. long long generation; // generation count
  47. long long histogram[GENOME_LEN+1]; // histogram of Hamming distance from the target
  48. double mutation_rate;
  49. double selection_coefficient;
  50. double floating_num_survivors;
  51. int num_survivors;
  52. char buffer[20];
  53.  
  54. // A genome consists of an array of loci, the associated fitness value, and a field indicating the number of matches
  55. // between the genome and the target phrase.
  56. typedef struct {
  57. char locus[GENOME_LEN];
  58. double fitness;
  59. int matches;
  60. } genome_t;
  61.  
  62. genome_t target; // target phrase
  63. genome_t genome_array[POPULATION_SIZE]; // genomes of all members of the population
  64.  
  65. char weasel_to_ascii(char c) {
  66. switch(c) {
  67. case 0: return ' ';
  68. case 1: return 'A';
  69. case 2: return 'B';
  70. case 3: return 'C';
  71. case 4: return 'D';
  72. case 5: return 'E';
  73. case 6: return 'F';
  74. case 7: return 'G';
  75. case 8: return 'H';
  76. case 9: return 'I';
  77. case 10: return 'J';
  78. case 11: return 'K';
  79. case 12: return 'L';
  80. case 13: return 'M';
  81. case 14: return 'N';
  82. case 15: return 'O';
  83. case 16: return 'P';
  84. case 17: return 'Q';
  85. case 18: return 'R';
  86. case 19: return 'S';
  87. case 20: return 'T';
  88. case 21: return 'U';
  89. case 22: return 'V';
  90. case 23: return 'W';
  91. case 24: return 'X';
  92. case 25: return 'Y';
  93. case 26: return 'Z';
  94. }
  95. }
  96.  
  97. char ascii_to_weasel(char c) {
  98. switch(c) {
  99. case ' ': return 0;
  100. case 'A': return 1;
  101. case 'B': return 2;
  102. case 'C': return 3;
  103. case 'D': return 4;
  104. case 'E': return 5;
  105. case 'F': return 6;
  106. case 'G': return 7;
  107. case 'H': return 8;
  108. case 'I': return 9;
  109. case 'J': return 10;
  110. case 'K': return 11;
  111. case 'L': return 12;
  112. case 'M': return 13;
  113. case 'N': return 14;
  114. case 'O': return 15;
  115. case 'P': return 16;
  116. case 'Q': return 17;
  117. case 'R': return 18;
  118. case 'S': return 19;
  119. case 'T': return 20;
  120. case 'U': return 21;
  121. case 'V': return 22;
  122. case 'W': return 23;
  123. case 'X': return 24;
  124. case 'Y': return 25;
  125. case 'Z': return 26;
  126. }
  127. }
  128.  
  129. // Swap two genomes given their indices.
  130. void swap_genomes(int g1, int g2) {
  131. genome_t temp_genome;
  132. if (g1 == g2) {
  133. return;
  134. } else {
  135. temp_genome = genome_array[g1];
  136. genome_array[g1] = genome_array[g2];
  137. genome_array[g2] = temp_genome;
  138. }
  139. }
  140.  
  141. // Thoroughly and randomly scramble the population within the genome array.
  142. void scramble_genomes() {
  143. for (int i=0; i<6; i++) {
  144. for (int j=0; j<POPULATION_SIZE; j++) {
  145. swap_genomes(j, rand()%POPULATION_SIZE);
  146. }
  147. }
  148. }
  149.  
  150. // Calculate a genome's fitness by determining the number of loci at which
  151. // it matches the target, then applying the formula (1+s)^m where s is the selection
  152. // coefficient and m is the number of matches.
  153. double fitness(genome_t *genome) {
  154. genome->matches = 0;
  155.  
  156. // scan the entire genome, counting the number of loci that match the target
  157. for (int i=0; i < GENOME_LEN; i++) {
  158. if (genome->locus[i] == target.locus[i]) ++genome->matches;
  159. }
  160.  
  161. return pow(1.0+selection_coefficient, genome->matches);
  162. }
  163.  
  164. // Display a genome as a character string.
  165. void display_genome(genome_t *genome) {
  166. for (int i=0; i < GENOME_LEN; i++) {
  167. putchar(weasel_to_ascii(genome->locus[i]));
  168. }
  169. }
  170.  
  171. // Convert an ascii string to a Weasel genome
  172. void weaselize (char *ascii_str, genome_t *genome) {
  173. for (int i=0; i<GENOME_LEN; i++) {
  174. genome->locus[i] = ascii_to_weasel(ascii_str[i]);
  175. }
  176. }
  177.  
  178. // Move cursor to the home position
  179. void cursor_home() {
  180. #ifdef WINDOWS
  181. gotoxy(0,0);
  182. #else
  183. printf("\x1b[H");
  184. #endif
  185. }
  186.  
  187. // Overwrite everything on the screen with spaces
  188. void clear_screen() {
  189. #ifdef WINDOWS
  190. clrscr();
  191. #else
  192. cursor_home();
  193. for (int i=0; i < NUM_LINES; i++) {
  194. for (int j=0; j < CHARS_PER_LINE; j++) {
  195. putchar(' ');
  196. }
  197. printf("\n\r");
  198. }
  199. #endif
  200. }
  201.  
  202. // Mutate 'old' genome, placing results in 'new'
  203. void mutate(genome_t *old, genome_t *new) {
  204. for (int i=0; i < GENOME_LEN; i++) {
  205. if ((double)rand()/(double)RAND_MAX < mutation_rate * 27.0/26.0 && (!ENABLE_LATCHING || old->locus[i] != target.locus[i])) {
  206. new->locus[i] = (char)(rand()%CHARSET_SIZE); // mutate this locus
  207. } else {
  208. new->locus[i] = old->locus[i]; // don't mutate
  209. }
  210. }
  211. new->fitness = fitness(new);
  212. }
  213.  
  214. // Initialize a genome with a random string
  215. void init_genome_rand(genome_t *genome) {
  216. for (int i=0; i < GENOME_LEN; i++) {
  217. genome->locus[i] = (char)(rand()%CHARSET_SIZE);
  218. }
  219. genome->fitness = fitness(genome);
  220. }
  221.  
  222. // Comparison function for use by qsort() library routine
  223. int compare_fitness (const void * elem1, const void * elem2) {
  224. if (((genome_t *)elem2)->fitness > ((genome_t *)elem1)->fitness) return 1;
  225. if (((genome_t *)elem2)->fitness == ((genome_t *)elem1)->fitness) return 0;
  226. return -1;
  227. }
  228.  
  229. // Use weighted random selection to choose offspring for the survivor slots
  230. void pick_survivors() {
  231.  
  232. double total_fitness = 0.0;
  233. double cumulative_fitness;
  234. double rand_num;
  235. int j;
  236.  
  237. // add up all fitnesses
  238. for (int i=0; i<POPULATION_SIZE; i++) {
  239. total_fitness += genome_array[i].fitness;
  240. }
  241.  
  242. // for each survivor slot
  243. for (int i=0; i<num_survivors; i++) {
  244.  
  245. cumulative_fitness = 0.0;
  246.  
  247. // generate a random number in the range 0-total_fitness
  248. rand_num = (double)rand()/(double)RAND_MAX * total_fitness;
  249.  
  250. // from the remaining offspring, select the one whose range is hit by the random number
  251. for (j=i; j<POPULATION_SIZE; j++) {
  252. cumulative_fitness += genome_array[j].fitness;
  253. if (cumulative_fitness >= rand_num) {
  254. break;
  255. }
  256. }
  257.  
  258. // swap it with the one in the survivor slot
  259. swap_genomes(i,j);
  260.  
  261. // adjust the total fitness to reflect only the remaining offspring
  262. total_fitness -= genome_array[i].fitness;
  263. }
  264. }
  265.  
  266. void display_histogram() {
  267. long long max = 0;
  268. long long width;
  269.  
  270. printf("Histogram of Hamming distances:\n\n\r");
  271.  
  272. for (int i=0; i<=GENOME_LEN; i++) {
  273. if (histogram[i] > max) {
  274. max = histogram[i];
  275. }
  276. }
  277.  
  278. if (max > 0) {
  279. for (int i=0; i<=GENOME_LEN; i++) {
  280. printf("%2d: %8lld ", i, histogram[i]);
  281. width = 60 * histogram[i]/max;
  282. for (int j=0; j<width; j++) {
  283. printf("X");
  284. }
  285. printf("\n\r");
  286. }
  287. }
  288. printf("\n\r");
  289. }
  290.  
  291. #ifndef WINDOWS
  292. //
  293. // From http://stackoverflow.com/questions/448944/c-non-blocking-keyboard-input
  294. //
  295.  
  296. struct termios orig_termios;
  297.  
  298. void reset_terminal_mode()
  299. {
  300. tcsetattr(0, TCSANOW, &orig_termios);
  301. }
  302.  
  303. void set_conio_terminal_mode()
  304. {
  305. struct termios new_termios;
  306.  
  307. /* take two copies - one for now, one for later */
  308. tcgetattr(0, &orig_termios);
  309. memcpy(&new_termios, &orig_termios, sizeof(new_termios));
  310.  
  311. /* register cleanup handler, and set the new terminal mode */
  312. atexit(reset_terminal_mode);
  313. cfmakeraw(&new_termios);
  314. tcsetattr(0, TCSANOW, &new_termios);
  315. }
  316.  
  317. int kbhit()
  318. {
  319. struct timeval tv = { 0L, 0L };
  320. fd_set fds;
  321. FD_ZERO(&fds);
  322. FD_SET(0, &fds);
  323. return select(1, &fds, NULL, NULL, &tv);
  324. }
  325.  
  326. int getch()
  327. {
  328. int r;
  329. unsigned char c;
  330. if ((r = read(0, &c, sizeof(c))) < 0) {
  331. return r;
  332. } else {
  333. return c;
  334. }
  335. }
  336. #endif
  337.  
  338. // Read a double-precision floating point number from keyboard input.
  339. // Can't use scanf here because of the non-blocking I/O.
  340. void get_double(double *ptr) {
  341. int i = 0;
  342. char c;
  343. do {
  344. while (!kbhit()) { // wait for keystroke
  345. ;
  346. }
  347. c = getch();
  348. putchar(c);
  349. fflush(stdout);
  350. buffer[i++] = c;
  351. } while (buffer[i-1] != '\r' && i < sizeof(buffer)); // accumulate characters in buffer until CR
  352. buffer[i-1] = '\n';
  353. i=0;
  354. fflush(stdout);
  355. sscanf(buffer, "%lf", ptr);
  356. }
  357.  
  358. int main() {
  359.  
  360. int i;
  361. int selection_enabled = ENABLE_SELECTION;
  362. char ascii_target[GENOME_LEN+1];
  363. char c;
  364.  
  365. #ifndef WINDOWS
  366. set_conio_terminal_mode();
  367. #endif
  368.  
  369. printf("\n\n\n\n\r");
  370.  
  371. // set cursor to home position, clear screen, and rehome
  372. clear_screen();
  373. cursor_home();
  374.  
  375. // seed the random number generator with the current time
  376. srand(time(NULL));
  377.  
  378. // set up ascii version of initial fitness target
  379. ascii_target[0] = 'M';
  380. ascii_target[1] = 'E';
  381. ascii_target[2] = 'T';
  382. ascii_target[3] = 'H';
  383. ascii_target[4] = 'I';
  384. ascii_target[5] = 'N';
  385. ascii_target[6] = 'K';
  386. ascii_target[7] = 'S';
  387. ascii_target[8] = ' ';
  388. ascii_target[9] = 'I';
  389. ascii_target[10] = 'T';
  390. ascii_target[11] = ' ';
  391. ascii_target[12] = 'I';
  392. ascii_target[13] = 'S';
  393. ascii_target[14] = ' ';
  394. ascii_target[15] = 'L';
  395. ascii_target[16] = 'I';
  396. ascii_target[17] = 'K';
  397. ascii_target[18] = 'E';
  398. ascii_target[19] = ' ';
  399. ascii_target[20] = 'A';
  400. ascii_target[21] = ' ';
  401. ascii_target[22] = 'W';
  402. ascii_target[23] = 'E';
  403. ascii_target[24] = 'A';
  404. ascii_target[25] = 'S';
  405. ascii_target[26] = 'E';
  406. ascii_target[27] = 'L';
  407.  
  408. // set up weasel version of the fitness target
  409. weaselize(ascii_target, &target);
  410.  
  411. // randomly initialize the genomes of the entire population
  412. for (i=0; i < POPULATION_SIZE; i++) {
  413. init_genome_rand(&genome_array[i]);
  414. }
  415. printf("\n\n\rInitial Organism 0 Genome:\n\n\r");
  416. display_genome(&genome_array[0]);
  417. printf("\n\n\n\r");
  418.  
  419. generation = 0;
  420. mutation_rate = MUTATION_RATE;
  421. selection_coefficient = SELECTION_COEFF;
  422. num_survivors = NUM_SURVIVORS;
  423.  
  424. // initialize the histogram to zeroes
  425. for (int i=0; i<=GENOME_LEN; i++) {
  426. histogram[i] = 0;
  427. }
  428.  
  429. while (1) {
  430.  
  431. // preserve the survivors, but convert the rest into mutated copies of the survivors
  432. for (i=0; i < POPULATION_SIZE - num_survivors; i++) {
  433. mutate(&genome_array[i%num_survivors], &genome_array[i+num_survivors]);
  434. }
  435.  
  436. // now mutate the survivors themselves
  437. for (i=0; i<num_survivors; i++) {
  438. mutate(&genome_array[i], &genome_array[i]);
  439. }
  440.  
  441. if (selection_enabled) {
  442. // use weighted random selection to choose offspring for the survivor slots
  443. pick_survivors();
  444. } else {
  445. // randomly scramble the genome array
  446. scramble_genomes();
  447. }
  448.  
  449. ++histogram[GENOME_LEN-genome_array[0].matches];
  450.  
  451. ++generation;
  452.  
  453. // in step mode, pause and display all genomes periodically
  454. if (STEP_MODE && generation % DISPLAY_INTERVAL == 0) {
  455. clear_screen();
  456. cursor_home();
  457. printf("Weasel version %s\n\n\r", VERSION);
  458. printf("Press 'h' at any time for help\n\n\n\n\r");
  459. printf("Target Phrase: ");
  460. printf("%s", ascii_target);
  461. printf("\n\n\rGeneration: %lld Number of survivors per generation: %d\n\r", generation, num_survivors);
  462. printf("\n\rSelection: %s Coefficient: %1.2lf Mutation rate: %1.2lf\n\n\r", selection_enabled ? "ON" : "OFF", selection_coefficient, mutation_rate);
  463. for (i = 0; i < GENOMES_TO_DISPLAY; i++) {
  464. printf("Organism %d Fitness %1.2le Hamming distance %d\n\n\r", i, genome_array[i].fitness, GENOME_LEN-genome_array[0].matches);
  465. display_genome(&genome_array[i]);
  466. printf("\r\n\n");
  467. }
  468. display_histogram();
  469. #ifdef WINDOWS
  470. sleep(PAUSE_TIME); // pause briefly before continuing
  471. #else
  472. usleep(1000*PAUSE_TIME); // pause briefly before continuing
  473. #endif
  474. }
  475.  
  476. if (!STEP_MODE && generation % DISPLAY_INTERVAL == 0) {
  477. printf("generation = %lld fitness = %1.2le Hamming distance = %d\n\r", generation, genome_array[0].fitness, GENOME_LEN-genome_array[0].matches);
  478. }
  479.  
  480. if (kbhit()) { // keystroke detected
  481. switch(getch()) {
  482.  
  483. case 'c': // clear the histogram
  484. for (int i=0; i<=GENOME_LEN; i++) {
  485. histogram[i] = 0;
  486. }
  487. break;
  488.  
  489. case 'f': // change the selection coefficient
  490. printf("Enter new selection coefficient: ");
  491. fflush(stdout);
  492. get_double(&selection_coefficient);
  493. break;
  494.  
  495. case 'h': // print help message
  496. printf("Commands:\r\n");
  497. printf(" c - clear the histogram data\r\n");
  498. printf(" f - change the selection coefficient\r\n");
  499. printf(" h - print this help message\r\n");
  500. printf(" m - change the mutation rate\r\n");
  501. printf(" p - pause until a key is pressed\r\n");
  502. printf(" q - quit the program\r\n");
  503. printf(" s - toggle selection on/off\r\n");
  504. printf(" t - change the target phrase\r\n");
  505. printf(" v - change the number of survivors per generation\r\n");
  506. printf("\r\nHit any key to continue\r\n");
  507. while(!kbhit()){
  508. ;
  509. }
  510. (void)getch(); // consume character
  511. break;
  512.  
  513. case 'm': // change the mutation rate
  514. printf("Enter new mutation rate: ");
  515. fflush(stdout);
  516. get_double(&mutation_rate);
  517. break;
  518.  
  519. case 'p': // pause until next keystroke
  520. printf("Paused: Hit any key to continue\r\n");
  521. while(!kbhit()){
  522. ;
  523. }
  524. (void)getch(); // consume character
  525. break;
  526.  
  527. case 'q': // quit the program
  528. exit(0);
  529.  
  530. case 's': // toggle selection on and off
  531. selection_enabled = !selection_enabled;
  532. break;
  533.  
  534. case 't': // set new target phrase
  535. printf("Enter new target phrase (spaces and letters only): ");
  536. fflush(stdout);
  537. i = 0;
  538. do {
  539. while (!kbhit()) { // wait for keystroke
  540. ;
  541. }
  542. c = toupper(getch());
  543. if (c == '\r') break;
  544. if (isalpha(c) || c == ' ') {
  545. putchar(c);
  546. fflush(stdout);
  547. ascii_target[i++] = c; // add character to target phrase
  548. }
  549. } while (ascii_target[i-1] != '\r' && i-1 < GENOME_LEN);
  550. for(; i<GENOME_LEN; i++) {
  551. ascii_target[i] = ' '; // pad remainder of target with spaces
  552. }
  553. weaselize(ascii_target, &target);
  554. // need to recompute fitness values because the target has changed
  555. for (i=0; i<POPULATION_SIZE; i++) {
  556. genome_array[i].fitness = fitness(&genome_array[i]);
  557. }
  558. break;
  559.  
  560. case 'v': // change the number of survivors per generation
  561. printf("Enter desired number of survivors: ");
  562. fflush(stdout);
  563. get_double(&floating_num_survivors);
  564. num_survivors = (int)floating_num_survivors;
  565. break;
  566.  
  567. }
  568. }
  569. }
  570.  
  571. printf("\n\n\rFinal generation = %lld fitness = %1.2le Hamming distance = %d\n\r", generation, genome_array[0].fitness, GENOME_LEN-genome_array[0].fitness);
  572. printf("\n\rFinal genome:\n\n\r");
  573. display_genome(&genome_array[0]);
  574. printf("\n\n\r");
  575. #ifndef WINDOWS
  576. tcsetattr(0, TCSANOW, &orig_termios);
  577. #endif
  578. }
Advertisement
RAW Paste Data Copied
Advertisement