Multi-Layered XOR Pattern Motion Generator

xorgen.py

#    Multi-Layered XOR Pattern Motion Generator
#    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.


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_rule225(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 225 
    (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 225;
        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).

    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
    """

    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_rule60(x):
    """
    Perform rule 60 ( http://mathworld.wolfram.com/Rule60.html ) 
    for every row using bitwise operators (only modulo 2.)
    """
    #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, dtype=x.dtype)
    return X

_last_mean = 0
def mean_split(x):
    """
    When modulo is large (like >3000), values of the state after 
    rule 60 are splitted into two 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)
    top = x.copy()
    bot = x.copy()
    #print top
    #print bot
    top[x<mean] = 0
    bot[x>=mean] = 0
    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 = 9

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

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

    # known state after applying rule 225
    br_pattern_state = np.zeros((size, size), int)
    for i in xrange(size):
        ri = bitreverse(i, power)
        br_pattern_state[ri, i] = 1

    for modulo in xrange(32, 1024*3+1): # image motion tiny cycle is 1024 (when power is 9)
        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_rule225(state, br_pattern_state, i, j, modulo)

        # the magic pass
        state2 = bitwise_rule60(state)

        # 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")