0
$\begingroup$

I'm trying to simulate a self-avoiding persistent random walk with periodic boundary conditions a large number of times. For some reason, the script I've written is extremely slow for self-avoiding walks, for overlapping walks it performs quite fast. Are there any suggestions to speed up the process, I might be coding very inefficient.

The process is as follows; the first two steps can be randomly placed on the lattice, the steps that follow are biased to go in the "forwards" direction via the parameters pf, where $0<pf<1$.

Nwalks = 10^5;
Np=250;
pf=0.5;
step = {{-1, 0, 0}, {1, 0, 0}, {0, -1, 0}, {0, 1, 0}, {0, 0, -1}, {0, 
    0, 1}}; (*All posibilities to walk on the cubic lattice*)
l = 40;
lz = 75;
Until[Length[allSAW] == Nwalks,
  pol = {{RandomInteger[{0, l}], RandomInteger[{0, l}], 
     RandomInteger[{0, lz}]}}; (*Generate first 2 steps, 
  can be chosen at random*)
  next = pol[[-1]] + step[[RandomInteger[{1, 6}]]];
  next1 = next;
  If[next1[[1]] < 0, next1 = {l, next1[[2]], next1[[3]]}, 
   If[next1[[1]] > l, 
    next1 = {0, next1[[2]], 
      next1[[3]]}]]; (*apply periodic boundary conditions x*)
  If[next1[[2]] < 0, next1 = {next1[[1]], l, next1[[3]]}, 
   If[next1[[2]] > l, 
    next1 = {next1[[1]], 0, 
      next1[[3]]}]]; (*apply periodic boundary conditions y*)
  If[next1[[3]] < 0, next1 = {next1[[1]], next1[[2]], lz}, 
   If[next1[[3]] > lz, 
    next1 = {next1[[1]], next1[[2]], 
      0}]]; (*apply periodic boundary conditions z*)
  AppendTo[pol, next1];
  Do[
   direction = 
    next - pol[[-2]];(*Check the direction of the previous step*)
   stepnew = 
    DeleteCases[
     step, -direction];(*Generate new list of step directions, 
   removing the backwards step*)
   stepnew = 
    DeleteCases[stepnew, 
     direction];(*Generate new list of step directions, 
   removing the forwards step*)
   stepnew = 
    AppendTo[stepnew, 
     direction]; (*Add the forwards step as last entry*)
   ;
   next = 
    pol[[-1]] + 
     RandomChoice[{(1 - pf)/4, (1 - pf)/4, (1 - pf)/4, (1 - pf)/4, 
        pf} -> stepnew];
   next1 = next;
   If[next1[[1]] < 0, next1 = {l, next1[[2]], next1[[3]]}, 
    If[next1[[1]] > l, 
     next1 = {0, next1[[2]], 
       next1[[3]]}]]; (*apply periodic boundary conditions x*)
   If[next1[[2]] < 0, next1 = {next1[[1]], l, next1[[3]]}, 
    If[next1[[2]] > l, 
     next1 = {next1[[1]], 0, 
       next1[[3]]}]]; (*apply periodic boundary conditions y*)
   If[next1[[3]] < 0, next1 = {next1[[1]], next1[[2]], lz}, 
    If[next1[[3]] > lz, 
     next1 = {next1[[1]], next1[[2]], 
       0}]]; (*apply periodic boundary conditions z*)
   If[MemberQ[pol,next1],Break[]];(*check if its self avoiding, 
   if not, retry*)
   AppendTo[pol, next1];
   , {j, 3, Np}];
  If[Length[pol] == Np, AppendTo[allSAW, pol]]; (*Save all walks*)];

Kind regards

$\endgroup$
3
  • 1
    $\begingroup$ There are some lines in here that make no sense: RandomChoice[{(1 - pf)/4, (1 - pf)/4, (1 - pf)/4, (1 - pf)/4, pf} -> stepnew]; - is this supposed to be the weighted random choice? Also you're using AppendTo and this is slow - better to use Reap+Sow or the linked list data structure. And you should break this down into functions because it's hard to follow with all the conditionals. $\endgroup$
    – flinty
    Dec 14, 2022 at 16:07
  • $\begingroup$ @flinty, yes this is the weighted random choice. Its very crude, I am aware, but I check for the previous direction, then remove the backwards step and add the forwards step as the last entry in the array stepnew. Then I choose the next step with bias using RandomChoice, where the bias is weighted by pf for the forwards step and (1-pf)/4 for the backwards step such that all weights add up to one. $\endgroup$
    – M . M
    Dec 15, 2022 at 6:34
  • $\begingroup$ (1-pf)/4 for the perpendicular steps ofcourse* $\endgroup$
    – M . M
    Dec 15, 2022 at 6:40

1 Answer 1

1
$\begingroup$

This isn't perfect by any means and you should check it's doing the right thing. It's about 30x faster than your implementation, but can still be slow because with long paths, it becomes more and more likely the walk must be 'reset' when surrounded by occupied voxels in the SparseArray and unable to continue.

wrap[pos_, lSize_] := Thread[Mod[pos, lSize]];

(* directions: last one is the initial forwards direction *)
globalDirections = {{-1,0,0},{1,0,0},{0,-1,0},{0,1,0},{0,0,-1},{0,0,1}};

(* Check if all the surrounding positions are occupied.
  Note:Mathematica is 1-indexed so add 1 when extracting *)
surrounded[lat_, pos_] := 
  Sum[Extract[lat, 1+wrap[pos+d, Dimensions@lat]],
     {d, globalDirections}] == Length[globalDirections];

generateWalk[latticeSize_, steps_, globalDirections_, pf_] :=
  Module[{lattice, position, initialPosition, testPosition, n, walk, 
    surroundings, localDirections, testDir,
    probabilities = {(1 - pf)/4, (1 - pf)/4, (1 - pf)/4, (1 - pf)/4, (1 - pf)/4, pf}},
   While[True,
    (* setup a walk attempt with empty lattice *)
    lattice = SparseArray[{}, latticeSize];
    initialPosition = RandomInteger[{0, #}] & /@ latticeSize;
    position = initialPosition;
    testPosition;
    (* Initially the local directions are the same as the global directions, so first move is most likely {0,0,1} *)
    localDirections = globalDirections;
    n = 0;

    walk = Reap[
       Sow[initialPosition];
       (* Finish if n==steps or reset if a position is surrounded *)
       While[n<steps && Not[surrounded[lattice, position]],
        (* Generate a candidate direction from the distribution, 
           wrap it with the boundary, and get next position *)
        testDir = RandomChoice[probabilities -> localDirections];
        testPosition = wrap[position + testDir, latticeSize];

        (* if it's occupied then try again, else update the walk,  (Mathematica is 1-indexed so add 1 when extracting) *)
        If[Extract[lattice, 1 + testPosition] > 0, Continue[],
         (* move the chosen direction (the new forwards) to the end of the localDirections *)
         
         localDirections = Append[DeleteCases[globalDirections,testDir],testDir];
         position = testPosition;
         (* mark the new position as occupied in the lattice *)
         lattice = ReplacePart[lattice, (1 + position) -> 1];
         Sow[position];
         ++n;
         ];
        ]
       ][[2, 1]];
    If[n == steps, Return[walk]];
    ]
   ];

allSAW = 
  ParallelTable[
   generateWalk[{40, 40, 75}, 250, globalDirections, 0.5], {10^4}];
$\endgroup$
4
  • $\begingroup$ This will take me a little while to digest, seeing what exactly is being done here. As far as I can see, the bias with regards to the previous direction is not being taken into account (so if the direction from step -2 to -1 was {0,0,1} the weight to go forwards in that direction is $pf$). This should be easy enough for me to add I think. $\endgroup$
    – M . M
    Dec 15, 2022 at 8:01
  • $\begingroup$ @M.M I see - I misread that part as moving in the z direction globally, not the local forwards direction. I'll adjust my answer. $\endgroup$
    – flinty
    Dec 15, 2022 at 13:21
  • $\begingroup$ @M.M you may want to look into the Pivot Algorithm if you want to generate these walks more efficiently. It looks like it works by taking several short walks you pre-generate, and you flip those segments of the walk to generate new walks. This should even preserve your biased-forward property. With a periodic boundary, I guess you'd flip in a normal non-periodic space, then apply the periodic wrapping afterwards to make the walk stay inside the periodic boundary again. Have a look at this webgl demo $\endgroup$
    – flinty
    Dec 15, 2022 at 14:15
  • $\begingroup$ Will do! Thank you very much! $\endgroup$
    – M . M
    Dec 15, 2022 at 16:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.