# 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 (failedattempts<maxfailedattempts) ):
            # this sets x,y,u values
            invalid = biasedrandominsert(bias)
            failedattempts  = failedattempts + invalid 
        # endwhile

        if (invalid):
            #  just tried to make 'k' (which would be the k+1st)
            K=k
            print "Difficult to create! - Giving up with K = ",
            print K  
            break
        else:
            occupy(x,y,u,k)
        #endif
        if(verbose): 
            gdisplay(0)
            print state
            raw_input('Please press return to start monte carlo...\n')
    #endfor
  
    # OK, we have a starting state.
    gdisplay(0) 
  
    print "=== K=",K,"     | ",up, "     - ",ho," ===" 
    for t in range(0,T):

        # select random individual to mess with
        k = random.randint(0,K-1);
##        print "updating k = ",k
        sk = state[k] 
##        print sk
        #    [oldx,oldy,oldu] = sk;
        oldx=sk[0];oldy=sk[1];oldu=sk[2];
        x=oldx;y=oldy;u=oldu; # RESTORE
        # remove the individual.
        vacate(x,y,u) ; 
        alreadyinserted = 0 ; 
        
        engine = random.randint(1,ENGINES); 

        if (engine==1) :
            # try flipping orientation
            u = 1-u ; # maps 1<->0
            if checka(x,y,u) :
                x=oldx;y=oldy;u=oldu; # RESTORE
            elif (verbose):
                print "flipped"
            #            endif
               
        elif (engine==2) :
            #  try moving in long direction by step chosen slicily
            # first make a step of -6, -2, 2, or 6
            dist = random.randint(0,3) *4 -  6  
            x=(x+dist*dx[u])%NX   
            y=(y+dist*dy[u])%NY 
            if checka(x,y,u) :
                x=oldx;y=oldy;u=oldu; # RESTORE
            elif (verbose):
                print dist ;

        elif (engine==3):
            #  try moving in long direction by smaller step 
            dist = random.randint(0,1) *2 -  1  
            x=(x+dist*dx[u])%NX   
            y=(y+dist*dy[u])%NY 
            if checka(x,y,u) :
                x=oldx;y=oldy;u=oldu; # RESTORE
            elif (verbose):
                print dist ;
            # endif
        elif (engine==4) :
            #  try moving in the short direction by +/-1 
            dist = random.randint(0,1) *2 -  1  
            x=(x+dist*dx[1-u])%NX   
            y=(y+dist*dy[1-u])%NY 
            if checka(x,y,u) :
                x=oldx;y=oldy;u=oldu; # RESTORE
            elif (verbose):
                print "sideways" 
            # endif
        else:
            # try moving him somewhere totally new
            if biasedrandominsert(0.5) :
                x=oldx;y=oldy;u=oldu; # restorexyu ;
            else:
                if (verbose):
                    print "teleported"
                if (keepteleporters and (K<keepteleporters) ) :
                    # try to increase the density
                    # by retaining the old guy too, if it is valid
                    occupy(x,y,u,k) ;
                    x=oldx;y=oldy;u=oldu; # RESTORE old guy
                    if checka(x,y,u) :   # see if the old guy can be inserted 
                        # he can't, so can't have both at once
                        alreadyinserted = 1 
                    else:
                        # x,y,u now contains the _old_ guy's coordinates; he'll be inserted (as 'k') in a moment
                        K=K+1 ;
                        k=K-1 ;
                    # endif
	        # endif
            # endif

        # OK, we have a new state
        if (alreadyinserted == 0 ) :
            occupy(x,y,u,k) ;
        #     endif
        if( (tperiod>0) 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<L2):
            tu=uu*10+1
        else:
            tu=uu*10+2
        # endif
        fillit(tx,ty,tu)
    # endfor

    state[k] = [ x,y,u ] ;
    if( u == 1 ) :
        ho += 1 ;
    else:
        up += 1 ;
    # endif
# enddef

def fillit(tx,ty,tu):
    # makes use of the global array a
    global a,NY,NX
    a[ty%NY][tx%NX] = tu
# enddef


def display(t):
    global a,NX,NY
    dict={0:".",
          10:"<",
          11:"-",
          12:">",
          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))