1

This is what I have imported:

import random
import matplotlib.pyplot as plt
from math import log, e, ceil, floor
import numpy as np
from numpy import arange,array
import pdb
from random import randint

Here I define the function matrix(p,m)

def matrix(p,m):            # A matrix with zeros everywhere, except in every entry in the middle of the row
    v = [0]*m
    v[(m+1)/2 - 1] = 1
    vv = array([v,]*p)
    return vv

ct = np.zeros(5)            # Here, I choose 5 cause I wanted to work with an example, but should be p in general  

Here I define MHops which basically takes the dimensions of the matrix, the matrix and the vector ct and gives me a new matrix mm and a new vector ct

def MHops(p,m,mm,ct):

k = 0
while k < p :                 # This 'spans' the rows
    i = 0
    while i < m :             # This 'spans' the columns
        if mm[k][i] == 0 :
            i+=1
        else:
            R = random.random()
            t = -log(1-R,e)              # Calculate time of the hopping
            ct[k] =  ct[k] + t
            r = random.random()
            if  0 <= r < 0.5 :              # particle hops right
                if 0 <= i < m-1:
                    mm[k][i] = 0
                    mm[k][i+1] = 1
                    break
                else: 
                    break           # Because it is at the boundary
            else:                            # particle hops left 
                if 0 < i <=m-1:    
                    mm[k][i] = 0
                    mm[k][i-1] = 1
                    break
                else:              # Because it is at the boundary
                    break
            break
    k+=1    
return (mm,ct)               # Gives me the new matrix showing the new position of the particles and a new vector of times, showing the times taken by each particle to hop

Now what I wanna do is iterating this process, but I wanna be able to visualize every step in a list. In short what I am doing is: 1. creating a matrix representing a lattice, where 0 means there is no particle in that slot and 1 means there is a particle there. 2. create a function MHops which simulate a random walk of one step and gives me the new matrix and a vector ct which shows the times at which the particles move.

Now I want to have a vector or an array where I have 2*n objects, i.e. the matrix mm and the vector ct for n iterations. I want the in a array, list or something like this cause I need to use them later on.

Here starts my problem:

I create an empty list, I use append to append items at every iteration of the while loop. However the result that I get is a list d with n equal objects coming from the last iteration!

Hence my function for the iteration is the following:

def rep_MHops(n,p,m,mm,ct):
    mat = mm
    cct = ct
    d = []
    i = 0
    while i < n :
        y = MHops(p,m,mat,cct)       # Calculate the hop, so y is a tuple y = (mm,ct)
        mat = y[0]                 # I reset mat and cct so that for the next iteration, I go further
        cct = y[1]
        d.append(mat)
        d.append(cct)
        i+=1
    return d


z = rep_MHops(3,5,5,matrix(5,5),ct)      #If you check this, it doesn't work
print z

However it doesn't work, I don't understand why. What I am doing is using MHops, then I want to set the new matrix and the new vector as those in the output of MHops and doing this again. However if you run this code, you will see that v works, i.e. the vector of the times increases and the matrix of the lattice change, however when I append this to d, d is basically a list of n equal objects, where the object are the last iteration.

What is my mistake? Furthermore if you have any coding advice for this code, they would be more than welcome, I am not sure this is an efficient way.

Just to let you understand better, I would like to use the final vector d in another function where first of all I pick a random time T, then I would basically check every odd entry (every ct) and hence check every entry of every ct and see if these numbers are less than or equal to T. If this happens, then the movement of the particle happened, otherwise it didn't. From this then I will try to visualize with matpotlibt the result with an histogram or something similar.

Is there anyone who knows how to run this kind of simulation in matlab? Do you think it would be easier?

Benjamin
  • 11,560
  • 13
  • 70
  • 119
Euler_Salter
  • 3,271
  • 8
  • 33
  • 74
  • I am unsure about the logic because I am not familiar with the problem but I have some tips for your code. lists can be appended to by using the `.append(object_to_append)` method as shown. In some cases it is useful to replace the items in the list `list[index] = some_value` however in your case there is no reason to instantiate a list with n indexes. just create an empty list and use the append method on each iteration rather than replacing a pre created value! – TheLazyScripter Jun 10 '16 at 15:38
  • Thank you for your comment, however I had already tried that. I used to have `d = []` and after `cct = y[1]` I had `d.append(mat)` and `d.append(cct)` but it didn't work!! – Euler_Salter Jun 10 '16 at 15:44
  • What I had was: `def rep_MHops(n,p,m,mm,ct):` `mat = mm` `cct = ct` `d = []` ` i = 0` `while i < n :` `y = MHops(p,m,mat,cct)` ` mat = y[0]` `cct = y[1]` `d.append(mat)` `d.append(cct)` ` i+=1` ` return d` – Euler_Salter Jun 10 '16 at 15:45
  • You are passing a list (actually a 2-level array) into `MHops`, which modifies it and returns it. Basically, it's a reference. So you get 5 copies of the same reference in your output list. A variation on http://stackoverflow.com/questions/1132941/least-astonishment-in-python-the-mutable-default-argument. – Corley Brigman Jun 10 '16 at 17:15
  • Thank you Corley Brigman. I have had a look at the question you suggested me to read, however I do not get how they are correlated. The guy in the question gets different answers everytime, while I get always the same. How should I change my code to make it work and have the desidered list d? And what is a reference? – Euler_Salter Jun 10 '16 at 17:22
  • @CorleyBrigman is right. You're passing and storing by references not copies, so on the next iteration of your loop `MHops` alters your previously stored version in `d`. Use `import copy; d.append(copy.deepcopy(mat))` to instead store a copy which won't be altered later. – RedCraig Jun 10 '16 at 17:43
  • @CorleyBrigman Oh okay I understand, so throughout the while loop in the definition what I'm doing is basically just modifying it and giving it a new name, so then when I append it to the list d, it's like I'm appending the same 'name' n times. But then how am I supposed to append these outputs to a list? – Euler_Salter Jun 10 '16 at 17:44
  • @RedCraig , it works, thank you a lot! How can I vote your comment? – Euler_Salter Jun 10 '16 at 17:47

1 Answers1

0

You're passing and storing by references not copies, so on the next iteration of your loop MHops alters your previously stored version in d. Use import copy; d.append(copy.deepcopy(mat)) to instead store a copy which won't be altered later.

Why?

Python is passing the list by reference, and every loop you're storing a reference to the same matrix object in d.

I had a look through python docs, and the only mention I can find is "how do i write a function with output parameters (call by reference)".

Here's a simpler example of your code:

def rep_MHops(mat_init):
    mat = mat_init
    d = []
    for i in range(5):
        mat = MHops(mat)
        d.append(mat)
    return d


def MHops(mat):
    mat[0] += 1
    return mat

mat_init = [10]
z = rep_MHops(mat_init)
print(z)

When run gives:

[[15], [15], [15], [15], [15]]

Python only passes mutable objects (such as lists) by reference. An integer isn't a mutable object, here's a slightly modified version of the above example which operates on a single integer:

def rep_MHops_simple(mat_init):
    mat = mat_init
    d = []
    for i in range(5):
        mat = MHops_simple(mat)
        d.append(mat)
    return d


def MHops_simple(mat):
    mat += 1
    return mat

z = rep_MHops_simple(mat_init=10)
print(z)

When run gives:

[11, 12, 13, 14, 15]

which is the behaviour you were expecting.

This SO answer How do I pass a variable by reference? explains it very well.

Community
  • 1
  • 1
RedCraig
  • 503
  • 1
  • 4
  • 15