Python to R translation-CodePudding

I want to implement the Ising model in R :

My attempt is :

L=4
N=L*L;N
A= matrix(nr= N, nc=L, 0);A
for (i in 0:N-1){
  A[i 1,1] = (i %/% L) * L   (i   1) %% L
  A[i 1,2] = (i   L) %% N
  A[i 1,3] = (i %/% L) * L   (i - 1) %% L
  A[i 1,4] = (i - L) %% N
}

Temperature = 100
Spin = sample(c(1,-1),N,replace=TRUE)
spin_table = matrix(NA,L,L);spin_table
nsteps = 100 *L*L
for(step in 1:nsteps){
  for(n in A){
    k = sample(1,N-1,replace=TRUE)
    delta_E = 2.0 * Spin[k] * sum(Spin[n])
    if(runif(1) < exp(-delta_E/Temperature)){
      Spin[k] = Spin[k]*-1
    }
    }}
  

for(k in 1:N-1){
  y = k %/% L
  x = k %% L
  spin_table[x,y] = Spin[k]}
spin_table

But I have several problems implement it.I have mess up R code with some pseudocode. Can someone help me out with this ?

I want the resulting plot to be with ggplot with scale_fill_viridis(option="D"). Any help ?

Python

in Python this model is :

L = 64
N = L*L

#   // floor division
#   %  modulus
neighbors = {i : ((i // L) * L   (i   1) % L,   # RIGHT
                 (i   L) % N,                   # DOWN
                 (i // L) * L   (i - 1) % L,    # LEFT
                 (i - L) % N)                   # UP
                 for i in range(N)}

Temperature = 1.0
Spin = [np.random.choice([1, -1]) for k in range(N)]
spin_table = [[None for x in range(L)] for y in range(L)]

nsteps = N * 100
for step in range(nsteps):
    k = np.random.randint(0, N - 1)
    delta_E = 2 * Spin[k] * sum(Spin[n] for n in neighbors[k])
    if np.random.uniform(0.0, 1.0) < np.exp(-delta_E/Temperature):
        Spin[k] *=  1

for k in range(N):
    y = k // L
    x = k % L
    spin_table[x][y] = Spin[k]

pylab.close()
pylab.imshow(spin_table,extent=[0,L,0,L],interpolation='nearest')
pylab.title("Temperature: " str(Temperature) ", Grid Size: " str(L))
pylab.show()

CodePudding user response：

Note that in

if(runif(1) < exp(-delta_E/Temperature))

The if condition is a logical vector of length N. This is ambiguous, since you might mean that all vector elements should be TRUE or any element should be TRUE. Depending on what applies, you have to write:

if(all(runif(1) < exp(-delta_E/Temperature)))

if(any(runif(1) < exp(-delta_E/Temperature)))

Assuming the former, you can write:

L=4
N=L*L;N
A= matrix(nr= N, nc=L, 0);A
for (i in 0:N-1){
  A[i 1,1] = (i %/% L) * L   (i   1) %% L
  A[i 1,2] = (i   L) %% N
  A[i 1,3] = (i %/% L) * L   (i - 1) %% L
  A[i 1,4] = (i - L) %% N
}
 
Temperature = 100
Spin = sample(c(1,-1),N,replace=TRUE)
spin_table = matrix(NA,L,L);spin_table
nsteps = 100 *L*L
for(step in 1:nsteps){
  for(n in A){
    k = sample(1,N-1,replace=TRUE)
    delta_E = 2.0 * Spin[k] * sum(Spin[n])
    if(all(runif(1) < exp(-delta_E/Temperature))){
      Spin[k] = Spin[k]*-1
    }
  }
}
  
 
for(k in 1:N-1){
  y = k %/% L
  x = k %% L
  spin_table[x,y] = Spin[k]}
spin_table

The code does not break now. Whether this is a correct implementation of your model is another story.

P.S. If you are new to R, it might be worth noting that

0:N-1

is equivalent to

(0:N)-1

and not

0:(N-1)

As a second note, in R, array indices are not zero-based.

If v is a vector, the assignment

v[0] <- 100

writes in a sort of empty location, which leaves v unchanged, but you do not get any error (<- is the stylistically correct way to make assignments).