Cellular Automata with Modular Arithmetic: Advanced Munching Squares

xorgen2.py

#    Cellular Automata with Modular Arithmetic: Advanced Munching Squares
#    Copyright (c) 2013, Sergey Shishmintzev
#    License: MPL/GPL
#    Interpreter: Python 2.7.3 + numpy 1.6.2 + PIL 1.1.7
#
#    P.S. Sorry for my English. It isn't my mother tongue.
#
#    Reference:
#        [Khan1997] A.R. Khan, et.al., "VLSI architecture of a cellular automata machine", 1997, 
#                   http://dx.doi.org/10.1016/S0898-1221(97)00021-7
    

import numpy as np


def bitreverse(n, l):
    # Reversing a binary number
    r = 0
    for _ in range(l):
        r = 2*r + (n&1)
        n /= 2
    return r
# Self-check
assert bin(bitreverse(1, 4))=='0b1000'
assert bin(bitreverse(2, 4))=='0b100'
assert bin(bitreverse(5, 4))=='0b1010'


def modulo_solve4(x, y, z, a, modulo):
    # Solve: x + y + z + ? == a (mod modulo).
    #return x^y^z^a # bitwise, modulo 2
    return (a-x-y-z) % modulo
# Self-check
assert modulo_solve4(1, 2, 3, (1+2+3+4)%10, modulo=10)==4
assert modulo_solve4(1, 2, 3, (1+2+3+9)%10, modulo=10)==9
assert modulo_solve4(9, 2, 3, (9+2+3+1)%10, modulo=10)==1
assert modulo_solve4(4, 5, 6, (4+5+6+7)%10, modulo=10)==7


def undo_rule15(unk_state, next_state, i0, j0, modulo):
    """
    Find the a fourth value, when three values are known,
    satisfying a condition:
    unk_state[i0, j0] + unk_state[i0, j0+1] + unk_state[i0+1, j0] + unk_state[i0+1, j0+1] == next_state[i0, j0] (mod modulo).

    In the other words, find the state of the 2D additive cellular 
    automata using a known next state after the applying rule 15
    (rule numbering from [Khan1997]). State to find is a partially known.


    Input values:
        unk_state - partially known state;
        next_state - known next state, after applying rule 15;
        i0, j0 - working position index;
        modulo - modulo of the operators.

    Changed values:
        unk_state - state with zero or one solution added (solution will be added 
                    at the some adjacent position near working position index).
    """

    assert unk_state.ndim == 2
    assert next_state.shape == unk_state.shape

    # additional indices for the adjacent values
    i1 = (i0 + 1) % unk_state.shape[0]
    j1 = (j0 + 1) % unk_state.shape[1]

    # prefetch a partially known state
    s00 = unk_state[i0, j0]
    s01 = unk_state[i0, j1]
    s10 = unk_state[i1, j0]
    s11 = unk_state[i1, j1]

    # prefetch a next state
    n00 = next_state[i0, j0]

    # solve (s00 + s01 + s10 + s11) % modulo == n00

    # Walk over adjacent values to find an unknown one
    if s00 is not None and s01 is not None and s10 is not None and s11 is None:
        unk_state[i1, j1] = modulo_solve4(s00, s01, s10, n00, modulo)
    elif s00 is not None and s01 is not None and s10 is not None and s11 is not None:
        pass # all adjacent values are known
    elif s00 is not None and s01 is not None and s10 is None and s11 is not None:
        unk_state[i1, j0] = modulo_solve4(s11, s00, s01, n00, modulo)
    elif s00 is not None and s01 is None and s10 is not None and s11 is not None:
        unk_state[i0, j1] = modulo_solve4(s10, s11, s00, n00, modulo)
    elif s00 is None and s01 is not None and s10 is not None and s11 is not None:
        unk_state[i0, j0] = modulo_solve4(s01, s10, s11, n00, modulo)
    else:
        assert False, "less than three adjacent values are known"


def bitwise_rule33(x):
    """
    Perform rule 33 using bitwise operator [Khan1997].
    (it is rule 60 in 1D for every row http://mathworld.wolfram.com/Rule60.html )
    """
    #return x^np.roll(x, 1, axis=1) # considerably slow

    assert x.ndim == 2

    X = np.empty(x.shape, dtype=int)
    I = x.shape[0]
    J = x.shape[1]
    for i in xrange(I):
        for j in xrange(J):
            X[i, j] = x[i, (j - 1) % J] ^ x[i, j]
    return X


def normalize(x):
    # Scale all array values into range [0, 255]
    mn = np.min(x)
    mx = np.max(x)
    if mx - mn != 0:
        X = 255*(x - mn) / (mx - mn)
    else:
        X = np.zeros(x.shape)
    return X


_last_mean = 0
def mean_split(x):
    """
    When modulo is large (like >10000), values of the state after 
    rule 33 are splitted into two or more noticeable close parts.
    When modulo is small - split by mean as failback.
    """
    global _last_mean
    l = list(set(x.flat))
    mean = max(np.mean(l), _last_mean)
    _last_mean = mean # hack to preventing the oscilation
    #l.sort(); print len(l), np.mean(l), (l)
    #l = np.asarray(l); m = np.argmax(l[1:]-l[:-1]); print m, l[m], l[m+1]
    top = x.copy()
    bot = x.copy()
    top[x<mean] = 0
    bot[x>=mean] = 0
    #print top
    #print bot
    return top, bot


def makecolor2(a, b):
    # Convert two number arrays into the array of the bi-color RGB tuples.
    assert a.ndim == 2
    assert a.shape == b.shape

    color = np.empty(a.shape, dtype=tuple)

    for i in xrange(a.shape[0]):
        for j in xrange(a.shape[1]):
            color[i, j] = (a[i, j], (a[i, j]+b[i, j])/2, b[i, j])

    return color


if __name__ == '__main__':
    from PIL import Image

    # Image size
    power = 10

    size = 2**power
    half = 2**(power-1)

    image = Image.new("RGB", (size, size), 0)
    pix = image.load()

    # known state after applying rule 15
    next_state = np.ones((size, size), int)

    for modulo in xrange(64, 6135):
        print "modulo =", modulo

        # partially known state
        state = np.empty((size, size), object)

        # fill initial (known) values
        for i in xrange(size):
            state[i, 0] = i % modulo
            state[0, i] = i % modulo

        # find rest values of the state
        for i in xrange(size):
            for j in xrange(size):
                undo_rule15(state, next_state, i, j, modulo)

        # the magic pass
        state2 = bitwise_rule33(np.asarray(state, dtype=int))

        # color magic
        a, b = mean_split(state2 % modulo)
        #a, b = mean_split(state2) # alternative view
        c = makecolor2(normalize(a), normalize(b))

        # produce image
        for x in xrange(size):
            for y in xrange(size):
                pix[x, y] = c[y, x]

        image.save("img-p%d-m%05d.png"%(power, modulo), "png")