#!/usr/bin/env python

"""
EvoLife Cellular Automaton implementation using CUDA.
Rules are:
- Each living cell has its own birth/sustain ruleset and an energy level;
- Cell is loosing all energy if number of neighbours is not in its sustain rule;
- Cell is born with max energy if there are exactly N neighbours with N in their birth rule;
  - Same is applied for living cells (re-occupation case), but only with different genomes;
- If there are several birth situations with different N possible, we choose one with larger N;
- Newly born cell's ruleset calculated as crossover between 'parent' cells rulesets;
- If cell is involved in breeding as a 'parent', it's loosing `BIRTH_COST` units of energy per each non-zero gene passed;
  - This doesn't apply in re-occupation case;
- Every turn, cell is loosing `DEATH_SPEED` units of energy;
- Cell with zero energy is dying;
- Cell cannot have more than `MAX_GENES` non-zero genes in ruleset.
Additional rule is: board has torus topology.

So, if all cells initially has B3/S23 ruleset, DEATH_SPEED = BIRTH_COST = 0, MAX_GENES >= 3, we have exact Conway rules.
But if there were more than one ruleset initially, evolution may begin.
There are 2^18 possible rulesets, only a small fraction of which have been
studied in any detail. So, who knows what we may discover with evolutionary rules :)

CONTROLS:
Arrows    move field
+/-       zoom in/out
]/[       speed up/down
F         toggle fullscreen
S         dump board state to a file
Q/ESC     quit

Prerequisites: pycuda, numpy, scipy, pygame, scikit-image
Debian: apt-get install python-pycuda python-numpy python-pygame python-scipy python-setuptools

Author: a5kin
Copyright: MIT License.

"""

import sys, time, math, colorsys, random, traceback
import pygame
from pygame.locals import *
import numpy as np
from scipy.misc import imsave
import scipy.ndimage.interpolation
from skimage import transform as tf
import importlib

import pycuda.driver as drv
import pycuda.tools
import pycuda.autoinit
from pycuda.compiler import SourceModule
import pycuda.gpuarray as gpuarray
from pycuda.elementwise import ElementwiseKernel

try:
    expmod = importlib.import_module('experiments2.' + sys.argv[1])
    DEATH_SPEED = expmod.DEATH_SPEED
    BIRTH_COST = expmod.BIRTH_COST
    MAX_GENES = expmod.MAX_GENES
    FIELD_WIDTH = expmod.FIELD_WIDTH
    FIELD_HEIGHT = expmod.FIELD_HEIGHT
    SAVE_FRAMES = expmod.SAVE_FRAMES
    DOWNSCALE_FACTOR = expmod.DOWNSCALE_FACTOR
    FRAME_SKIP = expmod.FRAME_SKIP
    RANDOM_SEED = expmod.RANDOM_SEED
    FADE_IN = expmod.FADE_IN
    FADE_OUT = expmod.FADE_OUT
    fld_init = expmod.fld_init
except (ImportError, IndexError):
    print "No experiment preset found, loading default (big_bang)."
    DEATH_SPEED = 0
    BIRTH_COST = 0
    MAX_GENES = 9
    FIELD_WIDTH = 1280
    FIELD_HEIGHT = 720
    SAVE_FRAMES = False
    DOWNSCALE_FACTOR = 1
    FRAME_SKIP = 1
    RANDOM_SEED = None
    FADE_IN = 6
    FADE_OUT = 6
    def fld_init(a):
        return np.asarray([[(random.choice([0, 1]) * random.randint(0, 256*512) if (i < 100 and j < 100) else 0) for j in range(a.height)] for i in range(a.width)]).astype(np.int32)
except:
    print traceback.format_exc()
    sys.exit(0)
    

step_gpu = ElementwiseKernel("unsigned int *fld, unsigned int *fld_new, unsigned int *seeds, unsigned int *bufs, unsigned int *img, int w, int h", """
    int x = i / h;
    int y = i % h;
    // torus topology emulation
    int xm1 = x - 1; if (xm1 < 0) xm1 = w + xm1;
    int xp1 = x + 1; if (xp1 >= w) xp1 = xp1 - w;
    int ym1 = y - 1; if (ym1 < 0) ym1 = h + ym1;
    int yp1 = y + 1; if (yp1 >= h) yp1 = yp1 - h;
    // cache neighbours values
    uint f0 = fld[i];
    uint f1 = fld[xm1 * h + ym1];
    uint f2 = fld[x * h + ym1];
    uint f3 = fld[xp1 * h + ym1];
    uint f4 = fld[xm1 * h + y];
    uint f5 = fld[xp1 * h + y];
    uint f6 = fld[xm1 * h + yp1];
    uint f7 = fld[x * h + yp1];
    uint f8 = fld[xp1 * h + yp1];
    uint energy = (f0 >> 17);
    // total number of neighbours
    int N = EXISTS(f1) + EXISTS(f2) + EXISTS(f3) + EXISTS(f4) +
            EXISTS(f5) + EXISTS(f6) + EXISTS(f7) + EXISTS(f8);
    if (energy >= 0xff || N == 0 || f0 > 0 && (((f0 >> 8) & (1 << N)) == 0)) {
        // cell is dying
        fld_new[i] = 0;
        //img[i] = fadeout(img0, 5);
    } else {
        uint f00 = f0;
        for (int ni = 8; ni > 0; ni--) {
            // no re-occupation rule, breeding in empty cells only
            //if (f0 > 0) break;
            // cache neighbours breeding fitnesses 
            int ff1 = FIT(f1, ni);
            int ff2 = FIT(f2, ni);
            int ff3 = FIT(f3, ni);
            int ff4 = FIT(f4, ni);
            int ff5 = FIT(f5, ni);
            int ff6 = FIT(f6, ni);
            int ff7 = FIT(f7, ni);
            int ff8 = FIT(f8, ni);
            if (ff1 + ff2 + ff3 + ff4 + ff5 + ff6 + ff7 + ff8 == ni) {
                // neighbours able to breed, cell is born
                f0 = 0;
                uint gene_num = 0;
                //int genes_count = {2};
                //int gene;
                uint nit = (int) (ni / 2);
                uint seed = seeds[i];
                uint nonzero_genes_num = 0;
                while (gene_num < 17) {
                    // pseudorandom cross breeding
                    uint rng = ((((seed + gene_num) * 58321) + 11113)) % 65535;
                    uint fg1 = (f1 >> gene_num) & ff1;
                    uint fg2 = (f2 >> gene_num) & ff2;
                    uint fg3 = (f3 >> gene_num) & ff3;
                    uint fg4 = (f4 >> gene_num) & ff4;
                    uint fg5 = (f5 >> gene_num) & ff5;
                    uint fg6 = (f6 >> gene_num) & ff6;
                    uint fg7 = (f7 >> gene_num) & ff7;
                    uint fg8 = (f8 >> gene_num) & ff8;
                    int n1 = fg1 + fg2 + fg3 + fg4 + fg5 + fg6 + fg7 + fg8;

                    //if ((int) (n1 * 65535 / ni) < 65535 && (int) (n1 * 65535 / ni) > 0)
                    //    printf("%d %d | ", rng, (int) (n1 * 65535 / ni));
                    //if (n1 > nit) {
                    if ((int) (n1 * 65535 / ni) > rng) {
                        f0 += 1 << gene_num;
                        nonzero_genes_num += 1;
                        if ({1}) {
                            if (fg1) atomicAdd(&bufs[xm1 * h + ym1], ({1} << 17));
                            if (fg2) atomicAdd(&bufs[x * h + ym1], ({1} << 17));
                            if (fg3) atomicAdd(&bufs[xp1 * h + ym1], ({1} << 17));
                            if (fg4) atomicAdd(&bufs[xm1 * h + y], ({1} << 17));
                            if (fg5) atomicAdd(&bufs[xp1 * h + y], ({1} << 17));
                            if (fg6) atomicAdd(&bufs[xm1 * h + yp1], ({1} << 17));