# 1D-randomwalk: A walker makes several steps,

# with a given number of walks pr. trial.

# It computes the entropy by filling in bins with counts

​

import numpy, sys, math

​

def mc_trial(number_walks,move_probability,walk_cum,walk2_cum, probability):

    """

    Do a MonteCarlo trial, that is,

    random-walk one particle.

​

    Input:

    - number_walks:     Number of steps to walk the particle

    - move_probability: Probability that the particle

                        will step right when doing a step

    - walk_cum:         Numpy-array of length number_walks + 1,

                        containing the sum of the position

                        of the particles as a function of time

                        (usefull to calculate mean pos. as a function

                        of time)

    - walk2_cum:        Same as walk_cum, but with the sum of the

                        positions squared

    - probability:      Number of times each gridpoint is hit

​

    Output: As walk_cum, walk2_cum, and probability are (pointers to)

    numpy arrays, they are altered also in the calling function.

    """

    #Initial pos. As walk_cum[0]=walk2_cum[0] = 0.0

    #by initialization, it is uneccessary to add this step to

    #the arrays...

    pos = 0;

    for walk in range(number_walks+1):

        if numpy.random.random() <= move_probability:

            pos += 1

        else:

            pos -= 1

        walk_cum[walk]   += pos

        walk2_cum[walk]  += pos**2

        #Zero-position of the array is the leftmost

        #end of the grid

        probability[pos+number_walks] += 1

​

def mc_sample(length,trials, number_walks, move_probability):

    """

    Generate the probability distribution for finding

    a walker at a gridpoint, after a number of walks on a

    1d lattice with wrap-around boundary conditions

​

    Input:

    - length: Lattice-points away from x=0

    - trials: Number of MonteCarlo trials (number of walkers)

    - move_probability: Probability of moving right

​

    Output:

    Normalized probability of finding a walker on a

    specific grid position

    """

​

    #Grid position of every walker

    x = numpy.zeros(trials,numpy.int)

​

    #Loop over timesteps and walkers,

    #and find the walkers "ending positions"

    for t in range(number_walks):

        for i in range(trials):

            if numpy.random.random() <= move_probability:

                x[i] += 1

                #Wraparound?

                if x[i] > length:

                    x[i] = -length

            else:

                x[i] -= 1

                if x[i] < -length:

                    x[i] = +length

    #Calculate the probability of finding a walker

    #each grid-position

    probability = numpy.zeros(2*length+1)

    for i in range(len(probability)):

        pos = i-length

        #count number of occurences of this pos i x array

        count = 0

        for j in range(len(x)):

            if x[j] == pos:

                count += 1

        #Normalize and save

        probability[i] = count/float(trials)

    return probability

​

#Main program

​

length               =   10

trials               =   100000

number_walks         =   100

move_probability     =  0.5

​

#Do the MC

probability = mc_sample(length,trials,number_walks,move_probability);

​

#Not reliable: ln(0)

#entropy = - numpy.sum(probability*numpy.log(probability))

​

entropy = 0.0

for i in range(len(probability)):

    if probability[i] > 0.0:

        entropy -= probability[i]*math.log(probability[i])

​

print("Timesteps             =",number_walks)

print("Walkers (num. trials) =",trials)

print("Entropy               =",entropy)

​

if len(probability) <= 101:

    print("Probability distribution (Flat => high entropy):")

    print(probability)

else:

    print("Probability distribution to big to print")