/* NAME * ca - simulate arbitrary one-dimensional cellular automata * NOTES * The main portion of this code is relatively straight forward, * whereas the code to compute lambda (or random rules of a * specified lambda) is a bit tricky and is solved by a method * that resembles dynamic programming. * * Also, the states are stored in a circular buffer such that * the sums for every cell are computed as fast as possible with * the number of additions on the order of two per cell. * MISCELLANY * When supplying a lambda value for a random rule, it may not * be possible to find a string with that lambda value because * one may not exist. In this case, the program will do its * best to find one as close as possible. In any event, the * algorithm for finding random rules strings for a specified * lambda value is non-deterministic and may not always find * a perfect match even if one exists. However, it will work * well with high probability, and even when it doesn't find * a perfect match it almost always gets close. * BUGS * No sanity checks are performed to make sure that any of the * options make sense. * AUTHOR * Copyright (c) 1997, Gary William Flake. * * Permission granted for any use according to the standard GNU * ``copyleft'' agreement provided that the author's comments are * neither modified nor removed. No warranty is given or implied. */ #include #include #include #include "misc.h" double lambda = -1.0; int width = 640, height = 480, states = 2, radius = 1, mag = 1; int wrap = 1, seed = 0, binary = 0, invert = 0, sq = 1; char *term = NULL, *rules = "0110", *init = "11"; char help_string[] = "\ Computes a one-dimensional cellular automata. The evolution of the \ CA is determined by the number of states, the radius size, the initial \ state, and the supplied rule. A rule is specified by a (states - 1) * \ (radius * 2 + 1) length string. At each time step a sum of each cell \ plus all of its neighbors within the radius is computed. That sum is \ used as an index into the rules string which determines the next step \ For example, with radius = 1 and states = 2 the rule \"0110\" specifies \ that sums of 0 and 3 map to the 0 state, and sums of 1 and 2 map to the \ 1 state. A negative init string randomly initializes the starting \ states. If the init string is \"-N\" then each cell has a 1 in N chance \ of being non-zero.\ "; OPTION options[] = { { "-width", OPT_INT, &width, "Width of the plot in pixels." }, { "-height", OPT_INT, &height, "Height of the plot in pixels." }, { "-states", OPT_INT, &states, "Number of CA states." }, { "-radius", OPT_INT, &radius, "Radius of CA neighborhood." }, { "-seed", OPT_INT, &seed, "Random seed." }, { "-wrap", OPT_SWITCH, &wrap, "Use a wrap-around space?" }, { "-rules", OPT_STRING, &rules, "CA rules to use." }, { "-init", OPT_STRING, &init, "Starting state (< 0 is random)." }, { "-lambda", OPT_DOUBLE, &lambda, "Lambda value for random rules." }, { "-sq", OPT_SWITCH, &sq, "Enforce strong quiescence?" }, { "-bin", OPT_SWITCH, &binary, "Binary colors?" }, { "-inv", OPT_SWITCH, &invert, "Invert all colors?" }, { "-mag", OPT_INT, &mag, "Magnification factor." }, { "-term", OPT_STRING, &term, "How to plot points." }, { NULL, OPT_NULL, NULL, NULL } }; /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /* This function has two distinct modes of operation; thus, there are two ways of calling it. The first way is with rules equal to a valid rules string (and, therefore, non-NULL). In this case, the passed lambda value is ignored and the rules string is analyzed to determine its lambda value. In the second mode, with rule equal to NULL, a random string is generated with a lambda value as close to the requested value as possible. In either case, both the states, radius, and sq arguments are used. */ char *dorules(int states, int radius, int sq, double *lambda, char *rules) { int area = 2 * radius + 1; int len = (states - 1) * area + 1; int N = (int) pow(states, area); int **table, *bits; int i, j, k, sum, a, b, aval, bval, noimprove; double *vals, oldlambda, newlambda; /* Simple sanity check. */ if(rules && strlen(rules) != len) { fprintf(stderr, "Rule length should be %d not %d\n", len, strlen(rules)); exit(1); } /* Build table that contains the number of transitions to a particular * state. More specifically, table[i][j] will contain the number of * possible ways in which a string of length (i + 1) can sum up to j * using any number for a state from 0 to (states - 1). Thus, the * final row of table[area - 1] will contain the number of ways we * can sum up to the column for this area size. * * We compute this by exploiting the recurrence relationship * table[i][j] = sum(k = 0 .. states - 1) table[i-1][j-k] * with any nonsense pair of [i-1][j-k] yielding zero. */ table = xmalloc(sizeof(int *) * area); for(i = 0; i < area; i++) table[i] = xmalloc(sizeof(int) * len); /* Start with base condition. */ for(i = 0; i < len; i++) table[0][i] = ((i < states) ? 1 : 0); /* Fill via recurrence relationship. */ for(i = 1; i < area; i++) for(j = 0; j < len; j++) { sum = 0; for(k = 0; k < states; k++) if(j - k >= 0) sum += table[i-1][j-k]; table[i][j] = sum; } /* Now turn the final row of the table into lambda units. */ vals = xmalloc(sizeof(double) * len); for(i = 0; i < len; i++) vals[i] = table[area-1][i] / (double) N; /* If we fall into this statement, then we are in a query mode. */ if(rules) { *lambda = 0; for(i = 0; i < len; i++) if(rules[i] != '0') *lambda += vals[i]; /* Clean up stuff that is no longer needed. */ for(i = 0; i < area; i++) free(table[i]); free(table); free(vals); /*MRM begin*/ #if __dest_os != __mac_os fprintf(stderr, "supplied rule = '%s'\n", rules); fprintf(stderr, "actual lambda = %f\n", *lambda); #endif /*MRM end*/ /* Return the original rule, but also set *lambda to point to * the real lambda value. */ return rules; } /* Since we've made it this far into the code, it means that we need * randomly generated rules. */ /* Make a random vector of bits. */ bits = xmalloc(sizeof(int) * len); for(i = 0; i < len; i++) bits[i] = random() % 2; /* Optionally enforce strong quiescence. */ if(sq) { bits[0] = 0; for(i = 1; i < states; i++) bits[i * area] = 1; } /* Calculate lambda for these bits. Note that we are only * distinguishing quiescent from non-quiescent states at this point. */ oldlambda = 0; for(i = 0; i < len; i++) oldlambda += vals[i] * bits[i]; /* Now, try to improve the lambda count for a bunch of steps. */ noimprove = 0; while(1) { /* Pick two random bits to change and two new values */ while(1) { a = random() % len; /* Make sure it doesn't violate strong quiescence. */ if(sq && (a % area == 0)) continue; aval = random() % 2; break; } while(1) { b = random() % len; /* Make sure it doesn't violate strong quiescence. */ if(sq && (b % area == 0)) continue; bval = random() % 2; break; } newlambda = oldlambda + (aval - bits[a]) * vals[a] + (bval - bits[b]) * vals[b]; /* If this is better, then accept it. */ if(fabs(newlambda - *lambda) < fabs(oldlambda - *lambda)) { oldlambda = newlambda; bits[a] = aval; bits[b] = bval; noimprove = 0; } else noimprove++; /* We've gone a long time without improving things, so quit. */ if(noimprove == 1000) break; } /* Now that we have the bits, we can pick state transitions * and return the new rule. */ rules = xmalloc(sizeof(char) * (len + 1)); rules[len] = 0; for(i = 0; i < len; i++) rules[i] = bits[i] ? ((random() % (states - 1)) + '1') : '0'; if(sq) for(i = 1; i < states; i++) rules[i * area] = i + '0'; for(i = 0; i < area; i++) free(table[i]); free(table); free(bits); free(vals); /*MRM begin*/ #if __dest_os != __mac_os fprintf(stdout, "generated rules = '%s'\n", rules); fprintf(stdout, "generated lambda = %f\n", oldlambda); #endif /*MRM end*/ return rules; } /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ int main(int argc, char **argv) { extern int plot_mag; extern int plot_inverse; char *old, *new, *swap; int i, j, len, sum, odds; get_options(argc, argv, options, help_string); plot_mag = mag; plot_inverse = invert; plot_init(width, height, binary ? 2 : states, term); srandom(seed); /* The extra 2 * radius + 2 cells are "buffers" that allow us to * compute the sums rapidly. */ old = xmalloc((width + 2 * radius + 2) * sizeof(char)); new = xmalloc((width + 2 * radius + 2) * sizeof(char)); for(i = 0; i < width + 2 * radius + 2; i++) old[i] = new[i] = 0; /* Initialize the first state. If negative, use random init state. */ len = strlen(init); if(init[0] == '-') { odds = -atoi(init); for(i = radius + 1; i < width + radius + 1; i++) if(random() % odds == 0) old[i] = random() % (states - 1) + 1; } else for(i = (width - len) / 2 + radius + 1, j = 0; j < len; i++, j++) old[i] = init[j] - '0'; /* Generate a random lambda rule if needed. */ len = (states - 1) * (radius * 2 + 1) + 1; if(lambda >= 0.0 && lambda <= 1.0) rules = dorules(states, radius, sq, &lambda, NULL); else rules = dorules(states, radius, sq, &lambda, rules); /* The main loop. All of the real work is done here. */ for(i = 0; i < height; i++) { /* If wrapping, then the CA is really a ring. Place the * right edge values on the buffer to the left, and the * left edge value on the buffer to the right. */ if(wrap) for(j = 0; j < radius; j++) { old[j + 1] = old[width + 1 + j]; old[width + radius + 1 + j] = old[radius + 1 + j]; } /* Initialize the running sum. */ for(j = 0, sum = 0; j < 2 * radius + 1; j++) sum += old[j]; /* If we know the sum of the cell to the left then we just subtract * its left-most neighbor within radius and add our right-most * neighbor within radius to compute this cell's sum. */ for(j = radius + 1; j < width + radius + 1; j++) { sum = sum + old[j + radius] - old[j - radius - 1]; new[j] = rules[sum] - '0'; plot_point(j - radius - 1, i, binary ? (old[j] ? 1 : 0) : old[j]); } swap = old; old = new; new = swap; } plot_finish(); exit(0); } /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */