Markov Chain Monte Carlo Simulation Prooblem

Question

I'm trying to run a MC simulator for a Markov Chain that is uniformly distributed among all NxN matrices that have no neighboring 1's. My algo is supposed to fill up the state space by running the chain a bunch of times. However there's something horribly wrong with my logic somewhere and the state space just isn't filling up. Any help would be greatly appreciated. Here is my code.

import random
import numpy

M=numpy.zeros((52,52),dtype=int)
z=0
State_Space=[]

for i in range(1,100):
    x=random.randint(1,50)
    y=random.randint(1,50)

    T=M
    if T[x][y]==1:
        T[x][y]=0
    if T[x][y]==0:
        T[x][y]=1


    if T not in State_Space:
    if T[x+1][y+1]==0 and T[x+1][y-1]==0 and T[x-1][y-1]==0 and T[x-1][y+1]==0:
        State_Space.append(T)
        M=T


    else:
         if T[x+1][y+1]==0 and T[x+1][y-1]==0 and T[x-1][y-1]==0 and T[x-1][y+1]==0:
            M=T
    print State_Space

Answer 1

I notice two things:

First in line 12 you have T=M and I assume you want T=M.copy(). Doing T=M makes T and M reference the same matrix, so changing a value in T will affect M also. If you assign a copy of M to T then this won't happen.

Second, T not in State_Space is not checking for T in the State_Space array. Because of how numpy indexing works, the in operator cannot be used for arrays. If you tried T in State_Space with a non-empty State_Space you would get a ValueError about truth value ambiguity. Instead you need to check if any element of State_Space is equal to T. We should use if any(numpy.array_equal(T, X) for X in State_Space):

In the end, my code looks like this:

import random
import numpy

M=numpy.zeros((52,52),dtype=int)
z=0
State_Space=[]

for i in range(1,100):
    x=random.randint(1,50)
    y=random.randint(1,50)

    T=M.copy()
    if T[x][y]==1:
        T[x][y]=0
    if T[x][y]==0:
        T[x][y]=1

    if not any(numpy.array_equal(T, X) for X in State_Space):
        if T[x+1][y+1]==0 and T[x+1][y-1]==0 and T[x-1][y-1]==0 and T[x-1][y+1]==0:
            State_Space.append(T)
            M=T
    else:
        if T[x+1][y+1]==0 and T[x+1][y-1]==0 and T[x-1][y-1]==0 and T[x-1][y+1]==0:
            M=T
    print len(State_Space)

After running, I have ~90 entries in State_Space.