Self Avoiding Walk in Julia

Question

I am trying to code the self avoiding walk in Julia, here is my code:
Here I defined the position function in the NxN matrix

function position(vertical, horizontal, N)
            
    ud = BitArray(a == vertical  for a = 1:N )              # Goes up or down
    lr = transpose(BitArray(a == horizontal  for a = 1:N )) # Goes left or right

    return ud * lr
end

This is the code of the walker that describes the movement inside a Matrix, where 1 means that the walker has been on that position and 0 where it has never been. (Furthermore, the walker is not allowed to move diagonally. This part I managed to code correctly)

using StatsBase, LinearAlgebra

function SAW(N, Iterations=5) # SAW, Self Avoiding Walk
    
    x = ceil(N/2) 
    y = ceil(N/2)
    Mat = zeros(Int8, N , N)
    pos = zeros(Int8, N , N)
       
    for i in 1:Iterations
        
        RandomOnes = sample(-1:2:1, 2, replace = false)[1]  #gives the values -1 or 1
        ZeroOne = sample(0:1, 1, replace = false)[1]   #gives the values 0 or 1
    
        
        # In order to avoid moving diagonally, only one of the two variables (x,y)
        # can have the values 1 or -1, the other one has to be 0 (e.g. if x = 1 then y = 0)       
        if ZeroOne == 1 
            x += RandomOnes
            y += 0
        else
            x += 0
            y += RandomOnes
        end
               
        Mat += position(y,x,N)
        pos  = position(y,x,N)
    end
    
    for i in 1:N
        for j in 1:N
            if Mat[i,j] > 1
                Mat[i,j] += 1 # Here the process should restart, because the walker ran into itself
            end
        end
    end
                
    return Mat, pos
end

SAW(7)[1]

I am not sure how to restart the process when the walker runs into itself

Answer 1

Here is something that you maybe will find useful. This is a solution that assumes that your lattice might be very large (so that storing a matrix as in your approach will not fit into memory). I have not optimized it for speed but for simplicity:

function SAW(N::Integer, Iterations::Integer)
    @assert N >= 1
    @assert Iterations >= 0
    # this will store the path we have traversed
    path = Tuple{Int, Int}[]
    cur = ((N+1) ÷ 2, (N+1) ÷ 2)
    push!(path, cur)
    for i in 1:Iterations
        x, y = cur
        # now check the four possible movement directions
        # store moves that are allowed in valid vector
        valid = Tuple{Int, Int}[]
        if x > 1 && !((x-1, y) in path)
            push!(valid, (x-1, y))
        end
        if x < N && !((x+1, y) in path)
            push!(valid, (x+1, y))
        end
        if y > 1 && !((x, y-1) in path)
            push!(valid, (x, y-1))
        end
        if y < N && !((x, y+1) in path)
            push!(valid, (x, y+1))
        end
        if isempty(valid)
            # we are stuck and failed to generate a walk of length N
            # no moves are possible and we still have not made N moves
            return nothing
        end
        # randomly pick a move from valid moves
        cur = rand(valid)
        push!(path, cur)
    end
    # we have successfully generated the walk - return it
    return path
end

Clearly this code can be easily made much faster but I did not want to complicate it. Also you might want to have different distribution of the generated walks (i.e. this code can potentially generate all possible walks but the simulation procedure assumes one specific distribution over them) - again this can be changed.

EDIT:

The simplest way to store the results in a matrix is:

julia> N = 10
10

julia> Iterations = 5
5

julia> path = SAW(N, Iterations)
6-element Vector{Tuple{Int64, Int64}}:
 (5, 5)
 (4, 5)
 (3, 5)
 (3, 6)
 (3, 7)
 (4, 7)

julia> m = zeros(Int, N, N)
10×10 Matrix{Int64}:
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0

julia> foreach(p -> m[p...] = 1, path)

julia> m
10×10 Matrix{Int64}:
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  1  1  1  0  0  0
 0  0  0  0  1  0  1  0  0  0
 0  0  0  0  1  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0

(assuming N is small)