% make the adjacency graph of a random walk on a square grid.
% the vertices of the graph are the vertices of the grid that
% have been visited and two of these are adjacent if the random
% walk stepped from one to the other at some time.
% The argument N is the number of steps taken by the walk.

function  [dir_edges,out_list,verts]=make_trace(N)
   % side length of box
   S = 2*ceil(sqrt(N)*log(log(N)));
   % visited records time of first visit to a grid position
   visited=zeros(2*S,2*S);
   right_edges=zeros(2*S,2*S);
   up_edges=zeros(2*S,2*S);
   dir_edges=zeros(N,2);
   % G is the adjacency matrix of visited vertices
   G=sparse(N,N);
   start=[S,S];
   % old_pos is current position of the walk
   old_pos=start;
   old_vert=1;
   % verts holds (x,y) coordinates of visited points
   verts=zeros(N,2);
   verts(1,:)=start;
   num_verts=1;
   visited(start(1),start(2))=1;
   old_V=1;
   % take N random steps
   edge_counter=0;
   for k=1:N
       % choose a random step among [0,1],[0,-1],[1,0], [-1,0]
       dir=randi(4)-1;
       z=i^dir;
       step=[real(z),imag(z)];
       % new_pos is the new coordinates of the walk
       new_pos=old_pos+step;
       % check if this vertex has been visited before
       old_V=visited(old_pos(1),old_pos(2));
       V=visited(new_pos(1),new_pos(2));
       if V==0
           % if not previously visited, add to verts
           num_verts=num_verts+1;
           visited(new_pos(1),new_pos(2))=num_verts;
           verts(num_verts,:)=new_pos;
           new_vert=num_verts;
           V=num_verts;
           % add edge to dir_edges
           dir_edges(edge_counter+1,:)=[old_V,V];
           dir_edges(edge_counter+2,:)=[V,old_V];
           edge_counter=edge_counter+2;
           switch dir
               case 0 
                   right_edges(old_pos(1),old_pos(2))=1;
               case 1
                   up_edges(old_pos(1),old_pos(2))=1;
               case 2
                   right_edges(new_pos(1),new_pos(2))=1;
               case 3
                   up_edges(new_pos(1),new_pos(2))=1;
           end % switch
       else
           % if it has been visited, check if edge has 
           % been used before; if not add it to dir_edges
               switch  dir
               case 0
                   if right_edges(old_pos(1),old_pos(2))==0
                      dir_edges(edge_counter+1,:)=[old_V,V];
                      dir_edges(edge_counter+2,:)=[V,old_V];
                      edge_counter=edge_counter+2;
                      right_edges(old_pos(1),old_pos(2))=1;
                   end % if 
               case 1
                   if up_edges(old_pos(1),old_pos(2))==0
                      dir_edges(edge_counter+1,:)=[old_V,V];
                      dir_edges(edge_counter+2,:)=[V,old_V];
                      edge_counter=edge_counter+2;
                      up_edges(old_pos(1),old_pos(2))=1;
                   end % if    
               case 2
                   if right_edges(new_pos(1),new_pos(2))==0
                      dir_edges(edge_counter+1,:)=[old_V,V];
                      dir_edges(edge_counter+2,:)=[V,old_V];
                      edge_counter=edge_counter+2;
                      right_edges(new_pos(1),new_pos(2))=1;
                   end % if
               case 3
                   if up_edges(new_pos(1),new_pos(2))==0
                      dir_edges(edge_counter+1,:)=[old_V,V];
                      dir_edges(edge_counter+2,:)=[V,old_V];
                      edge_counter=edge_counter+2;
                      up_edges(new_pos(1),new_pos(2))=1;
                   end % if
           end %switch
           new_vert=V;   
       end % if 
       old_pos=new_pos;
       old_vert=new_vert; 
       old_V=V;
   end % for k
   verts=verts(1:num_verts,:);
   dir_edges=dir_edges(1:edge_counter,:);
   
   % list directed edges leaving each vertex; 1st column is degree
   out_list=zeros(num_verts,5);
   for k=1:edge_counter
       v=dir_edges(k,1);
       c=out_list(v,1);
       out_list(v,c+2)=k;
       out_list(v,1)=c+1;
   end % for k
   %figure 
   %spy(visited)
   num_verts;
   %figure;hold on; axis equal; grid on;
   %plot_edge_list(dir_edges,verts)
return
