# http://www.python.org/doc/current/lib/lib.html
#
#
# find that for 200 particles on 42x42, the autocorrelation time for the number of
# ups is about 60,000
# So to get good histogram, need to run for thousands of times that.
from Numeric import *
import random
import Gnuplot, Gnuplot.funcutils
import math, time, tempfile
# define some global variables
a=arange(0)
## a =array[0]
NY=5
NX=5
L2=5
x=10; y=10; u=0
# define two orientations
dy=[0,1] ; dx=[1,0]
# create some interesting state variables
ho = 0
up = 0
state = {}
# histogram of value of K and ho
histo = {}
# log of value of "up"
uplog = {}
def simulate(Y,X,K,T,tperiod,keepteleporters,verbose,bias,L): # simulate monte carlo moves of needles of length L on grid
global a,NX,NY,L2 ,x,y,u , uplog
make_u_log = 1
ENGINES=5
L2 = (L-1)/2 # intend L to be odd so that L2 is an integer.
# valid values for x,y are from 0 to NX-1
# a={} # a is going to be a global dictionary
NX=X;NY=Y
# create an empty grid
a = zeros( [NX,NY] )
# state is going be a list of K lists of object coordinates [ [x,y,u], [x,y,u], [x,y,u] ]
# where K is possibly going to vary
# instantiate K objects
print "trying to create ", K*L, "in ", NY*NX
if( K*L > NY*NX ):
print("impossible!")
exit(0)
# endif
failedattempts=0
maxfailedattempts=10000
## declare x,y,u so that they are global from the point of view of the insert that's coming.
for k in range(0,K):
if(verbose):
print "Making number ",k
invalid=1
while ( (invalid) and (failedattempts0) and (not(t%tperiod)) ) :
gdisplay(t) ;
## print "=== K=",K," | ",up, " - ",ho," ==="
if (verbose):
print state
raw_input('Please press return to continue...\n')
# endif
if ( make_u_log ) :
uplog[t] = up
# endif
if (K,up) in histo.keys() :
histo[(K,up)] += 1 ;
else:
histo[(K,up)] = 1 ;
# endif
# endfor
if ( make_u_log ) :
# print it out to file
f = open("uplog", 'w')
for t in sort(uplog.keys()):
f.write('%s %s\n' % (t,uplog[t]))
f.close()
# endif
return histo
# ENDFUNCTION
def vacate(x,y,u):
global up , ho
for l in range (-L2,L2+1):
tx=x+l*dx[u]
ty=y+l*dy[u]
fillit(tx,ty,0)
if( u == 1 ):
ho -= 1 ;
else:
up -= 1 ;
#endif
def checka(x,y,u):
global a, L2, dx, dy, NX, NY
invalid=0
for l in range (-L2,L2+1):
tx=x+l*dx[u]
ty=y+l*dy[u]
if( a[ty%NY][tx%NX] ):
invalid=1
break ;
# endif
# endfor
return invalid
# enddef
## I want x,y,u to be overwritten
def biasedrandominsert(bias):
global x , y , u
x = random.randint(0,NX-1)
y = random.randint(0,NY-1)
ur = random.random()
if(ur>bias):
u = 1
else:
u = 0
## endif
# check if all L spots around (x,y) are empty. if ok, return invalid = 0 , else set invalid=1
return checka(x,y,u)
# enddef
def occupy(x,y,u,k): # position(x,y) and orientation u
global dx, dy, state, ho, up, L2
for l in range (-L2,L2+1):
tx=x+l*dx[u]
ty=y+l*dy[u]
uu = 2-u # map u to 1,2 instead of 0,1
if ( l==-L2) :
tu=uu*10
elif ( l",
20:"^",
21:"|",
22:"v"}
empty = 0
for x in range(0,NX):
for y in range(0,NY):
ta = a[y][x]
if (ta==0):
empty += 1
if (ta in dict.keys() ):
print dict[ta],
else:
print "?",
# endif
# endfor
print
# endfor
print "=== ",empty," vacancies"
filename1 = tempfile.mktemp()
def gdisplay(t):
global filename1 , state , L2 , dx , g , NX , NY , up , ho
d=0.4 ## half-thickness of rectangle
K = max ( state.keys() ) + 1
g.title('%s: K = %s | %s - %s' % (t,K,up,ho) )
g('set xrange [-0.5:%s] ' % (0.5+NX) )
g('set yrange [-0.5:%s] ' % (0.5+NY) )
f = open(filename1, 'w')
for k in state.keys() :
sk = state[k]
x=sk[0];y=sk[1];u=sk[2]
if (0): ## simple line
f.write('%s %s\n' % (x+dx[u]*L2,y+dy[u]*L2) )
f.write('%s %s\n' % (x-dx[u]*L2,y-dy[u]*L2) )
else: ## rectangles
f.write('%s %s\n' % (x+dx[u]*L2+d,y+dy[u]*L2+d) )
f.write('%s %s\n' % (x+dx[u]*L2+d,y-dy[u]*L2-d) )
f.write('%s %s\n' % (x-dx[u]*L2-d,y-dy[u]*L2-d) )
f.write('%s %s\n' % (x-dx[u]*L2-d,y+dy[u]*L2+d) )
f.write('%s %s\n' % (x+dx[u]*L2+d,y+dy[u]*L2+d) )
f.write('\n' )
f.close()
## print "plotting ", filename1
g.plot(Gnuplot.File(filename1))
## raw_input('Please press return\n')
# debug=1 prints gnuplot commands as they are run
g = Gnuplot.Gnuplot(debug=0)
g.title('A simple example') # (optional)
g('set data style lines') # give gnuplot an arbitrary command
g('set nokey') # give gnuplot an arbitrary command
# Make a temporary file:
f = open(filename1, 'w')
for x in arange(100)/5. - 10.:
f.write('%s %s %s\n' % (x, math.cos(x), math.sin(x)))
f.close()
# ensure that file will be deleted upon exit:
# file1 = Gnuplot.PlotItems.TempFile(filename1)
print '############### test File ########################################'
#raw_input('Please press return to start monte carlo...\n')
g.plot(Gnuplot.File(filename1))
#raw_input('Please press return to start monte carlo...\n')
# g.plot(Gnuplot.File(file1))