Carpet Genetic Algorithm

genetic.c

#include <float.h>
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>

#define TARGET  81
#define COLORS  32
#define RULE(r, n, p) (0x1f & (r[n] >> ((p) * 5)))

#define DX(i) ((int)((0x0489a621UL >> (4 * (i) + 0)) & 3) - 1)
#define DY(i) ((int)((0x0489a621UL >> (4 * (i) + 2)) & 3) - 1)

#define COUNTOF(a) (sizeof(a) / sizeof(*a))

static uint32_t
pcg32(uint64_t *s)
{
    *s = *s * UINT64_C(0x9b60933458e17d7d) + UINT64_C(0xd737232eeccdf7ed);
    int shift = 29 - (*s >> 61);
    return *s >> shift;
}

static int
pgm_load(unsigned char *buf)
{
    int w, h, d;
    if (scanf("P2 %d %d %d", &w, &h, &d) != 3)
        return 0;
    if (w != TARGET || h != TARGET || d != 31)
        return 0;
    for (int i = 0; i < TARGET * TARGET; i++) {
        int v;
        scanf("%d", &v);
        if (v < 0 || v > d)
            return 0;
        buf[i] = v;
    }
    return 1;
}

static void
pgm_write(unsigned char *buf, int scale, FILE *o)
{
    fprintf(o, "P2\n%d %d\n%d\n", TARGET * scale, TARGET * scale, COLORS);
    for (int y = 0; y < TARGET * scale; y++)
        for (int x = 0; x < TARGET * scale; x++)
            fprintf(o, "%d\n", buf[(y / scale) * TARGET + (x / scale)]);
}

static void
ppm_write(unsigned char *buf, int scale, FILE *o)
{
    fprintf(o, "P6\n%d %d\n255\n", TARGET * scale, TARGET * scale);
    for (int y = 0; y < TARGET * scale; y++) {
        for (int x = 0; x < TARGET * scale; x++) {
            int v = buf[(y / scale) * TARGET + (x / scale)];
            v = v * 255 / (COLORS - 1);
            fputc(v, o);
            fputc(v, o);
            fputc(v, o);
        }
    }
}

static double
score(unsigned char *a, unsigned char *b)
{
    double error = 0.0;
    for (int y = 0; y < TARGET; y++) {
        for (int x = 0; x < TARGET; x++) {
            int c = 1;
            int amean = a[y * TARGET + x];
            int bmean = b[y * TARGET + x];
            for (int i = 0; i < 8; i++) {
                int xx = x + DX(i);
                int yy = y + DY(i);
                if (xx >= 0 && xx < TARGET && yy >= 0 && yy < TARGET) {
                    amean += a[yy * TARGET + xx];
                    bmean += b[yy * TARGET + xx];
                    c++;
                }
            }
            double delta = amean / (double)c - bmean / (double)c;
            error += delta * delta;
        }
    }
    return error;
}

static void
render(unsigned char *buf, unsigned long long *rules)
{
    for (int y = 0; y < TARGET; y++) {
        for (int x = 0; x < TARGET; x++) {
            int c = 0;
            int px = x;
            int py = y;
            int div = TARGET / 3;
            for (int i = 1; div > 0; i++) {
                c = RULE(rules, c, py / div * 3 + px / div);
                px %= div;
                py %= div;
                div /= 3;
            }
            buf[y * TARGET + x] = c;
        }
    }
}

static void
random_rules(unsigned long long *rules, uint64_t *rng)
{
    for (int i = 0; i < COLORS; i++) {
        uint64_t high = pcg32(rng);
        uint64_t low  = pcg32(rng);
        rules[i] = (high << 32) | low;
    }
}

static void
breed(unsigned long long *o,
      unsigned long long *a,
      unsigned long long *b,
      uint64_t *rng)
{
    uint32_t select = pcg32(rng);
    for (int i = 0; i < COLORS; i++)
        o[i] = (select >> i) & 1 ? a[i] : b[i];

    int mutations = pcg32(rng) % 32;
    for (int i = 0; i < mutations; i++) {
        select = pcg32(rng);
        int r = (select >>  0) % COLORS;
        int c = (select >>  5) %  9;
        int t = (select >>  9) % COLORS;
        uint64_t mask = UINT64_C(0x1f) << (c * 5);
        o[r] = (o[r] & ~mask) | ((uint64_t)t << (c * 5));
    }
}

static int
dblcmp(const void *pa, const void *pb)
{
    double a = *(double *)pa;
    double b = *(double *)pb;
    if (a < b)
        return -1;
    else if (a > b)
        return 1;
    return 0;
}

static void
rule_write(unsigned long long *rule, FILE *o)
{
    fputs("static const unsigned long long rules[] = {\n", o);
    for (int i = 0; i < COLORS; i++) {
        if (i % 4 == 0)
            fputs("    ", o);
        unsigned long long v = rule[i] & 0x1fffffffffff;
        fprintf(o, "0x%012llx,%c", v, " \n"[i % 4 == 3]);
    }
    fputs("};\n", o);
    for (int i = 0; i < COLORS; i++) {
        fprintf(o, "// %d ->", i);
        for (int c = 0; c < 9; c++) {
            int p = RULE(rule, i, c);
            fprintf(o, " %d%c", p, ",\n"[c == 8]);
        }
    }
}

int main(void)
{
    uint64_t rng[] = {0x9e8480dd162324e1};
    unsigned long long rules[200][COLORS];
    unsigned char input[TARGET * TARGET];

    if (!pgm_load(input)) {
        fputs("Invalid input\n", stderr);
        return 1;
    }

    for (size_t i = 0; i < COUNTOF(rules); i++)
        random_rules(rules[i], rng);

    FILE *video = fopen("video.ppm", "wb");
    double global_best = DBL_MAX;
    for (long long g = 0; ; g++) {
        struct {
            double score;
            size_t i;
        } scores[COUNTOF(rules)];
        unsigned long long best[10][COLORS];

        /* Find the 10 best rulesets */
        for (size_t i = 0; i < COUNTOF(rules); i++) {
            unsigned char buf[TARGET * TARGET];
            render(buf, rules[i]);
            scores[i].score = score(input, buf);
            scores[i].i = i;
        }
        qsort(scores, COUNTOF(scores), sizeof(*scores), dblcmp);

        /* Breed next generation */
        for (int i = 0; i < COUNTOF(best); i++)
            memcpy(best[i], rules[scores[i].i], sizeof(best[i]));
        for (int i = 0; i < COUNTOF(rules) - 5; i++) {
            uint32_t select = pcg32(rng);
            int a = (select >>  0) % 10;
            int b = (select >> 16) % 10;
            breed(rules[i], best[a], best[b], rng);
        }
        for (int i = COUNTOF(rules) - 5; i < COUNTOF(rules); i++)
            random_rules(rules[i], rng);

        /* Report on progress */
        if (scores[0].score < global_best) {
            global_best = scores[0].score;

            /* Write out the current best image */
            unsigned char buf[TARGET * TARGET];
            FILE *image = fopen("_progress.pgm", "w");
            render(buf, best[0]);
            pgm_write(buf, 10, image);
            fclose(image);
            rename("_progress.pgm", "progress.pgm");

            /* Write out best image as video frame */
            ppm_write(buf, 10, video);
            fflush(video);

            /* Write out the ruleset */
            FILE *save = fopen("best.c", "w");
            rule_write(best[0], save);
            fclose(save);
        }
        printf("%12lld %f %f\n", g, global_best, scores[0].score);
    }
}

Makefile

CC     = c99
CFLAGS = -Wall -Wextra -Ofast
all: genetic render
genetic: genetic.c
	$(CC) $(CFLAGS) $(LDFLAGS) -o $@ genetic.c $(LDLIBS)
render: render.c
	$(CC) $(CFLAGS) $(LDFLAGS) -o $@ render.c $(LDLIBS)

render.c

#include <stdio.h>

#define TARGET 81
#define SCALE 10
#define RULE(n, p) (0x1f & (rules[n] >> ((p) * 5)))
static const unsigned long long rules[] = {
    0x12cff8e3eefd, 0x0240c6930008, 0x1ff73beef75b, 0x16dfefffebfd,
    0x1e6ee3c67375, 0x1d2b8e100422, 0x1d17fbee6b0a, 0x1aed8664928b,
    0x0e4918f5fc67, 0x0509ec4e0c41, 0x0167184a847f, 0x02824300785a,
    0x0002a000a2e7, 0x17ef6e67ce76, 0x1af987ff3def, 0x1af8fbef94a5,
    0x14281121307f, 0x0cb535de8846, 0x1feffa65987f, 0x002c400099a8,
    0x0508bab67842, 0x1dd9b7457416, 0x080801111522, 0x02875e3f8bff,
    0x02f17a988ff0, 0x0e7ba2f8c451, 0x1dff79ff739f, 0x1f7ff0e6fe3d,
    0x07e995a15160, 0x1dfffbfffffd, 0x1be73fee7bff, 0x1dff7daff7da,
};

int main(void)
{
    printf("P2\n%d %d\n29\n", TARGET * SCALE, TARGET * SCALE);
    for (int y = 0; y < TARGET * SCALE; y++) {
        for (int x = 0; x < TARGET * SCALE; x++) {
            int c = 0;
            int px = x / SCALE;
            int py = y / SCALE;
            int div = TARGET / 3;
            for (int i = 1; div > 0; i++) {
                c = RULE(c, py / div * 3 + px / div);
                px %= div;
                py %= div;
                div /= 3;
            }
            printf("%d\n", c);
        }
    }
}